I have written the code below to select specific polygons and determine which grid tiles they intersect via a lookup table. Then depending on the grid tiles I want to merge the associated raster tiles from a single band chunked master raster file.

DROP TABLE dem_stats;
--  Select Features using lookup table and determine the raster tiles said features are intersecting
    feat AS 
        (SELECT title_no,
                grid_tile_name || '.asc' AS tile_name,
                FROM parcels_all_shapefile_building a, parcels_all_shapefile_lookup_osgb_grid_5km b
                    WHERE a.title_no = b.olp_title_no AND
                          a.title_no = 'AGL253216'
-- Merge rasters tiles from main raster file that intersect features        
    merged_rast AS
        (SELECT filename, ST_Union(rast,1) AS rast
            FROM dem, feat
                WHERE filename
                    IN (tile_name)
                        GROUP BY filename
-- As the tiles are now merged duplicates are not required
    feat_temp AS 
        (SELECT DISTINCT ON (title_no) * FROM feat
-- Clip merged raster and obtain pixel statistics
    b_stats AS
        (SELECT title_no, (stats).*
            FROM (SELECT title_no, ST_SummaryStats(ST_Clip(a.rast,1,b.wkb_geometry,-9999,true)) AS stats
                FROM merged_rast a
                    INNER JOIN feat_temp b
                ON ST_Intersects(b.wkb_geometry,a.rast)
            ) AS foo
--  Summarise statistics for each title number
    SELECT title_no,
           count As pixel_val_count,
           min AS pixel_val_min,
           max AS pixel_val_max,
           mean AS pixel_val_mean,
           stddev AS pixel_val_stddev
            FROM b_stats
                WHERE count > 0;

Everything works except that instead of a single row result I am getting two as seen below:

enter image description here

So it looks like the raster tiles are not being merged correctly and the pixel statistics are being calculated for each raster tile. Does anyone know what I am doing incorrectly?


I figured it out, I shouldn't have been using GROUP BY, here is the replacement code for anyone wanting to do this in the future:

merged_rast AS
    (SELECT ST_Union(rast,1) AS rast
        FROM dem, feat
            WHERE filename
                IN (tile_name)*

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.