[postgis-users] Raster Points in Polygon (was Re: PostGIS Raster for Met Data)

Pierre Racine Pierre.Racine at sbf.ulaval.ca
Fri Aug 26 06:08:00 PDT 2011


> st_intersect(raster, raster) sounds very interesting, but what would be the result
> from the function?

To extract summary statistics of raster values included inside a polygon we need ST_Intersection(raster, raster) so that you could make

SELECT ST_SummaryStats(ST_Intersection(rast, ST_AsRaster(geom,rast))) stats
FROM yourrastertable rast, yourpolygon geom
WHERE ST_Intersects(rast, geom)

ST_Intersects(raster, raster) is not necessary for this.

ST_AsRaster(geom,rast) rasterize your geometry to a raster having the same grid definition, srid, pixeltype as rast.

ST_Intersection(rast, ST_AsRaster(geom,rast)) return a two band raster with values only in the withvalue shared area.

ST_SummaryStats(ST_Intersection(rast, ST_AsRaster(geom,rast))) compute summary statistics (count, sum, mean, min, max, sdtdev) of the resulting raster.

This should be faster than the technique I showed you where we vectorize the raster to points and do a vectorial ST_Intersection() because:

-ST_Intersection(raster, raster) works in raster mode which should be faster than ST_Intersection(raster, geometry)

-ST_AsRaster(geom, raster) should theoretically be faster than ST_PixelAsPoints(), ST_PixelAsPolygon() and ST_DumpAsPolygon() which are or can be used to vectorize a raster to a geometry table on the fly (any comparative tests done with ST_DumpAsPolygon() ?).

The missing key in PostGIS Raster to do that is ST_MapAlgebra(raster, raster) which enable mapalgebra on the intersection of both raster. All this is already possible in pl/pgsql (see http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra.sql) but need C implementation and optimization (http://trac.osgeo.org/postgis/browser/trunk/raster/scripts/plpgsql/st_mapalgebra_optimized.sql).

For now that would only work if your raster are not tiled. But with ST_SummaryStats() working as an aggregate function (ticket #1048) that would work on tiled raster as well much more efficiently.

Pierre





More information about the postgis-users mailing list