PostGIS polygon clip inaccurate using ST_Clip and ST_Union from tiled raster

Question

I have a RGB raster tiled at 128x128 in PostGIS enabled database. I am trying to clip from the raster using a polygon as follows:

    SELECT ST_AsGDALRaster(rast,'GTiff') 
        FROM (
            SELECT ST_Union(ST_Clip(rast,geom)) AS rast
                FROM 
                    scenes.rgb
                    CROSS JOIN
                    ST_Transform(ST_GeomFromText('POLYGON((
                        148.296964462160986 -28.914219626449281, 148.300054703309883  -28.892809637280056, 148.288369276810897 -28.891374636835419,
                        148.290750573935583 -28.87714520356819, 148.305307129311188 -28.878270936973706, 148.305323220779684 -28.892336715938328, 
                        148.307563219233941 -28.892525982844003, 148.30733671875592 -28.893583124483161, 148.323105875277065 -28.895609384106994, 
                        148.319860482823429 -28.917021311138011, 148.319860482823429 -28.917021311138011, 148.296964462160986 -28.914219626449281
                    ))',4326),32655) AS geom
            WHERE st_intersects(rast,geom)) 
    AS rast;

See graphic, paddock bounds is black outline - image is meant to clip to this area.

My question is how can I change my query so the clipped raster matches the polygon?

Answer 1

You need to reproject your polygon to your raster projection, then you can clip the raster.

This is how you can perform it :

R = raster

u_R = union of raster

P = Polygon

1-Select tiles which intersects with your polygon

Create table tiles_R as SELECT * 
from R where
ST_Intersects(u_R.raster,st_transform(P.geom,st_srid(tiles_R.raster)))

2-Perform the union of your raster like :

Create table u_R as SELECT st_union(rast) as raster from tiles_R

3-Then clip your raster

SELECT ST_Clip(u_R.raster, 1, st_transform(P.geom,st_srid(u_R.raster)),false) as  raster 
FROM u_R, P
where ST_Intersects(u_R.raster,st_transform(P.geom,st_srid(u_R.raster)))

Answer 2

Seems like this is a bug in newer versions of postgis. Check out the ticket @ https://trac.osgeo.org/postgis/ticket/3457

Answer 3

I am also facing similar problems. It seems like a bug resulting from enhancement as suggested by @Koshy Thomas in https://gis.stackexchange.com/a/188515/75756. I am sending some self-explanatory screenshots, how different calls of ST_Clip behave:

• area to clip (gisdata.f_a_sampling_stratum_testing)

correct results with union of tiles made first

not possible when union of tiles is too large (not fit memory)

-- union_variant drop table if exists gisdata.f_a_sampling_stratum_testing_clip; create table gisdata.f_a_sampling_stratum_testing_clip as select 1 as rid, st_clip(st_union(rast), geom, false) from gisdata.fty_2012_020m_eu_03035_d01_full, gisdata.f_a_sampling_stratum_testing where st_intersects(rast, geom) group by geom ;

incorrect result when clipping tiles

when tile extent doesn't intersect bbox of geom, ST_Clip seems not performed

-- tiled_variant_I drop table if exists gisdata.f_a_sampling_stratum_testing_clip; create table gisdata.f_a_sampling_stratum_testing_clip as select rid, st_clip(rast, geom, false) as rast from gisdata.fty_2012_020m_eu_03035_d01_full, gisdata.f_a_sampling_stratum_testing where st_intersects(rast, geom);

incorrect result when clipping raster tiles with intersections

-- tiled_variant_II drop table if exists gisdata.f_a_sampling_stratum_testing_clip; create table gisdata.f_a_sampling_stratum_testing_clip as select rid, case when st_within(st_convexhull(rast), geom) then rast else st_clip(rast, st_intersection(st_convexhull(rast), geom), false) end as rast, st_within(st_convexhull(rast), geom) as is_within, st_intersection(st_convexhull(rast), geom)::geometry(polygon, 3035) as geom from gisdata.fty_2012_020m_eu_03035_d01_full, gisdata.f_a_sampling_stratum_testing where st_intersects(rast, geom);

Answer 4

The code in my above question does work for some queries, but is unreliable. I have ended up settling with for now this query which takes in elmo's process but without creating any new tables. I am not sold on it as it requires more work than I would like - an extra half a second to return.

SELECT ST_AsGDALRaster(rast,'GTiff') 
    FROM ( 
        SELECT ST_Clip(rast,ST_Transform(ST_GeomFromText('POLYGON((
            148.296964462160986 -28.914219626449281, 148.300054703309883  -28.892809637280056, 148.288369276810897 -28.891374636835419,
            148.290750573935583 -28.87714520356819, 148.305307129311188 -28.878270936973706, 148.305323220779684 -28.892336715938328, 
            148.307563219233941 -28.892525982844003, 148.30733671875592 -28.893583124483161, 148.323105875277065 -28.895609384106994, 
            148.319860482823429 -28.917021311138011, 148.319860482823429 -28.917021311138011, 148.296964462160986 -28.914219626449281
            ))',4326),32655)) AS rast 
        FROM ( 
            SELECT ST_Union(rast) as rast from (
                SELECT rast FROM scenes.rgb WHERE ST_Intersects(scenes.rgb.rast,
                    ST_Transform(ST_GeomFromText('POLYGON((
                        148.296964462160986 -28.914219626449281, 148.300054703309883  -28.892809637280056, 148.288369276810897 -28.891374636835419,
                        148.290750573935583 -28.87714520356819, 148.305307129311188 -28.878270936973706, 148.305323220779684 -28.892336715938328, 
                        148.307563219233941 -28.892525982844003, 148.30733671875592 -28.893583124483161, 148.323105875277065 -28.895609384106994, 
                        148.319860482823429 -28.917021311138011, 148.319860482823429 -28.917021311138011, 148.296964462160986 -28.914219626449281
                        ))',4326),32655))
            ) AS rast
        ) AS rast               
    ) AS rast;