﻿ Building Blocks of Spatial Analysis > Distance Operations > Cost distance

# Cost distance

The term cost distance is widely used within GIS and in this context has a dual meaning: (i) the notion of an alternative family of distance metrics; and (ii) a procedure for determining least cost paths across continuous surfaces, typically using grid representations

Within GIS the term is distinct from the computation of varying costs or times of travel across a network. The term cost in this connection is also a generic idea, not necessarily implying financial cost, but some composite measure that varies (i.e. is not constant) across the study region, and which needs to be taken into account when computing paths and distances.

An example serves to clarify this notion. Suppose that we wish to construct a section of road between two points, A(0,2) and B(4,2), where costs are constant over the region of interest. In this case the distance is simply 4 units and the path involved is a straight line in the plane (Figure 4‑57, bold line). Construction costs will simply be 4k units, where k is some constant (the cost per unit length of path construction). Now suppose that construction costs also depend on the cost of land acquisition, and that land costs vary in our study region in a simple manner — greater to the “north” of our region than to the “south”. In this situation it is cheaper to construct our road with a curve in it, taking advantage of the fact that land costs are lower to the south (nearer the x-axis in our model).

Figure 4‑57 Cost distance model

If costs increase in a precisely linear manner, the optimal (least cost) path will have a smooth curve form of the kind shown by the thinner line (a shape like a chain hanging between two posts). This case turns out to be solvable mathematically, as are cases where costs vary in radially symmetric manner (e.g. around a city center), whereas almost every other example of varying costs cannot be solved precisely and we need to resort to computational procedures to find the least cost path.

Accumulated cost surfaces and least cost paths

Within GIS packages solutions to such problems are found using a family of algorithms known as Accumulated Cost Surface (ACS) methods — for a fuller discussion of ACS methods and representational accuracy see Douglas (1994) and Eastman (1989). The ACS procedure is applied to raster datasets, in which the primary input surface is a complete grid of generalized costs, i.e. every cell is assigned an absolute or relative cost measure, where costs must be (positive) ratio-scaled variables. In Figure 4‑58 we show a hypothetical cost grid for the same example as shown in Figure 4‑57. Values increase linearly away from the x-axis, from 0.5 to 2.5. This essentially a much simplified representation of the function f(x,y)=y, so f(x,y)=0 when y=0 for all x and increases steadily so at y=2 f(x,y)=2. In practice cost functions should rarely if ever have zero-valued sections, so this example should be considered illustrative!

Providing a gridded representation of the cost surface highlights the importance of using a sufficiently fine representation of the surface, since the solution path must be computed by adding up the costs along alternative paths. If the grid resolution is quite coarse, as in Figure 4‑58, the solution path will not have the curved form of the optimal path, but will simply be a very blocky stepped line of cells approximating the curve. Amending the grid resolution (halving the cell size to 0.25 units) as shown in Figure 4‑59A improves both the cost surface representation and the possible solution paths that are considered.

Figure 4‑58 Cost surface as grid

Figure 4‑59 Grid resolution and cost distance

If we look at a path in Figure 4‑58 that starts at A(0,2) and goes diagonally down towards C(2,0) and then back up to B(4,2) we have cell costs of 2+1.5+1+.5+.5+1+1.5+2=10 units (times cell size of 0.5 gives a cost distance of 5 units). But this path is longer than the sum of the steps — all the diagonal steps are √2 (1.414… units) longer than horizontal and vertical steps, so the total cost must be increased for such steps by this factor, giving a final cost distance of 6.85 units. This is less than the cost of 8 units that would be incurred for a straight line path from A to B. But can we be sure this is the best route possible for this problem and is the calculation “correct”?

