Sunday 2 February 2020

PostGIS polygon clip inaccurate using ST_Clip and ST_Union from tiled raster


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?




Answer



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;

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...