Intersección de un raster con un polígono usando PostGIS - error de artefacto

15

Estoy usando PostGIS2.0 para hacer algunas intersecciones de trama / polígono. Tengo dificultades para comprender qué operación debo usar y cuál es la forma más rápida de realizar esta operación. Mi problema es el siguiente:

  • Tengo un polígono y un raster
  • Quiero encontrar todos los píxeles que caen dentro del polígono y obtener la suma del valor en píxeles
  • Y (problema actualizado): estoy obteniendo 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 valores de rast , con lo que no estoy seguro de qué hacer. Intenté calcular las estadísticas de resumen utilizando 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 asumo 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 sea raster, raster ?

EDITAR: Actualicé el enfoque # 2 con algunos detalles del enlace de R.K. .

Al utilizar el método n.º 2, he notado el siguiente error en los resultados, lo cual es parte de la razón por la que no entendí el resultado. Aquí está la imagen de mi ráster original y un contorno del polígono que se está utilizando para intersectarlo, superpuesto en la parte superior:

YaquíestáelresultadodelaintersecciónutilizandoPostGIS:

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é esto está ocurriendo. Sospecho que está relacionado con números pequeños en algún lugar (dividido por casi 0), pero devuelve el resultado incorrecto.

    
pregunta djq 06.02.2012 - 15:38

3 respuestas

6

Encontré un tutorial en intersección de polígonos vectoriales con una gran cobertura de trama utilizando Post Rast WKT Raster . Podría tener las respuestas que estás buscando.

El tutorial usó dos conjuntos de datos, un archivo de forma de punto que se almacenó en búfer para producir polígonos y una serie de 13 rásteres de elevación SRTM. Hubo muchos pasos intermedios, pero la consulta utilizada para cruzar el ráster y el vector tenía este aspecto:

 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 utilizando 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 conozco suficiente PostGIS para explicar esto, pero seguramente suena como lo que intentabas lograr. El tutorial debería arrojar luz sobre los pasos intermedios. Buena suerte :)

    
respondido por el R.K. 06.02.2012 - 16:27
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 convierte flots a ints. Si su ráster tiene valores flotantes en él, y estaba usando una versión de GDAL inferior a 1.9.0, o una versión de PostGIS 2.0 que no se llamó correctamente a GDALFPolygonize (), es posible que se encuentre con este error. Entradas en el PostGIS y GDAL se registraron y cerraron los rastreadores de errores. Este error estaba activo en el momento de la pregunta original.

En términos de rendimiento, encontrarás 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 de espacio ráster. La segunda versión vectoriza la geometría y hace una intersección de espacio vectorial, que puede ser un proceso mucho más costoso.

    
respondido por el Pete Clark 06.06.2012 - 22:15
0

También tuve problemas extraños al utilizar ST_SummaryStats con ST_Clip . Al consultar los datos de manera diferente, me dijeron que el valor mínimo de mi ráster era 32 y luego un máximo de 300, pero ST_SummaryStats estaba devolviendo -32700 para los valores de píxeles en mi polígono objetivo.

Terminé hackeando el problema como tal:

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
    
respondido por el jczaplew 16.11.2015 - 23:18

Lea otras preguntas en las etiquetas