Navigating Wilderness Areas with GRASS (Where 2.0 Presentation)

Submitted by dylan on Wed, 2006-06-14 00:06.
Example DRG graphic
Figure 1: area of interest
Features extracted from DRG: lakes
Figure 2: lake features
Features extracted from DRG: trees
Figure 3: wooded areas
Example travel cost map
Figure 4: composite friction map
Graphical Example of least-cost path
Figure 5: least-cost path
vectorized trails
Figure 6: vectorized
least-cost trails

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.

 
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

Better check google maps

Better check google maps it's cool and very informative.

Hikeability of final trail

It has been my experience that USGS forested areas are unreliable for forest density. This means that an attempt to stay in wooded areas may cause you to enter impassable areas of very dense woods.

Did you experience this with your final trail?

final trail

DRG vs. DOQQ: forrest density: simple visual comparisonDRG vs. DOQQ: forrest density: simple visual comparison
Hi Robert,

I am planning on hiking the final trail next summer, and will post back then. I do not imagine that one would encounter very dense stands of trees in this area, as it is very near the tree line. The biggest obstacles would therefore be impassible stands of manzanita or shear rock faces. Thanks for the comment on the calculation of tree density with respect to estimates derived from USGS topo maps. Out of curiosity i prepared a visual comparison of DRG and DOQQ images: which illustrates your point.



Download DOQ data from TerraServer

g.region res=1 -a
r.in.wms -c -o -k output=doq mapserver=http://terraserver.microsoft.com/ogcmap6.ashx layers=DOQ format=jpeg srs=EPSG:4326
d.his i=doq h=drg
d.vect trails width=4 col=white
d.vect my_trails width=2 col=white
d.barscale -m bcol=none tcol=white
d.out.png res=1 out=drg_vs_doq-forrest

So, how did it go?

It's been a while since I visited here, but see you never posted the results of your attempt to hike the path that GRASS produced. Did you ever do it? How did it go? Have you given any thought to doing this using r.walk, too?

Comparison

Simple comparison between r.cost and r.walk posted here.

Not unitl next year...

Thanks for checking in Tom. Unfortunately I wasn't able to get a group together to hike those trails this last summer. I spent most of the summer about 2 miles from the trailhead, but I didn't have the critical mass of hikers-- thesis work was also in the way. I'll see about it again next summer. Hopefully my hiking partners will be available earlier in the season.

I haven't had a chance to try the r.walk module yet. It would be nice to have a side-by-side comparison of these two methods, and the type of results each returns.