This tutorial demonstrates simple cases of $\chi$ calculation with the gtbox
package.
The first thing we have to do is to load the gtbox
package (once it has been installed).
It is required to have GRASS GIS (version 7) installed, as some of the functions used here will rely heavily on GRASS raster processing modules (through the rgrass7
package).
It will also be necessary to install the following Addons (through the command g.extension
in a GRASS session) :
r.stream.distance
r.stream.basins
r.stream.order
Note that we will use the terra
package for most geospatial operations in R. We will use the standard formats from this package : SpatRaster
for raster layers (eventually with multiple layers) and SpatVector
for vector layers.
More information on the terra package can be found here :
library("gtbox") library("terra")
We first load a DEM.
data("dem_cevennes",package = "gtbox") dem = terra::rast(dem_cevennes) rm(dem_cevennes)
We can compute a hillshade surface and use some transparency on the DEM for a nicer display.
slope <- terra::terrain(dem, "slope", unit="radians") aspect <- terra::terrain(dem, "aspect", unit="radians") hill <- terra::shade(slope, aspect, 40, 270) plot(hill, col=grey(0:100/100), legend=FALSE) plot(dem, col=rainbow(25, alpha=0.35), add=TRUE)
We then process the DEM, as illustrated in the Simple Digital Elevation Model operations vignette.
gisbase = "/usr/lib/grass78/" th_px = 2000 # stream initiation threshold in pixels st = process_dem(dem,th_px,gisbase=gisbase)
We are going to compute a chi map over the DEM, starting the integration at a specified elevation. We first get the corresponding outlets and basins.
outlets = get_outlets(st,elevation=250,gisbase=gisbase) basins = get_basins(st,outlets,gisbase=gisbase)
We then compute the $\chi$ values along the network.
chi = compute_chi(st,outlets,a0=1000,mn=0.5,gisbase)
We transfer the $\chi$ values defined along streams to adjacent areas to obtain a continuous map, easier to visualize.
chi_map = streams2basins(st,chi,FUN="mean") chi_map = mask(chi_map,basins)
We plot everything
plot(hill, col=grey(0:100/100), legend=FALSE) plot(chi_map, col=heat.colors(25, alpha=0.35), add=TRUE) lines(basins) points(outlets)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.