Generalidades

En esta entrega se hace un diagnóstico preliminar del estado actual en cuaunto a la producción de sedimentos y fósforo en la cuenca.

Sedimentos

La siguiente figura muestra un resumen de la producción anual de sedimentos ordenado según el uso de la tierra. Los usos EUCA y MONT tienen una producción de sedimentos bastante inferior AGRP, AGRC, LECH y GRAS. Adicionalmente AGRP, LECH, y GRAS tienen un comportamiento bastante similar. AGRC está ligeramente por encima y su variabilidad practicamente duplica la variabilidad del grupo AGRP, LECH y GRAS.

SYLD <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/SYLD_promedio_10anos.RDS")
hru_info <- readRDS("~/R_projects/SWAT_calibration_ETmodis/shapefiles_RDS/hru_info.RDS")

SYLD_GRAS = unlist(SYLD[(hru_info$LANDUSE=="GRAS")])
SYLD_EUCA = unlist(SYLD[(hru_info$LANDUSE=="EUCA")])
SYLD_MONT = unlist(SYLD[(hru_info$LANDUSE=="MONT")])
SYLD_AGRC = unlist(SYLD[(hru_info$LANDUSE=="AGRC")])
SYLD_AGRP = unlist(SYLD[(hru_info$LANDUSE=="AGRP")])
SYLD_LECH = unlist(SYLD[(hru_info$LANDUSE=="LECH")])

library(plotly)
fig <- plot_ly(y = ~SYLD_AGRC, type = "box", name="AGRC")
fig <- fig %>% add_trace(y = ~SYLD_AGRP, name="AGRP")
fig <- fig %>% add_trace(y = ~SYLD_LECH, name="LECH")
fig <- fig %>% add_trace(y = ~SYLD_GRAS, name="GRAS")
fig <- fig %>% add_trace(y = ~SYLD_EUCA, name="EUCA")
fig <- fig %>% add_trace(y = ~SYLD_MONT, name="MONT")
fig <- fig %>% layout(yaxis = list(title="Producción anual de sedimentos [ton/ha]"), 
                      showlegend = FALSE,
    autosize=FALSE,
    width=500,
    height=350)
fig

A continuación se muestra la producción de sedimentos en funcion de las 6 clases de pendientes modeladas en SWAT:

Sband = paste(unique(hru_info$SLOPE_BAND))
paste0("S", 1:6,"(",Sband,")")

## [1] "S1(0-3.0)"     "S2(3.0-6.0)"   "S3(6.0-9.0)"   "S4(9.0-12.0)" 
## [5] "S5(12.0-15.0)" "S6(15.0-9999)"

