Penjumlahan raster PostGIS (aljabar peta)

Saya memiliki tabel poligon yang menunjukkan waktu tempuh isochrone pada hari-hari tertentu. Untuk setiap titik asal, ada lima geometri isochrone (disimpan dalam baris terpisah). Untuk setiap titik asal, saya ingin rasterise kelima isochrones (NULL biner atau 1), dan kemudian menggabungkannya menjadi satu layer raster. Lapisan raster ini membutuhkan aljabar peta sederhana: jumlah / 5, sehingga setiap asal pada akhirnya akan dikaitkan dengan lapisan raster tunggal yang memiliki nilai dalam [NULL, 0,2, 0,4, 0,6, 0,8, 1] tergantung pada berapa banyak lapisan konstituen tumpang tindih. Ini adalah permukaan probabilitas.

Semua data saya disimpan di Postgres 9.3 (dengan PostGIS). Masalah saya adalah bahwa sementara saya ingin belajar menggunakan raster PostGIS, tampaknya memiliki kurva belajar yang benar-benar curam, dan semua contoh saya dapat menemukan berurusan dengan satu layer raster. Dalam contoh-contoh, lapisan ini digunakan sebagai bagian dari overlay poligon, mungkin rata-rata nilai raster untuk setiap poligon. Saya belum menemukan contoh yang bisa ditiru untuk menggabungkan: a) vektor -> raster b) peta aljabar; dan c) atribut GROUP BY sesuai paragraf pertama saya.

Saya baik-baik saja menggunakan GDAL atau GRASS jika saya harus melakukan tugas ini, tetapi ini sepertinya sesuatu yang PostGIS harus dapat tangani; akan lebih mudah untuk melakukannya mengingat data input saya sudah geometri PostGIS; dan saya benar - benar ingin berdamai dengan raster PostGIS.

Beberapa struktur data sampel:

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

Saya ingin rasterise, kelompokkan berdasarkan areaid, dan kemudian lakukan aljabar peta untuk datang ke:

areaid    isorast (raster)
1000      aaa
1006      bbb

Saya belum berhasil memuat ini ke PostGIS. Pendekatan saya adalah untuk mengubah vektor menjadi raster, membuang raster ke array, dan melakukan kombinasi dengan array numpy via psycopg2, sebelum menulisnya ke GeoTIFF (untuk dimasukkan kembali ke PostGIS). Tidak ideal, tetapi bisa dilakukan.

Anda harus menulis fungsi agregat Anda sendiri:

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

setelah itu Anda bisa menyebutnya seperti ini

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

Ini memberi Anda jumlah dari semua raster Anda dengan ID area yang sama. Anda masih perlu membagi nilai raster dengan jumlah pengamatan untuk setiap ID area. (Saya tidak memasukkannya dalam fungsi agregat. Anda bisa melakukannya di sini atau sesudahnya menggunakan MapAlegbra)

Pastikan raster input Anda semua sejajar, jika tidak ini tidak akan berfungsi.