| Title: | Algorithms for Helping Build Dynamic TOPMODEL Implementations from Spatial Data |
|---|---|
| Description: | A set of algorithms based on Quinn et al. (1991) <doi:10.1002/hyp.3360050106> for processing river network and digital elevation data to build implementations of Dynamic TOPMODEL, a semi-distributed hydrological model proposed in Beven and Freer (2001) <doi:10.1002/hyp.252>. The 'dynatop' package implements simulation code for Dynamic TOPMODEL based on the output of 'dynatopGIS'. |
| Authors: | Paul Smith [aut, cre] |
| Maintainer: | Paul Smith <[email protected]> |
| License: | GPL-2 |
| Version: | 0.3.0.1010 |
| Built: | 2024-09-04 03:06:02 UTC |
| Source: | https://github.com/waternumbers/dynatopgis |
Converts SpatialLinesDataFrame or SpatialPolygonsDataFrame to the correct format of SpatialPolygonsDataFrame for dynatopGIS.
convert_channel( vect_object, property_names = c(name = "DRN_ID", length = "length", startNode = "startNode", endNode = "endNode", width = "width", slope = "slope"), default_width = 2, default_slope = 0.001 )convert_channel( vect_object, property_names = c(name = "DRN_ID", length = "length", startNode = "startNode", endNode = "endNode", width = "width", slope = "slope"), default_width = 2, default_slope = 0.001 )
vect_object |
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 |
default_width |
the width in m to be used for buffering lines to produce polygons |
default_slope |
the slope in m/m to be used when none is provided |
If the property_names vector contains a width this is used for buffering lines to produce polygons, otherwise the default_width value is used.
channel_file <- system.file("extdata", "SwindaleRiverNetwork.shp", package="dynatopGIS", mustWork = TRUE) vect_lines <- terra::vect(channel_file) property_names <- c(name="identifier",endNode="endNode",startNode="startNode",length="length") chn <- convert_channel(vect_lines,property_names)channel_file <- system.file("extdata", "SwindaleRiverNetwork.shp", package="dynatopGIS", mustWork = TRUE) vect_lines <- terra::vect(channel_file) property_names <- c(name="identifier",endNode="endNode",startNode="startNode",length="length") chn <- convert_channel(vect_lines,property_names)
R6 Class for processing a catchment to make a Dynamic TOPMODEL
R6 Class for processing a catchment to make a Dynamic TOPMODEL
new()
Initialise a project, or reopen an existing project
dynatopGIS$new(projectFolder)
projectFolderfolder for data files
This loads the project data files found in the projectFolder if present. If not the folder is created. The project data files are given by projectFolder/<filename>.<tif,shp>
A new 'dynatopGIS' object
add_catchment()
Add a catchment outline to the 'dynatopGIS' project
dynatopGIS$add_catchment(catchment)
catchmenta SpatRaster object or the path to file containing one which contains a rasterised catchment map.
If not a SpatRaster object the the catchment is read in using the terra package. Finite values in the raster indicate that the area is part of the catchment; with each subcatchment taking a unique finite value. Note that in the later processing it is assumed that outflow from the subcatchments can occur only through the channel network. The resolution and projection of the project is taken from the provided catchment
invisible(self)
add_dem()
Import a dem to the 'dynatopGIS' object
dynatopGIS$add_dem(dem, fill_na = -9999)
dema raster layer object or the path to file containing one which is the DEM
fill_nashould NA values in dem be filled. See details
verboseShould additional progress information be printed
If not a raster the DEM is read in using the terra package. If fill_na is TRUE all NA values other then those that link to the edge of the dem are filled so they can be identified as sinks.
suitable for chaining
add_channel()
Import channel data to the 'dynatopGIS' object
dynatopGIS$add_channel(channel, verbose = FALSE)
channela SpatVect object 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, length, area, startNode, endNode and overlaying it on the DEM. In doing this a variable called id is created (or overwritten) other variables in the data frame are passed through unaltered.
suitable for chaining
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 by the raster 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 tif file.
suitable for chaining
get_layer()
Get a layer of geographical information or a list of layer names
dynatopGIS$get_layer(layer_name = character(0))
layer_namename of the layer give to the layer
a 'raster' layer of the requested information if layer_name is given else a vector of layer names
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
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,
flow_type = c("quinn", "d8")
)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
flow_typeThe type of flow routing to apply see details
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 adaptations made are to ensure that all cells drain only within the subcatchments if provided.
The flow_type can be either - "quinn" where flow is split across all downslope directions or - "d8" where all flow follows the steepest between cell gradient
compute_band()
Computes the computational band of each cell
dynatopGIS$compute_band(type = c("strict"), verbose = FALSE)typetype of banding
verboseprint out additional diagnostic information
Banding is used within the model to define the HRUs and control the order of the flow between them; HRUs can only pass flow to HRUs in a lower numbered band. Currently only a strict ordering of river channels and cells in the DEM is implemented. To compute this the algorithm passes first up the channel network (with outlets being in band 1) then through the cells of the DEM in increasing height.
compute_properties()
Computes statistics e.g. gradient, log(upslope area / gradient) for raster cells
dynatopGIS$compute_properties(min_grad = 1e-04, verbose = FALSE)
min_gradgradient that can be assigned to a pixel if it can't be computed
verboseprint out additional diagnostic information
The algorithm passed through the cells in decreasing height. Min grad is applied to all cells. It is also used for missing gradients in pixels which are partially channel but have no upslope neighbours.
compute_flow_lengths()
Computes flow length for each pixel to the channel
dynatopGIS$compute_flow_lengths(
flow_routing = c("expected", "dominant", "shortest"),
verbose = FALSE
)flow_routingTODO
verboseprint out additional diagnostic information
The algorithm passes through the cells in the DEM in increasing height. Three measures of flow length to the channel are computed. The shortest length (minimum length to channel through any flow path), the dominant length (the length taking the flow direction with the highest fraction for each pixel on the path) and expected flow length (flow length based on sum of downslope flow lengths based on fraction of flow to each cell). By definition cells in the channel that have no land area have a length of NA.
classify()
Create a catchment classification based cutting an existing layer into classes
dynatopGIS$classify(layer_name, base_layer, cuts)
layer_namename of the new layer to create
base_layername of the layer to be cut into classes
cutsvalues on which to cut into classes. These should be numeric and define either the number of bands (single value) or breaks between band (multiple values).
This applies the given cuts to the supplied landscape layer to produce areal groupings of the catchment. Cuts are implement using terra::cut with include.lowest = TRUE. Note that is specifying a vector of cuts values outside the limits will be set to NA.
combine_classes()
Combine any number of classifications based on unique combinations and burns
dynatopGIS$combine_classes(layer_name, pairs, burns = NULL)
layer_namename of the new layer to create
pairsa vector of layer names to combine into new classes through unique combinations. Names should correspond to raster layers in the project directory.
burnsa vector of layer names which are to be burnt on
This applies the given cuts to the supplied landscape layers to produce areal groupings of the catchment. Burns are added directly in the order they are given. Cuts are implement using terra::cut with include.lowest = TRUE. Note that is specifying a vector of cuts values outside the limits will be set to NA.
create_model()
Compute a Dynamic TOPMODEL
dynatopGIS$create_model(
layer_name,
class_layer,
sf_opt = c("cnst", "kin"),
sz_opt = c("exp", "bexp", "cnst", "dexp"),
rain_layer = NULL,
rain_label = character(0),
pet_layer = NULL,
pet_label = character(0),
verbose = FALSE
)layer_namename for the new model and layers
class_layerthe layer defining the topographic classes
sf_optSurface solution to use
sz_opttransmissivity profile to use
rain_layerthe layer defining the rainfall inputs
rain_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 define the HRUs. Flow between HRUs is based on the ordering of the catchment (see the compute_band method). Flow from a HRU can only go to a HRU with a lower band.
Setting the sf_opt and sz_opt options ensures the model is set up with the correct parameters present.
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 weights used to compute an averaged input value for each HRU are computed, otherwise an input table for the models generated with the value "missing" used in place of the series name.
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
get_method()
get the cuts and burns used to classify
dynatopGIS$get_method(layer_name)
layer_namethe name of layer whose classification method is returned
a list with two elements, cuts and burns
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) ## initialise processing ctch <- dynatopGIS$new(file.path(demo_dir,"test")) ## add a catchment outline based on the digital elevation model dem_file <- system.file("extdata", "SwindaleDTM40m.tif", package="dynatopGIS", mustWork = TRUE) dem <- terra::rast(dem_file) dem <- terra::extend(dem,1) catchment_outline <- terra::ifel(is.finite(dem),1,NA) ctch$add_catchment(catchment_outline) ## add digital elevation and channel data ctch$add_dem(dem) channel_file <- system.file("extdata", "SwindaleRiverNetwork.shp", package="dynatopGIS", mustWork = TRUE) sp_lines <- terra::vect(channel_file) property_names <- c(name="identifier",endNode="endNode",startNode="startNode",length="length") chn <- convert_channel(sp_lines,property_names) ctch$add_channel(chn) ## compute properties ctch$sink_fill() ## fill sinks in the catchment and computes dem flow directions ctch$compute_properties() # like topograpihc index and contour length ctch$compute_band() ctch$compute_flow_lengths() ## classify and create a model ctch$classify("atb_20","atb",cuts=20) # classify using the topographic index ctch$get_method("atb_20") ## see the details of the classification ctch$combine_classes("atb_20_band",c("atb_20","band")) ## combine classes ctch$create_model(file.path(demo_dir,"new_model"),"atb_20_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) ## initialise processing ctch <- dynatopGIS$new(file.path(demo_dir,"test")) ## add a catchment outline based on the digital elevation model dem_file <- system.file("extdata", "SwindaleDTM40m.tif", package="dynatopGIS", mustWork = TRUE) dem <- terra::rast(dem_file) dem <- terra::extend(dem,1) catchment_outline <- terra::ifel(is.finite(dem),1,NA) ctch$add_catchment(catchment_outline) ## add digital elevation and channel data ctch$add_dem(dem) channel_file <- system.file("extdata", "SwindaleRiverNetwork.shp", package="dynatopGIS", mustWork = TRUE) sp_lines <- terra::vect(channel_file) property_names <- c(name="identifier",endNode="endNode",startNode="startNode",length="length") chn <- convert_channel(sp_lines,property_names) ctch$add_channel(chn) ## compute properties ctch$sink_fill() ## fill sinks in the catchment and computes dem flow directions ctch$compute_properties() # like topograpihc index and contour length ctch$compute_band() ctch$compute_flow_lengths() ## classify and create a model ctch$classify("atb_20","atb",cuts=20) # classify using the topographic index ctch$get_method("atb_20") ## see the details of the classification ctch$combine_classes("atb_20_band",c("atb_20","band")) ## combine classes ctch$create_model(file.path(demo_dir,"new_model"),"atb_20_band") ## create a model list.files(demo_dir,pattern="new_model*") ## look at the output files for the model ## tidy up unlink(demo_dir)