DWBmodelUN examples in a tropical basin

DWBCalculator function

DWBCalculator is the function responsible for incorporating the DWB model calculations into the DWBmodelUN package. In this connection, the following example corresponds to an ungauged modeling exercise developed for the Sogamoso River basin (SRB) by (Duque, 2018).

As a first step, the example starts with loading datasets generated for SRB that are required by DWB model and subsequently, the coordinates system of the input data are verified. The last step is to run DWBCalculator function in order to obtain a monthly streamflow simulated.

library(DWBmodelUN)
library(raster)
#> Warning: package 'raster' was built under R version 3.6.3
#> Loading required package: sp
#> Warning: package 'sp' was built under R version 3.6.3
# Verify that the coordinates of the databases match
Coord_comparison(P_sogamoso, PET_sogamoso)
#> First data file is a data frame - Converting to raster
#> Second data file is a data frame - Converting to raster
#> Two data frames - Comparing headers
#> First date does not match according to the header
#> Warning in Coord_comparison(P_sogamoso, PET_sogamoso): The model will run, but
#> WARNING - Please verify dates
#> Extent verified
#> Resolution verified
# Load geographic info of GRU and parameters per cell
data(GRU, param)
# Construction of parameter maps from values by GRU
GRU.maps <- buildGRUmaps(GRU, param)
alpha1_v <- GRU.maps$alpha1
alpha2_v <- GRU.maps$alpha2
smax_v <- GRU.maps$smax
d_v <- GRU.maps$d

# Establish the initial modeling conditions
init <- init_state(GRU.maps$smaxR)
#> Warning in init_state(GRU.maps$smaxR): Strange number of initial state files
#>  Review files of initial states 
#>  Creation by default from first raster
g_v <- init$In_ground
s_v <- init$In_storage
rm(init)

# Load general characteristics of modeling
setup_data <- readSetup(Read = TRUE)
Dates <- seq(as.Date( gsub('[^0-9.]','',colnames(P_sogamoso)[3]), format = "%Y.%m.%d"), 
             as.Date(gsub('[^0-9.]','',tail(colnames(P_sogamoso),1)) , format = "%Y.%m.%d"), by = "month")
Start.sim <- which(Dates == setup_data[8,1]); End.sim <- which(Dates == setup_data[10,1])
# Sim.Period: the 1st two columns of the P and PET are the coordinates of the cells
Sim.Period <- c(Start.sim:End.sim)+2  

# Run DWB model
DWB.sogamoso <- DWBCalculator(P_sogamoso[ ,Sim.Period], 
                              PET_sogamoso[ ,Sim.Period],
                              g_v, s_v, alpha1_v, alpha2_v, smax_v, d_v)

Calibration DWB model with Dynamical Dimension Search algorithm

dds is the function that implements the calibration algorithm Dynamical Dimension Search by (Tolson and Shoemaker, 2007). In the example, the DWB model parameters are calibrated in the SRB with 10 Group Response Unit (GRU). The Nash-Sutcliffe efficiency (nse.cof) is created as an objective function between simulated and observed streamflow data, with the treatment of missing values and its coupled with the DWBCalculator function. Finally, the calibration is done with 2 iterations using the databases of the sample basin.

library(DWBmodelUN)
library(raster)
# Load P and PET databases
data(P_sogamoso, PET_sogamoso)

# Verify that the coordinates of the databases match
Coord_comparison(P_sogamoso, PET_sogamoso)
#> First data file is a data frame - Converting to raster
#> Second data file is a data frame - Converting to raster
#> Two data frames - Comparing headers
#> First date does not match according to the header
#> Warning in Coord_comparison(P_sogamoso, PET_sogamoso): The model will run, but
#> WARNING - Please verify dates
#> Extent verified
#> Resolution verified

# Load geographic info of GRU and basins where calibration will be performed
data(GRU,basins)
cellBasins <- cellBasins(GRU, basins)
#> Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
#> Raster

# Establish the initial modeling conditions
GRU.maps <- buildGRUmaps(GRU, param)
init <- init_state(GRU.maps$smaxR)
#> Warning in init_state(GRU.maps$smaxR): Strange number of initial state files
#>  Review files of initial states 
#>  Creation by default from first raster
g_v <- init$In_ground
s_v <- init$In_storage
rm(init)

# Load general characteristics of modeling
setup_data <- readSetup(Read = TRUE)
Dates <- seq(as.Date( gsub('[^0-9.]','',colnames(P_sogamoso)[3]), format = "%Y.%m.%d"), 
             as.Date(gsub('[^0-9.]','',tail(colnames(P_sogamoso),1)) , format = "%Y.%m.%d"), by = "month")

# For this calibration exercise, the last date of simulation is 
# the same as the final date of calibration
Start.sim <- which(Dates == setup_data[8,1])
End.sim <- which(Dates == setup_data[11,1])
# the first two columns of the P and PET are the coordinates of the cells
Sim.Period <- c(Start.sim:End.sim)+2 
Start.cal <- which(Dates == setup_data[9,1])
End.cal <- which(Dates == as.Date("2004-12-01"))
# the first two columns of the P and PET are the coordinates of the cells
Cal.Period <- c(Start.cal:End.cal)+2  

#Load observed runoff
data(EscSogObs)

