A tutorial on quantitative terrain analysis in GRASS GIS.
Contents
Geomorphometry is the quantitative analysis of topography. Geomorphometric analyses include slope, aspect, curvature, topographic indices, and landforms. GRASS GIS includes many modules and addons for geomorphometric analysis including:
In this tutorial we will use the modules
r.param.scale
and
r.geomorphon
to automatically classify the landforms on Governor’s Island.
We will also use the addon module
r.convergence
to compute the topographic convergence index and classify landforms.
This tutorial uses the
Governor’s Island Dataset for GRASS GIS.
Download, extract, and move this geospatial dataset
for Governor’s Island in New York City
to your grassdata
directory.
Start GRASS GIS,
set the GRASS GIS database directory to grassdata
directory,
select nyspf_governors_island
as your location,
and create a new mapset called geomorphometry
.
Zoom in on the landforms in the southwest of the island.
Either set the computation from the display
using the various zoom options dropdown or
run g.region
and set the boundaries for the region. Save the region.
Then set a mask to the vector map shoreline
with the module
r.mask.
g.region n=189850 s=189100 e=978550 w=976850 save=landforms
r.mask vector=shoreline
Digital Elevation Model |
---|
The module r.param.scale calculates the morphometric parameters of topography using differential calculus. Terrain is parameterized based on the first and second order partial derivatives of the quadratic surfaces fit to a moving window using least squares. The scale of the morphometric parameterization depends on the size of the moving window. Parameters include slope, aspect, curvature, and morphometric features. Morphometric features, i.e. landforms, are defined by the relationship of a cell to its neighbors in terms of convexity and concavity, the second derivatives of the topographic surface. r.param.scale can identify six general landforms - peaks, ridges, passes, planes, channels, and pits. Use the module r.param.scale to automatically classify landforms. Set the moving window size to an odd number.
r.param.scale input=elevation_2017 output=landforms size=33 method=feature --overwrite
Landforms |
---|
Try running
r.param.scale
with different moving window sizes.
Either change the name of the output map or set the --overwrite
flag.
Add a legend.
Note how this module characterizes most of the landforms here
as either ridges or channels.
The module r.geomorphon automatically recognizes and classifies landforms using machine vision. Landforms are classified as either flats, peaks, ridges, shoulders, spurs, slopes, hollows, footslope, valleys, or pits based on their visibility from 8 cardinal and ordinal directions. Geomorphons has been used for diverse application such as characterizing submarine dunes on the ocean floor and creating global landform maps.
Landforms |
---|
Classify landforms using
r.geomorphon.
Experiment with the search
, skip
, and flat
parameters.
The search radius determines the scale of the landform features,
the flatness threshold set the angle at which ground is considered flat,
and the skip radius eliminates small landforms, reducing noise.
To better visualize the landforms
compute shaded relief from the digital elevation model using
r.relief.
Set the units to US survey feet
and optionally the vertical scale to 2 or higher.
Then drape a map of landforms over the shaded relief using
r.shade.
g.region region=landforms
r.geomorphon elevation=elevation_2017 forms=geomorphons search=36 skip=6 flat=12 --overwrite
r.relief input=elevation_2017 output=relief_2017 zscale=2 units=survey
r.shade shade=relief_2017 color=geomorphons output=shaded_geomorphons brighten=45
Geomorphons with Shaded Relief |
---|
Note how r.geomorphon has clearly identified the ridge lines and their peaks and has classified the pathways as either valleys or footslopes.
Use map algebra to extract the ridges from either map of landforms.
In the raster map calculator, use an if, then, else
statement.
If cells in the geomorphons raster equal 3,
then write a value of 1 in the new raster,
else write null values.
r.mapcalc expression="ridges = if(geomorphons==3,1,null())"
Set a color for the maps of ridges.
Right click on ridges layer in the layer manager
and choose set color interactively.
Assign a category value with
r.category.
In the Define
tab enter category value directly as
1|ridge
r.category map=ridges separator=pipe
Ridges |
---|
To make a simpler, cleaner vector map of the ridges,
first convert the raster map to a vector using the module
r.to.vect.
Then remove small areas from the vector map with the module
v.clean.
All areas smaller than the threshold parameter will be removed.
Then simplify the boundaries of the areas using
v.generalize
with the reumann
method.
Smooth the boundaries of the areas using
v.generalize
with the snakes
or hermite
method.
Remove the intermediate cleaned and simplified maps with
g.remove.
Add a vector legend for the ridges with
d.legend.vect.
r.to.vect -s input=ridges output=ridges type=area
v.clean input=ridges output=ridges_cleaned type=point,line,area tool=rmarea thres=2
v.generalize input=ridges_cleaned type=area output=ridges_generalized method=reumann threshold=2
v.generalize input=ridges_generalized type=area output=ridges method=snakes threshold=2 alpha=1 beta=1 --overwrite
g.remove -f type=vector name=ridges_cleaned,ridges_generalized
d.legend.vect at=2,95 font=Lato-Regular fontsize=14
Vector Ridges |
---|
Valleys and ridges can be identified from the convergence and divergence of the topography. Use g.extension to install the addon module r.convergence. Compute the convergence index of the terrain with r.convergence with a moving window size of 15 cells. In the resulting raster values from 1 to 100 are convergent, 0 is planar, and values from -1 to -100 are divergent.
g.extension extension=r.convergence
r.convergence -c input=elevation_2017 output=convergence window=15 weights=standard
Topographic Convergence |
---|
Use map algebra to extract ridges from areas of topographic divergence.
In the raster map calculator, use an if, then, else
statement.
If cells in the convergence raster are less than or equal to a threshold,
then write a value of 1 in the new raster,
else write null values.
The threshold for ridges should be a negative value between -1 and -100.
Experiment to find the right threshold.
This example uses a threshold of -15.
Then convert the raster map to a vector using the module
r.to.vect,
remove small areas from the vector map with the module
v.clean,
simplify the boundaries of the areas using
v.generalize,
and then smooth the boundaries of the areas using
v.generalize.
Remove the intermediate maps with
g.remove.
Add a vector legend for the ridgelines with
d.legend.vect.
r.mapcalc expression="ridgelines = if(convergence <= -15, 1, null())"
r.to.vect -s input=ridgelines output=ridgelines type=area
v.clean input=ridgelines output=ridges_cleaned type=point,line,area tool=rmarea thres=75
v.generalize input=ridges_cleaned type=area output=ridges_generalized method=reumann threshold=2
v.generalize input=ridges_generalized type=area output=ridgelines method=snakes threshold=2 alpha=1 beta=1 --overwrite
g.remove -f type=vector name=ridges_cleaned,ridges_generalized
d.legend.vect at=2,95 font=Lato-Regular fontsize=14
Ridges Derived from Topographic Divergence |
---|