<jupyter_start><jupyter_text>Geospatial Analysis II:  buffers, cost surfaces, least cost pathResources:* [GRASS GIS overview and manual](https://grass.osgeo.org/grass76/manuals/index.html)*  [Recommendations](data_acquisition.htmlcommands)and [tutorial](./grass_intro.html)how to use GUI from the first assignmentText files with recode rules, color rules and site locations:* [noise_cats.txt](data/noise_cats.txt)* [friction_rules.txt](data/friction_rules.txt)* [friction_color.txt](data/friction_color.txt)* [fire_pt.json](data/fire_pt.json)* [lostperson.txt](data/lostperson.txt) Start GRASS GISIn 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) andfor GRASS Mapset create a new mapset (called e.g. HW_buffers_cost) andclick Start GRASS session.<jupyter_code># using Python to initialize GRASS GIS
# This is a quick introduction into Jupyter Notebook.
# Python code can be executed like this:
a = 6
b = 7
c = a * b
print("Answer is", c)
# Python code can be mixed with command line code (Bash).
# It is enough just to prefix the command line with an exclamation mark:
!echo "Answer is $c"
# Use Shift+Enter to execute this cell. The result is below.
# using Python to initialize GRASS GIS
import os
import sys
import subprocess
from IPython.display import Image

# create GRASS GIS runtime environment
gisbase = subprocess.check_output(["grass", "--config", "path"], text=True).strip()
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))

# do GRASS GIS imports
import grass.script as gs
import grass.script.setup as gsetup

