[postgis-users] Announcement: PostGIS Raster now support map algebra on two rasters!

[Pierre Racine]

Hi all PostGIS Raster fans,

Thanks to Bborie Park form the Center for Vectorborne Diseases at the University of California, PostGIS Raster now supports map algebra operations on two rasters. The new function, along with its sister working on one raster, grants PostGIS with one of the most flexible map algebra operator available today in any GIS package.

ST_MapAlgebra(raster, raster, expression) -> raster...

-compute a new raster from the expression of two input raster,
-the expression is evaluated by the PostgreSQL parser so that any valid SQL expression can be used including conditional expressions. an upcoming variant will allow you to refer to your own custom function instead of using an expression,
-can compute the new raster over the extent of the first input raster, the second inout raster, the intersection of the two rasters or their union,
-accept alternative expressions to deal with nodata value occurring in the first, the second or both raster.

This major milestone which will be part of PostGIS 2.0 is the foundation for many future cool developments:

-ST_Intersection(raster, raster) and all other overlay operation like ST_Union(), ST_Difference(), ST_SymDifference()
-ST_Clip(raster, geometry) (which is essentially the same as ST_Intersection(raster, ST_AsRaster(geometry)))
-ST_BurnToRaster(raster, geometry)

You can obviously build your own raster model from an expression of many raster layers by imbricating many ST_Mapalgebra calls. e.g to implement r1 = (r2 + r3) / r4 you can do 

SELECT ST_MapAlgebra(rasttemp.rast, r4.rast, "rast1/rast2") AS r1
FROM (
             SELECT ST_MapAlgebra(r2.rast, r3.rast, "rast1 + rast2") AS rast
             FROM r2, r3
             WHERE ST_Intersects(r2.rast, r3.rast)) rasttemp,
             r4
WHERE ST_Intersects(rasttemp.rast, r4.rast))

Feel free to be creative and suggest your own plpgsql or C code to get map algebra to work!

Thanks,

Pierre