In-Class Exercise 11.2 Cost-Distance
NOTE: Due to requests to cite this method in publications I have written a short report based on this workshop and published it on eScholarship along with sample data at this stable URL.
Please cite the report available here
Tripcevich, N. 2009 "Tobler's walking model using Optimal Path as Line or Pathdistance in ArcGIS", ARF Practical Workshop series, Archaeological Research Facility, UC Berkeley.
Links to an external site.https://escholarship.org/uc/item/30c414w4
Links to an external site.
This workshop follows on Ex. 11 part 1 (Viewshed Analysis). Please download the data for that exercise.
Most cartographic distance calculations are performed in coordinate systems such as the metric Universal Transverse Mercator (UTM) where space is evenly divided into equal-sized units. In practice, however, people often conceive of distances in terms of time or effort. In mountainous settings the isotropic assumptions of metric map-distance calculations are particularly misleading because altitude gained and lost, and other obstacles such as cliff faces, are typically not factored into the travel estimate. Cost-distance functions are a means of representing such distances in terms of units of time or caloric effort.
In this exercise we will use an algorithm, one based primarily on slope from a topographic surface, to derive both an optimal route across that surface (Least-cost Path) as well as estimated time to cross that surface. Please see the general Cost-distance explanation Links to an external site. from ESRI as well as the manual on Path-distance Links to an external site. for more specifics on the feature we'll be using.
The ESRI Pathdistance function calculates movement with an 8 directional (Queen's move) algorithm and has the ability to incorporate customized anisotropy. This is a relatively sophisticated cost function as compared with other GIS environments (see I. Herzog's 2014 paper Links to an external site.comparing various platforms).
Here we are performing a Path Distance calculation where the cost function is modified by supplying a custom Vertical Factor Table Links to an external site. that is built from an empirically derived hiking function published by Tobler, 1993. Links to an external site. Tobler's function estimates time (hours) to cross a raster cell based on slope (rise/run) as calculated on a DEM.
I have converted data provided by Tobler into a structure that ArcGIS will accept as a Vertical Factor Table Links to an external site.for the Path Distance function. In other words, ArcGIS makes these sorts of calculations easy by allowing you to enter customized factors to the cost function.
2015 Update: the PathDistance custom vertical table function is currently broken in Arcmap 10.3 however it works as expected in the new ArcGIS Pro v1x application.
2017 Update: In ArcGIS Pro 2.0 the custom VF Table seems to be broken.
A. ISOTROPIC TIME ESTIMATE
1. Distance and time estimates based on 4 km/hr speed on an isotropic surface
This exercise begins by calculating distance and estimating time done in the old fashioned way, assuming a flat surface, as you would do on a map with a ruler.
After opening ArcGIS pro, you can import Ex11.1 map from the previous exercise, by selecting 'insert' in the upper portion of the screen, 'import map' and 'recent maps' or browse for the map you need.
We will begin by creating an isotropic cost surface and then generating 15 min isolines on this surface.
- If you don't already have the data used in the Viewshed exercise, download it here (10.6 Mb), unzip into Documents, and load the files into Arcmap. Also download, unzip, and load this Chivay Obsidian Source Point shapefile Download Chivay Obsidian Source Point shapefile showing the actual location of the Chivay Obsidian Source.
- Go to Toolbox > Spatial Analyst Tools > Distance > Euclidean Distance [Make sure Spatial Analyst extension is turned on under Tools > Extensions]
Input Feature: ChivaySource
Output Dist Raster: EucDist_Chivay
[click folder and save it in the same DEM folder]
Output cell size: 30
Leave the other fields blank.
- Before clicking OK click the Environments… button
- Zoom in your map pane so that the Chivay Source point and the Sites fill the display.
- Change Extent: <Default> to read "Same as Current Display"
2. Create raster showing isotropic distance in minutes
The output is a raster showing the metric distance from the site.
Let's assume for the purposes of this example that hiking speed is 4 km/hr regardless of the terrain, which works out to 1000m each 15 min.
- To transform the EucDist_shp grid into units of time (minutes) we can perform a Spatial Analyst menu > Raster Calculator expression using this information:
15 min / 1000 meters = 1 min / 66.66 meters
To construct that equation in Raster Calculator to get the raster to display in units of minutes we would divide the distance raster by 66.66. Built up the equation like so:
Double click EucDist_shp layer, click / button, then type in 66.66
If you've done it this way the output raster (Calculation) will range in value from
0 to 902.268
3. Convert to isochrones
- Create contours using the output raster with
Spatial Analyst toolbar > Surface Analysis > Contour…
- Then use the good Calculation raster as your input surface, choose 30 as your interval and use the defaults for the rest. Save the file out as "isochrone30" (click the folder to specify the target directory).
What do you think of the results of these lines?
The land is very steep in this area, especially in the cliff bands near the obsidian source. There must be a better way to estimate travel times!
B. ANISOTROPIC TIME ESTIMATE USING PATH DISTANCE
1. An algorithm is needed to convert from Slope to Time
In ArcGIS the PathDistance function can be used to conduct anisotropic distance calculations if we supply the function with a customized Vertical Factor table Links to an external site. with the degree slope in Column 1 and the appropriate vertical factor in Column 2. The vertical factor acts as a multiplier, and multiplying by a Vertical Factor of 1 has no effect. The vertical factor table we will use was generated in Excel and it reflects Tobler's Hiking Function by using the reciprocal of Tobler's function. This function was used in Excel:
TIME (HOURS) TO CROSS 1 METER or the reciprocal of Meters per Hour =0.000166666*(EXP(3.5*(ABS(TAN(RADIANS(slope_deg))+0.05))))
Thus for 0 deg slope which Tobler's function computes to 5.037km/hr, the result should be 1/5037.8 = 0.000198541
You can download the file that we will use for these values here.
Click ToblerAway.txt Links to an external site.and click the Download arrow, then move it into the same /cost folder. The reverse ToblerTowards.txt Links to an external site. is used for travel in the opposite direction, towards the source.
An abbreviated version of this ToblerAway file looks like this
|
Slope (deg) |
Vertical Factor |
|
-90 |
-1 |
|
-80 |
-1 |
|
-70 |
2.099409721 |
|
-50 |
0.009064613 |
|
-30 |
0.001055449 |
|
-10 |
0.00025934 |
|
-5 |
0.000190035 |
|
-4 |
0.000178706 |
|
-3 |
0.000168077 |
|
-2 |
0.000175699 |
|
-1 |
0.000186775 |
|
0 |
0.000198541 |
|
5 |
0.000269672 |
|
10 |
0.000368021 |
|
30 |
0.001497754 |
|
50 |
0.012863298 |
|
70 |
2.979204206 |
|
80 |
-1 |
|
90 |
-1 |
2. Pathdistance Calculation
Zoom your Data extent to 250,000 and pan the screen around until you can see both the Chivay Source point and the sites in the Study Area. If you are working on a laptop you may need to use a smaller scale like 1:300,000.
Go to ArcToolbox > Spatial Analyst Tools > Distance > Path Distance...
- Input: ChivaySource
- Output: <default> PathDis_shp1
- Input cost:<blank>
- Input surface: srtm30_utm (our 30m DEM)
- Output Backlink raster: backlink1
- Open the "Vertical factor parameters" field
- Input Vertical Raster: srtm30_utm
- Vertical Factor: Table (scroll down to bottom of menu list)
- If you haven't already stored the ToblerAway file click ToblerAway.txt Links to an external site.and click the Download arrow, then in Arcmap click Browse folder and browse find "Tobler_away.txt" on your hard drive
One more thing: Once again we'll restrict the processing extent
-
- click the Environments… tab
- Extent: Set it to "Same as Display"
- the first results include the Backlink Grid
- Turn off the Backlink Grid and reposition Hillshade.
Do you miss Arcmap?
*Note* that the previous step (use of the PathDistance function with a custom Vertical Factor table) is the reason we need to use ArcGIS Pro. If you prefer to return to Arcmap you may do so at this point. All the data layers generated into the default Geodatabase GDB by ArcGIS Pro should be readable by Arcmap.
To move back to Arcmap efficiently you might run both programs concurrently while you load the files into Arcmap, then close ArcGIS Pro so there aren't permissions errors.
3. Inspect your results carefully and convert units if necessary
Turn off the Backlink1 file so that you're viewing PathDist1_shp.
Why aren't the values larger? Recall what units the Tobler_Away function calculates?
We need to multiply this anisotropic result by 60 because…. It's in hours.
- Spatial Analyst Tools > Map Algebra > Raster Calculator…
PathDist1_shp * 60 ( note that you must have spaces between the variables and the operators in Raster Calculator.When in doubt, add a space ) - Right-click the Output Calculation … Call it "aniso_min" and save it into
\Ex11_surfaces\cost\
- Now generate 30 min contours on this layer as you did earlier and call the output "aniso30".
- Turn off all layers except "srtm30_utm", "Hydro", "isochrone30" and "aniso30".
- Do you see the difference in the two cost surfaces?
In this image the contours are Isochrones incrementing by 30 minutes, while the background raster shows the Anisotropic distance calculation with red lines showing 30 min increments away from the source.
Using the Info tool (the blue 'i') you can point-sample distances along the way against the layer aniso1_min and isodist_min.
The differences between the two surfaces can be interesting
- Spatial Analyst > Raster Calculator...
- [aniso_min] - [isodist_min]
Symbolize the output raster using a Red to Blue gradient fill and look at the result. What are the high and low values showing?
PART C. GENERATING LEAST-COST PATHS BETWEEN SITES
Recall that you generated a Backlink grid called aniso1_bklk as part of the PathDistance calculation. We will now use that to calculate least-cost paths between the two sites and other sites in the region.
First, we need to generate a
- Add Backlink Grid to Arcmap called "aniso1_bklk".
Look at the results and consider how Backlink Grids work (read about the Direction raster in the Help on Cost Distances Links to an external site.), do the values 0 through 8 make sense?
Input Feature: Sites
Destination Field: FID
Input Cost Distance Raster: aniso_hours (use the hours not the minutes one)
Input Cost Backlink Raster: aniso_bklk
Note that the output is a raster GRID file. Wouldn't these be more useful as Polylines?
- Use Conversion > From Raster > Raster to Polyline...
Input Raster: <the name of your cost paths file>
Field: Value
Geometry type: Polyline
Zoom in and explore the results with only the following layers visible:
Callalli_DEM, Hydro_lines100k, aniso_minutes isolines, and the cost-paths you just generated. Note how the paths avoid steep hills and choose to follow valleys, as you might do if you were planning to walk that distance.
Produce a map with the cost isochrones and the paths showing. You'll also want some text explaining the isochrone intervals and the fact that you used an anisotropic cost function based on Tobler 1993.
Click "Share" tab and choose Export Map... to produce a PDF.