# GDAL - Perform Simple Least Cost Path Analysis

I am investigating methods to perform a simple least cost path analysis with gdal. By simple, I mean using the slope of a dem as the only cost factor.

I would prefer to do using the python or .net bindings, but will take anything. Can anyone suggest any good tutorials or the like?

-
For analytical questions, perhaps better use a GIS rather than a data format abstraction library... – markusN Jul 1 '12 at 8:12
Out of curiosity, what is the application? It is difficult to think of anything for which the DEM's slope would be a realistic proxy for the cost of travel. Are you sure this is what you need? It would be a pity if, after going to the effort of writing this code, you discovered it did not actually solve your problem! – whuber Jul 1 '12 at 17:42
Slope could be useful as a travel cost if you're modelling a gravity-dependent dispersal model of some sort, though I'd expect some other factors too rather than just slope. – MappaGnosis Jul 2 '12 at 8:45
Also, slope usually shows the maximum slope at each cell, even if the route is not travelling directly downhill or uphill. – Matthew Snape Jul 25 '12 at 20:05

You can use the A* search algorithm using slope as the cost between generated nodes. To see a quick visualization of what that looks like:

to understand A*.

For a slope map there are options out there - Here is one.

With a slope map (raster) you can get cost values out of it with GDAL.

-
Can you explain how to make the slope raster to a graph so it can be used in the Python A* Search Algorithm Code you pointed out? I know how to get the value out of the raster with GDAL. but as what shall I store it in order to use it as a graph (e.g. Dictionary?)? – ustroetz Oct 1 '13 at 21:44

The following script performs a least cost path analysis. Input parameters are a cost surface raster (e.g. slope) and start and stop coordinates. A raster with the created path is returned. It requires the skimage library and GDAL.

For example the least cost path between point 1 and point 2 is created based on a slope raster:

``````import gdal, osr
from skimage.graph import route_through_array
import numpy as np

def raster2array(rasterfn):
raster = gdal.Open(rasterfn)
band = raster.GetRasterBand(1)
return array

def coord2pixelOffset(rasterfn,x,y):
raster = gdal.Open(rasterfn)
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
xOffset = int((x - originX)/pixelWidth)
yOffset = int((y - originY)/pixelHeight)
return xOffset,yOffset

def createPath(CostSurfacefn,costSurfaceArray,startCoord,stopCoord):

# coordinates to array index
startCoordX = startCoord[0]
startCoordY = startCoord[1]
startIndexX,startIndexY = coord2pixelOffset(CostSurfacefn,startCoordX,startCoordY)

stopCoordX = stopCoord[0]
stopCoordY = stopCoord[1]
stopIndexX,stopIndexY = coord2pixelOffset(CostSurfacefn,stopCoordX,stopCoordY)

# create path
indices, weight = route_through_array(costSurfaceArray, (startIndexY,startIndexX), (stopIndexY,stopIndexX),geometric=True,fully_connected=True)
indices = np.array(indices).T
path = np.zeros_like(costSurfaceArray)
path[indices[0], indices[1]] = 1
return path

def array2raster(newRasterfn,rasterfn,array):
raster = gdal.Open(rasterfn)
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = array.shape[1]
rows = array.shape[0]

driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, gdal.GDT_Byte)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()

def main(CostSurfacefn,outputPathfn,startCoord,stopCoord):

costSurfaceArray = raster2array(CostSurfacefn) # creates array from cost surface raster

pathArray = createPath(CostSurfacefn,costSurfaceArray,startCoord,stopCoord) # creates path array

array2raster(outputPathfn,CostSurfacefn,pathArray) # converts path array to raster

if __name__ == "__main__":
CostSurfacefn = 'CostSurface.tif'
startCoord = (345387.871,1267855.277)
stopCoord = (345479.425,1267799.626)
outputPathfn = 'Path.tif'
main(CostSurfacefn,outputPathfn,startCoord,stopCoord)
``````
-