Cálculo del valor medio del polígono a partir de la trama en PostGIS?

8

Comencé con un archivo raster Netgrf .gri y .grd del Reino Unido proporcionado por un colega. Lo recorté en R para que fuera solo Londres, lo exporté y lo convertí en un archivo ASC, y luego lo importé a PostGIS usando los siguientes comandos en R:

library(raster)
uk_raster <- raster("AnnMean2011.grd")
london_area <- extent(-720000.0,-630000.0,-50000.0,25000)
london_raster  <- crop(uk_raster, london_area)
writeRaster(london_raster, filename="AnnMean2011.asc", format="ascii")

Y luego en la línea de comandos de Ubuntu:

raster2pgsql -I -C -s 10001 -t 20x20 AnnMean2011.asc annualmean | psql -d james_traffic

Ahora tengo una tabla ráster en PostGIS. El SRID de 10001 es el siguiente, por cierto, proporcionado por un colega:

INSERT INTO spatial_ref_sys(srid, auth_name, auth_srid, proj4text)
VALUES (10001,'CMAQ_Urban',10001,'+proj=lcc +a=6370000 +b=6370000 +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=000000 +y_0=00000');

En la misma base de datos tengo un archivo de polígono, SRID 27700, que cubre Londres. Me gustaría calcular el valor promedio dentro de cada polígono, a partir del ráster.

Estoy intentando algo como esto, pero no está bien:

select polygons.postcode, avg(st_value(joined_data.rast))
    from (
       select (ST_Intersection(raster.rast, 1, polygons.geom)).*
    from raster, polygons 
       where ST_Intersects(raster.rast, 1, polygons.geom)
  ) joined_data
group by polygons.postcode

¿Cómo haría esto, por favor?

Parece que el polígono y el ráster se alinean correctamente, creo que tengo que convertirlos a WGS84.

Ráster con polígonos, visto en QGIS

TheRealJimShady
fuente
Tenga en cuenta un duplicado bastante, pero su respuesta probablemente esté aquí: stackoverflow.com/questions/24083732/…
GIS-Jonathan
Hmmm Gracias GIS-Jonathan, pero estoy luchando por traducir eso en mi conjunto de datos / situación. Estoy intentando algo como esto pero no está bien (pregunta editada arriba para incluirlo)
TheRealJimShady
Si aún no tiene una solución, puede valer la pena preguntar en la lista de PostGIS.
GIS-Jonathan
1
Creo que esto podría ser interesante para usted: gis.stackexchange.com/questions/76522/… . Una consulta exacta pero lenta en la pregunta y una consulta rápida y menos exacta en mi respuesta. Más información también aquí: postgis.net/docs/RT_ST_SummaryStats.html (PostGIS Doc !!!). Literatura: PostGIS Cookbook. Paolo Corti et al. !!!
Stefan

Respuestas:

6

Gracias al comentario de Stefan ayer, creo que puedo juntar algo de las preguntas relacionadas y la documentación oficial y ofrecer una gama de soluciones.

Documentación de PostGIS ( ST_SummaryStats)

Resumir píxeles que se cruzan con edificios de interés.

Este ejemplo tomó 574 ms en ventanas PostGIS de 64 bits con todos los edificios de Boston y mosaicos aéreos (mosaicos cada 150x150 píxeles ~ 134,000 mosaicos), ~ 102,000 registros de construcción

WITH 
-- our features of interest
   feat AS (SELECT gid As building_id, geom_26986 As geom FROM buildings AS b 
    WHERE gid IN(100, 103,150)
   ),
-- clip band 2 of raster tiles to boundaries of builds
-- then get stats for these clipped regions
   b_stats AS
    (SELECT  building_id, (stats).*
FROM (SELECT building_id, ST_SummaryStats(ST_Clip(rast,2,geom)) As stats
    FROM aerials.boston
        INNER JOIN feat
    ON ST_Intersects(feat.geom,rast) 
 ) As foo
 )
-- finally summarize stats
SELECT building_id, SUM(count) As num_pixels
  , MIN(min) As min_pval
  , MAX(max) As max_pval
  , SUM(mean*count)/SUM(count) As avg_pval
    FROM b_stats
 WHERE count > 0
    GROUP BY building_id
    ORDER BY building_id;

 building_id | num_pixels | min_pval | max_pval |     avg_pval
-------------+------------+----------+----------+------------------
         100 |       1090 |        1 |      255 | 61.0697247706422
         103 |        655 |        7 |      182 | 70.5038167938931
         150 |        895 |        2 |      252 | 185.642458100559

Evitar ST_Intersectionpor rendimiento

Tenga en cuenta que esto es menos exacto, ya que se ignoran los píxeles de intersección que cubren menos del 50% de una geometría de intersección.

Stefan tiene una respuesta aquí que evita el ST_Intersection.

Notas

  • También puede encontrar algunos consejos útiles en esta pregunta y en el tutorial de trama PostGIS WKT .
  • Si su ráster está en mosaico, una buena regla general es usar mosaicos más pequeños, tratando de hacerlos un poco más grandes que una característica vectorial típica que está intersectando.
alphabetasoup
fuente