Saya mulai dengan file raster dari NetCDF .gri dan .grd dari Inggris yang disediakan oleh seorang rekan. Saya memotongnya di R menjadi London, diekspor dan dikonversi ke file ASC, dan kemudian mengimpornya ke PostGIS menggunakan perintah berikut di 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")
Dan kemudian pada baris perintah Ubuntu:
raster2pgsql -I -C -s 10001 -t 20x20 AnnMean2011.asc annualmean | psql -d james_traffic
Saya sekarang memiliki tabel raster di PostGIS. SRID 10001 adalah yang berikut, sekali lagi, disediakan oleh seorang rekan:
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');
Dalam database yang sama saya punya file poligon, SRID 27700, yang mencakup London. Saya ingin menghitung nilai rata-rata dalam setiap poligon, dari raster.
Saya mencoba sesuatu seperti ini, tetapi itu tidak benar:
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
Bagaimana saya bisa melakukan ini?
Tampaknya mendapatkan poligon dan raster selaras dengan benar, saya harus mengubah keduanya menjadi WGS84 saya pikir.
sumber
Jawaban:
Berkat komentar dari Stefan kemarin, saya pikir saya bisa menyatukan sesuatu dari pertanyaan terkait dan dokumentasi resmi dan menawarkan berbagai solusi.
Dokumentasi PostGIS (
ST_SummaryStats
)Ringkas piksel yang memotong bangunan yang diminati
Contoh ini mengambil 574ms pada PostGIS windows 64-bit dengan semua Bangunan Boston dan Ubin Udara (setiap ubin 150x150 piksel ~ 134.000 ubin), ~ 102.000 catatan bangunan
Menghindari
ST_Intersection
kinerjaPerhatikan bahwa ini kurang tepat, dengan piksel berpotongan yang menutupi kurang dari 50% geometri berpotongan diabaikan.
Stefan punya jawaban di sini yang menghindari
ST_Intersection
.Catatan
sumber