GDAL - Lakukan Analisis Jalur Minimal Biaya Terkecil

Saya sedang menyelidiki metode untuk melakukan analisis jalur biaya paling sederhana dengan gdal. Secara sederhana, maksud saya menggunakan kemiringan dem sebagai satu-satunya faktor biaya.

Saya lebih suka melakukan menggunakan binding python atau .net, tetapi akan mengambil apa pun. Adakah yang bisa menyarankan tutorial yang baik atau sejenisnya?

Script berikut melakukan analisis jalur biaya terendah. Parameter input adalah raster permukaan biaya (misalnya kemiringan) dan mulai dan berhenti koordinat. Raster dengan path yang dibuat dikembalikan. Ini membutuhkan perpustakaan skimage dan GDAL.

Misalnya jalur biaya terendah antara titik 1 dan titik 2 dibuat berdasarkan slaster raster:

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 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 = 'CostSurface.tif'
    startCoord = (345387.871,1267855.277)
    stopCoord = (345479.425,1267799.626)
    outputPathfn = 'Path.tif'
    main(CostSurfacefn,outputPathfn,startCoord,stopCoord)

Anda dapat menggunakan algoritma pencarian A * menggunakan slope sebagai biaya antara node yang dihasilkan. Untuk melihat visualisasi cepat seperti apa itu:

Lihat A * Cari Algoritma (Wiki) dan Python A * Algoritma Pencarian (SO)

untuk memahami A *.

Untuk peta kemiringan ada opsi di luar sana - Ini ada satu.

Dengan peta kemiringan (raster) Anda bisa mendapatkan nilai biaya dengan GDAL.