,
Peter Metcalfe [aut]| Title: | An Implementation of Dynamic TOPMODEL Hydrological Model in R |
|---|---|
| Description: | An R implementation and enhancement of the Dynamic TOPMODEL semi-distributed hydrological model originally proposed by Beven and Freer (2001) <doi:10.1002/hyp.252>. The 'dynatop' package implements code for simulating models which can be created using the 'dynatopGIS' package. |
| Authors: | Paul Smith [aut, cre] ,
Peter Metcalfe [aut] |
| Maintainer: | Paul Smith <[email protected]> |
| License: | GPL-2 |
| Version: | 0.3.0.1010 |
| Built: | 2024-12-03 06:13:16 UTC |
| Source: | https://github.com/waternumbers/dynatop |
R6 Class for Dynamic TOPMODEL
R6 Class for Dynamic TOPMODEL
new()
Creates a dynatop class object from the a list based model description as generated by dynatopGIS.
dynatop$new(model, map = NULL, use_states = FALSE, delta = 1e-13)
modela dynamic TOPMODEL list object
mapfile name of the map layers for the model
use_stateslogical if states should be imported
deltaerror term in checking redistribution sums
drop_maplogical if the map should be dropped
This function makes some basic consistency checks on a list representing a dynamic TOPMODEL model. The checks performed and basic 'sanity' checks. They do not check for the logic of the parameter values nor the consistency of states and parameters. Sums of the redistribution matrices are checked to be in the range 1 +/- delta.
invisible(self) suitable for chaining
add_data()
Adds observed data to a dynatop object
dynatop$add_data(obs_data)
obs_dataan xts object of observed data
This function makes some basic consistency checks on the observations to ensure they have uniform timestep and all required series are present.
invisible(self) suitable for chaining
clear_data()
Clears all forcing and simulation data except current states
dynatop$clear_data()
invisible(self) suitable for chaining
initialise()
Initialises a dynatop object in the simplest way possible.
dynatop$initialise( vtol = sqrt(.Machine$double.eps), ftol = sqrt(.Machine$double.eps), max_it = 1000 )
vtoltolerance for the solution for the saturated zone storage (as volume)
ftoltolerance for the solution of the saturated zone storage (as difference of function from 0)
max_itmaximum number of iterations to use in the solution of the saturated zone
invisible(self) suitable for chaining
sim()
Simulate the hillslope and channel components of a dynatop object
dynatop$sim( output_defn, keep_states = NULL, sub_step = NULL, vtol = 0.001, ftol = sqrt(.Machine$double.eps), max_it = 1000 )
output_defna description of the output series
keep_statesa vector of POSIXct objects (e.g. from xts) giving the time stamp at which the states should be kept
sub_stepsimulation timestep in seconds, default value of NULL results in data time step
vtoltolerance on width of bounds in the numeric search for surface and saturated zone solutions (as volume)
ftol- not currently used
max_itmaximum number of iterations to use in the solution of the saturated zone
Saving the states at every timestep and keeping the mass balance can generate very large data sets!!
invisible(self) for chaining
get_output()
Return channel inflow as an xts series or list of xts series
dynatop$get_output(name = colnames(private$time_series$output))
nameone or more output series to return
plot_output()
Plot the channel inflow
dynatop$plot_output(name = colnames(private$time_series$output))
nameof series to plot
get_obs_data()
Get the observed data
dynatop$get_obs_data()
get_model()
Return the model
dynatop$get_model()
get_mass_errors()
Return the model
dynatop$get_mass_errors()
get_states()
Return states
dynatop$get_states(record = FALSE)
recordlogical TRUE if the record should be returned. Otherwise the current states returned
plot_state()
Plot a current state of the system
dynatop$plot_state(state = c("s_sf", "s_rz", "s_uz", "s_sz"))statethe name of the state to be plotted
clone()
The objects of this class are cloneable with this method.
dynatop$clone(deep = FALSE)
deepWhether to make a deep clone.
## the vignettes contains further details of the method calls. data("Swindale") ## example data mdl <- Swindale$model mdl$map <- system.file("extdata","Swindale.tif",package="dynatop",mustWork=TRUE) ctch_mdl <- dynatop$new(mdl$hru,map=mdl$map) ## create with model ctch_mdl$add_data(Swindale$obs) ## add observations ctch_mdl$initialise() ## initialise model ctch_mdl$sim(Swindale$model$output_flux) ## simulate model## the vignettes contains further details of the method calls. data("Swindale") ## example data mdl <- Swindale$model mdl$map <- system.file("extdata","Swindale.tif",package="dynatop",mustWork=TRUE) ctch_mdl <- dynatop$new(mdl$hru,map=mdl$map) ## create with model ctch_mdl$add_data(Swindale$obs) ## add observations ctch_mdl$initialise() ## initialise model ctch_mdl$sim(Swindale$model$output_flux) ## simulate model
Generate series of potential evapotranspiration
evap_est(ts, eMin = 0, eMax = 0)evap_est(ts, eMin = 0, eMax = 0)
ts |
as vector of POSIXct data/times |
eMin |
Minimum daily PE total (m or mm) |
eMax |
Maximum daily PE total (m or mm) |
Dynamic TOPMODEL requires a time series of potential evapotranspiration in order to calculate and remove actual evapotranspiration from the root zone during a run. Many sophisticated physical models have been developed for estimating potential and actual evapotranspiration, including the Priestly-Taylor (Priestley and Taylor, 1972) and Penman-Monteith (Montieth, 1965) methods. These, however, require detailed meteorological data such as radiation input and relative humidities that are, in general, difficult to obtain. Calder (1983) demonstrated that a simple approximation using a sinusoidal variation in potential evapotranspiration to be a good approximation to more complex schemes.
If the insolation is also taken to vary sinusoidally through the daylight hours then, ignoring diurnal meteorological variations, the potential evapotranspiration during daylight hours for each year day number can be calculated (for the catchment's latitude). Integration over the daylight hours allows the daily maximum to be calculated and thus a sub-daily series generated.
Time series (xts) of potential evapotranspiration totals for the time steps given in same units as eMin and eMax
Beven, K. J. (2012). Rainfall-runoff modelling : the primer. Chichester, UK, Wiley-Blackwell.
Calder, I. R. (1986). A stochastic model of rainfall interception. Journal of Hydrology, 89(1), 65-71.
## Generating daily PET data for 1970 ## the values of eMin and eMax may not by not be realistic st <- as.POSIXct("1970-01-02 00:00:00",tz='GMT') fn <- as.POSIXct("1971-01-01 00:00:00",tz='GMT') daily_ts <- seq(st,fn,by=24*60*60) dpet <- evap_est(daily_ts,0,1) ## create hourly data for the same period st <- as.POSIXct("1970-01-01 01:00:00",tz='GMT') fn <- as.POSIXct("1971-01-01 00:00:00",tz='GMT') hour_ts <- seq(st,fn,by=1*60*60) hpet <- evap_est(hour_ts,0,1) ## the totals should eb the same... stopifnot(all.equal(sum(hpet), sum(dpet)))## Generating daily PET data for 1970 ## the values of eMin and eMax may not by not be realistic st <- as.POSIXct("1970-01-02 00:00:00",tz='GMT') fn <- as.POSIXct("1971-01-01 00:00:00",tz='GMT') daily_ts <- seq(st,fn,by=24*60*60) dpet <- evap_est(daily_ts,0,1) ## create hourly data for the same period st <- as.POSIXct("1970-01-01 01:00:00",tz='GMT') fn <- as.POSIXct("1971-01-01 00:00:00",tz='GMT') hour_ts <- seq(st,fn,by=1*60*60) hpet <- evap_est(hour_ts,0,1) ## the totals should eb the same... stopifnot(all.equal(sum(hpet), sum(dpet)))
Takes an xts time series object and resamples then to a new time step.
resample_xts(obs, dt, is.rate = FALSE)resample_xts(obs, dt, is.rate = FALSE)
obs |
A times series (xts) object with a POSIXct index. |
dt |
New time interval in seconds |
is.rate |
If TRUE then these are rates i.e m/h. Otherwise they are absolute values accumulated within the preceding time interval. Values are scaled before returning so resampling is conservative. |
Time series of observation data are often of different temporal resolutions, however the input to most hydrological models, as is the case with the Dynamic TOPMODEL, requires those data at the same interval. This provides a method to resample a collection of such data to a single interval.
Because of the methods used the results:
- are not accurate when the input data does not have a constant timestep. The code issues a warning and proceeds assuming the data are equally spaced with the modal timestep. - do not guarantee the requested time step but returns a series with the timestep computed from an integer rounding the ratio of the current and requested time step.
An xts object with the new timestep
# Resample Swindale Rainfall to hourly intervals require(dynatop) data("Swindale") obs <- Swindale$obs cobs <- resample_xts(obs, dt=60*60) # hourly data dobs <- resample_xts(cobs,dt=15*60) # back to 15 minute data cdobs <- resample_xts(dobs,dt=60*60) # back to hourly data - checks time stamp conversion obs <- obs[zoo::index(obs)<=max(zoo::index(cobs)),] # check totals stopifnot( all.equal(sum(obs),sum(cobs)) ) stopifnot( all.equal(sum(obs),sum(dobs)) ) stopifnot( all.equal(cobs,cdobs) )# Resample Swindale Rainfall to hourly intervals require(dynatop) data("Swindale") obs <- Swindale$obs cobs <- resample_xts(obs, dt=60*60) # hourly data dobs <- resample_xts(cobs,dt=15*60) # back to 15 minute data cdobs <- resample_xts(dobs,dt=60*60) # back to hourly data - checks time stamp conversion obs <- obs[zoo::index(obs)<=max(zoo::index(cobs)),] # check totals stopifnot( all.equal(sum(obs),sum(cobs)) ) stopifnot( all.equal(sum(obs),sum(dobs)) ) stopifnot( all.equal(cobs,cdobs) )
This data set contains a processed model and observation data for Swindale.
data(Swindale)data(Swindale)
An object of class list of length 2.
require(dynatop) data(Swindale) # Show it # plot(obs)require(dynatop) data(Swindale) # Show it # plot(obs)