Python/GDAL: Best way to read in raster for use in A* search algortihm

Question

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()

Answer 1

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 ;) )

Answer 2

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.