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.
fuente
Respuestas:
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
Evitar
ST_Intersection
por rendimientoTenga 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
fuente