R6 Class for processing a catchment to make a Dynamic TOPMODEL

Source: R/dygis.R, R/dynatopGIS-package.R

This package contains the code for setting up a dynamic TOPMODEL implementation

Methods

Public methods

Method new()

Initialise a project, or reopen an existing project

Usage

dynatopGIS$new(projectFolder)

Arguments

projectFolder

folder for data files

Details

This loads the project data files found at projFile if present. If not the value is stored for later use. The project data files are given by projfile,<tif,shp>

Returns

A new `dynatopGIS` object

Method add_dem()

Import a dem to the `dynatopGIS` object

Usage

dynatopGIS$add_dem(dem, fill_na = TRUE)

Arguments

dem

a raster layer object or the path to file containing one which is the DEM

fill_na

should NA values in dem be filled. See details

verbose

Should additional progress information be printed

Details

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.

Returns

suitable for chaining

Method add_channel()

Import channel data to the `dynatopGIS` object

Usage

dynatopGIS$add_channel(channel)

Arguments

channel

a SpatialLinesDataFrame, SpatialPolygonsDataFrame or file path containing the channel information

Details

Takes the representation of the channel network as a SpatialPolygonsDataFrame 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.

Returns

suitable for chaining

Method add_layer()

Add a layer of geographical information

Usage

dynatopGIS$add_layer(layer, layer_name = names(layer))

Arguments

layer

the raster layer to add (see details)

layer_name

name to give to the layer

Details

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.

Returns

suitable for chaining

Method get_layer()

Get a layer of geographical information or a list of layer names

Usage

dynatopGIS$get_layer(layer_name = character(0))

Arguments

layer_name

name of the layer give to the layer

Returns

a `raster` layer of the requested information if layer_name is given else a vector of layer names

Method plot_layer()

Plot a layer

Usage

dynatopGIS$plot_layer(layer_name, add_channel = TRUE)

Arguments

layer_name

the name of layer to plot

add_channel

should the channel be added to the plot

Returns

a plot

Method sink_fill()

The sink filling algorithm of Planchona and Darboux (2001)

Usage

dynatopGIS$sink_fill(
  min_grad = 1e-04,
  max_it = 1e+06,
  verbose = FALSE,
  hot_start = FALSE
)

Arguments

min_grad

Minimum gradient between cell centres

max_it

maximum number of replacement cycles

verbose

print out additional diagnostic information

hot_start

start from filled_dem if it exists

Details

The algorithm implemented is 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>).

Method compute_properties()

Computes statistics e.g. gradient, log(upslope area / gradient) for raster cells

Usage

dynatopGIS$compute_properties(min_grad = 1e-04, verbose = FALSE)

Arguments

min_grad

gradient that can be assigned to a pixel if it can't be computed

verbose

print out additional diagnostic information

Details

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.

Method compute_flow_lengths()

Computes flow length for each pixel to the channel

Usage

dynatopGIS$compute_flow_lengths(verbose = FALSE)

Arguments

verbose

print out additional diagnostic information

Details

The algorithm passed through the cells in increasing height. For 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) and band (strict sequence to ensure that all contributing cell have a higher band value). By definition cells in the channel that have no land area have a length (or band) of NA.

Method classify()

Create a catchment classification based cutting an existing layer into classes

Usage

dynatopGIS$classify(layer_name, base_layer, cuts)

Arguments

layer_name

name of the new layer to create

base_layer

name of the layer to be cut into classes

cuts

values 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).

Details

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.

Method combine_classes()

Combine any number of classifications based on unique combinations and burns

Usage

dynatopGIS$combine_classes(layer_name, pairs, burns = NULL)

Arguments

layer_name

name of the new layer to create

pairs

a vector of layer names to combine into new classes through unique combinations. Names should correspond to raster layers in the project directory.

burns

a vector of layer names which are to be burnt on

Details

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.

Method create_model()

Compute a Dynamic TOPMODEL

Usage

dynatopGIS$create_model(
  layer_name,
  class_layer,
  dist_layer,
  sf_opt = c("cnst", "kin"),
  sz_opt = c("exp", "bexp", "cnst", "dexp"),
  dist_delta = 0,
  rain_layer = NULL,
  rain_label = character(0),
  pet_layer = NULL,
  pet_label = character(0),
  verbose = FALSE
)

Arguments

layer_name

name for the new model and layers

class_layer

the layer defining the topographic classes

dist_layer

the layer defining the distances to the channel

sf_opt

Surface solution to use

sz_opt

transmissivity transmissivity profile to use

dist_delta

TODO

rain_layer

the layer defining the rainfall inputs

rain_label

Prepended to rain_layer values to give rainfall series name

pet_layer

the layer defining the pet inputs

pet_label

Prepended to pet_layer values to give pet series name

verbose

print more details of progress

Details

The class_layer is used to define the HRUs. Flow between HRUs is based on the distance to a channel. For each HRU the shortest distance to a channel is computed. Flow from a HRU can only go to a HRU with a lower shortest distance to the channel. Flow from a HRU can occur from any raster cell within the HRU whose distance to the channel is within dist_delta of the shortest distance within the HRU. 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.

Method get_version()

get the version number

Usage

dynatopGIS$get_version()

Details

the version number indicates the version of the algorithms within the object

Returns

a numeric version number

Method get_method()

get the cuts and burns used to classify

Usage

dynatopGIS$get_method(layer_name)

Arguments

layer_name

the name of layer whose classification method is returned

Returns

a list with two elements, cuts and burns

Method clone()

The objects of this class are cloneable with this method.

Usage

dynatopGIS$clone(deep = FALSE)

Arguments

deep

Whether to make a deep clone.

Examples

## The vignettes contains more examples of the method calls.

## create temport 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/Rtmp1UiXow/dygis1e25478b6141/test

## add digital elevation and channel data
dem_file <- system.file("extdata", "SwindaleDTM40m.tif", package="dynatopGIS", mustWork = TRUE)
dem <- terra::rast(dem_file)
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
ctch$add_channel(chn)

## compute properties 
ctch$sink_fill() ## fill sinks in the catchment
# \donttest{
ctch$compute_properties() # like topograpihc index and contour length
ctch$compute_flow_lengths()
#> Warning: NAs introduced by coercion to integer range
#> Warning: NAs introduced by coercion to integer range
# }
## 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.3714 10.1066 10.8418 11.5771 12.3123 13.0476 13.7828 14.5181
#> [10] 15.2533 15.9885 16.7238 17.4590 18.1943 18.9295 19.6647 20.4000 21.1352
#> [19] 21.8705 22.6057 23.3410
#> 
ctch$combine_classes("atb_20_band",c("atb_20","band")) ## combine classes
ctch$create_model(file.path(demo_dir,"new_model"),"atb_20_band","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)