Average of multiple rasters in PostGIS raster?

Question

I am trying to calculate the average of multiple rasters stored as raster datatype in PostGIS enabled database. The rasters are stored in three different tables (spi3, sri3, and smdi) and have date as common column. I am trying to use ST_Union function of PostGIS, where it says we can specify setof raster.

This is what I have done so far.

select st_union(SETOF[spi3_raster, sri3_raster, smdi_raster], 'MEAN')
from
(
select spi3.raster as spi3_raster, sri3.raster as sri3_raster, smdi.raster as smdi_raster
from spi3 as spi3
join sri3 as sri3
  on spi3.date = sri3.date
  and sri3.date = '2013-01-01'
join smdi as smdi
  on sri3.date = smdi.date
where spi3.date = '2013-01-01'
) as foo;

However, I am not getting the result as expected.

Answer 1

There may be another way, but you can use UNION SELECT to stack rasters, as you would to combine any other type into a single field containing various rows. So, including the where condition in each select, rather than as a join, you would have

SELECT ST_Union (foo.raster, 'mean') 
  FROM (
     SELECT raster FROM spi3 WHERE date='2013-01-01'
      UNION 
     SELECT raster FROM sri3 WHERE date='2013-01-01'
      UNION 
     SELECT raster FROM smdi WHERE date='2013-01-01'
    ) foo;

I suspect that there might be a more elegant, but not necessarily more efficient, way to do this.

Answer 2

I had a similar problem, trying to get the standard deviation of all of the values across a selected set of raster layers. Most existing functions seem to be designed for a single layer. I couldn't use @John Powell's solution because ST_Union() doesn't have a standard deviation function. My solution uses ST_DumpValues() and it should work for the mean (or any other function) as well:

WITH rasts as (SELECT rast from mytable WHERE ...),
vals as (SELECT unnest(ST_DumpValues(rast,1,TRUE)) as val from rasts rs)
SELECT count(val), stddev(val) from vals WHERE val IS NOT NULL

ST_DumpValues() returns a record that is a composite type unless you specify the band, as I have done here, and you can alternatively use composite type grammar to get at the elements: unnest((st_dumpvalues(rs.rast)).valarray)

Also, I learned from another answer by John Powell that exclude_nodata_value=true, which is the default in ST_DumpValues(), "doesn't mean return no data, it means return null instead of a numeric value representing no data, such as -9999;" which is why I added WHERE val IS NOT NULL

I haven't tried this solution on a big problem so I don't know how efficient it is.