Modelo GR4J en una cuenca sintética

Iniciamos limpiando las variables del entorno de trabajo y definiendo las librerías a ser utilizadas

rm(list=ls())

library(airGRteaching)
library(tcltk)
library(Hmisc)
library(SoilHyP)
library(lubridate)
library(zoo)
library(airGR)

Defina la ruta de trabajo

path <- "Z:/My Drive/PERSONALES/TRABAJOS/DOCENCIA/2023 CURSO MODELACION UNICAUCA/2- MATERIAL/CLASES/SESION 4/TALLER GR4J"
setwd(path)

Defina el área de la cuenca.

Area_cuenca_km2 <- 156.6
Area_cuenca_m2 <- Area_cuenca_km2*1000000

Cargue el archivo con las variables de entrada al modelo

Entradas <- read.csv("DatosDiarios_CuencaSintetica.csv")
Entradas <- data.frame("DatesR"=Entradas$Fecha, "P"=Entradas$P_mm, "E"=Entradas$Etp_mm, "Qmm"=Entradas$Q_m3.s)
Entradas$Qmm <- Entradas$Qmm/Area_cuenca_m2*(60*60*24)*(1000)

Entradas$DatesR <- as.POSIXct(Entradas$DatesR, format="%d/%m/%Y", tz="UTC")

Defina el número de corridas

n <- 100 #numero de corridas

Defina los periodos de modelación

per.cal_warmup <- c("2003-01-01", "2003-12-31")
per.cal_warmup_extend <- seq(as.Date(per.cal_warmup[1]), as.Date(per.cal_warmup[2]), by="day")

per.cal <- c("2004-01-01", "2011-12-31")
per.cal_extend <- seq(as.Date(per.cal[1]), as.Date(per.cal[2]), by="day")

per.val_warmup <- c("2011-01-01", "2011-12-31")
per.val_warmup_extend <- seq(as.Date(per.val_warmup[1]), as.Date(per.val_warmup[2]), by="day")

per.val <- c("2012-01-01", "2015-12-31")
per.val_extend <- seq(as.Date(per.val[1]), as.Date(per.val[2]), by="day")

Defina el límite de los parámetros

## Simulation step using model parameters set by the user
par_names <- c("X1", "X2", "X3", "X4")
lim_inf <- c(20, -5,20,1)
lim_sup <- c(1200, 3, 500, 5)

Genere los conjuntos de parámetros a través del método Monte Carlo

# Montecarlo Parameters set - uniform distribution
parameters_set <- matrix(ncol = 4, nrow = n)
colnames(parameters_set) <- par_names

for(i in 1:length(par_names)){
  parameters_set[,i] <- runif(n, min=lim_inf[i], max=lim_sup[i])
}

Realice la calibración con las simulaciones Monte Carlo.

#aqui vamos a guardar los resultados de las n simulaciones y las n FO
Qsim <- list()
FO <- list()

PREP <- PrepGR(ObsDF = Entradas, HydroModel = "GR4J", CemaNeige = FALSE)

#barra de avance
pb <- tkProgressBar(title = "Calibrando GR4J", min = 0,
                    max = n, width = 300)

for(j in 1:n){
  
  SIM <- SimGR(PrepGR = PREP, Param = parameters_set[j,], EffCrit = "NSE",
               WupPer = per.cal_warmup, SimPer = per.cal )
  Qsim[[j]] <- SIM$OutputsModel$Qsim
  FO[[j]] <- SIM$EffCrit$CritValue
  
    #Darle paso a la barra de avance
  setTkProgressBar(pb, j, label=paste( round(j/n*100, 0),
                                       "% done"))
  
}

#Cerrar la barra de avance
close(pb)

Organice los resultados de las simulaciones

#Save Results

Qsim <- as.data.frame(do.call(cbind, Qsim))
Qsim <- cbind(seq(as.Date(per.cal[1]), as.Date(per.cal[2]), by="day"), Qsim)

FO <- do.call(rbind, FO)

Evalue la incertidumbre con la metodología GLUE

threshold <- 0.4 #limite para definir corridas comportamentales
parameters_beh <- parameters_set[FO > threshold,]
NSE_beh <- FO[FO > threshold]
Qsim_cal_beh <- Qsim[, FO > threshold]

Qsim_best <- Qsim[,(which.max(FO)+1)]
Qsim_best_m3s <- Qsim_best/(1000*60*60*24)*Area_cuenca_m2

weights <- NSE_beh - threshold
weights <- weights / sum(weights)

#con esto construyo la franja de 0.90 
limits <- apply(Qsim_cal_beh, 1, "wtd.quantile", weights = weights, probs = c(0.05,0.95), normwt=T)

limits_m3s <- limits/(1000*60*60*24)*Area_cuenca_m2

Qobs_cal <- as.Date(Entradas$DatesR)
Qobs_cal <- which(( Qobs_cal %in% per.cal_extend)==T)
Qobs_cal <- Entradas$Qmm[Qobs_cal]
Qobs_cal_m3s <- Qobs_cal/(1000*60*60*24)*Area_cuenca_m2

Grafico el intervalo de confianza

jpeg(filename = paste0( path,"/GLUE.jpeg"), width = 15, height = 10, units = "cm", res = 300) 
par(par(mfrow=c(2,2), mar=c(1, 5, 4, 1)+0.1, oma=c(0,0,0,0)+0.1))
plot(Qobs_cal_m3s, type="l", ylab="Q (m3/s)", xlab="", xaxt="n",
     main = "Banda de confianza")
