NCSU GIS/MEA582:
Geospatial Modeling and Analysis

Geomorphometry II: Spatial and Temporal Terrain Analysis

Resources: GRASS GIS overview and manual, GRASSbook.

Start GRASS GIS

Create a new mapset (called e.g. HW_terrain_analysis) in nc_spm_08_grass7 project and 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: 
cd
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

To explore computation of basic topographic parameters we set our computational region to NW section of our elevation DEM and compute steepest slope and aspect maps (make sure to zoom to computational region) 
g.region raster=elevation -p
g.region s=s+7000 e=e-7000 -p
r.slope.aspect elevation=elevation slope=myslope aspect=myaspect
r.colors myaspect co=aspectcolr

DEMs are sometimes distributed with elevations stored as integer values, changing the pattern of topographic parameters. To show impact of integer elevation values in meters on slope and aspect pattern we first compute an integer DEM, its slope and aspect and then compare the resulting maps with the original floating point DEM results 

r.mapcalc "elev_int = int(elevation)"
r.slope.aspect elevation=elev_int aspect=aspect_int_10m slope=slope_int_10m
r.colors aspect_int_10m co=aspectcolr
d.erase
d.rast myslope
d.legend myslope at=2,40,2,6
d.out.file myslopemap
d.rast slope_int_10m
d.out.file slope_intmap
d.rast myaspect
d.legend myaspect at=2,40,2,6
d.out.file myaspect
d.rast aspect_int_10m
d.out.file aspect_intmap

It is also helpful to compare the histograms (you can switch off the orange rectangle showing computational region by clicking on the Map Display Settings icon on top of Map Display and switching off "Show computational extent") 

d.erase
d.histogram myslope
d.out.file histmyslope
d.histogram slope_int_10m
d.out.file histslope_int_10m
d.histogram myaspect
d.out.file myaspectmap
d.histogram aspect_int_10m
d.out.file aspect_intmap

Can you explain the difference in slope and aspect maps and their histograms derived from integer (m vertical resolution) and floating point DEMs?

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
v.to.rast input=busroute11 type=line output=busroute11_dir use=dir
d.rast busroute11_dir
d.legend busroute11_dir
d.out.file bus11direction
Compute the steepest slope of the topography around the route by assigning the values of slope derived from a DEM 
r.slope.aspect elevation=elevation slope=slope_dem aspect=aspect_dem 
r.mapcalc "route_slope_dem = if(busroute11_dir, slope_dem)"
d.rast route_slope_dem
Then compute the slope in the direction of the route using the difference between aspect of the topography and the route direction angles. Display the results along with the elevation contours and compute univariate statistics. 
r.mapcalc "route_slope_dir = abs(atan(tan(slope_dem) * cos(aspect_dem - busroute11_dir)))"
r.colors map=route_slope_dem,route_slope_dir color=bcyr -e
r.contour input=elevation output=contours step=3
d.erase
d.vect contours color=220:200:200
d.rast route_slope_dem
d.legend route_slope_dem
d.out.file route_slope_dem
d.rast route_slope_dir
d.out.file route_slope_dir
r.univar route_slope_dem
r.univar route_slope_dir
What is the difference between the two slope maps? Which slope is steeper on average and why?

Curvatures

Compute curvatures simultaneously with interpolation of DEM from lidar data using regularized spline with tension. Evaluate the influence of the spline parameters on the spatial pattern of the curvature maps.

g.region rural_1m  -p
v.surf.rst input=elev_lid792_bepts elevation=elev_lid_1m pcurvature=pc_lid_1m tcurvature=tc_lid_1m npmin=120

Display the profile and tangential curvatures as 2D images.

d.erase
d.rast pc_lid_1m
d.legend pc_lid_1m 
d.out.file profcurvature_dtm
d.erase
d.rast tc_lid_1m
d.legend tc_lid_1m
d.out.file tangcurvature_dtm

The curvature maps reflect the lidar survey pattern rather than topographic features.
We can lower the tension and increase the smoothing to reduce the noise and map the curvatures of topography.

v.surf.rst input=elev_lid792_bepts elevation=elev_lidt15_1m aspect=asp_lidt15_1m pcurvature=pc_lidt15_1m tcurvature=tc_lidt15_1m npmin=120 tension=15 smooth=1.
d.erase
d.rast pc_lidt15_1m
d.legend pc_lidt15_1m
d.out.file profcurvature_dtm2
d.rast tc_lidt15_1m
d.legend tc_lidt15_1m
d.out.file tangcurvature_dtm2

Optionally, you can display the results in 3D view: switch off all layers except for elevation surface that you want to view and then drape the curvature maps over DEMs as color. How does changing the spline parameters influence the spatial pattern of curvatures?

Landforms

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.erase
d.rast feature9c_1m
d.legend -c feature9c_1m at=2,20,2,6 font=Arial
d.out.file landfroms9
d.rast feature45c_1m
d.vect elev_lid792_cont1m color=grey
d.out.file landfroms45

Map basic landforms using geomorphons DSM, read the manual page and the associated paper to learn how it works 

r.geomorphon elevation=elev_lid792_1m forms=geomorphon_s50d10 search=50 dist=10 
d.erase
d.rast geomorphon_s50d10
d.vect elev_lid792_cont1m color=grey
d.legend -c geomorphon_s50d10 at=2,30,2,8
d.out.file geomorphon_s50d10
Apply geomorphon to the 10m resolution DEM "elevation" (set region to raster elevation) and explain what is shown in the resulting map. Include the output map in your report.

Optional: Raster time series analysis

We will explore elevation time series analysis using a series of coastal DEMs provided in NagsHead_series Mapset

Download Mapset and color tables:

  • Download NagsHead time series and copy it into your nc_spm_08_grass7 project (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
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.erase
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
r.colors -n NH_9908_range color=inferno

d.erase
d.rast NH_9908_stddev
d.legend NH_9908_stddev
d.out.file series_stddev
d.rast NH_9908_range
d.legend NH_9908_range
d.out.file series_range

You can use cutting plane in 3D view to show the core and envelope.
Add constant elevation plane at -1m for reference, set zexag to 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. Include the output maps and the 3D view with a cutting plane in your report to illustrate the answers to the two questions.

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.erase
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