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

Tengo una tabla de polígonos que representan las isócronas del tiempo de viaje en días particulares. Para cada punto de origen, hay cinco geometrías de isócronas (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 trama. Esta capa ráster requiere un álgebra de mapa simple: suma / 5, de modo que cada origen se asociará al final con una sola capa ráster que tiene valores en [NULL, 0.2, 0.4, 0.6, 0.8, 1] dependiendo de cuántos de Las capas constituyentes se superponen. Es una superficie de probabilidad.

Todos mis datos están almacenados en Postgres 9.3 (con PostGIS). Mi problema es que, si bien quiero aprender a usar el ráster PostGIS, parece tener una curva de aprendizaje realmente empinada, y todos los ejemplos que puedo encontrar tratan de una sola capa ráster. En los ejemplos, esta capa se usa 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) atributo GROUP BY según mi primer párrafo.

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

Alguna estructura de datos de muestra:

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

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

areaid    isorast (raster)
1000      aaa
1006      bbb

No he tenido éxito al contener esto en PostGIS. Mi enfoque ha sido convertir el vector en ráster, volcar los rásteres en matrices y realizar la combinación con matrices numpy a través de psycopg2, antes de escribirlas en un GeoTIFF (para quizás volver a colocarlo en PostGIS). No es ideal, pero factible.

Necesitará escribir su propia función agregada:

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

luego puedes llamarlo así

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

Esto le da la suma de todos sus rásteres con la misma ID de área. Aún deberá dividir los valores ráster por el número de observaciones para cada ID de área. (No lo incluí en la función agregada. Puede hacerlo aquí o después usando MapAlegbra)

Asegúrese de que todos los rásteres de entrada estén alineados; de lo contrario, esto no funcionará.