---
title: "Examples DWB model"
author: "Camila Garcia and David Zamora"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Examples DWB model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

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

```{r echo=TRUE, results='hide'}
library(DWBmodelUN)
# Load P and PET databases
data(P_sogamoso, PET_sogamoso)
# Verify that the coordinates of the databases match
Coord_comparison(P_sogamoso, PET_sogamoso)
# 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)
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.

```{r echo=TRUE, results='hide'}
library(DWBmodelUN)
# Load P and PET databases
data(P_sogamoso, PET_sogamoso)

# Verify that the coordinates of the databases match
Coord_comparison(P_sogamoso, PET_sogamoso)

# Load geographic info of GRU and basins where calibration will be performed
data(GRU,basins)
cellBasins <- cellBasins(GRU, basins)

# Establish the initial modeling conditions
GRU.maps <- buildGRUmaps(GRU, param)
init <- init_state(GRU.maps$smaxR)
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
  nGRU <- max(terra::values(terra::rast(GRU)), na.rm = TRUE)
  param <- matrix(parameters, nrow = nGRU)
  
  # 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 assignment to a very high value
  sim <- Esc.Sogamoso$varAverage[Cal.Period - 2, ]
  # align the observations with the simulated basins (named by gauge code)
  obs <- EscSogObs[Cal.Period - 2, paste0("X", colnames(sim))]
  
  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 four types of graphs to be generated:

1. The first one (`tp = 1`) corresponds to a hydrograph of any continuous variable, such as the precipitation entries.
2. The second one (`tp = 2`) compares the simulated runoff with the observations.
3. The third one (`tp = 3`) generates a combined chart with a rain hyetograph and the simulated and observed streamflow hydrographs.
4. The fourth one (`tp = 4`) adds a hydrograph of evapotranspiration to the previous comparison.

The present example shows how to generate three of the four chart types, using the datasets included in **DWBmodelUN** for the SRB. To draw one of the 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

```{r echo=TRUE, fig.align= "center", fig.height= 4, fig.width= 8}
library(DWBmodelUN)
library(dygraphs)
data(P_sogamoso)
P.est <- ts(c(t(P_sogamoso[1, -2:-1])), start = c(2001, 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

```{r echo=TRUE, fig.align= "center", fig.height=5, fig.width=7, dpi = 600}
library(DWBmodelUN)
library(dygraphs)
data(P_sogamoso, simDWB.sogamoso, EscSogObs)
P.est <- ts(c(t(P_sogamoso[1, -2:-1])), start = c(2001, 1), frequency = 12)
runoff.sim <- ts(simDWB.sogamoso[, 1], start = c(2001, 1), frequency = 12)
runoff.obs <- ts(EscSogObs[, 1], start = c(2001, 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

```{r echo=TRUE, fig.align= "center", fig.height=5, fig.width=7, dpi = 600}
library(DWBmodelUN)
library(dygraphs)
data(P_sogamoso, PET_sogamoso, simDWB.sogamoso)
P <- ts(c(t(P_sogamoso[1, -2:-1])), start = c(2001, 1), frequency = 12)
PET <- ts(c(t(PET_sogamoso[1, -2:-1])), start = c(2001, 1), frequency = 12)
runoff.sim <- ts(simDWB.sogamoso[, 1], start = c(2001, 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á. https://repositorio.unal.edu.co/handle/unal/71663.

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