# Examples DWB model

#### Camila Garcia and David Zamora

#### 2020-08-10

# 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:

1. The first one corresponds to a hydrograph of the precipitation entries.
2. The second one generates a combined chart of a rain hydrograph and streamflow hydrograph.
3. 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.