Calculating mean value of polygon from raster in PostGIS?

Question

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:

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

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.

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

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

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.

Notes

• You might also find some useful tips in this question and PostGIS WKT raster tutorial.
• If your raster is tiled, a good rule of thumb is to use smaller tiles, trying to get them a little bit bigger than a typical vector feature that you're intersecting.