Alat model gravitasi / Huff

Saya mencari cara untuk mensimulasikan model gravitasi menggunakan lapisan berbasis titik.

Semua poin saya diberi nilai z dan semakin tinggi nilai ini, semakin besar "lingkup pengaruh" mereka. Pengaruh ini berbanding terbalik dengan jarak ke pusat.

Ini adalah model khas Huff, setiap titik adalah maksimum lokal dan lembah di antara mereka menunjukkan batas zona pengaruh di antara mereka.

Saya mencoba beberapa algoritma dari Arcgis (IDW, alokasi biaya, interpolasi polinomial) dan QGIS (plugin heatmap), tetapi saya tidak menemukan apa pun yang dapat membantu saya. Saya juga menemukan utas ini , tetapi tidak terlalu membantu bagi saya.

Sebagai alternatif, saya juga bisa puas dengan cara menghasilkan diagram Voronoi jika ada cara untuk mempengaruhi ukuran setiap sel dengan nilai-z dari titik yang sesuai.

Inilah sedikit fungsi python QGIS yang mengimplementasikan ini. Ini membutuhkan plugin rasterlang (repositori harus ditambahkan ke QGIS secara manual).

Ia mengharapkan tiga parameter wajib: Lapisan poin, lapisan raster (untuk menentukan ukuran dan resolusi output), dan nama file untuk lapisan output. Anda juga dapat memberikan argumen opsional untuk menentukan eksponen fungsi peluruhan jarak.

Bobot untuk poin harus di kolom atribut pertama dari layer poin.

Raster yang dihasilkan secara otomatis ditambahkan ke kanvas.

Berikut ini contoh cara menjalankan skrip. Poin memiliki bobot antara 20 dan 90, dan grid berukuran 60 kali 50 unit peta.

points = qgis.utils.iface.mapCanvas().layer(0)
raster = qgis.utils.iface.mapCanvas().layer(1)
huff(points,raster,"output.tiff",2)

from rasterlang.layers import layerAsArray
from rasterlang.layers import writeGeoTiff
import numpy as np

def huff(points, raster, outputfile, decay=1):
    if points.type() != QgsMapLayer.VectorLayer:
        print "Error: First argument is not a vector layer (but it has to be)"
        return
    if raster.type() != QgsMapLayer.RasterLayer:
        print "Error: Second argument is not a raster layer (but it has to be)"
        return
    b = layerAsArray(raster)
    e = raster.extent()
    provider = points.dataProvider()
    extent = [e.xMinimum(),e.yMinimum(),e.xMaximum(),e.yMaximum()]
    xcols = np.size(layerAsArray(raster),1)
    ycols = np.size(layerAsArray(raster),0)
    xvec = np.linspace(extent[0], extent[2], xcols, endpoint=False)
    xvec = xvec + (xvec[1]-xvec[0])/2
    yvec = np.linspace(extent[3], extent[1], ycols, endpoint=False)
    yvec = yvec + (yvec[1]-yvec[0])/2
    coordArray = np.meshgrid(xvec,yvec)
    gravity = b
    point = QgsFeature()
    provider.select( provider.attributeIndexes() )
    while provider.nextFeature(point):
      coord = point.geometry().asPoint()
      weight = point.attributeMap()[0].toFloat()[0]
      curGravity = weight * ( (coordArray[0]-coord[0])**2 + (coordArray[1]-coord[1])**2)**(-decay/2)
      gravity = np.dstack((gravity, curGravity))
    gravitySum = np.sum(gravity,2)
    huff = np.max(gravity,2)/gravitySum
    np.shape(huff) 
    writeGeoTiff(huff,extent,outputfile)
    rlayer = QgsRasterLayer(outputfile)
    QgsMapLayerRegistry.instance().addMapLayer(rlayer)