Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
2  
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

2 Answers 2

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.

share|improve this answer
1  
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: 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)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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