Rendimiento en el cálculo de estadísticas raster en PostGIS

Question

Estoy tratando de calcular las estadísticas de rasterización (mínimo, máximo, medio) para cada polígono en una capa vectorial utilizando PostgreSQL/PostGIS.

Esta respuesta de GIS.SE describe cómo hacerlo, calculando la intersección entre el polígono y el raster y luego calculando un promedio ponderado: https://gis.stackexchange.com/a/19858/12420

Estoy utilizando la siguiente consulta (donde dem es mi trama, topo_area_su_region es mi vector, y toid es un identificador único:

SELECT toid, Min((gv).val) As MinElevation, Max((gv).val) As MaxElevation, Sum(ST_Area((gv).geom) * (gv).val) / Sum(ST_Area((gv).geom)) as MeanElevation FROM (SELECT toid, ST_Intersection(rast, geom) AS gv FROM topo_area_su_region,dem WHERE ST_Intersects(rast, geom)) foo GROUP BY toid ORDER BY toid;

Esto funciona, pero es demasiado lento. Mi capa vectorial tiene 2489k características, y cada una toma alrededor de 90ms para procesar - se necesitaría días para procesar toda la capa. La velocidad de cálculo no parece mejorar significativamente si sólo calculo el mínimo y el máximo (lo que evita las llamadas a ST_Area).

Si hago un cálculo similar usando Python (GDAL, NumPy y PIL) puedo reducir significativamente la cantidad de tiempo que toma procesar los datos, si en lugar de vectorizar el raster (usando ST_Intersection) rasterizo el vector. Vea el código aquí: https://gist.github.com/snorfalorpagus/7320167

Realmente no necesito una media ponderada - un enfoque de "si toca, está dentro" es suficiente - y estoy razonablemente seguro de que esto es lo que está ralentizando las cosas.

Pregunta : ¿Existe alguna forma de conseguir que PostGIS se comporte así? Es decir, que devuelva los valores de todas las celdas del raster que toca un polígono, en lugar de la intersección exacta.

Soy muy nuevo en PostgreSQL/PostGIS, así que tal vez hay algo que no estoy haciendo bien. Estoy corriendo PostgreSQL 9.3.1 y PostGIS 2.1 en Windows 7 (2.9GHz i7, 8GB RAM) y he ajustado la configuración de la base de datos como se sugiere aquí: http://postgis.net/workshops/postgis-intro/tuning.html

Answer 1

Tienes razón, usar ST_Intersection ralentiza notablemente su consulta.

En lugar de utilizar ST_Intersection es mejor recortar ( ST_Clip ) su trama con los polígonos (sus campos) y volcar el resultado como polígonos ( ST_DumpAsPolygons ). Así, cada celda rasterizada se convertirá en un pequeño rectángulo poligonal con valores distintos.

Para recibir el mínimo, el máximo o la media de los volcados se pueden utilizar las mismas sentencias.

Esta consulta debería servir:

SELECT 
    toid,
    Min((gv).val) As MinElevation,
    Max((gv).val) As MaxElevation,
    Sum(ST_Area((gv).geom) * (gv).val) / Sum(ST_Area((gv).geom)) as MeanElevation
FROM (
    SELECT 
        toid,
        ST_DumpAsPolygons(ST_Clip(rast, 1, geom, true)) AS gv
    FROM topo_area_su_region,dem 
        WHERE ST_Intersects(rast, geom)) AS foo 
            GROUP BY toid 
            ORDER BY toid;

En la declaración ST_Clip se define la trama, la banda de trama (=1), el polígono y si el cultivo debe ser TRUE o FALSE.

Además puedes utilizar avg((gv).val) para calcular el valor medio.

EDITAR

El resultado de su enfoque es el más exacto, pero el más lento. Los resultados de la combinación de ST_Clip et ST_DumpAsPolygons están ignorando las celdas ráster que se cruzan con menos del 50% (o 51%) de su tamaño.

Estas dos capturas de pantalla de una intersección de CORINE Land Use muestran la diferencia. La primera imagen con ST_Intersection , el segundo con ST_Clip et ST_DumpAsPolygons .