Última actualización: 27 mayo 2020

En esta entrega se muestra la producción de sedimentos en la cuenca simulada con el esquema2.

La siguiente figura resume 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_27/SDYL_hru.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 función 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)"

Para AGRC, AGRP, LECH no hay clases de pendiente S3, S4 y S5. Sus valores aumentan según la pendiente.

Para GRAS, existen todas las clases de pendiente, llama la atención que la producción de sediemntos disminuye en las clases S3, S4 y S5.

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"))

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)

A work by Rafael Navas