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)add_catchment()Add a catchment outline to the `dynatopGIS` project
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
add_dem()Import a dem to the `dynatopGIS` object
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
add_channel()Import channel data to the `dynatopGIS` object
channela SpatVect object or file path that can be loaded as one containing the channel information
verboseShould additional progress information be printed
add_layer()Add a layer of geographical information
dynatopGIS$add_layer(layer, layer_name = names(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.
get_layer()Get a layer of geographical information or a list of layer names
dynatopGIS$get_layer(layer_name = character(0))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)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
compute_flow_lengths()Computes flow length for each pixel to the channel
dynatopGIS$compute_flow_lengths(
flow_routing = c("expected", "dominant", "shortest"),
verbose = FALSE
)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
combine_classes()Combine any number of classifications based on unique combinations and burns
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
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.
## 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"))
#> Creating new folder
#> Starting new project at /tmp/RtmpeAcr7b/dygis1b5748337241/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)
#> Warning: [vect] Z coordinates ignored
property_names <- c(name="identifier",endNode="endNode",startNode="startNode",length="length")
chn <- convert_channel(sp_lines,property_names)
#> Warning: Modifying to spatial polygons using default width
#> Warning: Adding default slope
ctch$add_channel(chn)
## compute properties
ctch$sink_fill() ## fill sinks in the catchment and computes dem flow directions
# \donttest{
ctch$compute_properties() # like topograpihc index and contour length
ctch$compute_band()
ctch$compute_flow_lengths()
# }
## classify and create a model
# \donttest{
ctch$classify("atb_20","atb",cuts=20) # classify using the topographic index
ctch$get_method("atb_20") ## see the details of the classification
#> $type
#> [1] "classification"
#>
#> $layer
#> [1] "atb"
#>
#> $cuts
#> [1] 8.6361 9.2646 9.8931 10.5217 11.1502 11.7787 12.4072 13.0357 13.6642
#> [10] 14.2927 14.9212 15.5498 16.1783 16.8068 17.4353 18.0638 18.6923 19.3208
#> [19] 19.9493 20.5779 21.2064
#>
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
#> [1] "new_model.rds" "new_model.tif"
# }
## tidy up
unlink(demo_dir)