# 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... Jul 1, 2012 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! Jul 1, 2012 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. Jul 2, 2012 at 8:45
• Also, slope usually shows the maximum slope at each cell, even if the route is not travelling directly downhill or uphill. Jul 25, 2012 at 20:05

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.

The following was copied from a very useful website for gdal related stuff.

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)
``````
• I like your answer. How do you deal with e.g. lakes where the cost value is the same for a larger area. My path goes through a lake and sort of meanders around like a snake until the area is covered before it continues as expected? Thanks. Jul 8, 2017 at 8:16
• I haven't worked on this for a long time. You probably thought of this already, but I would just set the cost for the lake really high. This way the path should avoid the lake, shouldn't it? Jul 10, 2017 at 9:48
• Yeah, I set the lake to be just a bit over 0, that way there is a cost and the meandering disappears. Jul 13, 2017 at 7:48
• In general, one should drop a reference to code they copy/paste pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html Oct 18, 2022 at 20:35

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?)? Oct 1, 2013 at 21:44