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..