Outils de modèle de gravité / Huff

Je cherche un moyen de simuler un modèle de gravité en utilisant une couche basée sur des points.

Tous mes points reçoivent une valeur z et plus cette valeur est élevée, plus leur "sphère d'influence" est grande. Cette influence est inversement proportionnelle à la distance au centre.

C'est un modèle typique de Huff, chaque point est un maximum local et les vallées entre elles indiquent les limites de la zone d'influence entre elles.

J'ai essayé plusieurs algorithmes d'Arcgis (IDW, allocation des coûts, interpolation polynomiale) et QGIS (plugin heatmap), mais je n'ai rien trouvé qui pourrait m'aider. J'ai également trouvé ce fil , mais il n'est pas très utile pour moi.

Comme alternative, je pourrais également être satisfait par un moyen de générer des diagrammes de Voronoi s'il existe un moyen d'influencer la taille de chaque cellule par la valeur z du point correspondant.

Voici une petite fonction python QGIS qui implémente cela. Il nécessite le plugin rasterlang (le référentiel doit être ajouté à QGIS manuellement).

Il attend trois paramètres obligatoires: la couche de points, une couche raster (pour déterminer la taille et la résolution de la sortie) et un nom de fichier pour la couche de sortie. Vous pouvez également fournir un argument facultatif pour déterminer l'exposant de la fonction de décroissance de la distance.

Les pondérations des points doivent figurer dans la première colonne d'attributs de la couche de points.

Le raster résultant est automatiquement ajouté au canevas.

Voici un exemple de la façon d'exécuter le script. Les points ont des poids compris entre 20 et 90, et la grille a une taille de 60 par 50 unités cartographiques.

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)