Herramientas del modelo Gravity / Huff

26

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 thread , pero no es ' t muy útil para mí.

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

    
pregunta Damien 29.01.2013 - 16:00

1 respuesta

13

Aquí hay una pequeña función python de QGIS que implementa esto. Requiere el plugin rasterlang (el repositorio debe ser añadido 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 la distancia.

Los pesos de 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 es de 60 por 50 unidades de mapa en tamaño.

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)     
respondido por el Jake 29.01.2013 - 18:51

Lea otras preguntas en las etiquetas