J'ai commencé avec un fichier raster NetCDF .gri et .grd du Royaume-Uni fourni par un collègue. Je l'ai découpé dans R pour être uniquement Londres, exporté et converti en fichier ASC, puis importé dans PostGIS à l'aide des commandes suivantes dans 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")

Et puis sur la ligne de commande Ubuntu:

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

J'ai maintenant une table raster dans PostGIS. Le SRID de 10001 est d'ailleurs le suivant, fourni par un collègue:

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');

Dans la même base de données, j'ai un fichier polygonal, SRID 27700, qui couvre Londres. Je voudrais calculer la valeur moyenne à l'intérieur de chaque polygone, à partir du raster.

J'essaie quelque chose comme ça, mais ce n'est pas 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

Comment pourrais-je m'y prendre s'il vous plaît?

Il semble que le polygone et le raster soient correctement alignés, je dois les convertir tous les deux en WGS84 je pense.

Grâce au commentaire de Stefan hier, je pense que je peux rassembler quelque chose à partir de questions connexes et de la documentation officielle et proposer une gamme de solutions.

Documentation PostGIS ( ST_SummaryStats)

Résumer les pixels qui croisent les bâtiments d'intérêt

Cet exemple a pris 574 ms sur les fenêtres PostGIS 64 bits avec tous les bâtiments de Boston et les tuiles aériennes (tuiles de 150x150 pixels environ ~ 134 000 tuiles), ~ 102 000 enregistrements de construction

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

Éviter les ST_Intersectionperformances

Notez que cela est moins exact, les pixels qui se croisent qui couvrent moins de 50% d'une géométrie qui se croisent sont ignorés.

Stefan a une réponse ici qui évite le ST_Intersection.

Remarques

• Vous trouverez également des conseils utiles dans cette question et dans le didacticiel raster PostGIS WKT .
• Si votre raster est en mosaïque, une bonne règle de base consiste à utiliser des tuiles plus petites, en essayant de les rendre un peu plus grandes qu'une entité vectorielle typique que vous intersectez.