TransWikia.com

Finding Least Cost Path from GeoTIFF file using Python

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?

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP