3

I am running the following code which I found on: GDAL - Perform Simple Least Cost Path Analysis which I adapted to my case, but I get the following error when running the code in the python console of QGIS. I am rather new to the programming & QGIS world so if somebody would be able to help me that would be wonderful!

import gdal, osr
import osgeo.gdal
osgeo.gdal.GetDriverByName
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] # East/West location of Upper Left corner
    originY = geotransform[3] # North/South location of Upper Left corner
    pixelWidth = geotransform[1] # X pixel size
    pixelHeight = geotransform[5] # Y pixel size
    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] # East/West location of Upper Left corner
    originY = geotransform[3] # North/South location of Upper Left corner
    pixelWidth = geotransform[1] # X pixel size
    pixelHeight = geotransform[5] # Y pixel size
    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 = "Administratief_perseel_1polygoon_PXLKOST_6000_6000.tif"
    startCoord = (174084.004,176969.786)
    stopCoord = (173697.916,175206.172)
    outputPathfn = 'Path.tif'
    main(CostSurfacefn,outputPathfn,startCoord,stopCoord)

Error Code:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/Users/PeterVanvoorden/Documents/GroepT/Thesis/c)QGIS/CLIP_AGIV_GRBgis_e40179/Projects/python_console_scripts/least_cost_path0.1.py", line 77, in <module>
    main(CostSurfacefn,outputPathfn,startCoord,stopCoord)
  File "/Users/PeterVanvoorden/Documents/GroepT/Thesis/c)QGIS/CLIP_AGIV_GRBgis_e40179/Projects/python_console_scripts/least_cost_path0.1.py", line 65, in main
    costSurfaceArray = raster2array(CostSurfacefn) # creates array from cost surface raster
  File "/Users/PeterVanvoorden/Documents/GroepT/Thesis/c)QGIS/CLIP_AGIV_GRBgis_e40179/Projects/python_console_scripts/least_cost_path0.1.py", line 9, in raster2array
    raster = gdal.Open(rasterfn)
AttributeError: 'module' object has no attribute 'Open'

I assume there is a problem with opening the file that I want to analyze "Administratief_perseel_1polygoon_PXLKOST_6000_6000.tif" but I just don't know how to resolve this..

2
  • I think that you need to import gdal inside your def.
    – radouxju
    Mar 14, 2014 at 10:23
  • @radouxju - no, importing gdal in the module (global) namespace is fine.
    – user2856
    Mar 15, 2014 at 6:09

1 Answer 1

5

import gdal is deprecated, gdal, ogr, osr etc are now part of the osgeo package. Try changing this:

import gdal, osr
import osgeo.gdal
osgeo.gdal.GetDriverByName

To this:

from osgeo import gdal, osr
4
  • Than I'm getting the following error (replaced path by 2 dots because it was to long): Traceback (most recent call last): File "<input>", line 1, in <module> File "/Users/../least_cost_path0.2.py", line 75, in <module> main(CostSurfacefn,outputPathfn,startCoord,stopCoord) File "/Users/../least_cost_path0.2.py", line 63, in main costSurfaceArray = raster2array(CostSurfacefn) # creates array from cost surface raster File "/Users/../least_cost_path0.2.py", line 8, in raster2array band = raster.GetRasterBand(1) AttributeError: 'NoneType' object has no attribute 'GetRasterBand' Mar 16, 2014 at 9:23
  • 1
    @PeterVanvoorden - That's because gdal couldn't open your raster and just returned None instead of raising an exception. Force gdal to fail by adding gdal.UseExceptions() at the beginning of the script (after the imports) so you can see what the error is. Also, make sure your raster path is correct, particularly if using relative paths. You should mark this question as answered and start a new one if you are still having trouble as your original issue is resolved
    – user2856
    Mar 17, 2014 at 4:22
  • Also use from osgeo import gdal, osr, ogr
    – pyaj
    Feb 19 at 14:31
  • Usually, yes. But the OP wasn't using ogr.
    – user2856
    Feb 19 at 19:55

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.