ACS provides a systematic procedure for answering this question. With this algorithm we start at A and sum up the costs in all directions, step by step, choosing the least cost increment at each step. Eventually, when we arrive at B, we will have the overall least cost of getting to B, and if we have kept track of our route (the steps we have taken) we will have the least cost path from A to B. Thus ACS is a form of spread algorithm, with values increasing in all directions depending on the grid resolution, cost surface and algorithmic implementation. Almost all GIS implementations of ACS use immediate neighbors in calculating values (i.e. queen’s moves, to the 8 adjacent cells) but some facilitate calculation to neighbors further away. GRASS, for example, offers the option in the r.cost function, of using knight’s move neighbors (one step up and one diagonal) in its calculations. Distance transform methods utilize a much wider range of neighbors, and hence are not restricted to 8-directional cost surfaces and path alignments, and thus do not exhibit the common directional bias and computational errors of ACS distance models.

As we start to consider the ACS idea more closely we see that a number of additional issues need to be addressed. The first is that in a GIS-enabled grid structure cells (or more specifically cell centers) are frequently used as references rather than grid lines. Thus our starting point A cannot be (0,2) exactly but must be a cell center and our end points will also all be cell centers. For example, see Figure 4‑59B where A(0,2) is approximated on the grid by the point A(0.125, 1.875). For the best results we should estimate the value of the cost function at cell centers if possible, rather than at cell edges or grid lines. In the example we are using we have ignored this latter factor, which will often be the case since finer detail for within cell estimation may not be available (although could be obtained by interpolation).

A second issue is that if we take grid cell centers as our reference point, the start point will be half way across the first cell, and the next point will be half way across an adjacent cell, so it can be argued that the cell values (costs) should be allocated 50:50 over this one step, and likewise for steps 2, 3 etc. This can be seen in Figure 4‑59B where the starting point for A is shown as a selected cell center and the step to A′ in the diagonally adjacent cell lies 50% in a cell with cost 2 and 50% in a cell with cost 1.75. This 50:50 division of costs is the normal approach used in GIS packages. The averaged cost from A to A is then (2+1.75)/2=3.75/2 and thus the incremental cost distance is √2*0.25*3.75/2 units (adjusting for the diagonal path using √2 and for the cell size of 0.25 units in this example).

We now have all the basic elements we need in order to compute accumulated costs and identify the least cost path from A to B. Using our example cost grid, with a cell size of 0.25 units, we need to carry out the following operations using a typical GIS:

 • Create a source grid that contains a 0 in the cell in which our starting point A lies, and some other value (e.g. ‑1, 1, 9999) for all other grid cells. This can be done with a text editor and converted to a GIS compatible grid format if necessary
 • Create a cost grid in the same manner, exactly the same size and the same resolution as the source grid, but containing the costs associated with each cell position (i.e. similar to the cells in Figure 4‑59a)
 • Run the cost distance facility within the GIS on these two input grids, requesting the output as a cost distance (ACS) grid, and also requesting (if necessary) that a direction or tracking dataset is generated that can be subsequently used to identify least cost paths
 • Finally, to create the least cost paths, a target grid is required, which contains the target or destination points (again as 0s or a similar cell identifier for the target points). The least cost path facility in the GIS is then run using the target grid, ACS grid and tracking grid

A sample result for our test problem is shown in Figure 4‑60. Each grid cell has been separately colored and a series of additional destination points (cells) have been included in the least cost path finding stage. The least cost path from A to B, corresponding to our original problem, hugs the bottom of the diagram as it seeks to take advantage of the low cost first row (cells with a value of only 0.25 units). Routes to other cells may share common paths for parts of their route and may not be unique, as indicated by the two routes to the cell diagonally below point A to the right. If the cost surface had been flatter, for example f(x,y)=y/4+1, the ACS model with the cell resolution illustrated would fail to select the correct path, but would show a horizontal path from A to B. To obtain a more accurate estimate of the correct path and cost distance a much finer grid would be required. For this latter example the true path is a gentle curve with a minimum value at around y=1.65 and total cost of 5.88 units.

Figure 4‑60 Accumulated cost surface and least cost paths

This simple example provides a flavor of cost distance computation, but this procedure (or family of procedures) has far more to offer. The first and perhaps most important point to note is that the source and target features do not need to be single points. They can be multiple distinct points (cells), lines (linear groups of cells), or regions (clusters of adjacent cells). In all such cases the ACS procedure computes the lowest total distance to a target cell from its nearest source cell. If requested the GIS may also produce an allocation raster, which identifies which source each cell in the grid is closest to (typically by giving each cell a unique integer identifier). Obstacles (e.g. lakes, cliffs etc.) may be included in these models by setting cells in the cost raster to very high values or by the GIS explicitly supporting unavailable (masked) cells (e.g. in Idrisi, coding barrier cells as -1).

The input cost raster may be, and often is, a composite grid created from a number of other (normalized) input grids by simple map algebraic operations. The final cost grid may also be normalized to provide an indication for each cell of the relative cost, difficulty or friction of crossing that cell. Figure 4‑61 illustrates the result of such a process, which has been applied to determining possible alternative routes for a main road in part of South-East England. The map shows the existing road (the A259T) from Hastings in the south-west to Brenzett in the north-east. Three alternative routes, identified here by arrows, were generated by ACS methods. These routes were generated from cost surfaces that included different weighted values for a series of inputs: land-use (grade and type of land); SSSIs (Sites of special scientific interest); areas of designated Ancient Woodland; the slope of the land across the study area (derived from a DEM); existing (buffered) built-up areas (many of which are conservation areas) and ancient monuments; National Trust land; and existing (buffered) road and rail routes. The final route alternatives were then computed from the composite ACS surface and smoothed to reduce the effects of the gridding of all the datasets and output paths. This information was then provided as input to a public consultation exercise.

This is a relatively simple example of a general class of problems categorized as corridor planning ― selection of alternative routes for roads, rail lines, power lines or other purposes frequently requires the involvement of many parties, and formal multi-criteria decision support systems, such as AHP and ANP, have been extensively used in this field ― see, for example, Azis (1990), Bailey (2003) and Piantanakulchai (2005).

Figure 4‑61 Alternative route selection by ACS

The above examples are based on the procedures utilized in ArcGIS, Spatial Analyst. Many other packages offer similar facilities. For example, Idrisi offers two alternative algorithms: one, called the pushbroom algorithm, is similar to the distance transform procedures described in Section 4.4.2, Cost distance, and is fast but not suitable for use with more complex problems (for example truly implementing barriers); a second Idrisi procedure, costgrow, is much slower but does facilitate complex friction patterns and fixed barriers. Idrisi also allows for anisotropy in the cost modeling process, with some directions (e.g. downhill) being preferable to others (e.g. uphill). SAGA, GRASS and PCRaster provide similar functions, with inputs optionally including a friction surface (with relative friction values and absolute barriers) and a physical surface, with movements restricted to uphill, downhill and either stopping at or crossing flat areas. Such procedures are clearly targeted at hydrological and related studies, so whilst being well suited to these environments are not readily generalized. Many such programs, particularly those with an environmental sciences and remote sensing background, perform the path finding stage using a form of steepest descent procedure (often described as a drainage procedure) applied to the ACS. This approach is not guaranteed to provide the correct shortest path — tracked paths should be used. This is illustrated in Figure 4‑62, where local Euclidean grid distances and 50:50 splitting of costs between the cells shaded in gray has been used to generate the ACS grid, Figure 4‑62B. The steepest path from the target cell back to the source cell (Figure 4‑62C) is not the least cost path, as its total accumulated cost is 13.726. The correct least cost path is shown in Figure 4‑62D. ACS procedures are subject to a number of other problems, for example: becoming trapped in localized plateaus or pits; inadequate handling of multiple equivalent paths; or failing when input and/or output datasets are specified to operate in integer rather than float mode.

A number of GIS packages incorporate extensions to the standard cost distance/ACS model. For example, ArcGIS provides a variant of the ACS procedure which it describes as Path Distance. This is essentially the same as ACS but includes options such as inclusion of surface distance (using a surface raster or DEM), and vertical and horizontal movement rasters that are utilized in adjusting the optimum path computation. TNTMips also includes a path analysis procedure that takes into account the physical surface form, but the authors are cautious about its optimality.

Figure 4‑62 Steepest path vs tracked path

 A. Cost surface B. Accumulated cost surface (ACS) 1 5 5 5 1 4.414 7.414 12.414 8.828 5.828 1 5 5 5 1 3.414 6.656 11.656 8.07 4.828 1 5 5 5 1 2.414 5.656 6.656 5.414 3.828 1 5 5 1 1 1.414 3 4.242 2.414 3.414 1 1 1 1 1 1 0 1 2 3 C. ACS Steepest descent path D. ACS Least cost path

Idrisi provides a similar facility, Disperse, designed to provide a simple force-based model of anisotropic object dispersal (SAGA provides comparable facilities). In this case the inputs are a set of four grid files: (i) the location of the target and source features (non-zero entries are targets, zeros are sources); (ii) and (iii) two grids, one containing the magnitude and the other containing the direction (azimuth) of the anisotropic force field, i.e. a vector field modeled using a pair of grids; and (iv) an output grid. The notion of the magnitude of the force field in this implementation is as the inverse of a friction field, F, force=1/F where F in this case is taken to be greater than 1. The anisotropic dispersal model is either user specified or utilizes an inbuilt cosine function that modifies the stated friction field:

F'=F*f, where

α is an angular difference between the angle being considered and the stated friction angle, and k is a parameter (k=2 is the default). The process may be modeled as a series of steps (time periods) by pre-specifying a maximum dispersal distance at each stage.

Distance transforms

In the image processing world a procedure similar to ACS has been developed in connection with work on pattern recognition. The procedure is known as a Distance Transform (DT), and is typically implemented on binary images in a scanning manner rather than the spreading approach of ACS. Extremely efficient algorithms have been developed to generate exact Euclidean DTs and excellent approximations to exact Euclidean DTs (EDTs) on large images. Some of these procedures have even been implemented within hardware chipsets for real-time applications. GIS vendors do not generally utilize DT methods, but they are included in many image processing toolsets and packages, including the MATLab Image Processing Toolkit: e.g. see the bwdist() function, which provides a variety of metrics for black and white images, including City Block and exact Euclidean. A second example with more direct links with GIS is the product suite, ERDAS.

Using DTs, distances over the image or grid are calculated in an incremental manner based entirely on the distance to nearby (not necessarily adjacent) cells. The standard algorithm involves a two-pass scan of a square or rectangular lattice: a forward scan from top left (NW corner) to bottom right (SE corner), and then a backwards scan from bottom right to top left (i.e. SENW corner). The scan commences by placing the forward scan mask over the first available cell in the grid to be transformed – in a 3x3 scan mask (MxM cells, where M=3) this will be grid cell (row 2, col 2). A basic DT 3x3 mask is illustrated in Figure 4‑63, where 5 values are used based on the two-value (3,4) integer mask. (3,4) is useful as it is extremely fast to process, and dividing through by 3 at the end yields direct steps of 1 and diagonal steps of 1.33... giving a reasonable approximation to local Euclidean measure (1.414...) and within 6% of exact global Euclidean measure.

Figure 4‑63 3x3 Distance transformation – scan elements

 4 3 4 4 3 4 3 0 0 3 3 0 3 4 3 4 4 3 4

Forwards     Backwards         Combined

Typically the underlying grid to be transformed is an array in which unreached cells are initially assigned a large number (e.g. 999) and target or source cells are assigned a value of 0. Cells that act as barriers or restricted areas (e.g. buildings) can be coded with even larger numbers (e.g. 9999), in which case they are never altered. On completion of the two-pass scan each cell in the resulting lattice will contain the distance to the nearest point in the set of source points. If barriers are included they must be in blocks of at least (M‑1)x(M‑1) cells, and accessible areas should consist of at least (M+1)x(M+1) cells. In addition, two scan pairs: (NWSE/SENW) and (SWNE/NESW), with one or more iterations, are required to ensure full scanning of the grid. On completion of the scan each cell in the resulting lattice will contain the distance to the nearest point in the set of source points. The pushbroom algorithm described by Eastman (1989) and used in Idrisi for their cost distance modeling is very similar to this approach.

The central function in this algorithm is of the form: d0=min{d+D(i),d0}, where d0 is the current value at the central point (0) of the template, D(i) is the local distance to the ith element of the template, and d is the current value at (r,c). The basic algorithm involves O(Mn2) computations where n is the maximum dimension of the grid. Details of DT procedures are provided in the references cited, and code samples in MATLab and Python are available from: http://www.desmith.com/MJdS/DT1.htm.

