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

[Tom van Tilburg]

Pierre,

Very nice tutorial you wrote, wish I had seen it earlier (for sure I 
didn't look well enough).
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.
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).

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



On 24-11-2011 16:42, Pierre Racine wrote:
> Tom,
>
>> An image of just over 4 million pixels with 11 discrete raster values
>> takes around 1 minute on my 2.3 Ghz, 6Gig. mem. laptop.
>> I can sent the corresponding SQL if you like.
> I would be glad to have a look at your SQL.
>
>> 2) Glueing together tiled rasters
>> Not sure if this is the right tool for it though, but we are calculating
>> average pixel values within an overlaying polygon. When a polygon
>> overlaps 2 or more raster tiles, we need to union the tiles first to
>> make a proper average.
>> The complete calculation works as follows:
>>
>> (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
>>
>> Doing this calculation with 600,000 polygons over a 5m grid for all of
>> the Netherlands takes over 8 hours.
> Did you try the vector approach to do that? Like in the tutorial:
>
> http://trac.osgeo.org/postgis/wiki/WKTRasterTutorial01
>
> The rules to choose one approach over the other is:
>
> -If your vector coverage is sparse the vector method is generally faster (your raster coverage has to be tiled small)
>
> -If your vector coverage has no gap (like a grid covering your study area), the technique you used is the right one.
>
> Today I will post a function helping to summarize the average and other weighted stats when applying the first technique.
>
> Pierre
> _______________________________________________
> postgis-devel mailing list
> postgis-devel at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-devel