GIS/MEA582: Geospatial Modeling and Analysis

Geospatial Analysis II: buffers, cost surfaces, least cost path

Resources:

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

Text files with recode rules, color rules and site locations:

Start GRASS GIS

In startup pannel set GRASS GIS Database Directory to path to datasets, for example on MS Windows, C:\Users\myname\grassdata. For GRASS Location select nc_spm_08_grass7 (North Carolina, State Plane, meters) and for GRASS Mapset create a new mapset (called e.g. HW_buffers_cost) and click Start GRASS session.

grass83

Remove MASK in case it was left over from previous assignment (use command r.mask -r).

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 all text files with recode rules, color rules and site locations (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.

Buffers

Find developed areas close to lakes

Set region, create buffers:

g.region swwake_30m -p
r.buffer lakes output=lakes_buff distances=60,120,240,500
d.rast lakes_buff

List categories in land use map to identify category numbers for developed areas. Then run r.mapcalc to extract the developed areas within the buffers and use r.support to assign the labels from the original buffer raster to the new developed buffer raster map (needed for legend):

r.category landuse96_28m
r.mapcalc "developed_lake = if(landuse96_28m==1 || landuse96_28m==2, lakes_buff, null())"
r.support developed_lake raster=lakes_buff
r.category developed_lake

Display results:
In GUI: switch off previous map, use Add raster, Add vector icons.

d.rast developed_lake
d.vect streets_wake color=grey
d.rast lakes
d.legend -c developed_lake use=2,3,4,5 at=5,25,2,5
d.out.file mylakesbuffmap

Find the total area within buffers and the developed area in ha:

r.report -n lakes_buff units=h
r.report -n developed_lake units=h

Find developed areas impacted by noise from highways

Set region and create buffers along major roads.
Intersect developed areas from landuse map with road buffers.
Transfer the category labels and compute the affected area.
Run r.category from GUI and browse to the txt file, or run r.category command (assuming noise_cats.txt is in your working directory).

g.region raster=landuse96_28m -p
r.buffer roadsmajor output=roads_buffers distances=250,500,2500
r.mapcalc "noise = if(landuse96_28m==1 || landuse96_28m==2, roads_buffers, null())"
r.colors noise color=ryg
r.category noise rules=noise_cats.txt separator=:
r.report -n noise units=p,h

Display the results. In GUI switch off previous maps.

d.rast noise
d.vect streets_wake
d.legend -c noise at=5,25,2,5

Find schools affected by high levels of noise

Convert schools to raster using CORECAPACI attribute (school capacity).
Use map algebra to overlay with noise impact buffers and compute the number of students exposed to noise (see result of r.univar).

v.to.rast schools_wake output=schoolscap_10m use=attr attrcolumn=CORECAPACI type=point
r.mapcalc "schools_noise = if(int(schoolscap_10m) && roads_buffers == 2, int(schoolscap_10m), null())"
r.univar schools_noise
r.to.vect schools_noise output=schools_noise type=point
d.vect schools_wake icon=basic/circle size=10 fill_color=blue
d.vect schools_noise size=14 fill_color=cyan icon=basic/circle color=black
d.out.file mynoisemap

Cost surfaces

Compute cumulative cost surface to a given accident site based on speed limits

Use 5mi/hr speed limit for off-road areas (nulls).
For legend set the discrete speed limit values in GUI legend dialog:
tab Required > streets_speed, Subset > List of discrete category numbers

g.region swwake_30m -p
v.info -c streets_wake
v.to.rast streets_wake output=streets_speedtmp use=attr attrcolumn=SPEED type=line
r.mapcalc "streets_speed = if(isnull(streets_speedtmp),5,streets_speedtmp)"
r.info streets_speed
r.colors streets_speed color=gyr
d.rast streets_speed
d.legend streets_speed at=5,40,2,5 use=5,25,35,45,65
d.out.file myspeedmap

Import the accident point location from GeoJSON file fire_pt.json
Run v.in.ogr from GUI, use file browser to provide full path to the file, or run the command assuming the file is in your working directory.

v.in.ogr input=fire_pt.json output=fire_pt

Assign travel time to cross a 30m grid cell in hours.
Note that cost in GRASS GIS is defined as travel time per cell.

Compute cumulative cost surface to the given point.
Where does the 0.018641 constant come from? You can modify the expression to get time in minutes.

r.mapcalc "streets_travtime = 0.018641/streets_speed"
r.cost -k streets_travtime output=streets_cost start_points=fire_pt

Compute isochrones and display the cumulative cost surface map.

r.contour streets_cost output=streets_cost_04 step=0.04
r.colors streets_cost color=elevation
d.rast streets_cost
d.vect fire_pt color=red icon=basic/marker size=20
d.vect streets_cost_04
d.legend streets_cost at=5,50,2,5
d.out.file mycostmap

Find cost (travel time) from selected firestations

First make your own copy of the firestations map and list attributes.
Then query the cumulative cost surface at the firestations location.
The travel time in hours will be stored in the attribute column CVLAG.

g.copy vector=firestations,myfirestations
v.info -c myfirestations
v.what.rast myfirestations raster=streets_cost column=CVLAG

Add myfirestations to Layer Manager and right click on Show attribute data.
Click on CVLAG to order them and find the lowest cost (shortest time) > 0 (firestations with 0 cost are outside the region).
You should see a table like this.
Export firestations with traveltime less than 0.1 hr:

v.out.ascii input=myfirestations separator=space precision=3 columns=ID,LOCATION,CVLAG where="CVLAG<0.1 AND CVLAG>0"

You should get something like the following:

635775.565 228121.693 19 19 4021 District Dr 0.076
635940.262 225912.796 20 0 5001 Western Blvd 0.037
637386.831 222569.152 21 0 1721 Trailwood Dr 0.071
633178.155 221353.037 52 27 6000 Holly Springs Rd 0.060

To get the computed time, you can also query the cumulative cost raster directly using coordinates (in this example it's Western Blvd firestation), or using GUI query tool by selecting the cost layer and clicking on the map:

r.what map=streets_cost coordinates=635940.262,225912.796 separator=space

Find the least cost path for the two closest stations:

r.drain -n input=streets_cost output=route_20Westernb start_coordinates=635940.3,225912.8
r.drain -n input=streets_cost output=route_52Hollyb start_coordinates=633178.2,221353.0
r.colors route_20Westernb color=grey
r.colors route_52Hollyb color=grey

Display the results in GUI.
You already may have some of these maps in GUI layer tree, so just switch them on and off and re-order as needed:

d.erase
d.rast streets_cost
d.vect fire_pt color=red icon=basic/marker size=20
d.rast route_20Westernb
d.rast route_52Hollyb
d.out.file mylcpmap

Print the length of the path in cells (multiply by 30m to get approx. m).
You should have the time in hr already from the cost map.
At what average speed [km/hr] needs the truck travel to get there in estimated time? Is the time, speed and distance realistic?

r.describe route_20Westernb
r.describe route_52Hollyb

If you would like to display streets_cost surface in 3D as shown in the lecture, you need to multiply the surface by 1000 using r.mapcalc because the numbers in hours are very low.
Then switch off all layers except for the streets_cost_1000 and change display to 3D.

Compute accessibility map for help in search for lost person

Create friction map based on land cover.
Reminder: Time to time, you should remove or at least uncheck previously used map layers in the Layer Manager, so the layers are not rendered when you don't need them anymore.
First display land cover classes:

g.region swwake_30m -p
r.category landclass96

Recode the landuse map to friction map using the rules in friction_rules.txt.

r.recode landclass96 out=friction rules=friction_rules.txt

Add the streets to friction map - they are missed by the landuse map.

r.mapcalc "friction2 = if(streets_speed > 6, 0.1, friction)"
r.colors friction2 rules=friction_color.txt
d.erase
d.rast friction2
d.legend friction2 at=5,50,2,5

Compute the cost map and generate isochrones:

r.walk -k elevation=elev_ned_30m friction=friction2 output=walkcost start_coordinates=635576,216485 lambda=0.5 max_cost=10000
r.contour walkcost output=walkcost step=1000

To display the input, import coordinates of the point where the lost person was last seen given in the file lostperson.txt.

v.in.ascii input=lostperson.txt output=lostperson separator=comma

Uncheck the friction map and display the results:

d.erase
r.colors -e map=walkcost color=viridis
d.rast lakes
d.rast walkcost
d.vect streets_wake
d.vect walkcost color=red
d.vect walkcost color=red where=level=6000 width=3
d.vect lostperson color=yellow fill_color=yellow icon=basic/marker size=30
d.legend walkcost at=5,50,2,5
d.out.file mywalkcostmap

Cost surfaces to line features

Compute the shortest distance to a major road.
Set region, convert vector road map to raster:

g.region swwake_30m -p
v.to.rast roadsmajor output=roadsmajor use=val type=line

Compute the distance map and cost surface to highways:

r.mapcalc "area_one = 1"
r.cost input=area_one output=dist_toroad start_rast=roadsmajor
r.mapcalc "dist_meters = dist_toroad * (ewres() + nsres())/2."
r.mapcalc "dist_class = int(dist_meters/500)"
d.erase
d.rast dist_class
d.vect roadsmajor
d.out.file dist_roadsmaj
r.cost -k input=streets_travtime output=cdist_toroadk start_rast=roadsmajor
r.colors cdist_toroadk color=bgyr
d.rast cdist_toroadk
d.vect streets_wake color=black
d.vect roadsmajor color=red
d.legend cdist_toroadk at=5,50,2,5
d.out.file cost_roadsmaj