In VincentGodard/gtbox: A Geomorphology Toolbox With Various Functions For The Analysis of Topography and Surface Processes
Objectives
This tutorial demonstrates the first steps of topographic analysis with the gtbox package.
- getting a Digital Elevation Model into the working environment
- computing various relevant raster layers
- extracting basins
- extracting river profiles
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.distancer.stream.basinsr.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")
Importing a DEM
The typical workflow consists in reading a Digital Elevation Model (DEM) and computing a number of basic rasters (flow accumulation and direction, etc ...) which will be used as a basis for further processing.
We start by reading and plotting a DEM using the.
This is a SRTM DEM (1'') over the Cévennes area in SE France (downloaded from opentopography.org), which is included in the gtbox package.
Alternatively we could import a DEM using the rast function of the terra package.
data("dem_cevennes",package = "gtbox")
dem = terra::rast(dem_cevennes)
rm(dem_cevennes)
#dem = terra::rast("srtm_cevennes.tif") # we could import our own DEM this way
plot(dem)
And here are the properties of this raster.
print(dem)
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)
Processing the DEM
Then we use the process_dem function to launch an analysis of the DEM and compute various associated rasters.
This part of the treatment relies on GRASS GIS, as this function will call some GRASS GIS modules, so we need to explain where to find them on the system with the gisbase argument.
The other important argument is the accumulation value (number of pixels) for channel initiation.
The output will be organized as a multi-layers SpatRaster.
gisbase = "/usr/lib/grass78/"
th_px = 2000 # stream initiation threshold in pixels
st = process_dem(dem,th_px,gisbase=gisbase)
print(st)
This SpatRaster the following layers (plus eventually some other optional layers) :
|Raster|Explanation|
|---|---|
|z|Elevation from input DEM |
|acc|Flow accumulation (in pixels) computed with GRASS module r.watershed |
|dir|Flow direction computed with GRASS module r.watershed |
|dist| Distance along network computed with GRASS module r.stream.distance (Addons) |
|st_id| Stream segments unique identifiers (integer) computed with GRASS module r.watershed |
|sto| Strahler orders of the stream segments|
|bs_id| Elementary sub-basins identifiers (integer) corresponding to st_id (r.watershed) |
We can export the stream raster (st_id) for visualization purposes in a GIS (e.g. QGIS).
We extract the corresponding SpatRaster layer prior to exporting into a GeoTif raster.
terra::writeRaster(st$st_id,"/tmp/streams.tif",overwrite=TRUE)
Extracting basins
Now that we have our raster layers computed and neatly organized, we want to extract some basins to process.
Using user-defined outlets
We can first import some corresponding outlets, which should be positioned on the river network pixels (that's why we exported the streams.tif raster above).
data("outlets_cevennes",package = "gtbox")
outlets = terra::vect(outlets_cevennes)
rm(outlets_cevennes)
#outlets = terra::vect("../data/outlets_cevennes.gpkg") # we could import our own outlets this way (EPSG:32631)
Once the outlets positions have been imported as SpatVector object, they can be used to compute the corresponding basins using the get_basins function (it again relies on GRASS modules), which require the flow direction raster (dir) from our SpatRaster as an input.
We can export the vector to the GeoPackage format for external visualization with QGIS.
basins = get_basins(st,outlets,gisbase)
writeVector(basins,"/tmp/basins.gpkg",overwrite=TRUE)
print(basins)
We can then plot these basins and their outlets on a map.
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
plot(st$z, col=rainbow(25, alpha=0.35), add=TRUE)
lines(basins)
points(outlets,pch=21,bg="white",cex=1.5)
text(outlets,basins$cat,cex=0.7)
The attribute table of the outlets vector is passed to the new basins vector.
We can now select a particular basin (using cat values for example) and subset the SpatRaster accordingly.
bas = basins[basins$cat==3,] # get basin of interest
tmp = terra::mask(terra::crop(st,bas),bas) # clip and crop
plot(tmp$z)
lines(bas)
According to their Strahler order
The get_outlets function allows also to define outlets for a given Strahler order.
outlets_sto = get_outlets(st,strahler=3)
basins_sto = get_basins(st,outlets_sto,gisbase)
We can plot these basins.
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
plot(st$z, col=rainbow(25, alpha=0.35), add=TRUE)
lines(basins_sto)
points(outlets_sto,pch=21,bg="white",cex=1.5)
Extracting the largest possible basins
Another approach is to try to extract the largest possible basins from the DEM.
Note that pixels with negative accumulation, where there is a possible contributions of areas located outside of the DEM, are not considered.
outlets_large = get_outlets(st,large=10e4) # we specify a minimum area in pixels
basins_large = get_basins(st,outlets_large,gisbase)
And we plot these basins.
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
plot(st$z, col=rainbow(25, alpha=0.35), add=TRUE)
lines(basins_large)
points(outlets_large,pch=21,bg="white",cex=1.5)
Extracting river network
For the whole network
We use the get_network function to process the SpatRaster object.
net = get_network(st,gisbase)
plot(net)
We obtain a SpatVector object, which contains the network as well as an attribute table with various information concerning the stream segments (the field stream correspond to st_id layer of the input SpatRaster).
We can plot this network, using the Strahler order of the streams segments for width.
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
plot(st$z, col=rainbow(25, alpha=0.35), add=TRUE)
lines(net,lwd=net$strahler,col="blue")
For a particular basin
We can also extract the network for a particular basin, after cropping and clipping the SpatRast to the contours of this basin.
bas0 = basins[basins$cat==3,] # get basin of interest
st0 = terra::mask(terra::crop(st,bas0),bas0) # clip and crop
ids = terra::unique(st0$st_id)$st_id # get the stream ids located in the clipped area
net0 = net[net$stream%in%ids,] # subset the network according to that
plot(st0$z)
lines(bas0)
lines(net0,lwd=net0$strahler,col="blue")
Extracting river profiles
We can also extract the river profile quite easily. Let's obtain the largest basin, using the same steps as above.
bas = basins_large[1,] # largest basin
st0 = terra::mask(terra::crop(st,bas),bas) # clip and crop
ids = terra::unique(st0$st_id)$st_id # get the stream ids located in the clipped area
net0 = net[net$stream%in%ids,] # subset the network according to that
We now identify the stream id corresponding to the stream with the largest accumulation inside the basin. We use the get_streams function to collect all the ids of stream segments moving upstream from this starting point.
id0 = net0[which.max(net0$flow_accum)]$stream # get id of stream from which we start
ids0 = get_streams(net0,id0,mode="up") # get upstream ids
We now convert the SpatRaster as a dataframe, subset to keep only network pixels with ids corresponding to our selection, and ordering according to distance along stream.
stream0 = terra::as.data.frame(st0) # convert to data frame
stream0 = stream0[stream0$st_id%in%ids0,] # select pixel along the profile
stream0 = stream0[order(stream0$dist),] # order
plot(stream0$dist,stream0$z,type="l",lwd=3,xlab="Distance along stream (m)",ylab="Elevation (m)")
VincentGodard/gtbox documentation built on Sept. 4, 2022, 3:46 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Objectives
This tutorial demonstrates the first steps of topographic analysis with the gtbox package.
- getting a Digital Elevation Model into the working environment
- computing various relevant raster layers
- extracting basins
- extracting river profiles
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.distancer.stream.basinsr.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")
Importing a DEM
The typical workflow consists in reading a Digital Elevation Model (DEM) and computing a number of basic rasters (flow accumulation and direction, etc ...) which will be used as a basis for further processing.
We start by reading and plotting a DEM using the.
This is a SRTM DEM (1'') over the Cévennes area in SE France (downloaded from opentopography.org), which is included in the gtbox package.
Alternatively we could import a DEM using the rast function of the terra package.
data("dem_cevennes",package = "gtbox") dem = terra::rast(dem_cevennes) rm(dem_cevennes) #dem = terra::rast("srtm_cevennes.tif") # we could import our own DEM this way plot(dem)
And here are the properties of this raster.
print(dem)
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)
Processing the DEM
Then we use the process_dem function to launch an analysis of the DEM and compute various associated rasters.
This part of the treatment relies on GRASS GIS, as this function will call some GRASS GIS modules, so we need to explain where to find them on the system with the gisbase argument.
The other important argument is the accumulation value (number of pixels) for channel initiation.
The output will be organized as a multi-layers SpatRaster.
gisbase = "/usr/lib/grass78/" th_px = 2000 # stream initiation threshold in pixels st = process_dem(dem,th_px,gisbase=gisbase) print(st)
This SpatRaster the following layers (plus eventually some other optional layers) :
|Raster|Explanation| |---|---| |z|Elevation from input DEM | |acc|Flow accumulation (in pixels) computed with GRASS module r.watershed | |dir|Flow direction computed with GRASS module r.watershed | |dist| Distance along network computed with GRASS module r.stream.distance (Addons) | |st_id| Stream segments unique identifiers (integer) computed with GRASS module r.watershed | |sto| Strahler orders of the stream segments| |bs_id| Elementary sub-basins identifiers (integer) corresponding to st_id (r.watershed) |
We can export the stream raster (st_id) for visualization purposes in a GIS (e.g. QGIS).
We extract the corresponding SpatRaster layer prior to exporting into a GeoTif raster.
terra::writeRaster(st$st_id,"/tmp/streams.tif",overwrite=TRUE)
Extracting basins
Now that we have our raster layers computed and neatly organized, we want to extract some basins to process.
Using user-defined outlets
We can first import some corresponding outlets, which should be positioned on the river network pixels (that's why we exported the streams.tif raster above).
data("outlets_cevennes",package = "gtbox") outlets = terra::vect(outlets_cevennes) rm(outlets_cevennes) #outlets = terra::vect("../data/outlets_cevennes.gpkg") # we could import our own outlets this way (EPSG:32631)
Once the outlets positions have been imported as SpatVector object, they can be used to compute the corresponding basins using the get_basins function (it again relies on GRASS modules), which require the flow direction raster (dir) from our SpatRaster as an input.
We can export the vector to the GeoPackage format for external visualization with QGIS.
basins = get_basins(st,outlets,gisbase) writeVector(basins,"/tmp/basins.gpkg",overwrite=TRUE) print(basins)
We can then plot these basins and their outlets on a map.
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4)) plot(st$z, col=rainbow(25, alpha=0.35), add=TRUE) lines(basins) points(outlets,pch=21,bg="white",cex=1.5) text(outlets,basins$cat,cex=0.7)
The attribute table of the outlets vector is passed to the new basins vector.
We can now select a particular basin (using cat values for example) and subset the SpatRaster accordingly.
bas = basins[basins$cat==3,] # get basin of interest tmp = terra::mask(terra::crop(st,bas),bas) # clip and crop plot(tmp$z) lines(bas)
According to their Strahler order
The get_outlets function allows also to define outlets for a given Strahler order.
outlets_sto = get_outlets(st,strahler=3) basins_sto = get_basins(st,outlets_sto,gisbase)
We can plot these basins.
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4)) plot(st$z, col=rainbow(25, alpha=0.35), add=TRUE) lines(basins_sto) points(outlets_sto,pch=21,bg="white",cex=1.5)
Extracting the largest possible basins
Another approach is to try to extract the largest possible basins from the DEM. Note that pixels with negative accumulation, where there is a possible contributions of areas located outside of the DEM, are not considered.
outlets_large = get_outlets(st,large=10e4) # we specify a minimum area in pixels basins_large = get_basins(st,outlets_large,gisbase)
And we plot these basins.
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4)) plot(st$z, col=rainbow(25, alpha=0.35), add=TRUE) lines(basins_large) points(outlets_large,pch=21,bg="white",cex=1.5)
Extracting river network
For the whole network
We use the get_network function to process the SpatRaster object.
net = get_network(st,gisbase) plot(net)
We obtain a SpatVector object, which contains the network as well as an attribute table with various information concerning the stream segments (the field stream correspond to st_id layer of the input SpatRaster).
We can plot this network, using the Strahler order of the streams segments for width.
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4)) plot(st$z, col=rainbow(25, alpha=0.35), add=TRUE) lines(net,lwd=net$strahler,col="blue")
For a particular basin
We can also extract the network for a particular basin, after cropping and clipping the SpatRast to the contours of this basin.
bas0 = basins[basins$cat==3,] # get basin of interest st0 = terra::mask(terra::crop(st,bas0),bas0) # clip and crop ids = terra::unique(st0$st_id)$st_id # get the stream ids located in the clipped area net0 = net[net$stream%in%ids,] # subset the network according to that plot(st0$z) lines(bas0) lines(net0,lwd=net0$strahler,col="blue")
Extracting river profiles
We can also extract the river profile quite easily. Let's obtain the largest basin, using the same steps as above.
bas = basins_large[1,] # largest basin st0 = terra::mask(terra::crop(st,bas),bas) # clip and crop ids = terra::unique(st0$st_id)$st_id # get the stream ids located in the clipped area net0 = net[net$stream%in%ids,] # subset the network according to that
We now identify the stream id corresponding to the stream with the largest accumulation inside the basin. We use the get_streams function to collect all the ids of stream segments moving upstream from this starting point.
id0 = net0[which.max(net0$flow_accum)]$stream # get id of stream from which we start ids0 = get_streams(net0,id0,mode="up") # get upstream ids
We now convert the SpatRaster as a dataframe, subset to keep only network pixels with ids corresponding to our selection, and ordering according to distance along stream.
stream0 = terra::as.data.frame(st0) # convert to data frame stream0 = stream0[stream0$st_id%in%ids0,] # select pixel along the profile stream0 = stream0[order(stream0$dist),] # order plot(stream0$dist,stream0$z,type="l",lwd=3,xlab="Distance along stream (m)",ylab="Elevation (m)")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.