Гравитация / Хафф модели инструментов

26

Я ищу способ имитации гравитационной модели с использованием точечного слоя.

Всем моим точкам присваивается z-значение, и чем выше это значение, тем больше их «сфера влияния». Это влияние обратно пропорционально расстоянию до центра.

Это типичная модель Хаффа, каждая точка является локальным максимумом, а долины между ними указывают границы зоны влияния между ними.

Я попробовал несколько алгоритмов от Arcgis (IDW, распределение затрат, полиномиальная интерполяция) и QGIS (плагин heatmap), но я не нашел ничего, что могло бы мне помочь. Я также нашел эту ветку , но она не очень полезна для меня.

В качестве альтернативы, я также мог бы быть удовлетворен способом генерирования диаграмм Вороного, если есть способ влиять на размер каждой ячейки по z-значению соответствующей точки.

Damien
источник

Ответы:

13

Вот небольшая функция Python QGIS, которая реализует это. Для этого требуется плагин rasterlang (репозиторий необходимо добавить в QGIS вручную).

Он ожидает три обязательных параметра: слой точек, растровый слой (для определения размера и разрешения выходных данных) и имя файла для выходного слоя. Вы также можете указать необязательный аргумент для определения степени функции затухания расстояния.

Веса для точек должны быть в первом столбце атрибутов слоя точек.

Полученный растр автоматически добавляется на холст.

Вот пример того, как запустить скрипт. Точки имеют вес от 20 до 90, а размер сетки составляет 60 на 50 единиц карты.

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)
Джейк
источник
3
(+1) Подход выглядит хорошо. Но почему вы берете квадратный корень, а затем заново вычисляете его в вычислениях curGravity? Это пустая трата вычислительного времени. Еще один потерянный набор вычислений включает в себя нормализацию всех «гравитационных» сеток перед нахождением максимума: вместо этого найдите их максимум и нормализуйте их по сумме.
whuber
Разве это не квадрат всей фракции?
lynxlynxlynx
1
Джейк, тебе все еще никогда не нужен квадратный корень: просто забудь об этом совсем и используй половину предполагаемого показателя. Другими словами, если z является суммой квадратов разностей координат, вместо вычисления (sqrt (z)) ^ p, то есть двух умеренно дорогих операций, просто вычислите z ^ (p / 2), которое (потому что p / 2 является предварительно вычисленным числом ) - это всего лишь одна растровая операция, что также приводит к более четкому коду. Эта идея выходит на передний план, когда вы применяете гравитационные модели так, как они изначально предназначались: для путешествий во времени. Больше нет формулы квадратного корня, поэтому вы увеличиваете время прохождения до степени -p / 2.
whuber
Большое спасибо, это выглядит как то, что мне нужно. Просто проблема, я не так привык к Python и никогда не использовал расширение Rasterlang. Я установил его на мою версию QGIS, но я застрял с «синтаксической ошибкой». Ваша функция уже реализована в расширении rasterlang? Если нет, как мне это сделать? Спасибо за помощь! http://i.imgur.com/NhiAe9p.png
Дэмиен
1
@ Джейк: Ладно, думаю, я начинаю понимать, как работает консоль. Я сделал, как вы сказали, и код, кажется, правильно поняли. Теперь у меня есть еще одна ошибка, связанная с пакетом Python "shape_base.py". В моей установке QGIS отсутствуют некоторые функции? http://i.imgur.com/TT0i2Cl.png
Дэмиен