2

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. enter image description here

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

2

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)))
  • Thanks, but it seems expensive to union the whole raster and create a new table especially if I am clipping <1% of area. In my original code I do transform polygon to raster projection, and union. Id like to try keep it in one query if possible. – Ben Aug 28 '15 at 20:52
  • what's dimensions of your raster and is your raster is tiled? – elmo Aug 28 '15 at 22:16
  • Pansharpened Landsat8 scene so X: 15681 Y: 15841 Bands: 4, tiled at 128x128 – Ben Aug 28 '15 at 22:18
  • so I think that you need to union tiles wich are intersected with your polygon, then clip your raster results on your polygon – elmo Aug 28 '15 at 22:31
  • Yes I was thinking about trying that. Cheers – Ben Aug 28 '15 at 23:03
1

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

  • Cheers - hope it is fixed soon. – Ben Apr 8 '16 at 3:19
1

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)

polygon to clip

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 ;

union_variant

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);

tiled_variant_I

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);

tiled_variant_II enter image description here

  • Looks like we stumbled upon the same problem - my post highlighting this here: gis.stackexchange.com/questions/294482/… (took some digging to realize it was the st_clip operation that was not working as expected) – mtreg Sep 5 '18 at 18:16
0

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;

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.