PostGIS Raster Summation (Kartenalgebra)

Ich habe eine Tabelle mit Polygonen, die die Isochronen der Reisezeit an bestimmten Tagen darstellen. Für jeden Ursprungspunkt gibt es fünf Isochronengeometrien (in separaten Zeilen gespeichert). Für jeden Ursprungspunkt möchte ich die fünf Isochronen (binär NULL oder 1) rastern und sie dann zu einer einzelnen Rasterebene kombinieren. Diese Rasterebene erfordert eine einfache Kartenalgebra: sum / 5, sodass jeder Ursprung am Ende einer einzelnen Rasterebene zugeordnet wird, die Werte in [NULL, 0,2, 0,4, 0,6, 0,8, 1] aufweist, je nachdem, wie viele davon vorhanden sind die konstituierenden Schichten überlappen sich. Es ist eine Wahrscheinlichkeitsfläche.

Meine Daten sind alle in Postgres 9.3 (mit PostGIS) gespeichert. Mein Problem ist, dass ich zwar lernen möchte, PostGIS-Raster zu verwenden, es aber eine sehr steile Lernkurve zu haben scheint und alle Beispiele, die ich finden kann, sich mit einer einzelnen Rasterebene befassen. In den Beispielen wird diese Ebene als Teil einer Polygon-Überlagerung verwendet, wobei möglicherweise der Rasterwert für jedes Polygon gemittelt wird. Ich habe kein reproduzierbares Beispiel für die Kombination gefunden: a) Vektor -> Raster b) Kartenalgebra; und c) GROUP BY-Attribut gemäß meinem ersten Absatz.

Ich kann GDAL oder GRASS verwenden, wenn ich diese Aufgabe ausführen muss, aber dies scheint etwas zu sein, mit dem PostGIS umgehen kann. Dies wäre praktisch, da meine Eingabedaten bereits PostGIS-Geometrie sind. und ich möchte mich unbedingt mit PostGIS-Raster auseinandersetzen.

Einige Beispieldatenstruktur:

areaid    time        date          isogeom (polygon)
1000      07:15:00    2014-05-05    xxx
1000      07:15:00    2014-05-06    xxy
...
1006      07:15:00    2014-05-05    zzz

Ich möchte rastern, nach Gebiet gruppieren und dann die Kartenalgebra ausführen, um Folgendes zu erreichen:

areaid    isorast (raster)
1000      aaa
1006      bbb

Ich bin nicht erfolgreich gewesen, dieses zu PostGIS zu enthalten. Mein Ansatz war es, den Vektor in Raster umzuwandeln, die Raster in Arrays zu speichern und die Kombination mit Numpy-Arrays über psycopg2 durchzuführen, bevor ich sie in ein GeoTIFF schreibe (um sie möglicherweise wieder in PostGIS einzufügen). Nicht ideal, aber machbar.

Sie müssen Ihre eigene Aggregatfunktion schreiben:

CREATE OR REPLACE function sum_raster_state(raster, raster)
returns raster
language sql
as $f$

SELECT
CASE $1
WHEN NULL THEN
      $2
ELSE
    ST_MapAlgebra(
        $1, 
        $2, 
        '([rast1] + [rast2])', 
        NULL, 
        'UNION', 
        '[rast2]', 
        '[rast1]', 
        NULL)
END;
$f$;

CREATE OR REPLACE FUNCTION sum_raster_final(raster)
returns raster
language sql
as $f$
SELECT 
   $1
$f$;

create aggregate sum_raster(raster) (
    SFUNC = sum_raster_state,
    STYPE = raster,
    FINALFUNC = sum_raster_final
);

danach kannst du es so nennen

SELECT areaid, sum_raster(st_asraster(isogeom, ...)) FROM isochrone GROUP BY areaid

Auf diese Weise erhalten Sie die Summe aller Raster mit derselben Gebiets-ID. Sie müssen die Rasterwerte weiterhin durch die Anzahl der Beobachtungen für jede Bereichs-ID dividieren. (Ich habe es nicht in die Aggregatfunktion aufgenommen. Sie können es entweder hier oder später mit MapAlegbra tun.)

Stellen Sie sicher, dass alle Eingabe-Raster ausgerichtet sind, da dies sonst nicht funktioniert.