[postgis-users] summarizing a polygon values in a raster

[Pierre Racine]

> I would like to create a raster of an area-weighted average value derived from a
> polygon layer. The raster template is a 2km grid and the polygon layer is also a
> grid but at 10km. If some one could recommend a (high level) strategy for this if
> be grateful.

1) Create a vector grid having a x and a y like this:

CREATE TABLE yourgrid AS
SELECT (gvxy).geom geom, (gvxy).x x,  (gvxy).y y
FROM (SELECT ST_PixelsAsPolygons(rast) gvxy
             FROM (SELECT ST_MakeEmptyRaster(your own parameters here covering your polygon extent)) foo1) foo2

2) Make sure your polygon layer is indexed

CREATE INDEX yourpolylayer_geom_idx  ON yourpolylayer USING gist  (geom );

3) Intersect mygrid with your polygon layer and generate summary stats at the same time:

CREATE TABLE weightedvectorgrid AS
SELECT g.x, g.y, g.geom, (aws).weightedmean mean
FROM (SELECT ST_AreaWeightedSummaryStats(intgeom, p.value) aws, g.x, g.y, g.geom
      FROM (SELECT ST_Intersection(g.geom, p.geom) intgeom, p.value, g.x, g.y, g.geom
            FROM yourgrid g, yourpolylayer p
            WHERE ST_Intersects(g.geom, p.geom)
           ) foo
      GROUP BY g.x, g.y, g.geom
     ) foo2

The ST_AreaWeightedSummaryStats() function is not part of the PostGIS release but is available as part of the source tree:

http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_areaweightedsummarystats.sql

4) Union all the geometry as a raster (this one might be a bit slow depending on the size of your grid)

CREATE TABLE weightedrastergrid AS
SELECT ST_Union(ST_AsRaster(geom, 'MEAN')) rast
FROM weightedvectorgrid

I wrote all these without testing so you might have some adjustment to do.

Let me know about your progress...

Pierre