Table 4‑10 3x3 Chamfer metrics

 Case Description Adjacency matrix 1 Distances are determined by the L1 or ‘city-block’ metric and paths correspond to the “rook’s move” in chess parlance. Maximum error: 29.29% 2 1 2 1 0 1 2 1 2 2 Distances are defined as case 1 for any immediate neighbor, i.e. horizontal or diagonal (“rook’s move” and “bishop’s move”). Maximum error: 41.41% 1 1 1 1 0 1 1 1 1 3 As per 2 but with diagonal distances determined by the Euclidean metric applied locally — sometimes referred to as Octagonal or Local Euclidean metric. Maximum error: 7.61% √2 1 √2 1 0 1 √2 1 √2 4 Chamfer (3,4)/3 metric. This uses integer values to provide a better estimate of Euclidean distance than Cases 1 or 2. Maximum error: 6.07% 4 3 4 3 0 3 4 3 4 5 Fractional chamfer (Borgefors, 1986) – non-integer values for both rook’s move and bishop’s move directions (revised after Butt and Maragos, 1998). Maximum error: 3.96% 1.36039 0.96194 1.36039 0.96194 0 0.96194 1.36039 0.96194 1.36039

A selection of 3x3 masks is shown in Table 4‑10. These are known as chamfer metrics because the locus of the metric generates a figure similar to a cross-section of a piece of wood with chamfered or beveled edges. The generic form of such masks is defined by the two unique mask values (a,b), e.g. (2,1) or (3,4). The maximum error percentages shown in Table 4‑10 provide the global maximum errors for the metric as compared with global Euclidean measure. The use of 3x3 local metrics is fast and simple, but results in directional bias in the resulting DT. These show the familiar 4- or 8-direction patterns seen in many GIS processes that rely on 3x3 neighborhoods (see, for example, Figure 4‑34).

This problem can be considerably reduced by extending the neighborhood to a 5x5 mask (M=5), in which case 3 values define the mask: (a,b,c) as shown in Figure 4‑64. Amongst the best integer values for this mask size are the triple (5,7,11), dividing all cell values by 5 at the end to give (1,1.4,2.2). The best fractional values are (0.9866, 1.4141, 2.2062). Use of the integer values above ensures that every grid cell distance computed (i.e. global values) is within 2% of the exact Euclidean distance and within 1.36% for the fractional masks. This compares with global errors of up to 7.6% with a pure local Euclidean metric.

 2b c 2a c 2b c b a b c 2a a 0 a 2a c b a b c 2b c 2a c 2b

Figure 4‑65A provides the simplest test case for the DT procedure. This shows a uniform region with a single target point in the center (originally coded as 0, with all other grid cells coded 999). This figure shows the best non-exact DT possible using a fractional-valued 5x5 scanning template. The exact Euclidean distance transform (EDT) of this uniform space yields a perfectly circular set of filled contours around the central point. The vector equivalent would be a multi-ring buffer operation. When applied to multiple points (e.g. as part of a point-pattern analysis exercise) the set of computed distances, F(), is known as the empty space function.

If a simple barrier is introduced, e.g. representing a river with a bridge point, shortest paths appear diffracted through the gap (Figure 4‑65B). In this example the target point has been set as (2,50). With a complex barrier, as in Figure 4‑65C, a single extra iteration has been used to ensure full convergence. In all cases shortest paths to target points can be drawn as orthogonal to the distance bands, but are best constructed using tracking of the selected steps, as described further in the examples below. Although the DT approach may seem somewhat divorced from practical GIS analysis, DT procedures have been extended and developed by: Ratti et al. (architecture and the built environment; 2001, 2004); de Smith (transport routes; 2004, 2006), Ikonen and Toivanen (image processing; 2005); and Kristinsson et al. (geothermal pipelines, 2005) to provide a new and fast family of ACS-like procedures for solving many grid-related distance problems.

Figure 4‑65A Distance transform, single point. 5x5 DT, no constraints, target (50,50)

