Soma de varredura PostGIS (álgebra de mapa)

Eu tenho uma tabela de polígonos representando isócronas no tempo de viagem em dias específicos. Para cada ponto de origem, existem cinco geometrias de isócronas (armazenadas em linhas separadas). Para cada ponto de origem, desejo rasterizar as cinco isócronas (um NULL binário ou 1) e depois combiná-las em uma única camada raster. Essa camada raster requer uma álgebra simples do mapa: soma / 5, para que cada origem seja associada a uma única camada raster que possua valores em [NULL, 0,2, 0,4, 0,6, 0,8, 1], dependendo de quantos as camadas constituintes se sobrepõem. É uma superfície de probabilidade.

Meus dados são todos armazenados no Postgres 9.3 (com PostGIS). Meu problema é que, embora eu queira aprender a usar a varredura PostGIS, ela parece ter uma curva de aprendizado muito acentuada, e todos os exemplos que posso encontrar lidam com uma única camada de varredura. Nos exemplos, essa camada é usada como parte de uma sobreposição de polígono, talvez com a média do valor da varredura para cada polígono. Não encontrei um exemplo replicável para combinar: a) vetor -> raster b) álgebra de mapa; e c) atributo GROUP BY conforme meu primeiro parágrafo.

Eu estou bem usando GDAL ou GRASS se for necessário para executar esta tarefa, mas isso parece algo que o PostGIS deve ser capaz de lidar; seria conveniente fazer isso, já que meus dados de entrada já são geometria PostGIS; e eu realmente quero aceitar o raster do PostGIS.

Alguma estrutura de dados de amostra:

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

Quero rasterizar, agrupar por areaid e, em seguida, executar a álgebra do mapa para chegar a:

areaid    isorast (raster)
1000      aaa
1006      bbb

Não consegui contê-lo no PostGIS. Minha abordagem foi converter o vetor em varredura, despejar os rasters em matrizes e executar a combinação com matrizes numpy via psycopg2, antes de gravá-las em um GeoTIFF (talvez para ser postado novamente no PostGIS). Não é o ideal, mas é capaz.

Você precisará escrever sua própria função 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
);

depois você pode chamar assim

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

Isso fornece a soma de todos os seus rasters com o mesmo ID de área. Você ainda precisará dividir os valores de varredura pelo número de observações para cada ID de área. (Eu não o incluí na função agregada. Você pode fazê-lo aqui ou depois usando o MapAlegbra)

Certifique-se de que todos os controles de entrada estejam alinhados, caso contrário, isso não funcionará.