Building a Dynamic TOPMODEL
The purpose of this vignette is to provide an outline of the steps
needed to build a dynamic TOPMODEL implementation using the dynatopGIS
package.
Implementation notes
The dynatopGIS package implements a structured, object orientated,
data flow. The steps outlined below create a dynatopGIS
catchment object to which actions are then applied to generate a
model.
The dynatopGIS package is written using the object
orientated framework provided by the R6 package. This means
that some aspects of working with the objects may appear idiosyncratic
for some R users. In using the package as outlined in this vignette
these problems are largely obscured, except for the call structure.
However, before adapting the code, or doing more complex analysis users
should read about R6 class objects (e.g. in the
R6 package vignettes or in the Advanced R book). One
particular gotcha is when copying an object. Using
my_new_object <- my_object
creates a pointer, that is altering my_new_object also
alters my_object. To create a new independent copy of
my_object use
my_new_object <- my_object$clone()
Getting started
The dynatopGIS packages works through a number of steps to generate a
Dynamic TOPMODEL object suitable for use in with the
dynatop package. Each step generates one or more layers
which are saved as raster or shape files into the projects working
directory (which is not necessarily the R working directory). A record
of these layers is kept in the json format meta data file.
This vignette demonstrates the use of the dynatopGIS
package using data from the Swindale catchment in the UK.
To start first load the library
For this vignette we will store the data into a temporary
directory
demo_dir <- tempfile("dygis")
dir.create(demo_dir)
and initialise the analysis by creating a new object specifying the
location of the meta data file, which will be created if it doesn’t
exist.
ctch <- dynatopGIS$new(file.path(demo_dir,"meta.json"))
#> Warning in initialize(...): Creating meta file
#> at/tmp/RtmpABiYZg/dygis2aa62ba12b33/meta.json
#> Warning in private$check_meta(verbose): No checks on the meta are currently
#> performed
Adding catchment data
The basis of the analysis is a rasterised Digital Elevation Model
(DEM) of the catchment and a vectorised representation of the river
network with attributes. Currently these can be in any format supported
by the raster and sp libraries
respectively.
However, within the calculations used for sink filling, flow routing
and topographic index calculations the raster DEM is presumed to be
projected so that is has square cells such that the difference between
the cell centres (in meters) does not alter.
For Swindale the suitable DEM and channel files can be found
using:
dem_file <- system.file("extdata", "SwindaleDTM40m.tif", package="dynatopGIS", mustWork = TRUE)
channel_file <- system.file("extdata", "SwindaleRiverNetwork.shp", package="dynatopGIS", mustWork = TRUE)
Either the DEM or channel files can be added to the project first. In
this case we add the DEM with
dem <- terra::rast(dem_file)
ctch$add_dem(dem)
Note: An additional row and column of NA is
added to each edge of the DEM. All raster layers in the project have the
same projection, resolution and extent of the extended DEM.
Adding river channel data is more complex. The
add_channel method requires a
SpatialLinesDataFrame or
SpatialPolygonsDataFrame as generated by the
sp package. Each entry in the data.frame is treated as a
length of river channel which requires the following properties
- endNode - a label for the downstream end of the river
length
- startNode - a label for the upstream end of the river
length
- length - the length in meters
Additional properties are currently kept but ignored with two
exceptions: - id - this is copied to original_id
with a warning since id is used internally - width
- if the channel is specified with a line sections then the
width property is used to buffer the lines to create channel
polygons.
Since it is possible that these properties are present in a data file
under different names the add_channel method allows for
renaming. To illustrate this let us examine the river network for
Swindale
sp_lines <- terra::vect(channel_file)
#> Warning: [vect] Z coordinates ignored
head(sp_lines)
#> class : SpatVector
#> geometry : lines
#> dimensions : 6, 11 (geometries, attributes)
#> extent : 349884.5, 351675.7, 508614.5, 513074.7 (xmin, xmax, ymin, ymax)
#> coord. ref. : OSGB36 / British National Grid (EPSG:27700)
#> names : name1 identifier startNode endNode
#> type : <chr> <chr> <chr> <chr>
#> values : Hawthorn Gill 16D0AC09-E0B6-~ 730F01D2-0F4F-~ 72B1A40B-E106-~
#> Little Mosedal~ 643017BB-01BC-~ 9B06188C-F4C7-~ A19E0D3A-2FA1-~
#> Mosedale Beck D2C2703E-A32F-~ B03DB3BD-0E33-~ D649B60C-E631-~
#> form flow fictitious length name2 sinkdepth Shape_Leng
#> <chr> <chr> <chr> <int> <chr> <num> <num>
#> inlandRiver in direction false 431 NA -1 431.3
#> inlandRiver in direction false 513 NA -1 513.1
#> inlandRiver in direction false 740 NA -1 739.8
The main properties are present under appropriate names and a call to
the add_channel method with no options would be successful.
However if we want to carry over the additional information in the
“identifier” as “channel_id” a named vector giving the variable names to
be used could be provided. In this case:
property_names <- c(channel_id="identifier",
endNode="endNode",
startNode="startNode",
length="length")
The river network can then be created in the correct format by
ctch$add_channel(sp_lines,property_names)
#> Warning in private$apply_add_channel(channel, property_names, default_width):
#> Modifying to spatial polygons using default width
Since the data set for Swindale does not contain a channel width the
default width of 2m is used.
Computing Basic Properties
So far, the DEM and channel data exist in isolation. Next, we compute
some basic summaries for each square in the DEM, specifically: - land
area - the area in the DEM cell covered by land - channel area - the
area in the DEM cell covered by channel - channel id - the id of the
channel within the cell, corresponding to the id in the channel
object.
If multiple river lengths intersect a DEM cell the properties of the
channel length with the largest area of intersection are used.
This computation is done by calling
Filling sinks
For the hill slope to be connected to the river network all DEM cells
must drain to those that intersect with the river network.
The algorithm of implemented in the sink_fill method
ensures this is the case. Since the algorithm is iterative the execution
time of the function is limited by capping the maximum number of
iterations. If this limit is reached without completion the method can
call again with the “hot start” option to continue from where it
finished.
For Swindale, where the example DEM is already partially filled the
algorithm only alters a small area near the foot of the catchment.
ctch$sink_fill()
terra::plot( ctch$get_layer('filled_dem') - ctch$get_layer('dem'),
main="Changes to height")
Computing properties
Two sets of properties are required for Dynamic TOPMODEL. The first
set is those required within the evaluation of the model; gradient and
contour length. The second set are those used for dividing the catchment
up into Hydrological Response Units (HRUs). Traditionally the summary
used for the separation of the HRUs is the topographic index, which is
the natural logarithm of the upslope area divided by gradient.
These are computed using the formulae in Quinn et al. 1991.
The upstream area is computed by routing down slope with the fraction
of the area being routed to the next downstream pixel being proportional
to the gradient times the contour length.
The local value of the gradient is computed using the average of a
subset of between pixel gradients. For a normal ‘hill slope’ cell these
are the gradients to downslope pixels weighted by contour length. In the
case of pixels which contain river channels the average of the gradients
from upslope pixels weighted by contour length us used.
These properties are computed in an algorithm that passes over the
data once in descending height. It is called as follows
ctch$compute_properties()
The plot of the topographic index shows a pattern of increasing
values closer to the river channels
## plot of topographic index (log(a/tan b))
ctch$plot_layer('atb')
Adding additional layer
Properties may come in additional GIS layers. To demonstrate the
addition of an additional layer we will extract the filled dem
tmp <- ctch$get_layer("filled_dem")
then separate it into a layers representing land above and below
500m.
tmp <- terra::classify( tmp,
matrix(c(0,500,NA,
500,1000,-999),
ncol=3, byrow=TRUE))
The resulting raster object can now be written to a file
terra::writeRaster(tmp,file.path(demo_dir,"greater_500.tif"))
which is added to the meta data with
ctch$add_layer("greater_500",file.path(demo_dir,"greater_500.tif"))
ctch$get_layer()
#> [1] "dem" "channel" "filled_dem" "land_area" "channel_area"
#> [6] "channel_id" "gradient" "upslope_area" "atb" "greater_500"
Flow distances and ordering
Since dynatop simulations make use of ordered HRUs to
work downslope, a metric is required to order the downslope sequencing.
The calculation of four such metrics is supported
- shortest flow length - the shortest length based on the
pixel flow paths to a channel
- Dominant flow length - the distance to a channel moving in
the dominant (largest fraction) flow direction from any grid cell
- Expected flow length - the distance to the channel based on
a weighted average of the downslope flow lengths. Weights are given by
the fraction of flow in each direction.
- Band - A strict computational order, starting at the
channel arranged such that if all pixels in band \(i+1,i+2,\ldots\) are evaluated the inflows
to band \(i\) are known.
The computation is initiated with
ctch$compute_flow_lengths()
The additional layers can be examined as expected
ctch$get_layer()
#> [1] "dem" "channel" "filled_dem"
#> [4] "land_area" "channel_area" "channel_id"
#> [7] "gradient" "upslope_area" "atb"
#> [10] "band" "shortest_flow_length" "dominant_flow_length"
#> [13] "expected_flow_length" "greater_500"
ctch$plot_layer("band")
Classifying into Hydrological Response Units
Methods are provided for the classification of the catchment areas of
similar hydrological response. The classifications generated in this
process are augmented with a further distance based separation when
generating a dynatop model (see following section).
By definition each channel length is treated as belonging to a single
class.
To classify the hillslope two methods can be used.
The classify method of a dynatopGIS allows
a landscape property to be cut into classes.
For example to cut the topographic index for Swindale into 21
classes:
ctch$classify("atb_20","atb",cuts=20)
ctch$plot_layer("atb_20")
Providing a single value to the cuts argument determines the number
of classes. The values used to cut the variable can be extracted from
the meta data with
ctch$get_class_method("atb_20")
#> $layer
#> [1] "atb"
#>
#> $cuts
#> [1] 8.636791 9.265304 9.893817 10.522330 11.150843 11.779356 12.407868
#> [8] 13.036381 13.664894 14.293407 14.921920 15.550433 16.178946 16.807458
#> [15] 17.435971 18.064484 18.692997 19.321510 19.950023 20.578536 21.207048
The combine_classes method of a dynatopGIS
allows classes to be combined in two ways, which are applied in the
order shown
- pairing - where unique combinations of classes create one
new class
- burning - where a single class is imposed upon an area
To demonstrate a pairing combination consider combining the atb
classes generated above with the classification provided by the distance
band
ctch$combine_classes("atb_20_band",c("atb_20","band"))
ctch$plot_layer("atb_20_band")
Additionally the land greater then 500 in altitude can be burnt in
with
ctch$combine_classes("atb_20_band_500",pairs=c("atb_20","band"),burns="greater_500")
ctch$plot_layer("atb_20_band_500")
The each class in the combined classification the values of the
classes used in the computations can be returned
head( ctch$get_class_method("atb_20_band_500") )
#> atb_20_band_500 atb_20 band greater_500
#> 1 1 1 1 NA
#> 2 2 2 1 NA
#> 3 3 3 1 NA
#> 4 4 4 1 NA
#> 5 5 5 1 NA
#> 6 6 4 2 NA
Note that by giving the area to be burnt in a negative value when it
was generated above we have ensured that the values do not clash with
those generated by the cuts which (except potentially when a cut is NA)
will always be positive.
Generating a dynamic TOPMODEL
A Dynamic TOPMODEL suitable for use with the dynatop
package can be generated using the create_model method.
This uses an existing classification to generate the model. The required
model structure is given in the vignettes of dynatop
package and is not described here in details.
Since dynatop simulations make use of ordered HRUs to
work downslope, a classification which used a distance layer (see
earlier section) which represents the ordered downslope sequencing of
the pixels is recommended. It is strongly recommended that the
‘band’ distance metric is used directly as shown below.
Even if a distance layer is not used in the classification one must
be given to the create_model method, so the resulting HRUs
can be ordered.
For example, in the case of the division of Swindale by topographic
index into 21 classes and the bands directly the resulting model can be
generated by
ctch$create_model("new_model","atb_20_band","band")
#> The following Channel HRUs are outflows: 1
Looking at the files within the demo_dir folder
list.files(demo_dir,pattern="new_model*")
#> [1] "new_model.rds" "new_model.tif"
shows that an addition raster map of the HRUs has been created in
new_model.tif along with a file new_model.rds
containing a model suitable for dynatop
The map of HRUs can be plotted as with any layer
ctch$plot_layer("new_model")
The values on the map correspond to the ìd column of the
hillslope table in the dynatop model.
