Suma de ráster PostGIS (álgebra de mapas)

Question

Tengo una tabla de polígonos que representan isocronas de tiempo de viaje en días particulares. Para cada punto de origen, hay cinco geometrías de isocrona (almacenadas en filas separadas). Para cada punto de origen, quiero rasterizar las cinco isocronas (un NULL binario o 1), y luego combinarlas en una sola capa de ráster. Esta capa de ráster requiere un álgebra de mapas simple: suma/5, de modo que cada origen esté asociado al final con una única capa de ráster que tiene valores en [NULL, 0.2, 0.4, 0.6, 0.8, 1] dependiendo de cuántas capas constituyentes se superponen. Es una superficie de probabilidad.

Mis datos están todos almacenados en Postgres 9.3 (con PostGIS). Mi problema es que aunque quiero aprender a usar el ráster de PostGIS, parece tener una curva de aprendizaje muy pronunciada, y todos los ejemplos que puedo encontrar tratan con una sola capa de ráster. En los ejemplos, esta capa se utiliza como parte de una superposición de polígonos, quizás promediando el valor del ráster para cada polígono. No he encontrado un ejemplo replicable para combinar: a) vector --> ráster b) álgebra de mapas; y c) AGRUPAR POR atributo según mi primer párrafo.

Estoy bien usando GDAL o GRASS si tengo que hacerlo para realizar esta tarea, pero parece algo que PostGIS debería poder manejar; sería conveniente hacerlo dado que mis datos de entrada ya son geometría de PostGIS; y realmente quiero dominar el ráster de PostGIS.

Algunos ejemplos de estructura de datos:

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

Quiero rasterizar, agrupar por areaid, y luego realizar el álgebra de mapas para llegar a:

areaid    isorast (ráster)
1000      aaa
1006      bbb

---

No he tenido éxito conteniendo esto en PostGIS. Mi enfoque ha sido convertir el vector a ráster, volcar los rásters en matrices, y realizar la combinación con matrices de numpy a través de psycopg2, antes de escribirlos en un GeoTIFF (para tal vez volver a ponerlo en PostGIS). No es ideal, pero se puede hacer.

Answer 1

Necesitarás escribir tu propia función de agregado:

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

después puedes llamarlo de esta manera

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

Esto te dará la suma de todos tus rasters con el mismo ID de área. Aún necesitarás dividir los valores de los rasters por el número de observaciones para cada ID de área. (No lo incluí en la función de agregado. Puedes hacerlo aquí o después usando MapAlegbra)

Asegúrate de que tus rasters de entrada estén alineados, de lo contrario esto no funcionará.