[postgis-devel] ST_MapAlgebraExpr(): Raster aggregates?

[Pierre Racine]

tom,

> Nevertheless, I think we got the correct way to do it. Like you describe it for
> the second variant, both the raster and the vectors cover the whole country
> and the vectors have no gap.

Yes. Up to now I guess this is the better way to do it. 

Do we agree that if ST_SummaryStats() would be an aggregate (computing global summary for a set of raster) the ST_Union() call would be useless?See http://trac.osgeo.org/postgis/ticket/1048

It seems to me that the ST_MapAlgebraExpr() alternate expressions should be left to null (ST_MapAlgebraExpr(ST_AsRaster(a.geom, b.rast, '32BF') ,b.rast ,'rast2', '32BF','INTERSECTION')) so that pixel outside the geometry get assigned nodata instead of 0...

Why do you surround ST_Union with parenthesis and take only the raster part like if the return value was a geomval? (ST_Union()).rast the return value of ST_Union is a raster so you should not have to do that.

> The thing I'm not sure about though: at the moment we made use of
> roughly 20 X 20 km raster tiles, each overlain by a couple of hundred
> polygons. Do you think it makes sense to make the raster tiles smaller like
> you recommend in the tutorial? It seems to me that I would have to do a
> lot more ST_Union on the rasters that share the same polygon. By the way:
> the polygons vary greatly in size and shape (postal codes).

I have the intuition that it should not make much difference. I haven't experimented with this technique yet so I might be wrong. Need testing...

Pierre

> See for the code below.
> 
> TIA.
>   Tom
> 
> 
> CODE:
> --------------------------------
> WITH
> postcodes As (
>      SELECT * FROM pc6_nl
> )
> 
> SELECT
> a.gid, a.pc6,a.geom,
> (ST_SummaryStats(
>      (ST_Union(
>          ST_MapAlgebraExpr(
>              ST_AsRaster(a.geom, b.rast, '32BF')
>              ,b.rast
>              ,'rast2', '32BF','INTERSECTION','0','0',0
>          )
>      )).rast
> )).mean As avg_height
> 
> FROM
> postcodes a LEFT JOIN
> rasters.ahn2 b
>      ON ST_Intersects(a.geom, b.rast)
> GROUP BY a.gid, a.pc6, a.geom