PostGIS raster summation (map algebra)

Question

I have a table of polygons representing travel time isochrones on particular days. For each origin point, there are five isochrone geometries (stored in separate rows). For each origin point, I want to rasterise the five isochrones (a binary NULL or 1), and then combine them into a single raster layer. This raster layer requires a simple map algebra: sum/5, so that each origin will in the end be associated with a single raster layer that has values in [NULL, 0.2, 0.4, 0.6, 0.8, 1] depending on how many of the constituent layers overlap. It is a probability surface.

My data is all stored in Postgres 9.3 (with PostGIS). My problem is that while I want to learn to use PostGIS raster, it seems to have a really steep learning curve, and all of the examples I can find deal with a single raster layer. In the examples, this layer is used as part of a polygon overlay, perhaps averaging the value of the raster for each polygon. I haven't found a replicable example for combining: a) vector --> raster b) map algebra; and c) GROUP BY attribute as per my first paragraph.

I'm OK using GDAL or GRASS if I have to in order to perform this task, but this seems like something PostGIS should be able to handle; it would be convenient to do so given my input data is already PostGIS geometry; and I really want to come to terms with PostGIS raster.

Some sample data structure:

areaid    time        date          isogeom (polygon)
1000      07:15:00    2014-05-05    xxx
1000      07:15:00    2014-05-06    xxy
...
1006      07:15:00    2014-05-05    zzz

I want to rasterise, group by areaid, and then perform the map algebra to come to:

areaid    isorast (raster)
1000      aaa
1006      bbb

---

I have not been successful containing this to PostGIS. My approach has been to convert the vector to raster, dump the rasters to arrays, and perform the combination with numpy arrays via psycopg2, before writing them to a GeoTIFF (to maybe be put back in PostGIS). Not ideal, but do-able.

Answer 1

You will need to write your own aggregate function:

CREATE OR REPLACE function sum_raster_state(raster, raster)
returns raster
language sql
as $f$

SELECT
CASE $1
WHEN NULL THEN
      $2
ELSE
    ST_MapAlgebra(
        $1, 
        $2, 
        '([rast1] + [rast2])', 
        NULL, 
        'UNION', 
        '[rast2]', 
        '[rast1]', 
        NULL)
END;
$f$;

CREATE OR REPLACE FUNCTION sum_raster_final(raster)
returns raster
language sql
as $f$
SELECT 
   $1
$f$;

create aggregate sum_raster(raster) (
    SFUNC = sum_raster_state,
    STYPE = raster,
    FINALFUNC = sum_raster_final
);

afterwards you can call it like this

SELECT areaid, sum_raster(st_asraster(isogeom, ...)) FROM isochrone GROUP BY areaid

This gives you the sum of all your rasters with the same area ID. You will still need to divide raster values by the number of observations for each area ID. (I didn't include it in the aggregate function. You can either do it here or afterwards using MapAlegbra)

Make sure your input rasters are all aligned, otherwise this won't work.