Pre-processing, and GRASS commands used to locate optimal routes as presented in the talk this afternoon. Presentation can be found here. Hints and comments are inline with commands. Looks like there was some coverage of this talk in the Wired Online Magazine .
General concept: want to plan a wilderness trip, deviating from the main trail, hoping to visit numerous lakes. In the notes, the lakes that will be visited are called points of interest (poi). We would like to automatically locate the paths between points of interest and the trail which require the least amount of effort: i.e. minimizing the total change in slope encountered while walking. This makes sense, as we all know that walking uphill or down steep hills is hard. Furthermore, walking a long way up or down a steep slope is even harder. Thus, a simple slope map can be used as a travel "friction" map for the input to the command r.cost. Additional sources of "travel friction" can be added to this slope map as follows:
- traversing lakes is not acceptable (extreme "travel friction" added to slope map)
- traversing wooded areas is faster due to shading (moderate "travel friction" subtracted from slope map)
Note that these features can be automatically extracted from our USGS topo map (Figures 2 and 3). Figure 4 depicts this composite travel "friction" map. By specifying start and stop points (one lake to another, or from the main trail to a lake), along with our composite travel friction map we can compute the least-cost path between these points via the GRASS commands r.cost and r.drain. Figure 5 demonstrates the output from r.cost and r.drain, in this 3D rendition it is possible to see how r.drain is actually simulating the flow of a drop of water down an elevation surface: an energy minimizing style operation. However in this case, the elevation surface is actually a travel cost map (computed by r.cost). The path of the theoretical drop of water represents the least-cost route from the starting point to the end point. Note that the starting point for r.cost (point of origin) is actually the ending point for r.drain. See code examples below for details. Notice that by specifying very large travel costs for traversing lakes (peaks in Figure 5), and slightly lower travel costs for wooded areas, the computed lease cost path in Figure 5 goes "around" lakes and "through" wooded areas. Figure 6 depicts the new trail network (red lines) for the wilderness area, connecting Park Service trails (black) and our points of interest (yellow boxes).
A final, press-ready map was produced by GMT suitable for printing and taking out into the field for assisting with navigation. UTM ticks are useful in associating realtime GPS positional data with the printed map. This map can be downloaded here as a PDF. The script used to produce the map can be found here. Note that this is a slightly modified r.out.gmt script. I will post these changes back to the GRASS wikki for integration into the main codebase.
Figure 1: area of interest
Figure 2: Lake Features
Figure 3: wooded areas
Figure 4: composite friction map
Figure 5: least-cost path
Figure 6: vectorized least-cost trails
Download and Process DRG data (USGS Topo Maps)
#kaiser peak wget http://gis.ca.gov/casil/gis.ca.gov/drg/7.5_minute_series_albers_nad27_trimmed/37119/o37119b2.tif wget http://gis.ca.gov/casil/gis.ca.gov/drg/7.5_minute_series_albers_nad27_trimmed/37119/o37119b2.tfw # #huntington lake wget http://gis.ca.gov/casil/gis.ca.gov/drg/7.5_minute_series_albers_nad27_trimmed/37119/o37119c2.tif wget http://gis.ca.gov/casil/gis.ca.gov/drg/7.5_minute_series_albers_nad27_trimmed/37119/o37119c2.tfw # #convert to UTM zone 11 NAD83: gdalwarp -tps -t_srs "+proj=utm +zone=11 +datum=NAD83" o37119c2.tif o37119c2-utm.tif gdalwarp -tps -t_srs "+proj=utm +zone=11 +datum=NAD83" o37119b2.tif o37119b2-utm.tif # #import into GRASS: r.in.gdal in=o37119b2-utm.tif out=o37119b2 r.in.gdal in=o37119c2-utm.tif out=o37119c2 # #fix null values in DRG files: #check to see what white is defined as d.rast o37119b2 r.what.rast # #set this value to null in both maps r.null map=o37119b2 setnull=0 r.null map=o37119c2 setnull=0 # #patch the two maps together g.region rast=o37119b2,o37119c2 # r.patch in=o37119b2,o37119c2 out=drg # #reset the output file's color table to that of one of the input files: r.colors map=drg rast=o37119c2
Extract wooded areas and lakes from topo maps:
#reclass drg to find water features: r.reclass in=drg out=l <<< EOF 9 = 1 lakes * = NULL EOF # #grow water features, to make them contiguous (remove some noise) r.grow in=l out=lg radius=3 # #make the grown water features blue: echo "1 blue "| r.colors map=lg color=rules # #find lakes by isolating areas greater than 0.3ha r.reclass.area in=lg out=lakes greater=0.3 echo "1 blue "| r.colors map=lakes color=rules # #optionally vectorize these features: r.to.vect -s in=lakes out=lakes feature=area # #locate wooded areas from topo map r.reclass in=drg out=trees g.region res=3 ; r.grow in=trees out=trees_grow radius=3 r.reclass.area in=trees_grow out=trees_final greater=5 d.rast trees_final
Create travel "friction" map:
### reset the res to that of the slope map:
g.region align=slope
#
#update our slope map, to include traversing water features, and prefering wooded areas
#add a "cost" of 1000 to lake areas
r.mapcalc "new_slope = if(isnull(lakes) == 0, 1000.0+slope, slope)"
#
#subtract a small amount of cost for wodded areas:
r.mapcalc "new_slope = if(isnull(trees_final) == 0, abs(new_slope - 10.0), new_slope)"
Compute least-cost paths:
# ------- trail --> 3
#
start_poi=trail_to_
stop_poi=3
#
start_xy=`d.where -1 | awk '{print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope from point
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
#
# ------- 3 --> 4
#
start_poi=3
stop_poi=4
#
start_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$start_poi"'$/ {print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope from point
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
# ------- 3 --> 10
#
start_poi=3
stop_poi=10
#
start_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$start_poi"'$/ {print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope from point
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
#
# ------- 4 --> 5
#
start_poi=4
stop_poi=5
#
start_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$start_poi"'$/ {print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
# ------- 5 --> trail
#
start_poi=5
stop_poi=_to_trail
#
start_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$start_poi"'$/ {print $1","$2}'`
stop_xy=`d.where -1 | awk '{print $1","$2}'`
#
#compute a cost surface based on slope
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
#
# ------- trail --> 8
#
start_poi=trail_to
stop_poi=8
#
start_xy=`d.where -1 | awk '{print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
#
# ------- 8 --> 9
#
start_poi=8
stop_poi=9
#
start_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$start_poi"'$/ {print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
#
# ------- 9 --> 11
#
start_poi=9
stop_poi=11
#
start_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$start_poi"'$/ {print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
# ------- 11 --> 13
#
start_poi=11
stop_poi=13
#
start_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$start_poi"'$/ {print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
# ------- trail --> 12
#
start_poi=trail_to
stop_poi=12
#
start_xy=`d.where -1 | awk '{print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
#
#
# ------- 12 --> 10
#
start_poi=12
stop_poi=10
#
start_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$start_poi"'$/ {print $1","$2}'`
stop_xy=`v.out.ascii poi | awk -F"|" '$3 ~ /^'"$stop_poi"'$/ {print $1","$2}'`
#
#compute a cost surface based on slope
r.cost -k in=new_slope out=slope_cost coordinate=$start_xy stop_coordinate=$stop_xy --o
#
#compute least cost path:
r.drain in=slope_cost out=d$start_poi$stop_poi coordinate=$stop_xy --o
Combine all least-cost paths and vectorize:
#setup the color table for all of the least cost paths: for x in `g.mlist type=rast pattern=d*[3,4,5,6,8,9,1,2,3,10,to_trail] sep=" "`; do echo "1 red" | r.colors map=$x color=rules ; done # ### patch together our least cost paths, and thin r.patch in=`g.mlist type=rast pattern=d*[3,4,5,6,8,9,1,2,3,10,to_trail] sep=","` out=my_trails --o r.thin in=my_trails out=my_trails_thin --o # #vectorize: r.to.vect -s in=my_trails_thin out=my_trails feature=line --o # # ### patch together my trails and the main trails v.patch in=my_trails,tf out=all_trails # ### use v.path to loctate small "breaks" v.path all_trails d.zoom # #then fix the breaks with "move vertex" tool v.digit all_trails # #clean the line work and start again: v.clean in=all_trails out=all_trails_clean tool=break thresh=2 --o # #convert lines to ploylines v.build.polylines in=all_trails_clean out=poly # #add cats v.category in=poly out=poly_cats option=add # #add attribute table v.db.addtable map=poly_cats # #populate table with line length echo "alter table poly_cats add column length double" | db.execute v.to.db map=poly_cats option=length units=mi column=length
