Week 9 - Anisotropic Cost surfaces

Introduction: Non-Cartesian Distance Units
Spatial analysis in anthropology is frequently based on Cartesian coordinates that describe space in terms of regular units on a grid coordinate system. The UTM system is a good example of a Cartesian coordinate system in metric units. GIS permits the calculation of distance costs in terms of non Cartesian units that, some have argued, are a better reflection of real perspectives and human (agent) decision-making.

Consider your commute to the University today. You may calculate it in terms of linear distance units like miles or kilometers, but you likely arrived here by any number of ways: bike, bus, walking, or private car. Furthermore, there are added complexities to some modes of travel, such as catching the proper bus or finding parking for your car.

In terms of time budgeting, you might be better off considering the commute to the university in terms of time units (hours) rather than linear distance units. In terms of a general model we could assign some costs to certain features: wading across the river perhaps takes 5x longer per meter than crossing land. Bushwacking through chaparral takes 10x longer. These types of costs can be implemented systematically in the GIS in order to convert distance measures in Cartesian space (e.g., UTM) into time. More explicit ecological approaches might convert costs into calories.
Clearly these kinds of cost-surfaces will involve many assumptions about transportation and impedance.

In this week's exercise we will perform a cost distance calculation using a customized cost function that reflects empirical data published by Tobler, 1993. Tobler's function estimates time (hours) to cross a raster cell based on slope (degrees) as calculated on a DEM.

I have converted data provided by Tobler into a structure that ArcGIS will accept as a "Vertical Factor Table" 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.

First, skim over the ArcGIS 9 description of the PathDistance function.

PathDistance is an elaboration of Cost Distance algoritm, as discussed in class.

PART A. ISOTROPIC TIME ESTIMATION

1. Distance and time estimates based on 4 km/hr speed on an isotropic surface

This exercise begins by calculating distance and estimating time the old fashioned way, as you would do on a map with a ruler.

We will begin by creating an isotropic cost surface and then generating 15 min isolines on this surface.

• From the Callalli.zip file load "All_ArchID", "Callalli_DEM", and "hydro_lines100k" into Arcmap

• Open the attribute table for All_ArchID and choose the rows with ArchID = 675 and ArchID = 985.

Hint: Sort Ascending ArchID column, and hold down ctrl button while you click the left-most column in the table for those two rows.

• Go to Toolbox > Spatial Analyst Tools > Distance > Euclidean Distance

Input Feature: All_ArchID
Output Dist Raster: Iso_dist1
<click folder and save it in the same DEM folder>
Leave the other fields blank.

• Before clicking OK click the Environments… button at the bottom of the window.

• Click the General Settings… tab and look at the different settings. Change the Output Extent field to read "Same as layer callalli_dem".

2. Create raster showing isotropic distance in minutes

The output is a raster showing the metric distance from the two sites.
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 isodist1 grid into units of time (minutes) we can perform a Spatial Analyst > Raster Calculator expression using this information:
  15 min / 1000 meters = 1 min / 66.66 meters
  

How would you construct that equation in Raster Calculator to get the raster to display in units of minutes?
If you've done it right the output raster (Calculation) will range in value from
0 to 112.614
Hint: divide the raster by 66.66

• When you have this raster right click and choose "Make Permanent…", and call it "isodist_min". Remove the Calculation grid and add this "isodist_min" grid to your project.

3. Convert to isochrones

• Create contours using the output raster with
  Spatial Analyst > Surface Analysis > Contour…

• Then use the good Calculation raster as your input surface, choose 15 as your interval and use the defaults for the rest. Save the file out as "isochrone15".

What do you think of the results of these lines?

This photo was taken facing south towards Taukamayo which lies in the background on the right side of the creek in the center-left of the photo. Taukamayo is the western-most of the two sites you just worked with in Arcmap. You can see how steep the hills are behind the site and on the left side of the photo you can see the town of Callalli for scale. Clearly the land is very steep around the site, especially in those cliff bands on the south end of the DEM. There are also cliffbands on the north side of the river but the river terrace in the foreground where the archaeologists are standing is quite flat.

PART B. ANISOTROPIC COST SURFACE

We will use ArcWorkstation GRID program for this calculation because the custom (directionally dependent) slope calculation is not supported directly in the Arcmap CostDistance or Pathdistance functions.
ArcGRID dates to before the Shapefile and Geodatabase formats were used, so we must work with GRID rasters in this program.

Step 1: Create the source_grid

• Once again, open the attribute table for All_ArchID and choose the rows with ArchID = 675 and ArchID = 985.
• Set the extent of the output grid
  Spatial Analyst > Options > Extent tab > Analysis extent: Same as layer "callalli_dem"

Convert the two points to a GRID layer

• Spatial Analyst > Convert… > Features to Raster…

Input: All_ArchID
Field: ARCHID
Output (default)
Output raster: click folder and save into /DEM/ folder. Call it "sites"

Note: For today's exercise you can not have spaces anywhere in the path to these files.

Turn off the other layers and look closely at the two sites in the raster you just created. Do you see a square at each site?

Step 2: Calculate PathDistance in ArcWorkstation GRID

• Start ArcGRID
• Start Menu > Programs > ArcGIS > ArcInfo Workstation > GRID
• Click on the GRID window with the prompt that reads "Grid: " to make it active.
• Type "Help" (leave out the quotes)
• Then in the Help window click the Index tab and type "Pathdistance"

You will see a description of the function with a breakdown of how it is issued. It looks like this with a variety of required and optional arguments.

PATHDISTANCE(<source_grid>, {cost_grid}, {surface_grid}, {horiz_factor_grid}, {horiz_factor_parm}, {vert_factor_grid}, {vert_factor_parm}, {o_backlink_grid}, {o_allocate_grid}, {max_distance}, {value_grid})

• Scroll down through this page. As you can see, it is a comprehensive description, and it has a lot of similarities to the help description webpage we saw earlier. Close the Help box.
  

The formula used by PATHDISTANCE to calculate the total cost from cell a to cell b is:
Cost_distance = Surface_distance * Vertical_factor * (Friction(a) * Horizontal_factor(a) + Friction(b) * Horizontal_factor(b)) /2)

In ArcInfo GRID the PATHDISTANCE function can be used to conduct anisotropic distance calculations if we supply the function with a customized Vertical Factor table with the degree slope in Column 1 and the appropriate vertical factor in Column 2. The vertical factor acts as a multiplier where 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.

Right-click ToblerAway.txt, choose"Save Link as..." and save it in the same folder as your GRID files (the reciprocal ToblerTowards.txt is used for travel in the opposite dirction, towards the source).

An abbreviated version of this ToblerAway file looks like this

• In Windows browse to the folder that holds the Callalli_DEM, Site GRIDs, and the ToblerAway file.
• Look at the PATH to this directory (confirm that there are no spaces and only alphanumeric characters)

• Set the Working directory in GRID

• At the GRID: prompt type
• &wo <complete path to DEM>

on my computer this command reads: &wo c:\gis_data\anth197\dem\

• If it is successful, then type "lg" at the prompt.
  You should see a list of the grid contents of the working directory.
  

• Issue this command at the GRID: prompt
  Pathdistance
• You'll see an example of how to strcuture this command
• Following that example, type the following at the GRID: prompt
  aniso_hours = Pathdistance (source_grid , # ,DEM_grid, #, #, DEM_grid, "TABLE ToblerAway.txt", aniso_bklk)
  

Substituting "sites" for source_grid and "callalli_dem" for DEM_grid in the above function. Make sure you type it exactly! If you have an error, use the up arrow key to re-type the previous command and the left/right arrow keys to move around. Include the # sign. # tells GRID to use the default values.

• When the function runs successfully (i.e., no error), return to Arcmap.

• Add the resulting grid to Arcmap (aniso_hours). Does the file range from 0 to 2.5?

Why aren't the values larger? Recall what units the Tobler_Away function calculates?

We need to multiply this result by 60 because…. It's in hours.

• Spatial Analyst > Raster Calculator…
  [aniso_hours] * 60
• Right-click the Output Calculation and Make Permanent… Call it "aniso_min"

• Now generate 15 min contours on this layer as you did earlier and call the output "aniso15".
• Turn off all layers except Callalli_DEM", "Hydro", "isochrone15" and "aniso15".
• Do you see the difference in the two cost surfaces?
  

Using the Info tool (the blue 'i') you can point-sample distances along the way.

The differences between the two surfaces are interesting

• Spatial Analyst > Raster Calculator...
• [aniso_min] - [isodist_min]

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 Aniso_bklk as part of your ArcGRID work. We will now use that to calculate least-cost paths between the two sites and other sites in the region.

• Add Backlink Grid to Arcmap called "aniso_bklk".

Look at the results. Recall how Backlink Grids work (read the Help on COST SURFACES), do the values 0 through 8 make sense?

Choose four sites from ArchID for this calculation
588, 611, 811, 835

• How would you construct this selection in ArchID from the Attribute Table using Options > Select by Attributes…?

It should look something like this
[ARCHID] = 588 OR [ARCHID] =611 OR [ARCHID] = 811 OR [ARCHID] = 835

• Now close the attribute table. You have four sites selected around the perimeter of the study area?

• Choose Toolbox > Spatial Analyst Tools > Distance > Cost Path…

Input Feature: All_ArchID
Destination Field: ArchID
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 Spatial Analyst > Convert > Raster to feature….

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.

The Lab 4 assignment (due next week at the beginning of class) involves conducting these analyses on your own dataset.

To quit GRID type "quit" at the prompt.