Geographic Information Systems Asked by Tamil selvan on April 4, 2021
I’am trying find least cost path from a GeoTIFF file using Python. I have found the following solution from Google:
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 pixel2coords(arrayData, rasterfn):
coords = []
raster = gdal.Open(rasterfn)
xoff, a, b, yoff, d, e = raster.GetGeoTransform()
array = np.array(arrayData)
Orginx = array[1]
Orginy = array[0]
for x,y in zip(Orginx,Orginy):
xp = a * x + b * y + a * 0.5 + b * 0.5 + xoff
yp = d * x + e * y + d * 0.5 + e * 0.5 + yoff
coords.append((xp,yp))
return coords
def utmToLatlng(utmCoordinates):
latLngs = []
p = pyproj.Proj(proj='utm', zone="39U", ellps='WGS84')
for i in utmCoordinates:
lat, lng = p(i[0],i[1],inverse=True)
latLngs.append([lat, lng])
return latLngs
def main(CostSurfacefn,startCoord,stopCoord):
costSurfaceArray = raster2array(CostSurfacefn) # creates array from cost surface raster
indicesArray = findPath(CostSurfacefn,costSurfaceArray,startCoord,stopCoord) # creates path array
utmCoordinates = pixel2coords(indicesArray, CostSurfacefn)
latlng = utmToLatlng(utmCoordinates)
print(latlng)
if __name__ == "__main__":
CostSurfacefn = 'CostSurface.tif'
startCoord = (345387.871,1267855.277)
stopCoord = (345479.425,1267799.626)
main(CostSurfacefn,startCoord,stopCoord)
But in that above code they converting pixel value from float to integer.
xOffset = int((x - originX)/pixelWidth)
yOffset = int((y - originY)/pixelHeight)
Because route through array function accepts integer’s only.So when i convert pixel to coordinate, it doesn’t showing correct lat long.
def pixel2coords(arrayData, rasterfn):
coords = []
raster = gdal.Open(rasterfn)
xoff, a, b, yoff, d, e = raster.GetGeoTransform()
array = np.array(arrayData)
Orginx = array[1]
Orginy = array[0]
for x,y in zip(Orginx,Orginy):
xp = a * x + b * y + a * 0.5 + b * 0.5 + xoff
yp = d * x + e * y + d * 0.5 + e * 0.5 + yoff
coords.append((xp,yp))
return coords
I need to find least cost path between two points. Can anyone suggest any better solution for this problem?
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP