Estoy usando PostGIS2.0 para hacer algunas intersecciones ráster / polígono. Tengo dificultades para comprender qué operación debo usar y cuál es la forma más rápida de realizarla. Mi problema es el siguiente:
- Tengo un polígono y una trama
- Quiero encontrar todos los píxeles que se encuentran dentro del polígono y obtener la suma del valor del píxel
- Y (problema actualizado): obtengo valores masivos para algunos píxeles que no existen en el ráster original cuando realizo la consulta
Tengo dificultades para entender si debo usar ST_Intersects()
o ST_Intersection()
. Tampoco sé cuál es el mejor enfoque para sumar mis píxeles. Aquí está el primer enfoque que he probado (# 1):
SELECT
r.rast
FROM
raster as r,
polygon as p
WHERE
ST_Intersects(r.rast, p.geom)
Esto devuelve una lista de rast
valores, con lo que no estoy seguro de qué hacer. Intenté calcular las estadísticas de resumen usando, ST_SummaryStats()
pero no estoy seguro de si esta es la suma ponderada de todos los píxeles que se encuentran dentro del polígono.
SELECT
(result).count,
(result).sum
FROM (
SELECT
ST_SummaryStats(r.rast) As result
FROM
raster As r,
polygon As p
WHERE
ST_Intersects(r.rast, p.geom)
) As tmp
El otro enfoque que he probado (# 2) usa ST_Intersection()
:
SELECT
(gv).geom,
(gv).val
FROM
(
SELECT
ST_Intersection(r.rast, p.geom) AS gv
FROM
raster as r,
polygon as p
WHERE
ST_Intersects(r.rast, p.geom)
) as foo;
Esto devuelve una lista de geometrías que analizo más a fondo, pero supongo que esto es menos eficiente.
No tengo claro cuál es el orden de operación más rápido también. ¿Debo elegir siempre raster, polygon
o polygon, raster
, o convertir el polígono en un ráster para que así sea raster, raster
?
EDITAR: Actualicé el enfoque # 2 con algunos detalles del R.K.
enlace.
Usando el enfoque # 2, he notado el siguiente error en los resultados, que es parte de la razón por la que no entendí el resultado. Aquí está la imagen de mi ráster original, y un esquema del polígono que se está utilizando para intersecarlo, superpuesto en la parte superior:
Y aquí está el resultado de la intersección usando PostGIS:
El problema con el resultado es que hay valores de 21474836 que se devuelven, que no están en el ráster original. No tengo idea de por qué ocurre esto. Sospecho que está relacionado con números pequeños en algún lugar (dividido por casi 0), pero devuelve el resultado incorrecto.
ST_SummaryStats()
para el # 1, pero no estoy seguro de cómo hacerlo para el # 2.Respuestas:
Encontré un tutorial sobre la intersección de polígonos vectoriales con una gran cobertura de trama usando PostGIS WKT Raster . Podría tener las respuestas que estás buscando.
El tutorial utilizó dos conjuntos de datos, un archivo de forma de punto que se almacenó para producir polígonos y una serie de 13 rásteres de elevación SRTM. Hubo muchos pasos intermedios, pero la consulta utilizada para intersectar el ráster y el vector se veía así:
Los valores se resumieron luego usando lo siguiente:
Realmente no sé suficiente PostGIS para explicar esto, pero seguro que suena como lo que estabas tratando de lograr. El tutorial debería arrojar luz sobre los pasos intermedios. Buena suerte :)
fuente
Con respecto al punto 2 de la pregunta original, varias de las versiones de desarrollo de postgis 2.0 utilizaron una versión de la biblioteca GDAL que convirtió los flotantes en ints. Si su ráster tiene valores flotantes y estaba usando una versión de GDAL inferior a 1.9.0, o una versión de prelanzamiento de PostGIS 2.0 que no llamó correctamente a GDALFPolygonize (), entonces podría estar encontrando este error. Las entradas en los rastreadores de errores PostGIS y GDAL se archivaron y cerraron. Este error estaba activo en el momento de la pregunta original.
En términos de rendimiento, encontrará que usar
ST_Intersects(raster, geom)
es mucho más rápido que usarST_Intersects(geom, raster)
. La primera versión rasteriza la geometría y hace una intersección espacio-ráster. La segunda versión vectoriza la geometría y realiza una intersección espacio-vector, que puede ser un proceso mucho más costoso.fuente
También estaba teniendo problemas extraños al usar
ST_SummaryStats
conST_Clip
. Consultar los datos de manera diferente me dijo que el valor mínimo de mi ráster era 32, y luego máximo 300, sin embargo,ST_SummaryStats
estaba devolviendo -32700 para los valores de píxeles en mi polígono objetivo.Terminé hackeando el tema así:
fuente