Narzędzia modeli grawitacyjnych / Huff

26

Szukam sposobu na symulację modelu grawitacyjnego przy użyciu warstwy punktowej.

Wszystkim moim punktom przypisano wartość Z i im wyższa ta wartość, tym większa jest ich „sfera wpływów”. Wpływ ten jest odwrotnie proporcjonalny do odległości od centrum.

Jest to typowy model Huffa, każdy punkt jest lokalnym maksimum, a doliny między nimi wskazują granice strefy wpływów między nimi.

Wypróbowałem kilka algorytmów Arcgis (IDW, alokacja kosztów, interpolacja wielomianowa) i QGIS (wtyczka Heatmap), ale nie znalazłem nic, co mogłoby mi pomóc. Znalazłem też ten wątek , ale nie jest to dla mnie zbyt pomocne.

Alternatywnie, mógłbym również zaspokoić sposób generowania diagramów Voronoi, jeśli istnieje sposób, aby wpłynąć na wielkość każdej komórki przez wartość z odpowiedniego punktu.

Damien
źródło

Odpowiedzi:

13

Oto mała funkcja python QGIS, która to implementuje. Wymaga wtyczki rasterlang (repozytorium należy dodać do QGIS ręcznie).

Oczekuje trzech obowiązkowych parametrów: warstwy punktowej, warstwy rastrowej (w celu określenia rozmiaru i rozdzielczości danych wyjściowych) oraz nazwy pliku dla warstwy wyjściowej. Możesz także podać opcjonalny argument, aby określić wykładnik funkcji zanikania odległości.

Wagi dla punktów muszą znajdować się w pierwszej kolumnie atrybutów warstwy punktów.

Powstały raster jest automatycznie dodawany do obszaru roboczego.

Oto przykład uruchomienia skryptu. Punkty mają wagę od 20 do 90, a siatka ma rozmiar 60 na 50 jednostek mapy.

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)
Jake
źródło
3
(+1) Podejście wygląda dobrze. Ale dlaczego bierzesz pierwiastek kwadratowy, a następnie przekształcasz go ponownie w obliczeniach curGravity? To strata czasu obliczeniowego. Kolejny zmarnowany zestaw obliczeń obejmuje normalizację wszystkich siatek „grawitacyjnych” przed znalezieniem maksimum: zamiast tego znajdź ich maksimum i znormalizuj je o sumę.
whuber
Czy to nie obciąża całej frakcji?
lynxlynxlynx
1
Jake, wciąż nie potrzebujesz pierwiastka kwadratowego: po prostu całkowicie o nim zapomnij i użyj połowy zamierzonego wykładnika. Innymi słowy, jeśli z jest sumą kwadratów różnic współrzędnych, zamiast obliczania (sqrt (z)) ^ p, czyli dwóch umiarkowanie kosztownych operacji, wystarczy obliczyć z ^ (p / 2), który (ponieważ p / 2 jest liczbą obliczoną wstępnie ) jest tylko jedną operacją rastrową - i prowadzi również do wyraźniejszego kodu. Pomysł ten pojawia się w przypadku zastosowania modeli grawitacyjnych w pierwotnym zamierzeniu: do czasów podróży. Nie ma już wzoru na pierwiastek kwadratowy, więc zwiększasz czas podróży do mocy -p / 2.
whuber
Wielkie dzięki, wygląda na to, czego potrzebuję. Po prostu problem, nie jestem przyzwyczajony do pythona i nigdy nie korzystałem z rozszerzenia Rasterlang. Zainstalowałem go w mojej wersji QGIS, ale utknąłem z „błędem składni”. Czy twoja funkcja jest już zaimplementowana w rozszerzeniu rasterlang? Jeśli nie, jak to zrobić? Dzięki za pomoc! http://i.imgur.com/NhiAe9p.png
Damien
1
@Jake: Ok, chyba zaczynam rozumieć, jak działa konsola. Zrobiłem tak, jak powiedziałeś, a kod wydaje się być właściwie zrozumiany. Teraz mam inny błąd związany z pakietem Pythona „shape_base.py”. Czy w mojej instalacji QGIS brakuje niektórych funkcji? http://i.imgur.com/TT0i2Cl.png
Damien