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)
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.