I started with a NetCDF .gri and .grd raster file of the UK provided by a colleague. I clipped it in R to be only London, exported and converted to a ASC file, and then imported it to PostGIS using the following commands in R:

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")

And then on the Ubuntu command line:

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

I now have a raster table in PostGIS. The SRID of 10001 is the following by the way, again, provided by a colleague:

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

In the same database I have a polygon file, SRID 27700, which covers London. I would like to calculate the average value within each polygon, from the raster.

I'm trying something like this, but it's not right:

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

How would I go about this please?

It seems to get the polygon and the raster aligned properly, I have to convert them both to WGS84 I think.

Raster with Polygons, viewed in QGIS

  • Note quite a duplicate, but your answer is probably here: stackoverflow.com/questions/24083732/… Commented Jul 27, 2015 at 15:05
  • Hmmm. Thanks GIS-Jonathan, but I'm struggling to translate that into my dataset/situation. I'm trying something like this but it's not right (question edited above to include it) Commented Jul 28, 2015 at 14:29
  • If you still haven't got a solution, it may be worth asking on the PostGIS list. Commented Jul 29, 2015 at 8:18
  • 1
    I think this could be interesting for you: gis.stackexchange.com/questions/76522/…. One exact but slow query in the question and a fast, less exact query in my answer. Further information also here: postgis.net/docs/RT_ST_SummaryStats.html (PostGIS Doc!!!). Literatur: PostGIS Cookbook. Paolo Corti et al.!!!
    – Stefan
    Commented Nov 9, 2015 at 20:08

1 Answer 1


Thanks to the comment from Stefan yesterday, I think I can put something together from related questions and the official documentation and offer a range of solutions.

PostGIS documentation (ST_SummaryStats)

Summarize pixels that intersect buildings of interest

This example took 574ms on PostGIS windows 64-bit with all of Boston Buildings and aerial Tiles (tiles each 150x150 pixels ~ 134,000 tiles), ~102,000 building records

-- 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

Avoiding ST_Intersection for performance

Note that this is less exact, with intersecting pixels that cover less than 50% of an intersecting geometry being ignored.

Stefan has an answer here that avoids the ST_Intersection.


  • Is it possible to calculate raster statistics (mean value of crossed raster cells) along linestrings in a similar way?
    – Mario
    Commented Feb 18, 2021 at 16:58
  • 1
    I can't see why not, ST_Intersects will work the same with a linestring as a polygon. Commented Feb 18, 2021 at 21:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.