2

I want to create a least cost path program with Python and GDAL since I couldn't find any open source cost path additions for use in Python.

As recommended here the best approach is to use A* search algorithm using slope as the cost between generated nodes.

I want to read in my Slope Raster with GDAL. What is the best way to read in the Raster in order to further process it in my (not written yet) python A* search algorithm?

So far I read it in as a array but I am not sure if this is the best:

import gdal

# Open data
input_raster = "slope2.tif"
raster = gdal.Open(input_raster)

# read raster as array
bandraster = raster.GetRasterBand(1)
datamask = bandraster.ReadAsArray()
4
  • SAGA (open source) within the QGIS 2.0 processing toolbox has a least-cost-path function. It also allows to create cost rasters
    – Curlew
    Oct 2, 2013 at 22:49
  • Yeah I know that. But isn't it rather difficult to access SAGA via Python?
    – ustroetz
    Oct 2, 2013 at 22:51
  • What do you need python for in the first place? Both complex models and batch processing can be done in the modeler. But well, you could call saga modules via the processing runalg command like that processing.runalg("saga:leastcostpath", ...)
    – Curlew
    Oct 2, 2013 at 23:04
  • Well my entire script is written in Python. The least-cost-path is just a minor part of it. I am not familiar with the processing runalg command you suggest. Can I use it also in a Python stand alone script?
    – ustroetz
    Oct 2, 2013 at 23:18

2 Answers 2

4

NetworkX provides a ready-to-use library for the A* Algorithm.

Basically the steps you want to take are:

  1. Read the slope (the slope numbers are the weight, the more weight the less optimal)
  2. Create the graph from the slope matrix. This is the hardest part.
  3. Feed the NetworkX lib the graph and according to docs it should do the rest.

This is a canned solution imo. If you want to learn how the A* algorithm works you can create a T/F matrix, where T is the 'walkable' cells (with a slope less than 4?) and F the not walkables.

There are many worth reading tutorials on the subject for the A* around the intertubes that operate with that kind of arrays. Good learning material! (That will be my weekend project ;) )

1
  • Any recommendations for step 2?
    – ustroetz
    Oct 3, 2013 at 18:05
1

I think the best way to read in a raster for any purpose with Python/GDAL is by using a scanline and the unpack struct function. The code is more compact, the control is more effective and the execution time is faster than the one with 'ReadAsArray'. The scanline/struct method depends on fmttypes and their values can be supplied in a dictionary. In the next code I include a complete example of use to determine, by using the Python Console Editor of QGIS, the total average and the average by columns (only the first average value is printed) for values of a raster loaded in the Map View.

from osgeo import gdal
import struct
layer = iface.activeLayer()
provider = layer.dataProvider()
fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}
path= provider.dataSourceUri()
dataset = gdal.Open(path)
band = dataset.GetRasterBand(1)
totHeight = 0
totColumns = 0
BandType = gdal.GetDataTypeName(band.DataType)
column_means = []

for x in range(band.XSize):
    scanline = band.ReadRaster(x, 0, 1, band.YSize,1, band.YSize, band.DataType)
    values = struct.unpack(fmttypes[BandType] * band.YSize, scanline)

    for value in values:
        totHeight += value
        totColumns += value

    column_means.append(totColumns/float(band.YSize))
    totColumns = 0

average = totHeight / float((band.XSize * band.YSize))
print "Average = %0.5f" % average
print "First mean = %0.5f" % column_means[0]  
dataset = None

The results (Average = 1824.71801, First mean = 1685.04298) were reached in only one second for my raster of 791 rows x 1680 columns (1,328,880 pixels) and 'Int16' band type.

Your Answer

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

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