GIS/MEA582: Geospatial Modeling and Analysis

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

Geospatial Analysis I: global, zonal and focal operations, map algebra

Resources:

GRASS GIS overview and manual
Recommendations and tutorial how to use wxGUI from the first assignment

To run mapcalc expressions, you can either run the entire command in the GUI Command console, or in case of any problems, type or copy the expression in the GRASS GIS Raster Map Calculator which can be launched from Layer Manager toolbar.

Text files with color rules:

srtmneddiff_color.txt

Start GRASS GIS

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_map_algebra) and click Start GRASS.

grass70

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 with color rules (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 summaries

Compute areas for each category at two different resolutions.

Are results equal? Explain (see also Lecture 1). Copy and paste the report from the output window or save the report in a text file: Output window > Save. Use fixed width font (e.g., Courier, Andale Mono in your report to preserve formatting).

g.region raster=landuse96_28m res=12 -ap
r.report landuse96_28m unit=c,h,p
g.region raster=landuse96_28m res=30 -ap
r.report landuse96_28m unit=c,h,p

Compute areas for each category of land use for each zipcode. Compare 27601 Raleigh with 27511 Cary. Include only the relevant part of the table in your report.

r.report zipcodes,landuse96_28m unit=h,p

Compute zonal statistics maps representing average slope for each basin.
Add legends using Add map elements in Map Display toolbar. Reminder: d.out.file means Save to graphics file for your report.

g.region raster=slope -p
r.stats.zonal base=basin_50K cover=slope method=average output=slope_avgbasin
r.colors slope_avgbasin color=gyr
d.rast slope_avgbasin
d.legend slope_avgbasin at=90,50,5,8
d.vect streams color=15:25:110
d.out.file slope_avgbasin

Compute zonal statistics maps representing most common land use for each basin.
When adding legend, make sure you define List of discrete category numbers.

g.region raster=landuse96_28m -p
r.mode base=basin_50K cover=landuse96_28m output=landuse96_modebasin
d.rast landuse96_modebasin
d.vect streams
d.legend landuse96_modebasin at=40,20,2,5 -n
d.out.file landuse96_modezip

Apply neighborhood operators

Use neighborhood operator to compute land use diversity map.
How diverse is land use in NCSU area?
Before you start to display the new results, remove all previously added map layers from the Layer Manager.

g.region raster=landuse96_28m -p
r.info -g landuse96_28m
r.neighbors landuse96_28m output=lu_divers method=diversity size=7

d.rast lu_divers
d.legend lu_divers at=70,15,5,10 -v
d.vect streets_wake
r.report lu_divers unit=p

Use neighborhood operator to smooth the SRTM elevation map and compare the global statistical measures for the original and smoothed DEM. How would size of the neighborhood influence the result?

g.region raster=elev_srtm_30m -p
r.neighbors elev_srtm_30m output=elev_srtm30m_sm5 method=average size=5
d.rast elev_srtm_30m
d.rast elev_srtm30m_sm5
r.univar elev_srtm_30m
r.univar elev_srtm30m_sm5

Patch multiple raster layers into a single raster

Patch raster tiles for lidar based, 6m res. DEM for Centennial Campus. Before displaying new data, remove older map layers from Layer Manager.

g.region raster=el_D793_6m,el_D783_6m,el_D782_6m,el_D792_6m -p
r.patch input=el_D793_6m,el_D783_6m,el_D782_6m,el_D792_6m output=elevlidD_6m
r.colors elevlidD_6m raster=elevation

d.rast elevlidD_6m

Map Algebra

See r.mapcalc manual page for syntax and functions. If you are getting en error when running r.mapcalc in command console, launch the GUI version from Layer Manager toolbar.

Compute Normalized Difference Vegetation Index (NDVI).
Explain the difference between floating point and integer handling in ndvi1, ndvi2 and ndvi3 result.

g.region raster=lsat7_2002_40 -p
r.mapcalc "ndvi1 = (lsat7_2002_40 - lsat7_2002_30) / (lsat7_2002_40 + lsat7_2002_30)"
r.mapcalc "ndvi2 = 1.0 * (lsat7_2002_40 - lsat7_2002_30) / (lsat7_2002_40 + lsat7_2002_30)"
r.mapcalc "ndvi3 = float(lsat7_2002_40 - lsat7_2002_30) / float(lsat7_2002_40 + lsat7_2002_30)"
r.info -r ndvi1
r.info -r ndvi2
r.info -r ndvi3
r.colors ndvi3 color=ndvi
d.rast ndvi3
d.out.file ndvi3

Explore the difference between the SRTM DSM and lidar-based NED DEM. Compute the map of elevation differences.

g.region raster=elev_ned_30m -p
r.mapcalc "srtm_ned30_dif = elev_srtm_30m - elev_ned_30m"

Create a custom color table to distinguish the negative and positive values:

r.info -r srtm_ned30_dif

fp: Data range is -142.24... to 86.19...
Assign custom color table srtmneddiff_color.txt.
GUI: Right click on layer > Properties > Set color table > Colors > Path to rules file.

r.colors srtm_ned30_dif rules=srtmneddiff_color.txt

Zoom to computational region and switch off previous map layers. Display the difference map layer:

d.rast srtm_ned30_dif
d.legend srtm_ned30_dif at=70,15,5,10
d.vect streets_wake
d.out.file srtm_ned30_dif
r.univar elev_srtm_30m
r.univar elev_ned_30m

Is the srtm mostly higher or lower than elev_ned? Which result will you use to answer the above question - the srtm_ned30_dif map or numbers provided by r.univar?

Working with if statements

Create map of urban areas (high density and low density class) with 0 class elsewhere.

g.region raster=landuse96_28m -p
r.mapcalc "urban1_30m = if(landuse96_28m == 1,1,0) + if(landuse96_28m == 2,2,0)"
d.rast urban1_30m
d.legend urban1_30m at=10,30,5,8

Alternatively with logical 'or' operator and null elsewhere:

r.mapcalc "urban2_30m = if(landuse96_28m == 1 || landuse96_28m == 2,landuse96_28m,null())"
d.rast urban2_30m

Handling null values

Replace lakes in the landclass96 map with more accurate representation.
From raster map lakes, we divide lakes by 1000 because the lake map has large numbers that would have required custom color table.

r.mapcalc "luselakes = if(isnull(lakes), landclass96, lakes/1000)"
d.rast luselakes
d.out.file luselakes

Create mask for low lying developed areas (e.g. vulnerable to flooding) with elevation between 60 and 100m and land use 1 or 2.
Set the region, display the input maps and create a MASK.
Before you start new computations, remove or switch off previous map layers in the Layer Manger. You may also zoom to computational region in Map Display once you set a new one.

g.region raster=elevation -p

d.rast elevation
d.rast landuse96_28m
r.mapcalc "low_elev_developed = if((elevation < 100 && elevation > 60) && (landuse96_28m == 1 || landuse96_28m == 2),1,null())"
r.mask raster=low_elev_developed

Command r.mask creates a raster map "MASK" in your mapset.

Remove "low_elev_developed" layer if it was added.

Display watersheds to see the mask effect:

d.rast basin_50K
d.out.file basin_masked

Disable mask, and display basin_50K again to show that the mask was removed.

r.mask -r
d.rast basin_50K

Using coordinates of moving window in map algebra

Replace section of SRTM DSM with NED DEM elevation.
Try to explain how this r.mapcalc expression works.

r.mapcalc "elev_combined = if(y() < 224274. && x() > 637455., elevation, elev_srtm_30m)"

Tilted plane

Create a raster map representing tilted plane (e.g., geologic fault):

g.region rural_1m -p
r.mapcalc "tiltplane = 0.2*(0.1*row()+col())+50"
r.mapcalc "tiltpl_under = if(tiltplane < elev_lid792_1m + 2,tiltplane,null())"

View the elevation surface and subsurface plane in 3D. Switch off all layers in the layer manager except for elev_lid792_1m and tiltpl_under. Change display to 3D view, adjust viewing position to a view from South. Save an image for your report.

Optional - various additional useful tasks

Use map algebra to create map subsets.
Change region to the airphoto tile 792 and resolution 10m. Since we will work in different are, it is a good idea to remove all previously used map layers from the Layer Manager.

g.region raster=ortho_2001_t792_1m res=10 -ap

d.rast ortho_2001_t792_1m

Create a subset of the map elevation for this subregion.

r.mapcalc "elevation_792_10m = elevation"
d.rast elevation_792_10m

Zoom out to see that it is a subset.

Work with NULL and MASK. Set the mask and check its effect.

d.rast elevation
d.vect streets_wake
r.mask urban maskcats=55
d.rast elevation

Remove mask:

r.mask -r

Compute % area for each category in each zipcode. Which zipcode has the most high density development?

r.stats -pl zipcodes,landuse96_28m

Working with relative coordinates.
Enter the expression on a single line without \
Again, it is a good idea to remove the previously used map layers before we start to work on a new task.

g.region raster=elev_srtm_30m -p

r.mapcalc "elev_srtm30m_smooth = (elev_srtm_30m[-1,-1]   \
           + elev_srtm_30m[-1,0] + elev_srtm_30m[-1,1] \
           + elev_srtm_30m[0,-1] + elev_srtm_30m[0,0]  \
           + elev_srtm_30m[0,1]  + elev_srtm_30m[1,-1] \
           + elev_srtm_30m[1,0]  + elev_srtm_30m[1,1])/9."

Assign the resulting map the same color table as the original. Compare global statistics

r.colors elev_srtm30m_smooth raster=elev_srtm_30m
r.univar elev_srtm_30m
r.univar elev_srtm30m_smooth
d.rast elev_srtm_30m
d.rast elev_srtm30m_smooth