lines(limits_m3s[1,], type="l", col="blue", lty=2)
lines(limits_m3s[2,], type="l", col="red", lty=3)
lines(Qsim_best_m3s, type="l", col="gray40")
lab <- seq(1,length(per.cal_extend),365)
axis(1,at=lab ,labels = year(per.cal_extend[lab]), las=2)
legend( "topleft",c("Qsim 95%","Qsim 5%", "Qsim","Qobs"), col = c("blue","red", "gray40", "black"), lty = c(2,3,1,1),
        cex=0.9, lwd=1)
dev.off()

Identifique la mejor FO y el mejor conjunto de parámetros

NSE_max <- max(FO)
NSE_max_pars <- parameters_set[which(FO==NSE_max),]

Haga sus dotty plots

#Dotty plots
jpeg(filename = paste0( path,"/DOTTY_PLOTS.jpeg"), width = 15, height = 10, units = "cm", res = 300) 
par(mfrow = c(2, 2), mar=rep(3,4), oma=c(0.5,1,0.5,0.5))
plot(parameters_set[,1], 1-FO, ylim=c(0,2), xlab="", ylab="", cex.lab=1, cex.axis=1) ; mtext("1-NSE",2, line = 2) ; mtext("X1",1, line = 2)
points(NSE_max_pars[1], 1-NSE_max, col="red", pch=15, cex=1)
plot(parameters_set[,2], 1-FO, ylim=c(0,2), xlab="", ylab="", cex.lab=1, cex.axis=1) ; mtext("1-NSE",2, line = 2) ; mtext("X2",1, line = 2)
points(NSE_max_pars[2], 1-NSE_max, col="red", pch=15, cex=1)
plot(parameters_set[,3], 1-FO, ylim=c(0,2), xlab="", ylab="", cex.lab=1, cex.axis=1) ; mtext("1-NSE",2, line = 2) ; mtext("X3",1, line = 2)
points(NSE_max_pars[3], 1-NSE_max, col="red", pch=15, cex=1)
plot(parameters_set[,4], 1-FO, ylim=c(0,2), xlab="", ylab="", cex.lab=1, cex.axis=1) ; mtext("1-NSE",2, line = 2) ; mtext("X4",1, line = 2)
points(NSE_max_pars[4], 1-NSE_max, col="red", pch=15, cex=1)
dev.off()

Grafique la mejor hidrógrafa de calibración

jpeg(filename = paste0( path,"/Hidro_Cal.jpeg"), width = 15, height = 10, units = "cm", res = 300)
par(par(mfrow=c(2,2),  mar=c(1, 5, 4, 1)+0.1, oma=c(0,0,0,0)+0.1))
plot(Qobs_cal_m3s, type="l", col="black", ylab="Q (m3/s)", xlab="", xaxt="n")
lines(Qsim_best_m3s, type="l", col ="red")
title(paste("NSE_Cal =",round(NSE_max,3)))
lab <- seq(1,length(per.cal_extend),365)
axis(1,at=lab ,labels = year(per.cal_extend[lab]), las=2)
legend( "topleft",c("Qsim", "Qobs"), col = c("red", "black"), lty = c(1,1))
dev.off()

Ahora se hace la validación

SIM <- SimGR(PrepGR = PREP, Param = NSE_max_pars, EffCrit = "NSE",
             WupPer = per.val_warmup, SimPer = per.val)
Qsim_val <- SIM$OutputsModel$Qsim
Qsim_val_m3s <- Qsim_val/(1000*60*60*24)*Area_cuenca_m2

FO_val <- SIM$EffCrit$CritValue

Qobs_val <- as.Date(Entradas$DatesR)
Qobs_val <- which(( Qobs_val %in% per.val_extend)==T)
Qobs_val <- Entradas$Qmm[Qobs_val]
Qobs_val_m3s <- Qobs_val/(1000*60*60*24)*Area_cuenca_m2

Grafique la validación

jpeg(filename = paste0( path,"/Hidro_Val.jpeg"), width = 15, height = 10, units = "cm", res = 300)
  plot(Qobs_val_m3s, type="l", col="black", ylab="Q (m3/s)", xlab="", xaxt="n")
lines(Qsim_val_m3s, type="l", col ="red")
title(paste("NSE_Val =",round(FO_val,3)))
lab <- seq(1,length(per.val_extend),366)
axis(1,at=lab ,labels = year(per.val_extend[lab]), las=2)
legend( "topright",c("Simulados", "Observados"), col = c("red", "black"), lty = c(1,1))
dev.off()

Guarde las series simuladas

#calibración
q_cal_save <- cbind(per.cal_extend,Qsim_best, Qsim_best_m3s)
q_cal_save <- as.data.frame(q_cal_save)
colnames(q_cal_save) <- c("Fecha", "Q(mm/dia)", "Q(m3/s)")
q_cal_save$Fecha <- as.Date(q_cal_save$Fecha)
write.csv(q_cal_save, "Resultados_calibración.csv", row.names = F)

#parámetros
write.csv(NSE_max_pars, "Parámetros.csv")

#validación
q_val_save <- cbind(per.val_extend,Qsim_val, Qsim_val_m3s)
q_val_save <- as.data.frame(q_val_save)
colnames(q_val_save) <- c("Fecha", "Q(mm/dia)", "Q(m3/s)")
q_val_save$Fecha <- as.Date(q_val_save$Fecha)
write.csv(q_val_save, "Resultados_validación.csv", row.names = F)