qgis - Gravity/Huff model tools

I am looking for a way to simulate a gravity model using a point-based layer.

All my points are assigned a z-value and the higher this value, the larger is their "sphere of influence". This influence is inversely proportional to the distance to the center.

It is a typical Huff model, each point is a local maximum and valleys between them indicate the limits of the zone of influence between them.

I tried several algorithms from Arcgis (IDW, cost allocation, polynomial interpolation) and QGIS (heatmap plugin), but I found nothing that could help me. I also found this thread, but it isn't very helpful for me.

As an alternative, I could also be satisfied by an way to generate Voronoi diagrams if there is a way to influence the size of each cell by the z-value of the corresponding point.

Answer

Here's a little QGIS python function that implements this. It requires the rasterlang plugin (the repository has to be added to QGIS manually).

It expects three mandatory parameters: The points layer, a raster layer (to determine the size and resolution of the output), and a filename for the output layer. You can also provide an optional argument to determine the exponent of the distance decay function.

The weights for the points need to be in the first attribute column of the points layer.

The resulting raster is automatically added to the canvas.

Here's an example of how to run the script. The points have weights between 20 and 90, and the grid is 60 by 50 map units in size.

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)