Estoy tratando de calcular estadísticas ráster (min, max, mean) para cada polígono en una capa vectorial usando PostgreSQL / PostGIS.
Esta respuesta GIS.SE describe cómo hacer esto, calculando la intersección entre el polígono y el ráster y luego calculando un promedio ponderado: https://gis.stackexchange.com/a/19858/12420
Estoy usando la siguiente consulta (donde dem
está mi ráster, topo_area_su_region
es mi vector y toid
es una ID única:
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 tarda alrededor de 90 ms en procesarse; tomaría días procesar la capa completa. La velocidad del cálculo no parece mejorar significativamente si solo 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 lleva procesar los datos, si en lugar de vectorizar el ráster (usando ST_Intersection) rasterizo el vector. Ver código aquí: https://gist.github.com/snorfalorpagus/7320167
Realmente no necesito un promedio ponderado, un enfoque de "si toca, está en" es lo suficientemente bueno, y estoy razonablemente seguro de que esto es lo que está frenando las cosas.
Pregunta : ¿Hay alguna forma de hacer que PostGIS se comporte así? es decir, devolver los valores de todas las celdas del ráster que toca un polígono, en lugar de la intersección exacta.
Soy muy nuevo en PostgreSQL / PostGIS, así que tal vez hay algo más que no estoy haciendo bien. Estoy ejecutando PostgreSQL 9.3.1 y PostGIS 2.1 en Windows 7 (2.9GHz i7, 8GB RAM) y modifiqué la configuración de la base de datos como se sugiere aquí: http://postgis.net/workshops/postgis-intro/tuning.html
fuente
Respuestas:
Tienes razón, el uso
ST_Intersection
ralentiza tu consulta notablemente.En lugar de usarlo
ST_Intersection
, es mejor recortar (ST_Clip
) su ráster con los polígonos (sus campos) y volcar el resultado como polígonos (ST_DumpAsPolygons
). Por lo tanto, cada celda ráster se convertirá en un pequeño rectángulo poligonal con valores distintos.Para recibir min, max o mean de los volcados, puede usar las mismas declaraciones.
Esta consulta debería hacer el truco:
En la declaración
ST_Clip
, define el ráster, la banda de ráster (= 1), el polígono y si el recorte debe ser VERDADERO o FALSO.Además puedes usar
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
yST_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 CORINE Land Use muestran la diferencia. Primera foto con
ST_Intersection
, segunda conST_Clip
yST_DumpAsPolygons
.fuente