Calcular el valor de la media de polígono de PostGIS raster en

Question

Empecé con un NetCDF .gri y .catálogo de ráster archivo del reino unido proporcionada por un colega. He recortado en R para ser sólo de Londres, exportado y se convierte en un archivo ASC, y luego importados 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, a continuación, en el Ubuntu de la línea de comandos:

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

Ahora tengo un ráster en la tabla PostGIS. El SRID de 10001 es la siguiente por el camino, de nuevo, proporcionada 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 polígono de archivo, SRID 27700, que cubre Londres. Quiero calcular el valor promedio dentro de cada polígono, de la trama.

Estoy tratando de hacer algo como esto, pero no es correcto:

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 se podría ir sobre esto, por favor?

Gracias JDS

PD: parece que el polígono y la trama alineados correctamente, tengo que convertirlos a ambos en WGS84 creo.

Answer 1

Gracias al comentario de Stefan ayer, creo que puedo hacer algo juntos de preguntas relacionadas y la documentación oficial y ofrecen una amplia gama de soluciones.

PostGIS documentación (ST_SummaryStats)

Resumen de los píxeles que se cruzan edificios de interés

En este ejemplo se tomó 574ms en PostGIS windows de 64 bits con todos los de Boston Edificios y aéreos Azulejos (azulejos cada uno de 150x150 píxeles ~ 134,000 de azulejos), ~102,000 la construcción de registros

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

Evitando ST_Intersects de rendimiento

Tenga en cuenta que esto es menos exacto, con la intersección de los píxeles que cubren menos del 50% de la intersección de la geometría de ser ignorados.

Stefan tiene una respuesta aquí que evita la ST_Intersects.

Notas

• Usted también puede encontrar algunos consejos útiles en esta pregunta y PostGIS WKT tutorial de ráster.
• Si su trama es de baldosas, una buena regla del pulgar es el uso de pequeños azulejos, tratando de conseguir un poco más grande que la de un vector de característica que estás de intersección.