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

13

Tengo una tabla de polígonos que representan los isócronos de tiempo de viaje en días específicos. Para cada punto de origen, hay cinco geometrías isócronas (almacenadas en filas separadas). Para cada punto de origen, quiero rasterizar los cinco isócronos (un NULL binario o 1) y luego combinarlos en una sola capa ráster. Esta capa ráster requiere un álgebra de mapa simple: suma / 5, para que cada origen se asocie al final con una sola capa ráster que tenga valores en [NULL, 0.2, 0.4, 0.6, 0.8, 1] dependiendo de la cantidad de las capas constituyentes se superponen. Es una superficie de probabilidad.

Todos mis datos se almacenan en Postgres 9.3 (con PostGIS). Mi problema es que, si bien quiero aprender a usar el ráster PostGIS, parece que tiene una curva de aprendizaje muy pronunciada, y todos los ejemplos que puedo encontrar tratan con una sola capa ráster. En los ejemplos, esta capa se usa como parte de una superposición de polígonos, tal vez 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; yc) atributo GROUP BY según mi primer párrafo.

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

Algunas estructuras 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 área y luego realizar el álgebra del mapa para:

areaid    isorast (raster)
1000      aaa
1006      bbb

No he tenido éxito conteniendo 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 volver a colocarlas en PostGIS). No es ideal, pero se puede hacer.

    
pregunta Richard Law 05.12.2014 - 13:33

1 respuesta

1

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

después 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 el mismo ID de área. Aún tendrá que 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á.

    
respondido por el Thomas 16.02.2016 - 01:21

Lea otras preguntas en las etiquetas