---
title: "Using Dynamic TOPMODEL"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using Dynamic TOPMODEL}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```
The purpose of this vignette is to provide an outline of the steps needed to
perform a Dynamic TOPMODEL simulation and introduce the formats of the data
input and returned. 

The data used in this example comes from Swindale and is contained within the package and can be loaded with
```{r, setup}
library(dynatop)
data("Swindale")
```
which returns a variable `Swindale` with the following components:
```{r data_loaded}
names(Swindale)
```

For better comparison with a likely analysis we separate these into a model
and observed data variables

```{r sep}
swindale_model <- Swindale$model
swindale_obs <- Swindale$obs
```

# The model structure
A dynamic TOPMODEL is described in a list object. The list has the following
elements
```{r model_parts}
names(swindale_model)
```
which are described in [associated
vignette](The_Model_Object.html). The
[dynatopGIS](https://waternumbers.github.io/dynatopGIS/) package
can be used for constructing models. 

# Setting map locations

While not required for simulations if the location of a file containing
the a raster map of the locations of the HRUs is provided the states can be visualised within
dynatop.

The locations of the files are set in the `map` element of the model. For this
example the maps are located within the `extdata` directory of the package and
can be set using commands

```{r set_map}
swindale_model$map <- system.file("extdata","Swindale.tif",package="dynatop",mustWork=TRUE)
```

# Preparing input data

The input to the model is expected to take the form of an ```xts``` object with
constant time step whose column names are found in the 'precip' and 'pet'
variables of the HRUs descriptions in the model. Helpful functions for creating and
manipulating ```xts``` objects can be found
[here](http://rstudio-pubs-static.s3.amazonaws.com/288218_117e183e74964557a5da4fc5902fc671.html),
see also the `resample_xts` function in this package.

The discharge, precipitation and potential evapotranspiration (PET) inputs for
Swindale are contained with `swindale_obs` on a 15 minute time step.
```{r, obs}
head(swindale_obs)
```
Note the discharge is in m$^{3}$/s while the precipitation and PET are in m
accumulated over the time step.

To use the data with the model we need to set the names of the time series
inputs within the model. In this case this is already done as can be seen by
inspecting the `precip` and `pet` values for the individual HRU's
```{r, set_obs_names}
head(swindale_model$hru[[1]]$precip)
head(swindale_model$hru[[1]]$pet)
```

# Altering parameters

The parameter values are stored within the HRU descriptions. 
Which parameters are present depends upon the options selected
for the transmissivity and channel solution. Details can be found in the
[HRU](HRU.html) Vignettes.

Altering parameter values requires changing their values in the HRU definitions. 
For this catchment all HRU have the same parameter values. For this simulation 
we change the parameter vectors to be more representative of the catchment
```{r, change_param}
hru <- swindale_model$hru
for(ii in 1:length(hru)){
    if(is.na(hru[[ii]]$class$endNode)){
        ## then HRU is not a channel
        ## saturated zone parameters
        hru[[ii]]$sz$parameters["m"] <- 0.0063
        hru[[ii]]$sz$parameters["t_0"] <- exp(7.46)
        ## unsaturated zone parameters
        hru[[ii]]$uz$parameters["t_d"] <- 8*60*60
        ## root zone parameters
        hru[[ii]]$rz$parameters["s_rzmax"] <- 0.1
        ## surface parameters
        hru[[ii]]$sf$parameters["c_sf"] <- 0.4
    }else{
        ## then HRU is a channel - set so no subsurface response
        ## saturated zone parameters
        hru[[ii]]$sz$parameters["t_0"] <- 0.001
        ## root zone parameters
        hru[[ii]]$rz$parameters["s_rzmax"] <- 0.001
        ## surface parameters
        hru[[ii]]$sf$parameters["c_sf"] <- 0.8
    }
    ## initialisation parameters
    hru[[ii]]$initialisation["s_rz_0"] <- 0.98
    hru[[ii]]$initialisation["r_uz_sz_0"] <- 1.755582e-07 ## initial outflow divided by catchment area
}
```

# Creating the dynatop Object

Simulations are performed by embedding the model and the observed data into
a `dynatop` object. First the object is created using the model in list form

```{r create_object}
ctch_mdl <- dynatop$new(hru,map=swindale_model$map)
```

This step performs some basic checks on the model for consistency.
The data can then be added

```{r add_data}
ctch_mdl$add_data(swindale_obs)
```

# Running dynamic TOPMODEL

The model can be run with the `sim` method. The methods check that suitable
input data is available, but not how it was generated.

The initial states of the simulations can be specified in the model
object. If, as in the case of this example, the states are not specified then
any attempt to perform a simulation will fail.

```{r sim_fail, error=TRUE, purl=FALSE}
ctch_mdl$sim(swindale_model$output_flux)
```

Rather then editing each HRU to provide states they can be initialised using the `initialise` method which requires
an initial recharge rate to be set for each HRU. In the following we initialise the states and plot
the initial saturated zone storage deficit, using the chaining of commands.

```{r initialise}
ctch_mdl$initialise()$plot_state("s_sz")
```

The simulation can now be performed and the flow at the gauge extracted with

```{r sim1}
sim1 <- ctch_mdl$sim(swindale_model$output_flux)$get_output()
```

Note that the states of the system are now those at the end of the simulation
for example:

```{r new_states}
ctch_mdl$plot_state("s_sz")
```

Rerunning the simulation with the new initial states will of course produce
different results. Output for the above examples can be plotted against observed discharge for comparison as follows:

```{r sim2}
sim2 <- ctch_mdl$sim(swindale_model$output_flux)$get_output()
out <- merge( merge(swindale_obs,sim1),sim2)
names(out) <- c(names(swindale_obs),'sim_1','sim_2')
plot(out[,c('flow','sim_1','sim_2')], main="Discharge",ylab="m3/s",legend.loc="topright")
```

# Mass balance

It is possible to output the mass balance check for each time step of the
simulation using the `get_mass_errors` method. The returned `xts` object gives 
the volumes in the states at the start and end of the time step along with the
other fluxes as volumes. This can easily be used to plot the errors as shown below.

```{r mass_check}
mb <- ctch_mdl$get_mass_errors()
plot( mb[,6] , main="Mass Error", ylab="[m^3]")
```