Opened 3 years ago

Closed 3 weeks ago

Last modified 3 weeks ago

#3457 closed defect (fixed)

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

Reported by: GMudambi Owned by: komzpa
Priority: critical Milestone: PostGIS 2.4.6
Component: raster Version: 2.1.x
Keywords: Cc:

Description

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

Change History (16)

comment:1 Changed 3 years ago by robe

Component: postgisraster
Owner: changed from pramsey to Bborie Park

comment:2 Changed 14 months ago by pramsey

Milestone: PostGIS 2.1.9PostGIS 2.2.6

comment:3 Changed 13 months ago by pramsey

Milestone: PostGIS 2.2.6PostGIS 2.2.7

comment:4 Changed 8 months ago by robe

Milestone: PostGIS 2.2.7PostGIS 2.5.0

comment:5 Changed 4 months ago by robe

Milestone: PostGIS 2.5.0PostGIS 2.4.5

comment:6 Changed 2 months ago by 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ří.

comment:7 Changed 2 months ago by pramsey

Milestone: PostGIS 2.4.5PostGIS 2.4.6

comment:8 Changed 3 weeks ago by 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

comment:9 Changed 3 weeks ago by komzpa

Owner: changed from Bborie Park to komzpa

comment:10 Changed 3 weeks ago by komzpa

Resolution: fixed
Status: newclosed

In 16940:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

comment:11 Changed 3 weeks ago by komzpa

In 16941:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

comment:12 Changed 3 weeks ago by komzpa

In 16942:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

comment:13 Changed 3 weeks ago by komzpa

In 16943:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

comment:14 Changed 3 weeks ago by komzpa

In 16944:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

comment:15 Changed 3 weeks ago by komzpa

In 16945:

[raster] Fix envelope Contains shortcut in ST_Clip

Patch by Sai-bot

Closes #3457

comment:16 Changed 3 weeks ago by komzpa

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

Note: See TracTickets for help on using tickets.