Figure 4‑65B Distance transform, single point 5x5 DT, barrier – target (2,50)

Figure 4‑65C Distance transform, single point. 5x5 DT, complex obstacle – target (50,90)

We limit discussion here to a small sample of applications. The first deals with radially symmetric models of traffic in cities ― see Angel and Hyman, Urban Fields (1976, pp 15-16 and 159) for a full discussion of the models involved. The second example examines pedestrian movement within urban areas, whilst the third set of examples address questions of optimal road and pipeline location under a range of constraints.

Figure 4‑66 shows the result of applying a least ‘cost’ distance transform (LCDT) to a city with traffic speeds symmetrically varying from almost 0kph in the center (where the crossed lines meet) to much higher levels (e.g. 50kph+) at the city edge. The cost field in this case is traffic speed (i.e. a velocity field) and the curved lines show contours of equal travel time (isochrones) to or from a point due South of the city center.

Figure 4‑66 Urban traffic modeling

The central function in the DT procedure in this instance is of the form:

d0 = min{d+D(i)*COST(r,c), d0}

where d0 is the current value at the central point (0) of the mask, D(i) is the local distance to the ith element of the mask, d is the current value at (r,c) and COST(r,c) is the (averaged) cost function at row r, column c. As can be seen, the velocity field distorts the distance bands. Shortest (quickest) paths across this transformed surface may be plotted by following routes to or from the source point in any given direction, remaining at right angles (orthogonal) to the equal time contours, or by tracking of the computed quickest paths. Hence the quickest route from the start point due South of the center to locations North of the city center would curve around (East or West) the center to avoid congestion (low traffic speeds. In this particular example these curve form log spirals (Wardrop, 1969).Angel and Hyman applied models of this type to the analysis of traffic flows in Manchester, UK. Distance transform methods enable such models to be calibrated against known analytical results and then applied to more realistic complex velocity and cost fields, for example to examine the possible effects of introducing central zone congestion charging.

A second example of this procedure is shown in Figure 4‑67. This shows a section of the Notting Hill area of central London, with the annual Carnival route shown in white and areas in dark gray being masked out (buildings etc.). Colors indicate increasing distances (in meters) from the Carnival route along nearby streets.

Figure 4‑67 Notting Hill Carnival routes

Although this example has been created by a DT scanning procedure the result is sometimes described as a “burning fire” analogy, with the fire spreading from the source locations (in this case a route) outwards along adjacent streets and squares (for true fire spread models see http://www.fire.org and Stratton, 2004). The same techniques can be applied in public open spaces and even within building complexes such as shopping malls, galleries and airports. This illustrates the ability of DT methods to operate in what might be described as complex friction environments.

Distance transforms are very flexible in their ability to incorporate problem-specific constraints. For example, in a typical ACS procedure variable slopes are modeled by regarding steeper slopes as higher cost, rather than relating the spreading process to slope directly (other than some uphill/downhill constrained implementations). The selection of a route for a road, railway or pipeline is affected by slope, but primarily in the direction of the path rather than the maximum slope of the physical surface (although this does have a bearing on construction cost).

Figure 4‑68 Alternative routes selected by gradient constrained DT

Figure 4‑68 illustrates the application of a gradient-constrained DT (GCDT) to the selection of road routes through a hilly landscape. The most northerly two routes (upper right in the diagram) were constrained to a 1:12 and 1:11 gradient (before engineering grading) whilst the most southerly two represent routes found using less restrictive 1:9 and 1:8 gradient constraints.

With minor modifications the GCDT procedure described above has been utilized with a high resolution (10meter) DEM to locate alternative pipeline routes for the Hellisheiði 80MW geothermal power station in Iceland (Kristinsson, Jonsson and Jonsdottir, 2005). In this case additional constraints were applied giving the maximum and minimum heights permitted for the route and differential directional constraints that permitted limited flow uphill, which is possible for these kinds of two-phase pipelines (Figure 4‑69A,B).

Figure 4‑69 Hellisheiði power plant pipeline route selection

 A. Overview of power plant area B. Alternative solution paths A, B1, B2 and C for Well 3