11

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?

4
  • 3
    For analytical questions, perhaps better use a GIS rather than a data format abstraction library...
    – markusN
    Jul 1, 2012 at 8:12
  • 2
    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, 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

2 Answers 2

14

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: enter image description here

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)
    array = band.ReadAsArray()
    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)
4
  • 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.
    – Michael
    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?
    – ustroetz
    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.
    – Michael
    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
3

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:

A Star Animated

See A* Search Algorithm (Wiki) and Python A* Search Algorithm (SO)

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.

1
  • 2
    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, 2013 at 21:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.