When a geometry is smaller than a raster but their extents are the same, one of the overloads of ST_Clip doesn't give expected output

[GMudambi]

Consider a square geometry

ST_GeomFromText('POLYGON((0 0,100 0,100 100,0 100,0 0))',4326)

Let us create a raster out of it. ST_ASRaster( ST_GeomFromText('POLYGON((0 0,100 0,100 100,0 100,0 0))',4326), ST_Addband(ST_MakeEmptyRaster(1,1,0,0,1,-1,0,0,4326),'32BF'::text,0,-1),'32BF'::text,1,-1) ——- {1}

Consider another geometry whose extent (Envelope) is same the raster but is actually smaller. For Ex:

ST_GeomFromText('POLYGON((0 0,100 100,100 0,0 0))',4326) ——- {2}

Clipping {1} and {2} and dumping the raster as polygons should give us a polygon closer to a triangle but gives us a square output.

This behaves right in 2.1.3 version but has got changed from 2.1.4

Query:

SELECT ST_AsText((ST_DumpAsPolygons(ST_Clip(ST_ASRaster(ST_GeomFromText('POLYGON((0 0,100 0,100 100,0 100,0 0))',4326),ST_Addband(ST_MakeEmptyRaster(1,1,0,0,1,-1,0,0,4326),'32BF'::text,0,-1),'32BF'::text,1,-1), ST_GeomFromText('POLYGON((0 0,100 100,100 0,0 0))',4326)))).geom)

PostGis v2.1.3 POLYGON((99 99,99 98,98 98,98 97,97 97,97 96,96 96,96 95,95 95,95 94,94 94,94 93,93 93,93 92,92 92,92 91,91 91,91 90,90 90,90 89,89 89,89 88,88 88,88 87,87 87,87 86,86 86,86 85,85 85,85 84,84 84,84 83,83 83,83 82,82 82,82 81,81 81,81 80,80 80,80 79,79 79,79 78,78 78,78 77,77 77,77 76,76 76,76 75,75 75,75 74,74 74,74 73,73 73,73 72,72 72,72 71,71 71,71 70,70 70,70 69,69 69,69 68,68 68,68 67,67 67,67 66,66 66,66 65,65 65,65 64,64 64,64 63,63 63,63 62,62 62,62 61,61 61,61 60,60 60,60 59,59 59,59 58,58 58,58 57,57 57,57 56,56 56,56 55,55 55,55 54,54 54,54 53,53 53,53 52,52 52,52 51,51 51,51 50,50 50,50 49,49 49,49 48,48 48,48 47,47 47,47 46,46 46,46 45,45 45,45 44,44 44,44 43,43 43,43 42,42 42,42 41,41 41,41 40,40 40,40 39,39 39,39 38,38 38,38 37,37 37,37 36,36 36,36 35,35 35,35 34,34 34,34 33,33 33,33 32,32 32,32 31,31 31,31 30,30 30,30 29,29 29,29 28,28 28,28 27,27 27,27 26,26 26,26 25,25 25,25 24,24 24,24 23,23 23,23 22,22 22,22 21,21 21,21 20,20 20,20 19,19 19,19 18,18 18,18 17,17 17,17 16,16 16,16 15,15 15,15 14,14 14,14 13,13 13,13 12,12 12,12 11,11 11,11 10,10 10,10 9,9 9,9 8,8 8,8 7,7 7,7 6,6 6,6 5,5 5,5 4,4 4,4 3,3 3,3 2,2 2,2 1,1 1,1 0,99 0,100 0,100 99,99 99))

PostGis v2.1.4+ POLYGON((0 100,0 0,100 0,100 100,0 100))

The issue seems to stem from the fact that in ST_Clip function, we have a short-circuit SQL.

— short-cut if geometry's extent fully contains raster's extent IF (nodataval IS NULL OR array_length(nodataval, 1) < 1) AND geom ~ ST_Envelope(rast) THEN

RETURN rast;

END IF;

geom ~ ST_Envelope(rast) is a necessary but not a sufficient condition for the short cut.

One of the sufficiency condition here could be ST_Contains(geom, ST_Envelope(rast))

which ensures that the rast is completely enveloped by the geometry (not just its bounding box).

The change occurred as part of this enhancement ticket ​https://trac.osgeo.org/postgis/ticket/2829

[jurafejfar]

Dear developers,

I created some SQL code with screenshots showing strange behavior of ST_Clip when processing tiles:

​https://gis.stackexchange.com/questions/159828/postgis-polygon-clip-inaccurate-using-st-clip-and-st-union-from-tiled-raster/294860#294860

It seems to me, that it is connected to this bug.

Best regards, Jiří.

[saibot]

Dear Jiří, users and developers,

the st_clip function works incorrectly. I strongly recommend the following correction to all users. I hope that the bug will be fixed soon.

CREATE OR REPLACE FUNCTION public.st_clip(
  rast      RASTER, nband INTEGER [], geom GEOMETRY,
  nodataval DOUBLE PRECISION [] DEFAULT NULL :: DOUBLE PRECISION [],
  crop      BOOLEAN DEFAULT TRUE)
  RETURNS RASTER AS
$BODY$
BEGIN

-- wrong version
  -- short-cut if geometry's extent fully contains raster's extent  
  -- IF (nodataval IS NULL OR array_length(nodataval, 1) < 1) AND geom ~ ST_Envelope(rast) 
-- corrected version (bug-fix)
  -- no more short-cut!!!  if geometry's extent fully contains raster's extent  
  IF (nodataval IS NULL OR array_length(nodataval, 1) < 1) AND st_contains(geom, ST_Envelope(rast))
  THEN RETURN rast;
  END IF;
  
  RETURN PUBLIC._ST_Clip($1, $2, $3, $4, $5);
END;
$BODY$
LANGUAGE plpgsql
IMMUTABLE
COST 100;

Thank you for all your work and best regards,

Sai-bot

[komzpa]

In 16940:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

[komzpa]

In 16941:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

[komzpa]

In 16942:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

[komzpa]

In 16943:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

[komzpa]

In 16944:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

[komzpa]

In 16945:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

[komzpa]

Thank you for the patch. It should be available with next release of PostGIS.