10

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.

Raster with Polygons, viewed in QGIS

4
8

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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