Generalidades

En esta entrega se hace una verificación complementaria del caudal. Se incluyen estaciones de caudal que fueron descartadas en el proceso de calibración/validación.

La verificación se hace a escala mensual en el período 2006-2015 (2005 de calentamiento). Esta verificacion muetra a grandes trazos que en los primeros años de simulacion (2006-2010) predomina la subestimación. Luego a partir del año 2010 se sobreestima el caudal.

A continuación leemos los resultados del modelo ejecutado en ese período así como los valores observados. Nótese que la simulación de caudal tiene de 100 ejecuciones del modelo, los estadísticos se dan en funcion de la mediana de las 100 simulaciones.

flow_m23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m23.RDS")
flow_m30 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m30.RDS")
flow_m36 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m36.RDS")
flow_m38 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m38.RDS")
flow_m41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m41.RDS")

flow_m23 = lapply(flow_m23[,-1],function(x){convertFlow(zoo(x, flow_m23$date),
                                                        from="cumecs", to="mm", area.km2 = 687,
                                                        timestep.default = "months")}) 
flow_m30 = lapply(flow_m30[,-1],function(x){convertFlow(zoo(x, flow_m30$date),
                                                        from="cumecs", to="mm", area.km2 = 1092,
                                                        timestep.default = "months")}) 
flow_m36 = lapply(flow_m36[,-1],function(x){convertFlow(zoo(x, flow_m36$date),
                                                        from="cumecs", to="mm", area.km2 = 2744,
                                                        timestep.default = "months")})
flow_m38 = lapply(flow_m38[,-1],function(x){convertFlow(zoo(x, flow_m38$date),
                                                        from="cumecs", to="mm", area.km2 = 3159,
                                                        timestep.default = "months")}) 
flow_m41 = lapply(flow_m41[,-1],function(x){convertFlow(zoo(x, flow_m41$date),
                                                        from="cumecs", to="mm", area.km2 = 4896,
                                                        timestep.default = "months")})

Qsim = cbind(Qsim_23 = apply(as.data.frame(flow_m23),1,median),
             Qsim_30 = apply(as.data.frame(flow_m30),1,median),
             Qsim_36 = apply(as.data.frame(flow_m36),1,median),
             Qsim_38 = apply(as.data.frame(flow_m38),1,median),
             Qsim_41 = apply(as.data.frame(flow_m41),1,median))


Qobs <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/runoff_mensual_2006_2015.RDS")

qo = cbind(Qobs, Qsim)

La siguiente tabla muestra los estadísticos KGE, NSE y Pbias a escala mensual en el período 2006-2015.

library(hydroGOF)

kge_stat = round(c(KGE(qo$Qsim_23,qo$Qobs_23),
             KGE(qo$Qsim_30,qo$Qobs_30),
             KGE(qo$Qsim_36,qo$Qobs_36),
             KGE(qo$Qsim_38,qo$Qobs_38),
             KGE(qo$Qsim_41,qo$Qobs_41)),2)

nse_stat = round(c(NSE(qo$Qsim_23,qo$Qobs_23),
             NSE(qo$Qsim_30,qo$Qobs_30),
             NSE(qo$Qsim_36,qo$Qobs_36),
             NSE(qo$Qsim_38,qo$Qobs_38),
             NSE(qo$Qsim_41,qo$Qobs_41)),2)

pbias_stat = round(c(pbias(qo$Qsim_23,qo$Qobs_23),
             pbias(qo$Qsim_30,qo$Qobs_30),
             pbias(qo$Qsim_36,qo$Qobs_36),
             pbias(qo$Qsim_38,qo$Qobs_38),
             pbias(qo$Qsim_41,qo$Qobs_41)),1)

stat = rbind(KGE=kge_stat,NSE=nse_stat, PBIAS=pbias_stat)

kable(stat, col.names = c("WL_23", "WL_30", "WL_36", "WL_38", "WL_41")) %>%
  kable_styling(c("striped", "bordered"))

res = cbind(r23 = (qo$Qsim_23 - qo$Qobs_23),
            r30 = (qo$Qsim_30 - qo$Qobs_30),
            r36 = (qo$Qsim_36 - qo$Qobs_36),
            r38 = (qo$Qsim_38 - qo$Qobs_38),
            r41 = (qo$Qsim_41 - qo$Qobs_41))

Res_plot = function(x){
  
  cc = substr(index(res),1,3) %in% "ene"
  nn = index(res)
  nn[!cc] = NA 
  nn = substr(nn,6,10)
  
  barplot(res[,x], las=1, col=ifelse(res[,x]>=0,4,2),
          main = x, ylab="Qsim-Qobs (mm)", names=nn, ylim = c(-100,100))
} 

par(mfrow=c(5,1), mar=c(4,4,1,1))
Res_plot("r23")
Res_plot("r30")
Res_plot("r36")
Res_plot("r38")
Res_plot("r41")