# set GRASS GIS session data
rcfile = gsetup.init(gisbase, "/home/jovyan/grassdata", "nc_spm_08_grass7", "user1")
# using Python to initialize GRASS GIS
# default font displays
os.environ['GRASS_FONT'] = 'sans'
# overwrite existing maps
os.environ['GRASS_OVERWRITE'] = '1'
gs.set_raise_on_error(True)
gs.set_capture_stderr(True)
# using Python to initialize GRASS GIS
# set display modules to render into a file (named map.png by default)
os.environ['GRASS_RENDER_IMMEDIATE'] = 'cairo'
os.environ['GRASS_RENDER_FILE_READ'] = 'TRUE'
os.environ['GRASS_LEGEND_FILE'] = 'legend.txt'<jupyter_output><empty_output><jupyter_text>Remove MASK in case it was left overfrom previous assignment (use command `r.mask -r`).Change working directory:_Settings_ > _GRASS working environment_ > _Change working directory_ > select/create any directoryor type `cd` (stands for change directory) into the GUI_Console_ and hit Enter:<jupyter_code># a proper directory is already set, download files
import urllib.request
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/noise_cats.txt", "noise_cats.txt")
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/friction_rules.txt", "friction_rules.txt")
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/friction_color.txt", "friction_color.txt")
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/fire_pt.json", "fire_pt.json")
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/lostperson.txt", "lostperson.txt")
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/noise_cats.txt", "noise_cats.txt")
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/fire_pt.json", "fire_pt.json")
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/friction_rules.txt", "friction_rules.txt")
urllib.request.urlretrieve("http://ncsu-geoforall-lab.github.io/geospatial-modeling-course/grass/data/lostperson.txt", "lostperson.txt")<jupyter_output><empty_output><jupyter_text>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 filewithout the need to specify the full path to the file. Buffers Find developed areas close to lakesSet region, create buffers:<jupyter_code>!g.region swwake_30m -p
!r.buffer lakes output=lakes_buff distances=60,120,240,500
!d.rast lakes_buff
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>List categories in land use map to identify category numbers for developed areas.Then run _r.mapcalc_ to extract the developed areas within the buffersand use _r.support_ to assign the labels from the original buffer raster to the new developed buffer raster map (needed for legend):<jupyter_code>!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<jupyter_output><empty_output><jupyter_text>Display results:In GUI: switch off previous map, use Add raster, Add vector icons.<jupyter_code>!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
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Find the total area within buffers and the developed area in ha:<jupyter_code>!r.report -n lakes_buff units=h
!r.report -n developed_lake units=h<jupyter_output><empty_output><jupyter_text>Find developed areas impacted by noise from highwaysSet 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](data/noise_cats.txt)is in your working directory).<jupyter_code>!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<jupyter_output><empty_output><jupyter_text>Display the results. In GUI switch off previous maps.<jupyter_code>!d.rast noise
!d.vect streets_wake
!d.legend -c noise at=5,25,2,5
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Find schools affected by high levels of noiseConvert schools to raster using CORECAPACI attribute (school capacity).Use map algebra to overlay with noise impact buffers andcompute the number of students exposed to noise (see result of _r.univar_).<jupyter_code>!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
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Cost surfaces Compute cumulative cost surface to a given accident site based on speed limitsUse 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_<jupyter_code>!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
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Import the accident point location from GeoJSON file [fire_pt.json](data/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.<jupyter_code>!v.in.ogr input=fire_pt.json output=fire_pt<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>!r.mapcalc "streets_travtime = 0.018641/streets_speed"
!r.cost -k streets_travtime output=streets_cost start_points=fire_pt<jupyter_output><empty_output><jupyter_text>Compute isochrones and display the cumulative cost surface map.<jupyter_code>!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
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Find cost (travel time) from selected firestationsFirst 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.<jupyter_code>!g.copy vector=firestations,myfirestations
!v.info -c myfirestations
!v.what.rast myfirestations raster=streets_cost column=CVLAG<jupyter_output><empty_output><jupyter_text>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](img/firestations_table.png).Export firestations with traveltime less than 0.1 hr:<jupyter_code>!v.out.ascii input=myfirestations separator=space precision=3 columns=ID,LOCATION,CVLAG where="CVLAG<0.1 AND CVLAG>0"<jupyter_output><empty_output><jupyter_text>You should get something like the following:```635775.565 228121.693 19 19 4021 District Dr 0.076635940.262 225912.796 20 0 5001 Western Blvd 0.037637386.831 222569.152 21 0 1721 Trailwood Dr 0.071633178.155 221353.037 52 27 6000 Holly Springs Rd 0.060```To get the computed time, you can also query the cumulative cost raster directlyusing coordinates (in this example it's Western Blvd firestation),or using GUI query tool by selecting the cost layer and clicking on the map:<jupyter_code>!r.what map=streets_cost coordinates=635940.262,225912.796 separator=space<jupyter_output><empty_output><jupyter_text>Find the least cost path for the two closest stations:<jupyter_code>!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<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>!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
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>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?<jupyter_code>!r.describe route_20Westernb
!r.describe route_52Hollyb<jupyter_output><empty_output><jupyter_text>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 numbersin hours are very low.Then switch off all layers except for the streets_cost_1000and change display to 3D. Compute accessibility map for help in search for lost personCreate friction map based on land cover.Reminder: Time to time, you should remove or at least uncheck previouslyused map layers in the _Layer Manager_, so the layers are not renderedwhen you don't need them anymore.First display land cover classes:<jupyter_code>!g.region swwake_30m -p
!r.category landclass96<jupyter_output><empty_output><jupyter_text>Recode the landuse map to friction map using the rules in [friction_rules.txt](data/friction_rules.txt).<jupyter_code>!r.recode landclass96 out=friction rules=friction_rules.txt<jupyter_output><empty_output><jupyter_text>Add the streets to friction map - they are missed by the landuse map.<jupyter_code>!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
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Compute the cost map and generate isochrones:<jupyter_code>!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<jupyter_output><empty_output><jupyter_text>To display the input, import coordinates of the point where the lost person was last seengiven in the file [lostperson.txt](data/lostperson.txt).<jupyter_code>!v.in.ascii input=lostperson.txt output=lostperson separator=comma<jupyter_output><empty_output><jupyter_text>Uncheck the friction map and display the results:<jupyter_code>!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
Image(filename="map.png")<jupyter_output><empty_output><jupyter_text>Cost surfaces to line featuresCompute the shortest distance to a major road.Set region, convert vector road map to raster:<jupyter_code>!g.region swwake_30m -p
!v.to.rast roadsmajor output=roadsmajor use=val type=line<jupyter_output><empty_output><jupyter_text>Compute the distance map and cost surface to highways:<jupyter_code>!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
Image(filename="map.png")
!d.vect streets_wake color=black
!d.vect roadsmajor color=red
Image(filename="map.png")
# end the GRASS session
os.remove(rcfile)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0958444581365099, 'fraction_numerical': 0.027452591223386048, 'mean_word_length': 6.921908893709327, 'pattern_counts': {'":': 0, '<': 110, '<?xml version=': 0, '>': 117, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '479380', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 5197, 'n_tokens_neox': 4809, 'n_words': 1344, 'path': 'notebooks/buffers_cost.ipynb', 'repo_name': 'wenzeslaus/geospatial-modeling-course-jupyter', 'size': '14607'}