[postgis-users] ST_SummaryStats

Simon SPDBA Greener simon at spdba.com.au
Sun Jul 2 16:26:17 PDT 2023


Following the documentation I have crafted some GDAL to load, and SQL to 
query, some rasters.

I note that :

1. If, when loading using GDAL, I don't tile the rasters I get a table 
with a single row (toddriver_2019_2p5m)

2. If, when loading using GDAL, I specify a tile size (eg -t 256*256) I 
get a table with multiple rows (toddriver_2019_3p0m)

Now if I execute ST_SummaryStats against 1 I get a single row which 
looks correct.:

WITH data AS (
   SELECT '1' as band,
          'SRID=4283;POLYGON ((133.84747311077376 -23.74668364533433, 
133.84747311077376 -23.7409001580403, 133.85607301067603 
-23.7409001580403, 133.85607301067603 -23.746683645334333, 
133.84747311077376 -23.74668364533433))'::geometry as geom
)
SELECT (stats).sum   as sum,
        (stats).mean  as mean,
        (stats).stddev as stddev,
        (stats).min   as min,
        (stats).max   as max,
        (stats).count as count
    FROM (SELECT 
ST_SummaryStats(ST_Clip(ST_Band(raster,d.band),geom,true)) As stats
            FROM gis.toddriver_2019_2p5m
                 INNER JOIN data as d
                 ON ST_Intersects(d.geom,ST_Band(raster,d.band))
         ) As foo;

If I execute the above against 2) I get multiple summary statistics 
rows, one for each of the tiles in the table.

To return the summary stats for 2 across all rows I did the following:

WITH data AS (
   SELECT '1' as band,
          'SRID=4283;POLYGON ((133.84747311077376 -23.74668364533433, 
133.84747311077376 -23.7409001580403, 133.85607301067603 
-23.7409001580403, 133.85607301067603 -23.746683645334333, 
133.84747311077376 -23.74668364533433))'::geometry as geom
)
SELECT sum((stats).sum)   as sum,
        avg((stats).mean)  as mean,
        sqrt(sum((stats).stddev*(stats).stddev)/min(groups)) as stddev,
        min((stats).min)   as min,
        max((stats).max)   as max,
        sum((stats).count) as count
   FROM (SELECT row_number() over () as groups,
ST_SummaryStats(ST_Clip(ST_Band(raster,d.band),geom,true)) As stats
           FROM gis.toddriver_2019_3p0m as t
                INNER JOIN data as d
                ON ST_Intersects(d.geom,ST_Band(raster,d.band))
        ) As f;

I have tested this SQL against tiled and untiled versions of the same raster and the result looks correct.
Even so, I could like confirmed that this approach is correct for tiled rasters?

regards
Simon

-- 
Simon Greener
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia
(m) +61 418 396 391
(w)www.spdba.com.au
(m)simon at spdba.com.au	
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20230703/1a109607/attachment.htm>


More information about the postgis-users mailing list