slope_plot = function(Luse){
c1 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[1]
S1 = unlist(SYLD[c1])
c2 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[2]
S2 = unlist(SYLD[c2])
c3 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[3]
S3 = unlist(SYLD[c3])
c4 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[4]
S4 = unlist(SYLD[c4])
c5 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[5]
S5 = unlist(SYLD[c5])
c6 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[6]
S6 = unlist(SYLD[c6])

if(is.null(S1)){S1=0}
if(is.null(S2)){S2=0}
if(is.null(S3)){S3=0}
if(is.null(S4)){S4=0}
if(is.null(S5)){S5=0}
if(is.null(S6)){S6=0}

fig <- plot_ly(y = ~S1, type = "box", name=paste0("s1 (",table(c1)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S2, name=paste0("s2 (",table(c2)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S3, name=paste0("s3 (",table(c3)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S4, name=paste0("s4 (",table(c4)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S5, name=paste0("s5 (",table(c5)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S6, name=paste0("s6 (",table(c6)["TRUE"],")" ))
fig <- fig %>% layout(title=paste(Luse),
                      yaxis = list(title="Producción anual de sedimentos [ton/ha]",
                                   range = c(0, 70)), 
                      showlegend = FALSE,
                      autosize=FALSE,
                      width=350,
                      height=300)
fig
}


leafsync::sync(slope_plot("AGRC"),slope_plot("AGRP"),
             slope_plot("LECH"),slope_plot("GRAS"))

Producción de sedimentos por subcuencas

La siguiente imagen muestra la producción de sedimentos por subcuencas (ton/ha). Se pueden observar 2-3 regiones.

La primera regíon está en la parte alta de la cuenca. La media de la prodducción de sedimentos es baja con respecto al resto de la cuenca (figura mean). Sin embargo, esta región tiene la variabilidad mas alta, alcanzando valores de hasta tres veces la media (figura sd/mean).

La segunda region se ubica en la parte baja de la cuenca, en esta región se observa una ligera diferencia entre la margen norte y sur del río, siendo la margen norte la que tiene los valores de producción de sedimentos.

subs1 = readRDS("~/R_projects/SWAT_calibration_ETmodis/shapefiles_RDS/subs1.RDS")

SYLD_sub = matrix(ncol=2, nrow = 0)
for(sb in 1:41){
  SYLD_sub = rbind(SYLD_sub, 
                   c(mean(unlist(SYLD[hru_info$Subbasin==sb])),
                     sd(unlist(SYLD[hru_info$Subbasin==sb]))))
}

SYLD_sub = round(SYLD_sub,2)
colnames(SYLD_sub) = c("mean", "sd")

subs1$SYLD_mean = SYLD_sub[,"mean"]
subs1$SYLD_Nsd = round(SYLD_sub[,"sd"]/SYLD_sub[,"mean"],2)


library(mapview)
map_SYLD_mean = mapview(subs1, zcol="SYLD_mean", layer.name ="mean")
map_SYLD_Nsd = mapview(subs1, zcol="SYLD_Nsd", layer.name ="sd/mean")

library(leafsync)
sync(map_SYLD_mean, map_SYLD_Nsd, sync.cursor=TRUE)

Producción de Fósforo

Se cuantificará la producción de fósforo sumando el fósforo orgánico (ORGP) y el fósforo mineral soluble (SOLP).

La siguiente imagen muestra la producción media de fósforo por tipo de uso de la tierra.

SYLD <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/SOLP_plus_ORGP_HRU_10anos.RDS")
hru_info <- readRDS("~/R_projects/SWAT_calibration_ETmodis/shapefiles_RDS/hru_info.RDS")

SYLD_GRAS = unlist(SYLD[(hru_info$LANDUSE=="GRAS")])
SYLD_EUCA = unlist(SYLD[(hru_info$LANDUSE=="EUCA")])
SYLD_MONT = unlist(SYLD[(hru_info$LANDUSE=="MONT")])
SYLD_AGRC = unlist(SYLD[(hru_info$LANDUSE=="AGRC")])
SYLD_AGRP = unlist(SYLD[(hru_info$LANDUSE=="AGRP")])
SYLD_LECH = unlist(SYLD[(hru_info$LANDUSE=="LECH")])

library(plotly)
fig <- plot_ly(y = ~SYLD_AGRC, type = "box", name="AGRC")
fig <- fig %>% add_trace(y = ~SYLD_AGRP, name="AGRP")
fig <- fig %>% add_trace(y = ~SYLD_LECH, name="LECH")
fig <- fig %>% add_trace(y = ~SYLD_GRAS, name="GRAS")
fig <- fig %>% add_trace(y = ~SYLD_EUCA, name="EUCA")
fig <- fig %>% add_trace(y = ~SYLD_MONT, name="MONT")
fig <- fig %>% layout(yaxis = list(title="Producción anual de fósforo (ORGP + SOLP) [kg/ha]"), 
                      showlegend = FALSE,
                      autosize=FALSE,
                      width=500,
                      height=350)
fig

A continuación se muestra la producción de fósforo en función del tipo de uso de la tierra y la pendiente

slope_plot = function(Luse){
c1 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[1]
S1 = unlist(SYLD[c1])
c2 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[2]
S2 = unlist(SYLD[c2])
c3 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[3]
S3 = unlist(SYLD[c3])
c4 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[4]
S4 = unlist(SYLD[c4])
c5 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[5]
S5 = unlist(SYLD[c5])
c6 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[6]
S6 = unlist(SYLD[c6])

if(is.null(S1)){S1=0}
if(is.null(S2)){S2=0}
if(is.null(S3)){S3=0}
if(is.null(S4)){S4=0}
if(is.null(S5)){S5=0}
if(is.null(S6)){S6=0}

fig <- plot_ly(y = ~S1, type = "box", name=paste0("s1 (",table(c1)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S2, name=paste0("s2 (",table(c2)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S3, name=paste0("s3 (",table(c3)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S4, name=paste0("s4 (",table(c4)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S5, name=paste0("s5 (",table(c5)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S6, name=paste0("s6 (",table(c6)["TRUE"],")" ))
fig <- fig %>% layout(title=paste(Luse),
                      yaxis = list(title="Producción anual de Fósforo [kg/ha]",
                                   range = c(0, 6)), 
                      showlegend = FALSE,
                      autosize=FALSE,
                      width=350,
                      height=300)
fig
}


leafsync::sync(slope_plot("AGRC"),slope_plot("AGRP"),
             slope_plot("LECH"),slope_plot("GRAS"))

La siguiente imagen muestra la producción de fósforo en la cuenca. A comentar…

subs1 = readRDS("~/R_projects/SWAT_calibration_ETmodis/shapefiles_RDS/subs1.RDS")

SYLD_sub = matrix(ncol=2, nrow = 0)
for(sb in 1:41){
  SYLD_sub = rbind(SYLD_sub, 
                   c(mean(unlist(SYLD[hru_info$Subbasin==sb])),
                     sd(unlist(SYLD[hru_info$Subbasin==sb]))))
}

SYLD_sub = round(SYLD_sub,2)
colnames(SYLD_sub) = c("mean", "sd")

subs1$SYLD_mean = SYLD_sub[,"mean"]
subs1$SYLD_Nsd = round(SYLD_sub[,"sd"]/SYLD_sub[,"mean"],2)


library(mapview)
map_SYLD_mean = mapview(subs1, zcol="SYLD_mean", layer.name ="mean")
map_SYLD_Nsd = mapview(subs1, zcol="SYLD_Nsd", layer.name ="sd/mean")

library(leafsync)
sync(map_SYLD_mean, map_SYLD_Nsd, sync.cursor=TRUE)