[GIS] Clipping tiled raster with polygon using PostGIS

clippostgisrasterunion

I have a table in postgis with a tiled raster and a table with polygons overlaying the raster. I want to get the number of pixels for each category (pixel value) that fall inside a specific polygon.

If I use the following query the result I am getting is wrong because, from what I understand, it hasn't unioned the tiles of the raster:

SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Clip(myraster.rast, myvector.geom)) As pvc
  FROM rastertable As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

With the following query the result I am getting is the right one, but it takes longer than it should, since it unions all the tiles before clipping

WITH
  unionedraster AS (SELECT ST_Union(rast) As rast FROM rastertable)   
SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Clip(myraster.rast, myvector.geom)) As pvc
  FROM unionedraster As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

I thought the following would work, but the results are the same with the first query (like no union took part).

SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Union(ST_Clip(myraster.rast, myvector.geom))) As pvc
  FROM rastertable As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

What is the correct way to do this?

Upon further investigation, it seems that the problem lies in the clipping of the raster tiles. I managed to export the result of the clipping as an image using this query:

COPY (
WITH clipraster As (
    SELECT ST_Clip(myraster.rast, myvector.geom) As pvc
    FROM rastertable As myraster, vectortable As myvector
    WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
)
SELECT encode(ST_AsPNG(
    clipraster.pvc
), 'hex') AS png FROM clipraster
) TO 'clip.hex';

and the resulting image was the following

enter image description here

The blue line is the polygon to which the raster tiles should be clipped.

However, it is obvious that some of the tiles where not clipped.

Any idea to what is the problem?

Best Answer

This thread is a little old. I thought it would be good to share my method for solving the issue presented. I have GEOS 3.8.0 installed so the ST_Union is pretty snappy. I haven't tried testing it on large datasets.

I use a function to determine the bounding tiles, merge those into one raster, and then apply a mask with a buffer returning the masked raster in the raster's original format. The code is specific to my current application and needs to be made more general, e.g. pulling rasters from a different table (not just ned_13), taking the raster table's SRID, and selecting the correct UTM zone's SRID based on the geometry's location.

CREATE OR REPLACE FUNCTION
ST_GeoMask ( geom GEOMETRY, buf DOUBLE PRECISION DEFAULT 60, 
                    srid INTEGER DEFAULT 6345 )
RETURNS RASTER AS
$$
    DECLARE
        mask GEOMETRY := ST_Transform (
            ST_Buffer (
                ST_Transform( geom, srid ), 
                buf -- use rounded edge buffer the default
            ),
            4269 -- mask in NAD83 coordinates EPSG::4269
        );
        rasterout RASTER;
    BEGIN
        SELECT
            ST_Clip ( ST_Union ( rast ), mask, TRUE ) INTO rasterout
        FROM public.ned_13
        WHERE ST_Intersects( rast, geom );
        RETURN rasterout;
    END;  
$$
LANGUAGE plpgsql RETURNS NULL ON NULL INPUT;
Related Question