Geospatial Modeling and Analysis

Geomorphometry II: Spatial and Temporal Terrain Analysis

Resources: GRASS GIS overview and manual, GRASSbook.

Download Mapset and color tables:

  • Download NagsHead time series and copy it into your nc_spm_08_grass7 Location (it should be the same level directory as PERMANENT). Do not let your unzipping program create additional level directory with the same name! If you are not sure about GRASS GIS Database structure read about it in the manual.
  • Custom color table for time series standard deviations map stddev_color.txt


Start GRASS - click on GRASS icon or type

In startup pannel set GIS Data Directory to path to datasets, for example on MS Windows, C:\Users\myname\grassdata.
For Project location select nc_spm_08_grass7 (North Carolina, State Plane, meters) and for Accessible mapset create a new mapset (called e.g. HW_terrain_analysis).
Click Start GRASS.

Change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI Console and hit Enter:

Download all text files (see above) to the selected directory. Now you can use the commands from the assignment requiring the text file without the need to specify the full path to the file.

Compute basic topographic parameters: slope and aspect

g.region raster=elevation -p
r.slope.aspect elevation=elevation slope=myslope aspect=myaspect

Display resulting maps with legend using GUI.

d.rast myslope
d.legend myslope at=2,40,2,6
d.out.file myslope
d.rast myaspect
d.legend myaspect at=2,40,2,6
d.out.file myaspect

Show impact of integer values in meters on slope and aspect pattern.
Compute integer DEM and derive its slope and aspect.
Use GUI to display the histogram: in Map Display > Analyze > Create histogram:

r.mapcalc "elev_int = int(elevation)"
r.slope.aspect elevation=elev_int aspect=aspect_int_10m slope=slope_int_10m

d.histogram slope_int_10m
d.histogram myslope
d.histogram aspect_int_10m
d.histogram myaspect

d.rast myslope

Zoom into NW area of the current region (relatively flat area near large interchange).
Can you explain the difference in slope maps derived from integer (m vertical resolution) and floating point (mm vertical resolution) DEMs?

d.rast slope_int_10m
d.out.file slope_int

Compute slope along road

First set the region to the extent of the bus route #11 and to 10m resolution. Then convert the vector line of the route to raster using the direction of the route.
g.region vect=busroute11 align=elevation res=10 -p input=busroute11 type=line output=busroute11_dir use=dir
r.colors map=busroute11_dir color=aspect
d.rast busroute11_dir
d.legend busroute11_dir
Compute the steepest slope of the topography along the route by assigning the values of slope derived from a DEM in the first part of this assignment to the grid cells along the route.
r.mapcalc "route_slope = if(busroute11_dir, myslope)"
Then compute the slope in the direction of the route using the difference between aspect of the topography and the route direction angles.
r.mapcalc "route_slope_dir = abs(atan(tan(myslope) * cos(myaspect - busroute11_dir)))"
r.colors map=route_slope,route_slope_dir color=slope
Display the results along with the elevation contours and compute univariate statistics. Comment on the difference of the two results.
r.contour input=elevation output=contours step=2
d.vect contours
d.rast route_slope
d.legend route_slope
d.out.file route_slope
d.rast route_slope_dir
d.out.file route_slope_dir

r.univar route_slope
r.univar route_slope_dir


Compute slope, aspect and curvatures simultaneously with interpolation.
You can do the examples below for the bare earth data only (first example),
multiple return example is optional (if you are curious how it differs from BE).

g.region rural_1m  -p input=elev_lid792_bepts elevation=elev_lid_1m aspect=asp_lid_1m pcurvature=pc_lid_1m tcurvature=tc_lid_1m npmin=120 segmax=25 input=elev_lidrural_mrpts elevation=elev_lidmr_1m aspect=asp_lidmr_1m pcurvature=pc_lidmr_1m tcurvature=tc_lidmr_1m npmin=120 segmax=25 tension=300 smooth=1.

Display the results as 2D images or in 3D view.
For 3D view, switch off everything except for elevation surface that you want to view.
Drape topographic parameters raster maps over DEMs as color.

d.rast elev_lid_1m
d.rast pc_lid_1m
d.out.file profcurvature1
d.rast elev_lidmr_1m
d.rast pc_lidmr_1m

The curvature maps reflect survey pattern rather than topographic features.
So we lower the tension and increase the smoothing.
You can use multiple displays to compare the results side-by-side.

