8

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?

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

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: 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)
  • 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 '17 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 '17 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 '17 at 7:48
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.

  • 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 '13 at 21:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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