Podsumowanie rastrów PostGIS (algebra map)

Mam tabelę wielokątów reprezentujących izochrony czasu podróży w określone dni. Dla każdego punktu początkowego istnieje pięć geometrii izochronicznych (przechowywanych w osobnych rzędach). Dla każdego punktu początkowego chcę zrasteryzować pięć izochronów (binarna wartość NULL lub 1), a następnie połączyć je w jedną warstwę rastrową. Ta warstwa rastrowa wymaga prostej algebry mapy: suma / 5, dzięki czemu każde źródło będzie w końcu powiązane z pojedynczą warstwą rastrową, która ma wartości w [NULL, 0,2, 0,4, 0,6, 0,8, 1] w zależności od liczby warstwy składowe zachodzą na siebie. Jest to powierzchnia prawdopodobieństwa.

Wszystkie moje dane są przechowywane w Postgres 9.3 (z PostGIS). Mój problem polega na tym, że chociaż chcę nauczyć się korzystać z rastra PostGIS, wydaje się, że ma on bardzo stromą krzywą uczenia się i wszystkie przykłady, które mogę znaleźć, dotyczą jednej warstwy rastrowej. W przykładach warstwa ta jest używana jako część nakładki wielokąta, być może uśredniając wartość rastra dla każdego wielokąta. Nie znalazłem powtarzalnego przykładu łączenia: a) wektor -> raster b) algebra mapy; oraz c) atrybut GROUP BY zgodnie z moim pierwszym akapitem.

Mogę używać GDAL lub GRASS, jeśli muszę, aby wykonać to zadanie, ale wydaje się, że PostGIS powinien być w stanie to zrobić; byłoby to wygodne, ponieważ moje dane wejściowe to już geometria PostGIS; i naprawdę chcę pogodzić się z rastrem PostGIS.

Niektóre przykładowe struktury danych:

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

Chcę zrasteryzować, pogrupować według areaid, a następnie wykonać algebrę mapy, aby przejść do:

areaid    isorast (raster)
1000      aaa
1006      bbb

Nie udało mi się tego zawrzeć w PostGIS. Moje podejście polegało na przekształceniu wektora w raster, zrzuceniu rastrów na tablice i wykonaniu kombinacji z tablicami numpy przez psycopg2, przed zapisaniem ich w GeoTIFF (może być ponownie wprowadzone w PostGIS). Nie idealny, ale możliwy do wykonania.

— zupa alfabetyczna
źródło

Fajne pytanie. Podzielam, że naprawdę powinieneś nauczyć się postgisowego rastra i jestem pewien, że to, czego chcesz, jest możliwe. Niestety dzisiaj jest zbyt zajęty, aby spróbować.

— John Powell,

Na blogu BostonGIS znajduje się dość ostry artykuł na temat algebry map . Autor tego bloga jest także autorem doskonałej książki Postgis in Action, która ma wiele możliwości rastrowych Postgis. Przepraszam, nie mogłem wymyślić bardziej bezpośredniego przykładu.

— John Powell,

Musisz napisać własną funkcję agregującą:

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
);

potem możesz to tak nazwać

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

Daje to sumę wszystkich rastrów o tym samym identyfikatorze obszaru. Nadal będziesz musiał podzielić wartości rastrowe przez liczbę obserwacji dla każdego identyfikatora obszaru. (Nie dodałem go do funkcji agregującej. Możesz to zrobić tutaj lub później za pomocą MapAlegbra)

Upewnij się, że wszystkie wejściowe rastry są wyrównane, w przeciwnym razie to nie zadziała.

— Tomasz
źródło