######################### HOMEWORK 1 CI71F: USEFUL COMMANDS ############################# # # Created by Pablo Mendoza, July 2017 # # Description: # # This script includes some commands to process and plot hydrometeorological data # required by GR models, and also to run GR models and display results. # # # Updates: # # July 30, 2019 (PM) # - Overall revision of script. Add a few lines to improve visualization # # # ##### Load libraries #################################################################### # If the package is not installed, run the following command # install.packages("airGR") require(airGR) # Define months Months = c('Jan','Feb','Mar','Apr','May','Jun', 'Jul','Aug','Sep','Oct','Nov','Dec') Months2 = c('J','F','M','A','M','J','J','A','S','O','N','D') ############ Load hydrometeorological datasets ### Basin 1 data(L0123001) # Time series with observations (Precip in mm/d, temp in C, # PET in mm/d, and streamflow) DataBas1TS <- BasinObs # Basin properties: Code, name, area (km2) and hypsometric curve # including min, quantiles from q1 to q99 and max (meters) DataBas1H <- BasinInfo ### Basin 2 data(L0123002) DataBas2TS <- BasinObs DataBas2H <- BasinInfo rm(BasinObs, BasinInfo) # Remove ancillary variables ############ Example plots for basin 2 par(mfrow=c(2,1), mar =c(4,4,1,1)) # plot with 2 rows and 1 column (i.e., 2 panels) ### Daily precipitation plot(DataBas2TS$DatesR, DataBas2TS$P, xlab = "Time", ylab = "Precipitation (mm/d)", type = "l", col = "blue", las = 1) ### Daily precipitation plot(DataBas2TS$DatesR, DataBas2TS$Qmm, xlab = "Time", ylab = "Runoff (mm/d)", type = "l", col = "red", las = 1) ### Monthly cycles # First, compute monthly amounts (P,PET,R) or averages (T) Tyears <- as.numeric(format(DataBas2TS$DatesR,"%Y")) # Time series with years only Tmonths <- as.numeric(format(DataBas2TS$DatesR,"%m")) # Time series with months only Tdays <- as.numeric(format(DataBas2TS$DatesR,"%d")) # Time series with days only # Define variables for number of years and number of months Nyears <- length(unique(Tyears)) Nmonths <- 12 # Initialize matrices with monthly values MTemp <- array(NA, dim=c(Nyears,Nmonths)) # Temperature MPrecip <- array(NA, dim=c(Nyears,Nmonths)) # Precipitation MPET <- array(NA, dim=c(Nyears,Nmonths)) # Potential ET MRobs <- array(NA, dim=c(Nyears,Nmonths)) # Monthly runoff # Before computing averages, examine that we don´t have NaN or NA which(is.na(DataBas2TS$T)) which(is.nan(DataBas2TS$T)) # Repeat these commands for all variables. Output should be zero! # Compute monthly values per year for (iyear in 1:Nyears){ # Start loop over years for (imonth in 1:Nmonths) { # Start loop over months mindex = intersect(which(Tyears == unique(Tyears)[iyear]), which(Tmonths == imonth)) MTemp[iyear,imonth] <- mean(DataBas2TS$T[mindex]) # Mean monthly temperature MPrecip[iyear,imonth] <- sum(DataBas2TS$P[mindex]) # Montly precipitation MPET[iyear,imonth] <- sum(DataBas2TS$E[mindex]) # Monthly PET MRobs[iyear,imonth] <- sum(DataBas2TS$Qmm[mindex]) # Monthly Runoff } # End loop over months } # End loop over years # Compute mean monthly values for the entire period Monthly_T <- apply(MTemp, 2, mean) Monthly_P <- apply(MPrecip, 2, mean) Monthly_PET <- apply(MPET, 2, mean) Monthly_Robs <- apply(MRobs, 2, mean) # Now plot the results! par(mfrow=c(2,1)) # plot with 2 rows and 1 column (i.e., 2 panels) # First we plot temperature plot(1:Nmonths,Monthly_T,type="l",xlab="",ylab="Temperature (ºC)", las=1,col="black",xaxt="n") axis(1,at=1:Nmonths,Months2) # Second: rest of variables plot(1:Nmonths,Monthly_P,type="l",xlab="",ylab="Variable (mm)", ylim=range(Monthly_P,Monthly_PET,Monthly_Robs),las=1,col="blue",xaxt="n") lines(1:Nmonths,Monthly_PET,xlab="",ylab="Variable (mm)", las=1,col="red",xaxt="n") lines(1:Nmonths,Monthly_Robs,xlab="",ylab="Variable (mm)", las=1,col="green",xaxt="n") axis(1,at=1:Nmonths,Months2) legend("topright",c("P","PET","R"),horiz="TRUE",bty="n",cex = 0.8, lty=rep(1,3),lwd=rep(2,3),col=c("blue","red","darkgreen") ) # Compute and plot a flow duration curve par(mfrow=c(1,1)) NdaysQ <- length(DataBas2TS$Qmm) # Number of days with flows ObsFDC <- sort(DataBas2TS$Qmm, decreasing = TRUE) # Sort flows from largest so smallest Pexc <- (1:NdaysQ)/(NdaysQ+1) # Probability of exceedance (Weibull formula) plot(Pexc, ObsFDC, type="l",log="y", xlab="Exceedance probability", ylab="Streamflow (mm/d)",las = 1) abline( v = c(0.02, 0.2, 0.7), lty = 2 ) # Draw vertical lines to separate segments ################### Examples on how to run GR4J ######################################### ## Create input forcing files # Basin 1 InputsModelBas14J <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = DataBas1TS$DatesR, Precip = DataBas1TS$P, PotEvap = DataBas1TS$E) # Basin 2 InputsModelBas24J <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = DataBas2TS$DatesR, Precip = DataBas2TS$P, PotEvap = DataBas2TS$E) ## Define simulation period Ind_Run <- seq(which(format(DataBas2TS$DatesR, format = "%d/%m/%Y %H:%M")=="01/01/1990 00:00"), which(format(DataBas2TS$DatesR, format = "%d/%m/%Y %H:%M")=="31/12/1999 00:00")) ## preparation of the RunOptions object for both basins, GR4J model # Basin 1 RunOptions4JBas1 <- CreateRunOptions(FUN_MOD = RunModel_GR4J, # Specify model InputsModel = InputsModelBas14J, # Forcing dataset IndPeriod_Run = Ind_Run) # Simulation period # Basin 2 RunOptions4JBas2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModelBas24J, IndPeriod_Run = Ind_Run) ## Run the model GR4J with a pre-defined set of parameter values Param <- c(257.238, 1.012, 88.235, 2.208) OutputsGR4JBas1 <- RunModel_GR4J(InputsModel = InputsModelBas14J, RunOptions = RunOptions4JBas1, Param = Param) # Run for basin 1 OutputsGR4JBas2 <- RunModel_GR4J(InputsModel = InputsModelBas24J, RunOptions = RunOptions4JBas2, Param = Param) # Run for basin 2 str(OutputsGR4JBas1) # Scatter plot of simulated versus observed flows par(mfrow=c(1,1)) plot(DataBas2TS$Qmm[Ind_Run],OutputsGR4JBas2$Qsim,las = 1, xlab= expression(paste("Q"["obs"]," (mm/d)")), ylab= expression(paste("Q"["sim"]," (mm/d)"))) lines(c(-10000,10000),c(-10000,10000)) # Add 1:1 line ######## Include CemaNeige in simulations! # Define vector with parameters. The last two values are Cemma Neige parameters ParamSnow <- c(257.238, 1.012, 88.235, 2.208, 0.962, 2.249) InputsModelBas24JCemma <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR4J, DatesR = DataBas2TS$DatesR, Precip = DataBas2TS$P, PotEvap = DataBas2TS$E, TempMean = DataBas2TS$T, ZInputs = median(DataBas2H$HypsoData), HypsoData = DataBas2H$HypsoData, NLayers = 5) RunOptions4JBas2 <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModelBas24JCemma, IndPeriod_Run = Ind_Run) OutputsGR4JBas2Cemma <- RunModel_CemaNeigeGR4J(InputsModel = InputsModelBas24JCemma, RunOptions = RunOptions4JBas2, Param = ParamSnow) str(OutputsGR4JBas2Cemma) # Scatter plot of simulated versus observed flows plot(DataBas2TS$Qmm[Ind_Run],OutputsGR4JBas2Cemma$Qsim,las = 1, xlab= expression(paste("Q"["obs"]," (mm/d)")), ylab= expression(paste("Q"["sim"]," (mm/d)"))) lines(c(-10000,10000),c(-10000,10000)) # Add 1:1 line ##################################### END OF SCRIPT #####################################