Herramientas de modelo Gravity / Huff

Estoy buscando una manera de simular un modelo de gravedad utilizando una capa basada en puntos.

A todos mis puntos se les asigna un valor z y cuanto mayor sea este valor, mayor será su "esfera de influencia". Esta influencia es inversamente proporcional a la distancia al centro.

Es un modelo típico de Huff, cada punto es un máximo local y los valles entre ellos indican los límites de la zona de influencia entre ellos.

Probé varios algoritmos de Arcgis (IDW, asignación de costos, interpolación polinómica) y QGIS (complemento de mapa de calor), pero no encontré nada que pudiera ayudarme. También encontré este hilo , pero no es muy útil para mí.

Como alternativa, también podría estar satisfecho con una forma de generar diagramas de Voronoi si hay una manera de influir en el tamaño de cada celda por el valor z del punto correspondiente.

Aquí hay una pequeña función de Python QGIS que implementa esto. Requiere el complemento rasterlang (el repositorio debe agregarse a QGIS manualmente).

Espera tres parámetros obligatorios: la capa de puntos, una capa ráster (para determinar el tamaño y la resolución de la salida) y un nombre de archivo para la capa de salida. También puede proporcionar un argumento opcional para determinar el exponente de la función de disminución de distancia.

Los pesos para los puntos deben estar en la primera columna de atributos de la capa de puntos.

El ráster resultante se agrega automáticamente al lienzo.

Aquí hay un ejemplo de cómo ejecutar el script. Los puntos tienen pesos entre 20 y 90, y la cuadrícula tiene un tamaño de 60 por 50 unidades de mapa.

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)