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

Create a new mapset (called e.g. HW_buffers_cost) in nc_spm_08_grass7 project. 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 potentially 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
v.to.rast roadsmajor out=roadsmajor_28m use=value
r.buffer roadsmajor_28m 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

What is the total developed area iin [ha] within the 250m from the major roads?

Find schools potentially 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 legend_label="schools"
d.vect schools_noise size=14 fill_color=cyan icon=basic/circle color=black legend_label="noise risk schools"
d.legend.vect
d.barscale
d.out.file mynoisemap

How many students are potentially affected by noise based on the school capacities and less than 250m distance to major roads?

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 -n streets_speed color=inferno
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 cut=100
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
d.vect myfirestations icon=basic/circle size=15 fill_color=red

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.vect myfirestations icon=basic/circle size=15 fill_color=red
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

Given a point where a lost person was last located we compute a map that represents how far the person could have walked within the surrounding environment. The environmental variables for accessibility are derived from elevation, land use and streets using the equation in the r.walk manual page and the friction_rules.txt file.

First, create friction map based on land cover.
Note: We remove previously used map layers in the Layer Manager or using the d.erase tool, so the layers are not rendered when we don't need them anymore.
Set the computational region and display land cover classes:

d.erase
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 impact of 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
d.out.file myfrictionmap

Compute the cumulative cost map from the given point 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 cut=100

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

How did we ensure that the person does not walk over the lakes?

Optional

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

Compute the shortest distance map and cost surface to highways

: Set region, convert vector road map to raster:

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

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