# Function that runs the DWB model
NSE_Sogamoso_DWB <- function(parameters, P, PET, g_v,s_v, Sim.Period, EscObs, Cal.Period){
  
  parameters <- as.vector(parameters)
  # Transform the parameters to the format that the model needs
  param <- matrix(parameters, nrow = raster::cellStats(GRU,stat="max"))  
  
  # Construction of parameter maps from values by GRU
  GRU.maps <- buildGRUmaps(GRU, param)
  alpha1_v <- GRU.maps$alpha1
  alpha2_v <- GRU.maps$alpha2
  smax_v <- GRU.maps$smax
  d_v <- GRU.maps$d
  DWB.sogamoso <- DWBCalculator(P_sogamoso[ ,Sim.Period], PET_sogamoso[ ,Sim.Period],
                                g_v,s_v, alpha1_v, alpha2_v, smax_v,d_v, calibration = TRUE)
  Esc.Sogamoso <- varBasins(DWB.sogamoso$q_total, cellBasins$cellBasins)
  
  # model evaluation; in case of possible NA results in the simulation, 
  # add a conditional assingment to a very high value
  sim <- Esc.Sogamoso$varAverage[Cal.Period - 2, ]
  obs <- EscSogObs[Cal.Period - 2, ]
  
  if (sum(!is.na(sim)) == prod(dim(sim))){
    numer <- apply((sim - obs)^2, 2, sum, na.rm = TRUE)
    demom <- apply((obs - apply(obs, 2, mean, na.rm = TRUE))^2, 2, sum, na.rm = TRUE)
    nse.cof <- 1 - numer / demom
  } else {
    nse.cof <- NA
  }
  
  Perf <- (-1)*nse.cof
  if(!is.na(mean(Perf))){ 
    Mean.Perf <- mean(Perf)
  } else {Mean.Perf <- 1e100}
  return(Mean.Perf)
}

# coupling with the DDS algorithm
xBounds.df <- data.frame(lower = rep(0, times = 40), upper = rep(c(1, 2000), times = c(30, 10)))
result <- dds(xBounds.df = xBounds.df, numIter=2, OBJFUN=NSE_Sogamoso_DWB,
              P = P_sogamoso, PET = PET_sogamoso, g_v = g_v, s_v = s_v, Sim.Period = Sim.Period, 
              EscObs = EscSogObs, Cal.Period = Cal.Period)

Interactive graphics to represent outputs from the DWB model

graphDWB is the function in charge of creating the graphical inputs and outputs of DWBmodelUN. The function has three types of graphs to be generated:

The first one corresponds to a hydrograph of the precipitation entries.
The second one generates a combined chart of a rain hydrograph and streamflow hydrograph.
The third one involves the graphic number 3 plus a hydrograph of evapotranspiration.

The present example shows how to generate each chart type, using the datasets included in DWBmodelUN for the SRB. To draw one of the four chart types it is necessary to load the variables involved in each one, to adjust the input series to a time vector and assign them to a list element. Finally to run graphDWB.

First: Example with precipitation time series

library(DWBmodelUN)
library(dygraphs)
#> Warning: package 'dygraphs' was built under R version 3.6.3
data(P_sogamoso)
P.est <- ts(c(t(P_sogamoso[1, -2:-1])), star = c(2012, 1), frequency = 12)
var <- list("Precipitation" = P.est)
 graphDWB(var, tp = 1, main = "Precipitation Lat:7.0 Lon:-72.94")

Second: Example with precipitation and runoff time series

library(DWBmodelUN)
library(dygraphs)
data(P_sogamoso, simDWB.sogamoso, EscSogObs)
P.est <- ts(c(t(P_sogamoso[1, -2:-1])), star = c(2012, 1), frequency = 12)
runoff.sim <- ts(simDWB.sogamoso[c(131:192) ,1], star = c(2012, 1), frequency = 12)
runoff.obs <- ts(EscSogObs[c(131:192) ,1] , star = c(2012, 1), frequency = 12)
var <- list("Precipitation" = P.est,"Runoff.sim" = runoff.sim, "Runoff.obs" = runoff.obs)
graphDWB(var, tp = 3, main = "DWB results at Sogamoso Basin")

Third: Example with precipitation and runoff time series

library(DWBmodelUN)
library(dygraphs)
data(P_sogamoso, PET_sogamoso, simDWB.sogamoso)
P <- ts(c(t(P_sogamoso[1, -2:-1])), star = c(2012, 1), frequency = 12)
PET <- ts(c(t(PET_sogamoso[1, -2:-1])), star = c(2012, 1), frequency = 12)
runoff.sim <- ts(simDWB.sogamoso[c(131:192), 1], star = c(2012, 1), frequency = 12)
var <- list("P" = P,"PET" = PET, "Runoff.sim" = runoff.sim)
graphDWB(var, tp = 4, main = "General Comparison Sogamoso Basin")

References

Duque, Nicolás. 2018. “Estimación de Campos de Precipitación En Cuencas Hidrográficas Colombianas Con Escasez de Datos, Combinando Datos Teledetectados Y de Estaciones En Tierra, Utilizando Funciones de Kernel.” Master’s thesis, Universidad Nacional de Colombia - Sede Bogotá. http://bdigital.unal.edu.co/71663/.

Tolson, Bryan A, and Christine A Shoemaker. 2007. “Dynamically dimensioned search algorithm for computationally efficient watershed model calibration” 43: 1–16.

Examples DWB model

Camila Garcia and David Zamora

2020-08-10