g.region rural_1m  -p input=elev_lid792_bepts elevation=elev_lidt15_1m aspect=asp_lidt15_1m pcurvature=pc_lidt15_1m tcurvature=tc_lidt15_1m npmin=120 segmax=25 tension=15 smooth=1. input=elev_lidrural_mrpts elevation=elev_lidmrt15_1m aspect=asp_lidmrt15_1m pcurvature=pc_lidmrt15_1m tcurvature=tc_lidmrt15_1m npmin=120 segmax=25 tension=15 smooth=1.

d.rast elev_lidt15_1m
d.rast pc_lidt15_1m
d.out.file profcurvature2
d.rast tc_lidt15_1m
d.rast elev_lidmrt15_1m
d.rast pc_lidmrt15_1m


Extract landforms at different levels of detail by adjusting the size of moving window.
Set rural subregion at 1m resolution, compute landforms using 9m and 45m neighborhood: read the manual to learn more.
Explain types of landforms and the role of the neighborhood size.
g.region rural_1m -p
r.param.scale elev_lid792_1m output=feature9c_1m size=9  method=feature
r.param.scale elev_lid792_1m output=feature45c_1m size=45 method=feature

Display with legend, save images for the report.
Optionally display the feature maps draped over elev_lid792_1m as color.

d.rast feature9c_1m
d.legend feature9c_1m at=2,20,2,6
d.rast feature45c_1m
d.legend feature45c_1m at=2,20,2,6
d.vect elev_lid792_cont1m color=brown

Raster time series analysis

For this exercise we will use NagsHead_series Mapset you downloaded.
You have to first make the mapset accessible.
In GUI: menu Settings -> GRASS working environment -> Mapset access or by using a command:
g.mapsets operation=add mapset=NagsHead_series
If you don't see the mapset, make sure you downloaded it and unzipped it correctly.

Run the series analysis and explain the results:
Which maps are core and envelope? Which landforms have high standard deviation and what does it mean?

g.region raster=NH_2008_1m -p

d.rast NH_2008_1m
r.series NH_1999_1m,NH_2001_1m,NH_2004_1m,NH_2005_1m,NH_2007_1m,NH_2008_1m out=NH_9908_min method=minimum
r.series NH_1999_1m,NH_2001_1m,NH_2004_1m,NH_2005_1m,NH_2007_1m,NH_2008_1m out=NH_9908_max method=maximum
r.series NH_1999_1m,NH_2001_1m,NH_2004_1m,NH_2005_1m,NH_2007_1m,NH_2008_1m out=NH_9908_mintime method=min_raster
r.series NH_1999_1m,NH_2001_1m,NH_2004_1m,NH_2005_1m,NH_2007_1m,NH_2008_1m out=NH_9908_maxtime method=max_raster
r.series NH_1999_1m,NH_2001_1m,NH_2004_1m,NH_2005_1m,NH_2007_1m,NH_2008_1m out=NH_9908_range method=range
r.series NH_1999_1m,NH_2001_1m,NH_2004_1m,NH_2005_1m,NH_2007_1m,NH_2008_1m out=NH_9908_avg method=average
r.series NH_1999_1m,NH_2001_1m,NH_2004_1m,NH_2005_1m,NH_2007_1m,NH_2008_1m out=NH_9908_stddev method=stddev

r.colors NH_9908_stddev rules=stddev_color.txt

d.rast NH_9908_stddev
d.rast NH_9908_range

Use cutting plane in 3D view to show the core and envelope.
Add constant elevation plane at -1m for reference, set zexag somewhere 3-5 (the default is too high).
Assign surfaces constant color, use top or bottom surface for crossection color.
When using top for color, lower the light source to make top surface dark and highlight the crossection.

Optional: Cut and fill and volume

Compute cut and fill for 4m deep excavation to build a facility. First, set the region and display facility on top of orthophoto.
g.region rural_1m -p

d.rast ortho_2001_t792_1m
d.rast facility
Then set (raster) mask to the facility map and find minimum elevation within the facility:
r.mask raster=facility
r.univar elevation
Minimum which you obtain should be 123.521m. Bottom of 4m excavation will be then
123.52 - 4 = 119.52
Use raster algebra to create the excavation:
r.mapcalc "excavation = elevation - 119.52"
r.univar excavation
d.rast excavation
Minimum you get should be 4.00057 and maximum 9.50554. Note that the excavation is limited by the mask we set earlier, so we can now do global operation to compute the volume which applies just the the facility.
r.volume excavation
Now remove mask. This is important so that your future work is not affected.
r.mask -r