Herramientas del modelo Gravity/Huff

Question

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 Huff típico, cada punto es un máximo local y los valles entre ellos indican los límites de la zona de influencia entre ellos.

He probado varios algoritmos de Arcgis (IDW, asignación de costes, interpolación polinómica) y de QGIS (plugin de mapas térmicos), pero no he encontrado nada que me ayude. También encontré esto rosca pero no es muy útil para mí.

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

Answer 1

Aquí hay una pequeña función de QGIS python que implementa esto. Requiere el plugin rasterlang (el repositorio tiene que ser añadido a QGIS manualmente).

Espera tres parámetros obligatorios: La capa de puntos, una capa rasterizada (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 descomposición de la distancia.

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

La trama resultante se añade 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)