GIS/MEA582: Geospatial Modeling and Analysis

This is an unmaintained course material, please see current material at:

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.
stddev_color.txt
regrslope_color.txt

Start GRASS GIS

Start GRASS - click on GRASS icon or type

grass70

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 command console and hit Enter:

cd

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 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 SW 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

Compute slope of a bus route. First extract values from slope map layer:

g.region raster=elevation res=20
v.to.rast input=busroute11 type=line output=busroute11 use=dir
r.mapcalc "route_slope = if(busroute11, myslope)"

Then compute slope in the direction of the route:

v.to.rast input=busroute11 type=line output=busroute11_dir use=dir
r.mapcalc "route_slope_dir = abs(atan(tan(myslope) * cos(myaspect - busroute11_dir)))"
r.colors map=route_slope,route_slope_dir color=slope

Display, zoom if necessary, and compute univariate statistics and comment on the difference of two results.

r.univar route_slope
r.univar route_slope_dir

d.shade shade=elevation_shade color=myslope brighten=20
r.contour input=elevation output=contours step=2
d.vect contours
d.rast route_slope
d.out.file route_slope
d.rast route_slope_dir
d.out.file route_slope_dir

Curvatures

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
v.surf.rst elev_lid792_bepts elevation=elev_lid_1m aspect=asp_lid_1m pcurvature=pc_lid_1m tcurvature=tc_lid_1m npmin=120 segmax=25
v.surf.rst 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
v.surf.rst 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.
v.surf.rst 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

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