| 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 building and simulating models. |
| Authors: | Paul Smith [aut, cre] (ORCID: <https://orcid.org/0000-0002-0034-3412>), Peter Metcalfe [aut] |
| Maintainer: | Paul Smith <[email protected]> |
| License: | AGPL-3 |
| Version: | 0.5.0.4 |
| Built: | 2026-10-01 11:31:33 UTC |
| Source: | https://codeberg.org/dynatop/dynatop |
Performs chacks on the values and properties of a channel then checks connnections to outlets if given
check_channel(chn, outlets = NULL, check_connectivity = TRUE)check_channel(chn, outlets = NULL, check_connectivity = TRUE)
chn |
a SpatVect object representing the channel newtork as returned by convert_channel |
outlets |
uid values of the outlets of the channel network |
check_connectivity |
check the connectivity of the channel network |
The processing follows the follwoing sequence: - Checks of the required properties of the channel netowrk for use in dynatopGIS - checks on the connectivity for loops - If outlets provided comparision of identified and known outlets
depth and width are only checked is present
TRUE if completes else failure message
chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth")) check_channel(chn)chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth")) check_channel(chn)
Performs chacks on the values and properties of a water body
check_water_bodies(wb)check_water_bodies(wb)
wb |
a SpatVect object representing the water bodies. |
The processing checks of the required properties of the water body for use in dynatopGIS: - uid: unique character uid - class: numeric - slope: positive numeric - width: positive numeric - depth: positive numeric
TRUE if completes else failure message
wbFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) wb <- terra::vect(wbFile) wb <- terra::buffer(wb,5) wb$uid <- paste0("wb",seq_len(nrow(wb))) wb$slope <- 1e-3 wb$class <- 1 wb$width <- 1 wb$depth <- 0.1 check_water_bodies(wb)wbFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) wb <- terra::vect(wbFile) wb <- terra::buffer(wb,5) wb$uid <- paste0("wb",seq_len(nrow(wb))) wb$slope <- 1e-3 wb$class <- 1 wb$width <- 1 wb$depth <- 0.1 check_water_bodies(wb)
Converts SpatVect of the channel network into SpatVect containing polgons suitable for dynatopGIS
convert_channel( chn, property_names = c(uid = "uid", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width = "width", depth = "depth"), min_slope = 0.001 )convert_channel( chn, property_names = c(uid = "uid", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width = "width", depth = "depth"), min_slope = 0.001 )
chn |
a SpatVect object or a file which can read by terra::vect to create one |
property_names |
a named vector of containing the columns of existing data properties required in the final SpatialPolygonsDataFrame |
min_slope |
the minimum slope |
The processing follows the follwoing sequence: - The 'property_names' input is used to rename the data associated with spatial objects. Remaining data is dropped - Varaibles are converted to the expected type then checked for missing values; in doing this - Missing or lower ‘slope’ values are populated by the min_slope and a warning issued
A SpatVect containing polygons of the channel network, with at least the following properties: uid, width, depth, slope, startNode and endNode.
chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth"))chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth"))
Two functions to generate daily profile properties then apply these to a time series
daily_evap_profile(eMin = 0, eMax = 3, max_day = 183) evap_est(ts, daily_profile, daily_data = NULL, pet_name = "pet")daily_evap_profile(eMin = 0, eMax = 3, max_day = 183) evap_est(ts, daily_profile, daily_data = NULL, pet_name = "pet")
eMin |
Minimum daily PE total (m or mm) |
eMax |
Maximum daily PE total (m or mm) |
max_day |
Year day of maximum PET |
ts |
as vector of POSIXct data/times |
daily_profile |
data.frame containing annual profile properties for PET calculation |
daily_data |
data frame containing daily PET values for disaggregation |
pet_name |
column of daily_data containing the PET values |
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.
The daily_evap_profile function computes a profile properties based on a
simple sinusoid giving year day values for PET, dawn and day length.
The evap_est function takes a time series of steps and computes the aggregated PET
over each time step baed on the profile and optionally daily PET values.
All inputs times are handled as though they are UTC!
daily_evap_profile returns a data.frame of year day, dawn (s since start of day), day length (s) and PET (units as eMin/eMax).
evap_est a data.frame with columns dateTime containing the time stamps of the data and pet containing the PET values.
Beven, K. J. (2012). Rainfall-runoff modelling : the primer. Chichester, UK, Wiley-Blackwell.
Calder, I. R., Harding, R. J. and Rosier, P. T. W. (1983) An objective assessment of soil-moisture deficit models. Journal of Hydrology, 60(1-4) 329-355
## Generating daily and hourly PET data for 1970 ## the values of eMin and eMax may not by not be realistic st_day <- as.POSIXct("1970-01-02 00:00:00",tz='GMT') st_hour <- as.POSIXct("1970-01-01 01:00:00",tz='GMT') fn <- as.POSIXct("1971-01-01 00:00:00",tz='GMT') daily_prof <- daily_evap_profile(0,1) daily_ts <- seq(st_day,fn,by=24*60*60) daily_pet <- evap_est(daily_ts,daily_prof) hourly_ts <- seq(st_hour,fn,by=1*60*60) hourly_pet <- evap_est(hourly_ts,daily_prof) ## the totals should be the same... stopifnot(all.equal(sum(hourly_pet$pet), sum(daily_pet$pet)))## Generating daily and hourly PET data for 1970 ## the values of eMin and eMax may not by not be realistic st_day <- as.POSIXct("1970-01-02 00:00:00",tz='GMT') st_hour <- as.POSIXct("1970-01-01 01:00:00",tz='GMT') fn <- as.POSIXct("1971-01-01 00:00:00",tz='GMT') daily_prof <- daily_evap_profile(0,1) daily_ts <- seq(st_day,fn,by=24*60*60) daily_pet <- evap_est(daily_ts,daily_prof) hourly_ts <- seq(st_hour,fn,by=1*60*60) hourly_pet <- evap_est(hourly_ts,daily_prof) ## the totals should be the same... stopifnot(all.equal(sum(hourly_pet$pet), sum(daily_pet$pet)))
R6 Class for Dynamic TOPMODEL
dynatop$new()Creates a dynatop class object from the a list based model description as generated by dynatopGIS.
dynatop$new(geo, param, delta = 1e-13)
geoa data.frame describing the HRUs geography or path to a gpkg file that can read to create one
parama data.frame describing the parameters for each HRU class or path to a jsonfile that can read to create one
deltaerror allowed in checking redistribution sums
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
dynatop$initialise()Initialises a dynatop object in the simplest way possible.
dynatop$initialise(initial_states = NULL)
initial_statesa date frame of states
is initial_states are present these are used for the HRUs given.
Otherwise a steady state initialisation is performed.
invisible(self) suitable for chaining
dynatop$sim()Simulate the dynatop model
dynatop$sim( obs_data, output_defn, keep_states = NULL, output_dir = NULL, sub_step = NULL, n_thread = 1 )
obs_datadata.frame of observations
output_defndefinition of the outputs from the model
keep_statesa vector of POSIXct objects (e.g. from xts) giving the time stamp at which the states should be kept
output_dirthe directory to write output to
sub_stepsimulation timestep in seconds, default value of NULL results in data time step
n_threadnumber of parallel threads to use in execution
Saving the states at every timestep and keeping the mass balance can generate very large data sets!!
invisible(self) for chaining
dynatop$get_states()Get the current states
dynatop$get_states()
data.frame of current statea
dynatop$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. mdlFile <- system.file("extdata","mdl","test_mdl.gpkg",package="dynatop",mustWork=TRUE) paramFile <- system.file("extdata","mdl","test_mdl_param.json",package="dynatop",mustWork=TRUE) obsFile <- system.file("extdata","mdl","SwindaleObs.csv",package="dynatop",mustWork=TRUE) ## create and initialise model dyn <- dynatop$new(mdlFile,paramFile) dyn$initialise() ## this initialisation isn't good - see the vignette ## define an output definition defn <- data.frame(name="outlet",hru=1,scale=1,flux="q_sf") ## read in observed obs <- read_obs(obsFile) ## simulate out <- dyn$sim(obs,defn)## the vignettes contains further details of the method calls. mdlFile <- system.file("extdata","mdl","test_mdl.gpkg",package="dynatop",mustWork=TRUE) paramFile <- system.file("extdata","mdl","test_mdl_param.json",package="dynatop",mustWork=TRUE) obsFile <- system.file("extdata","mdl","SwindaleObs.csv",package="dynatop",mustWork=TRUE) ## create and initialise model dyn <- dynatop$new(mdlFile,paramFile) dyn$initialise() ## this initialisation isn't good - see the vignette ## define an output definition defn <- data.frame(name="outlet",hru=1,scale=1,flux="q_sf") ## read in observed obs <- read_obs(obsFile) ## simulate out <- dyn$sim(obs,defn)
R6 Class for processing a catchment to make a dynatop TOPMODEL
dynatopGIS$new()Initialise a project, or reopen an existing project
dynatopGIS$new(projectDir)
projectDiris the directroy of project data
If the directory exists the existing tif files are loaded as the project.
A new 'dynatopGIS' object
dynatopGIS$add_catchment()Add a catchment outline and dem to the 'dynatopGIS' project
dynatopGIS$add_catchment(catchment, dem, fill_na = -9999)
catchmenta SpatVect object or the path to file containing one which contains the catchment outline
dema SpatRast object or the path to file containing one which is the DEM
fill_navalue used to replace missing dem vales within the catchment
If the inputs are not spatial objects they are read in using the terra package.
The resolution and projection of the project is taken from the provided dem. All NA values in the dem within the catchment are replaced by the fill_na values so they can be identified as sinks.
invisible(self)
dynatopGIS$add_channel()Import channel data to the 'dynatopGIS' object
dynatopGIS$add_channel(channel, water_bodies = NULL, verbose = FALSE)
channela SpatVect object of lines or file path that can be loaded as one containing the channel information
water_bodiesa SpatVect object of polygons or file path that can be loaded as one containing the channel information
verboseShould additional progress information be printed
Takes the representation of the channel network as a SpatVect with properties name, startNode, endNode, slope, depth and width. This is rasterised. The polygons in water_bodies are then burnt in. it is presumed that each polygon has only one outlet.
suitable for chaining
dynatopGIS$add_layer()Add a layer of geographical information
dynatopGIS$add_layer(layer, layer_name = names(layer))
layerthe raster layer to add (see details)
layer_namename to give to the layer
The layer should either be a raster layer or a file that can be read into a SpatRast by the terra package. The projection, resolution and extent are checked against the existing project data. Only layer names not already in use (or reserved) are allowed. If successful the layer is added to the project.
suitable for chaining
dynatopGIS$get_layer()Get a layer of geographical information or a list of layer names
dynatopGIS$get_layer(layer_name = NULL)
layer_namename of the layer give to the layer or NULL to return list of layers
a SpatRaster containing the requested layer or a vector of layer names
dynatopGIS$plot_layer()Plot a layer
dynatopGIS$plot_layer(layer_name, add_channel = TRUE)
layer_namethe name of layer to plot
add_channelshould the channel be added to the plot
a plot
dynatopGIS$sink_fill()The sink filling algorithm of Planchona and Darboux (2001)
dynatopGIS$sink_fill( min_grad = 1e-04, max_it = 1e+06, verbose = FALSE, hot_start = FALSE )
min_gradMinimum gradient between cell centres
max_itmaximum number of replacement cycles
verboseprint out additional diagnostic information
hot_startstart from filled_dem if it exists
The algorithm implemented is based on that described in Planchona and Darboux, "A fast, simple and versatile algorithm to fill the depressions in digital elevation models" Catena 46 (2001). A pdf can be found at (<https://horizon.documentation.ird.fr/exl-doc/pleins_textes/pleins_textes_7/sous_copyright/010031925.pdf>).
The flow_type can be either MFD algorithm in TODO
dynatopGIS$order()Computes the order (banding) to evaluate cells in
dynatopGIS$order(verbose = FALSE)
verboseprint out additional diagnostic information
The algorithm passed through the cells to compute the band starting at the river outlets hru numbers are assigned in approxamtly band order
dynatopGIS$combine_layers()Combine any number of layers into a unique classification
dynatopGIS$combine_layers(class_name, layers)
class_namename of the new layer to create
layersa vector of layer names to combine into new classes through unique combinations. Names should correspond to raster layers in the project directory.
This combines the layers into unique classess and stores a raster of outcomes along with a json file of the classifications
dynatopGIS$create_model()Create a Dynamic TOPMODEL
dynatopGIS$create_model( model_name, class_layer = NULL, precip_layer = NULL, precip_label = "precip", pet_layer = NULL, pet_label = "pet", verbose = FALSE )
model_namebase of the name for the new model files
class_layerthe layer defining the landscape classes
precip_layerthe layer defining the rainfall inputs
precip_labelPrepended to rain_layer values to give rainfall series name
pet_layerthe layer defining the pet inputs
pet_labelPrepended to pet_layer values to give pet series name
verboseprint more details of progress
The class_layer is used to assign classes to the HRUs. Flow between HRUs are based on the ordering of the catchment (see the order method) with low from a HRU can only go to a HRU with a lower band.
The rain_layer (pet_layer) can contain the numeric id values of different rainfall (pet) series. If the value of rain_layer (pet_layer) is not NULL the value for the precipitation and PET inputs for each HRU are taken as e.g. paste0(precip_label,precip) otherwise the value of e.g. precip_label is used. The model_name can contain the path to the files.
Two files model_name.gpkg which contains the spatial description of the model and model_name_param.json an example parameter file for the model.
dynatopGIS$get_version()get the version number
dynatopGIS$get_version()
the version number indicates the version of the algorithms within the object
a numeric version number
dynatopGIS$clone()The objects of this class are cloneable with this method.
dynatopGIS$clone(deep = FALSE)
deepWhether to make a deep clone.
## The vignettes contains more examples of the method calls. ## create temporary directory for output demo_dir <- tempfile("dygis") dir.create(demo_dir) ## locate the dem, channel and catchment outline files chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) demFile <- system.file("extdata", "gis", "SwindaleDTM.tif", package="dynatop", mustWork = TRUE) ctchFile <- system.file("extdata", "gis", "SwindaleBoundary.gpkg", package="dynatop", mustWork = TRUE) ## Preprocess channel data chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth")) chn <- simplify_channel(chn) ## create a project dy <- dynatopGIS$new(demo_dir) ## add dem and catchment data dy$add_catchment(ctchFile, demFile) ## add the channel dy$add_channel(chn) ## fill the sinks dy$sink_fill() ## compute the order dy$order() ## add a new layer tmp <- dy$get_layer("filled_dem") tmp <- terra::classify(tmp,5) names(tmp) <- "alt_band" dy$add_layer(tmp) ## classify the layer dy$combine_layers("class_alt_band","alt_band") ## build a model dy$create_model(file.path(demo_dir,"new_model"),"class_alt_band") ## create a model list.files(demo_dir,pattern="new_model*") ## look at the output files for the model ## tidy up unlink(demo_dir)## The vignettes contains more examples of the method calls. ## create temporary directory for output demo_dir <- tempfile("dygis") dir.create(demo_dir) ## locate the dem, channel and catchment outline files chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) demFile <- system.file("extdata", "gis", "SwindaleDTM.tif", package="dynatop", mustWork = TRUE) ctchFile <- system.file("extdata", "gis", "SwindaleBoundary.gpkg", package="dynatop", mustWork = TRUE) ## Preprocess channel data chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth")) chn <- simplify_channel(chn) ## create a project dy <- dynatopGIS$new(demo_dir) ## add dem and catchment data dy$add_catchment(ctchFile, demFile) ## add the channel dy$add_channel(chn) ## fill the sinks dy$sink_fill() ## compute the order dy$order() ## add a new layer tmp <- dy$get_layer("filled_dem") tmp <- terra::classify(tmp,5) names(tmp) <- "alt_band" dy$add_layer(tmp) ## classify the layer dy$combine_layers("class_alt_band","alt_band") ## build a model dy$create_model(file.path(demo_dir,"new_model"),"class_alt_band") ## create a model list.files(demo_dir,pattern="new_model*") ## look at the output files for the model ## tidy up unlink(demo_dir)
Method for reading both model maps and parameter files.
read_param(file_name) write_param(param, file_name) read_geo(file_name, verbose = FALSE) write_geo(geo, file_name) read_obs(file_name, time_zone = "UTC")read_param(file_name) write_param(param, file_name) read_geo(file_name, verbose = FALSE) write_geo(geo, file_name) read_obs(file_name, time_zone = "UTC")
file_name |
name of file to read or write into |
param |
parameter file containing classes |
verbose |
provide details of progress |
geo |
geography object describing the catchment with locations, hru numbers,connectivity, classes etc. |
time_zone |
time zone of data |
WARNING: There is no checking of format
Fill NA values in a SpatRast object
fill_na(rst, mask = NULL, layer_names = names(rst))fill_na(rst, mask = NULL, layer_names = names(rst))
rst |
SpatRast (or path to file) to infill |
mask |
SpatRast (or path to file) defining non-NA cells (omit for all cells) |
layer_names |
names of layers to infill (omit for all layers) |
rst and mask should have the same extent and resolution. The NA values in rst
are infilled using the most common neighbouring value through multiple passes of focal.
Simplify a channel object by removing bifurcations and merging spatial objects and removing bifurcations
simplify_channel(chn, simplify_length = 100, verbose = FALSE)simplify_channel(chn, simplify_length = 100, verbose = FALSE)
chn |
the channel object |
simplify_length |
minimum length of channel sections to remove |
verbose |
print a higher level of output |
Analysis is carriend out in the folowing order - Bifurcations are removed by leaving only the longest link from the node - Removing nodes nodes with a single input and output by merging adjacent river sections to get length > simplify_length
chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth")) tchn <- simplify_channel(chn) nrow(chn) nrow(tchn)chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth")) tchn <- simplify_channel(chn) nrow(chn) nrow(tchn)
Trim channel objects to those draining to named outlets
trim_channel(chn, outlets, removed = FALSE)trim_channel(chn, outlets, removed = FALSE)
chn |
the channel object |
outlets |
names of the outlets to consider |
removed |
should the removed objects be returned |
Starting with the named outlet objects the channel network is passed up using the connectivity given by the start and end node, only those flaged as connected are returned.
chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth")) tchn <- trim_channel(chn,chn$uid[1]) nrow(chn) nrow(tchn)chnFile <- system.file("extdata", "gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE) chn <- terra::vect(chnFile) chn$slope <- 1e-3 chn$class <- 1 chn$width <- 1 chn$depth <- 0.1 chn <- convert_channel(chn, property_names=c(uid = "identifier", startNode = "startNode", endNode = "endNode", slope = "slope", class = "class", width="width", depth="depth")) tchn <- trim_channel(chn,chn$uid[1]) nrow(chn) nrow(tchn)