Gravity / Huff-Modellwerkzeuge

26

Ich suche nach einer Möglichkeit, ein Gravitationsmodell mithilfe einer punktbasierten Ebene zu simulieren.

Allen meinen Punkten ist ein Z-Wert zugeordnet, und je höher dieser Wert ist, desto größer ist ihr "Einflussbereich". Dieser Einfluss ist umgekehrt proportional zum Abstand zum Zentrum.

Es ist ein typisches Huff-Modell, jeder Punkt ist ein lokales Maximum und Täler zwischen ihnen geben die Grenzen der Einflusszone zwischen ihnen an.

Ich habe verschiedene Algorithmen von Arcgis (IDW, Kostenverteilung, Polynominterpolation) und QGIS (Heatmap-Plugin) ausprobiert, aber ich habe nichts gefunden, was mir helfen könnte. Ich habe diesen Thread auch gefunden , aber er ist für mich nicht sehr hilfreich.

Alternativ könnte ich mich auch damit zufrieden geben, Voronoi-Diagramme zu erstellen, wenn es eine Möglichkeit gibt, die Größe jeder Zelle durch den z-Wert des entsprechenden Punkts zu beeinflussen.

Damien
quelle

Antworten:

13

Hier ist eine kleine QGIS-Python-Funktion, die dies implementiert. Es benötigt das Rasterlang Plugin (das Repository muss manuell zu QGIS hinzugefügt werden).

Es werden drei obligatorische Parameter erwartet: Der Punkt-Layer, ein Raster-Layer (um die Größe und Auflösung der Ausgabe zu bestimmen) und ein Dateiname für den Ausgabe-Layer. Sie können auch ein optionales Argument angeben, um den Exponenten der Abstandsfunktion zu bestimmen.

Die Gewichte für die Punkte müssen sich in der ersten Attributspalte des Punkte-Layers befinden.

Das resultierende Raster wird automatisch zur Zeichenfläche hinzugefügt.

Hier ist ein Beispiel, wie das Skript ausgeführt wird. Die Punkte haben eine Gewichtung zwischen 20 und 90 und das Raster ist 60 mal 50 Karteneinheiten groß.

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)
Jake
quelle
3
(+1) Der Ansatz sieht gut aus. Aber warum nimmst du die Quadratwurzel und quadrierst sie dann beim Rechnen neu curGravity? Das ist eine Verschwendung von Rechenzeit. Eine weitere Verschwendung von Berechnungen besteht darin, alle "Schwerkraft" -Gitter zu normalisieren, bevor das Maximum ermittelt wird: Ermitteln Sie stattdessen ihr Maximum und normalisieren Sie das durch die Summe.
whuber
Quadriert das nicht die ganze Fraktion?
Lynxlynxlynx
1
Jake, du brauchst die Quadratwurzel immer noch nicht: Vergiss sie einfach ganz und benutze die Hälfte des vorgesehenen Exponenten. Mit anderen Worten, wenn z die Summe der Quadrate der Koordinatendifferenzen ist, anstatt (sqrt (z)) ^ p zu berechnen, was zwei mäßig teure Operationen sind, berechne einfach z ^ (p / 2), was (weil p / 2) ist eine vorberechnete Zahl ) ist nur eine Rasteroperation - und führt auch zu einem klareren Code. Diese Idee kommt zum Ausdruck, wenn Sie Gravitationsmodelle anwenden, wie sie ursprünglich vorgesehen waren: um die Reisezeit zu bestimmen. Da es keine Quadratwurzelformel mehr gibt, erhöhen Sie die Laufzeit auf die Potenz -p / 2.
whuber
Vielen Dank, das sieht so aus, als ob ich es brauche. Nur ein Problem, ich bin nicht so an Python gewöhnt und ich habe nie die Rasterlang-Erweiterung verwendet. Ich habe es auf meiner QGIS-Version installiert, aber ich habe einen "Syntaxfehler". Ist Ihre Funktion bereits in der rasterlang-Erweiterung implementiert? Wenn nein, wie mache ich das? Danke für deine Hilfe! http://i.imgur.com/NhiAe9p.png
Damien
1
@Jake: Ok, ich glaube ich fange an zu verstehen, wie die Konsole funktioniert. Ich habe getan, wie Sie gesagt haben, und der Code scheint richtig verstanden zu werden. Jetzt habe ich einen anderen Fehler, der sich auf ein Python-Paket "shape_base.py" bezieht. Fehlen meiner QGIS-Installation einige Funktionen? http://i.imgur.com/TT0i2Cl.png
Damien