Estadísticas zonales de superposición de polígonos en PostGIS

Question

Estoy intentando realizar una tarea aparentemente sencilla: extraer los valores medios de los píxeles de un raster basado en una superposición de polígonos en PostGIS. He consultado varias fuentes, incluyendo esto blog , este blog y esto Correo electrónico: ). He cargado un mosaico Landsat en mi base de datos (in-db) y una superposición de polígonos que consta de tres polígonos (ID 1,2,3). El mosaico Landsat es una imagen panorámica con una resolución espacial de 15 m y 4 bandas.

Basándome en los archivos de ayuda y en los posts/preguntas mencionados, llegué a la siguiente consulta de PostGIS:

SELECT id, (SUM((ST_SummaryStats(a.rast, 1, true)).sum)/SUM((ST_SummaryStats(a.rast, 1, true)).count)) as mean FROM imagery.l8_2015_09_12 AS a, analysis_results.zonal_stats_test as b WHERE ST_Intersects(b.geom,a.rast) GROUP BY id;

Los resultados son:

1: 9902.49
2: 10079.68
3: 12355.90

Sin embargo, cuando pruebo los resultados con la herramienta "estadísticas zonales a tabla" de ArcGIS, obtengo los siguientes resultados:

1: 8089.61
2: 7527.62
3: 12290.05

Cerca, pero no idéntico. Ya casi lo he conseguido, pero me preguntaba si alguien podría ayudarme a averiguar el origen de la diferencia entre ambas técnicas. Sospecho que hay algo mal en mi método ST_summaryStats.

Answer 1

1) Me sorprende que no "corte" la trama con ST_Clip(rast, geom) o ST_Intersection(rast, geom) antes de calcular sus estadísticas. Lo que quieres son las estadísticas de los valores de los píxeles que intersectan los polígonos, ¿no? Esa podría ser una de las razones por las que terminas con un resultado diferente, pero eso sólo sería si sólo se cortan pequeñas porciones de los rásteres, ya que tu resultado es bastante parecido al de ArcGIS.

2) Incluso si se "corta" el ráster correctamente antes de calcular las estadísticas, la forma en que PostGIS recorta o interseca el ráster puede ser diferente a la de ArcGIS: ¿Cómo se eligen los píxeles de intersección? ¿Por su centroide? ¿Por el área total de intersección? ¿Cómo se tratan los valores de los nodos? ¿Cómo se convierten los polígonos a raster antes de recortar o intersecar, si es el caso? ¿Cómo se convierten los rásteres en polígonos si es el caso? Todo esto depende de las herramientas que se seleccionen y, a veces, del orden en que se realicen las operaciones.

Te sugiero que trabajes con un subconjunto de tus datos y visualices/analices los resultados paso a paso (en OpenJump por ejemplo) para entender la diferencia con ArcGIS.

Ver también este tutorial .

Answer 2

Ese podría ser el caso de que cuando importaste el raster de Landsat a tu base de datos, con suerte pusiste un parámetro -t 300x300, que divide la imagen en mosaicos de 300x300 píxeles. Entonces, imagina que tienes tus polígonos en un shapefile que importaste a tu base de datos.

Es posible que necesite más de un mosaico rasterizado para cubrir una característica poligonal, por lo que, para evitar este problema, debe añadir la función St_union a su consulta, así:

SELECT 
id,
(St_SummaryStats(St_Union(ST_Clip(rast,1,geom, true)))).*
FROM  imagery.l8_2015_09_12, analysis_results.zonal_stats_test 
WHERE st_intersects(rast,geom)
GROUP BY id;

Además, al usar .* obtienes la suma, el recuento, el máximo, el mínimo y el stddev como bonus. Si quieres sólo la media, cambia * por media, esa línea terminaría así:

(St_SummaryStats(St_Union(ST_Clip(rast,1,geom, true)))).mean