Intersección de un ráster con un polígono usando PostGIS - error de artefacto

15

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 rastvalores, 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, polygono 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:

ingrese la descripción de la imagen aquí

Y aquí está el resultado de la intersección usando PostGIS:

ingrese la descripción de la imagen aquí

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.

djq
fuente
Con respecto al punto dos, ¿desea obtener la suma de los valores de los píxeles que intersecan el polígono?
RK
Sí, lo he usado ST_SummaryStats()para el # 1, pero no estoy seguro de cómo hacerlo para el # 2.
djq
Ha publicado un enlace a una referencia. Espero que eso ayude.
RK
2
El valor máximo de la escala en su mapa es el máximo de un entero con signo de 32 bits. No sé qué significa eso en su caso, pero podría tener que ver con los valores de NA. El rango en su consulta puede tener valores nulos que no se manejan correctamente. en.wikipedia.org/wiki/2147483647#2147483647_in_computing
yellowcap
66
FWIW, 21474836 no es el valor máximo de un int con signo de 32 bits. Sin embargo, 2 ^ 31-1 = 2147483647 es el máximo, y observe que 21474836 = 2147483647/100 (división entera). Esto sugiere que internamente se generan algunos NA (quizás a lo largo de las celdas fronterizas), se representan como 2 ^ 31-1, y luego el código "olvida" que estos son NA y (¿tal vez en un proceso de remuestreo?) Los divide por 100.
whuber

Respuestas:

6

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í:

 CREATE TABLE caribou_srtm_inter AS
 SELECT id, 
        (gv).geom AS the_geom, 
        (gv).val
 FROM (SELECT id, 
              ST_Intersection(rast, the_geom) AS gv
       FROM srtm_tiled,
            cariboupoint_buffers_wgs
       WHERE ST_Intersects(rast, the_geom)
      ) foo;

Los valores se resumieron luego usando lo siguiente:

 CREATE TABLE result01 AS
 SELECT id, 
        sum(ST_Area(ST_Transform(the_geom, 32198)) * val) / 
        sum(ST_Area(ST_Transform(the_geom, 32198))) AS meanelev
 FROM caribou_srtm_inter
 GROUP BY id
 ORDER BY id;

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 :)

RK
fuente
Gracias @RK leí el tutorial. Creo que mi cálculo es más básico, ¡pero todavía estoy descubriendo este paso básico!
djq
2

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 usar ST_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.

Pete Clark
fuente
0

También estaba teniendo problemas extraños al usar ST_SummaryStatscon ST_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_SummaryStatsestaba devolviendo -32700 para los valores de píxeles en mi polígono objetivo.

Terminé hackeando el tema así:

WITH first AS (
   SELECT id, (ST_Intersection(geom, rast)).val
   FROM raster_table
   INNER JOIN vector_table ON ST_Intersects(rast, geom)
)
SELECT id, COUNT(val), SUM(val), AVG(val), stddev(val), MIN(val), MAX(val)
FROM first
GROUP BY id
jczaplew
fuente