NCSU GIS/MEA582:
Geospatial Modeling and Analysis

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

Resources:

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

Text file with color rules:

Start GRASS GIS

Create a new mapset (called e.g. assignment3 or HW_map_algebra) 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

Download text file with color rules "srtmneddiff_color.txt" (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 and map summarry statistics

Compute zonal statistics maps representing average slope for each basin, often needed for hydrologic modeling.
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.vect streets_wake color=grey
d.out.file slope_avgbasin

Compute zonal statistics maps representing most common land use for each basin.
To get the best result, make sure you use all the flags from the example below. 

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 -f -c
d.out.file landuse96_modebasin
What is the most common land use in the two basins with the steepest average slope? You can use map query or map algebra to find out the naswer.

Apply neighborhood operators and map algebra

Map land use diversity

Use neighborhood operator to compute land use diversity map and create a map of locations with the most homogeneous (single category) landuse where a land use diversification could be beneficial (e.g. by planting trees). Use 7x7 cells window to find areas which have single landuse within 200x200m.

First, remove all previously added map layers from the Layer Manager, then compute and display the land use diversity map. Adjust the legend size and placement as needed using right click and mouse. 

d.erase
g.region raster=landuse96_28m -p
r.neighbors landuse96_28m output=lu_divers method=diversity size=7
d.rast lu_divers
d.legend lu_divers at=90,70,5,8 -v
d.vect streets_wake co=white
r.report lu_divers unit=p
d.out.file lu_diversity_map

Use map algebra to extract the single category areas and find out the area totals for each category using the report tool. 

r.mapcalc "landuse_1cat = if(lu_divers == 1, landuse96_28m, null())"
r.colors landuse_1cat raster=landuse96_28m
r.category landuse_1cat raster=landuse96_28m
r.report landuse_1cat unit=h,p -n
d.rast landuse_1cat
d.legend -c landuse_1cat
d.out.file lu_signle_map
Which land use covers the largest area with a single landuse category within 7x7 moving windows (most homogeneous areas)?
Use map algebra to extract a map layer that represents homogeneous "Mixed hardwoods / Conifers" areas (category 18) from the "landuse_1cat" layer that we may target for protection. Include the command and a resulting map with extracted map displayed together with "streets_wake" in your report

Smoothing elevation raster and analyzing differences between DEMs

Use neighborhood operator to smooth the SRTM elevation map "elev_srtm_30m" and compare the summary statistics for the original and smoothed SRTM DEM. First, display the lidar-based "elev_ned_30m" map, assign the same colors to the "elev_srtm_30m" and compare their values using legends. Then use neighborhood operator applied to 5x5 window to smooth the "elev_srtm_30m" and compare the results using the r.univar tool. How does the size of the neighborhood influence the resulting DEM? (You can test different sizes yourself.) 

d.erase
g.region raster=elev_srtm_30m -p
d.rast elev_ned_30m
d.legend elev_ned_30m
r.colors elev_srtm_30m raster=elevation
d.rast elev_srtm_30m
d.legend elev_srtm_30m
d.out.file srtm_dem_original
r.neighbors elev_srtm_30m output=elev_srtm30m_sm5 method=average size=5
d.rast elev_srtm30m_sm5
d.out.file srtm_dem_smoothed
r.univar elev_ned_30m
r.univar elev_srtm_30m
r.univar elev_srtm30m_sm5
Further explore the difference between the SRTM DEM "elev_srtm_30m" and lidar-based NED DEM "elev_ned_30m". First, compute the map of elevation differences using map algebra and find the range of differences: 
r.mapcalc "srtm_ned30_dif = elev_srtm_30m - elev_ned_30m"
r.info -r srtm_ned30_dif

Assign a divergent color table to distinguish the negative and positive values, 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.erase
d.rast srtm_ned30_dif
d.legend -d srtm_ned30_dif at=90,40,5,7
d.out.file srtm_ned30_dif

Are the elevations in "elev_srtm_30m" higher in most areas or lower than in "elev_ned_30m"? Which result will you use to answer the above question - the "srtm_ned30_dif" map or the numbers provided by r.univar? Are there any values in the elev_srtm_30m map that are not realistic? Where are they located?

Patch multiple raster layers into a single raster

Patch raster tiles for lidar based, 6m res. DEM for Centennial Campus. 
d.erase
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
d.vect roadsmajor
d.out.file patched_dem
Using map algebra, Compute difference between the "elevlidD_6m" and "elevation" maps, assign it a diverging color ramp (e.g. "differences" provided by r.colors) and compute a summary statistics of the differences. Include the map and the statistics into your report and comment on the resulting pattern

More Map Algebra

See r.mapcalc manual page for syntax and functions. You can run r.mapcalc from console using command line or GUI (Raster map calculator).

NDVI

Compute Normalized Difference Vegetation Index (NDVI).
Explain the difference between floating point and integer handling in ndvi1, ndvi2 and ndvi3 result. Make sure to zoom into computational area after changing the region. 
d.erase
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.legend ndvi3
d.out.file ndvi3
Note that this is a simplified, map algebra example, for computing various vegetation indices in GRASS GIS, we would use the i.vi tool after performing atmospheric corrections.

Working with if statements and null()

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
d.out.file urban12_nonurbanzero

Alternatively with logical 'or' operator and null elsewhere: 

d.erase
r.mapcalc "urban2_30m = if(landuse96_28m == 1 || landuse96_28m == 2,landuse96_28m,null())"
d.rast urban2_30m
d.legend urban2_30m at=10,30,5,8
d.out.file urban12_nonurbannull
What is the difference between the two urban raster maps?

Optional - various additional useful tasks

Handling null values and creating a mask

Create mask for low lying developed areas (e.g. vulnerable to flooding) with elevation between 60 and 100m and land use 1 or 2. We could use this mask to exclude this areas from future high density urbanization in an urban growth simulation.
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.erase
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.

Map subsets

Use map algebra to create map subsets.
Change region to the airphoto tile 792 and resolution 10m. Since we will work in different area, it is a good idea to remove all previously used map layers from Layers in the Layer Manager. 
g.region raster=ortho_2001_t792_1m res=10 -ap
d.erase
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.

Working with relative coordinates

Compute smoothed SRTM DEM using 3x3 moving window average. 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
d.erase
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