Invalid output of ST_ClipByBox2D with an invalid polygon

[Algunenano]

I start with a valid polygon:

Select St_IsValid('SRID=3857;MULTIPOLYGON(((-8231365.02893734 4980355.83678553,-8231366.94866817 4980350.2976052,-8231368.82883367 4980344.52250402,-8231370.96255986 4980333.47364302,-8231372.09200747 4980327.69798564,-8231372.84505222 4980323.68014243,-8231374.2255656 4980316.77455645,-8231375.10554947 4980309.86865218,-8231380.6311111 4980275.9676398,-8231394.82332406 4980186.31880185,-8231394.8569833 4980133.57928934,-8231367.43081065 4979982.17443372,-8231372.85765532 4979985.17751813,-8231395.25669799 4980132.93828551,-8231396.69199339 4980227.59327083,-8231377.30092207 4980340.77515211,-8231368.35061694 4980357.62467295,-8231365.02893734 4980355.83678553)))'::geometry);
 st_isvalid 
------------
 t

I pass this polygon through ST_Simplify and ST_RemoveRepeatedPoints and it comes up as invalid:

Select St_IsValid(ST_Simplify(ST_RemoveRepeatedPoints('SRID=3857;MULTIPOL...'), 19.109257221221924), 19.109257221221924, true));
 st_isvalid 
------------
 f
(1 row)

Then I pass it through ST_ClipByBox2D:

Select St_AsEWKT(ST_ClipByBox2D(
    'MULTIPOLYGON(((-8231365.02893734 4980355.83678553,-8231394.82332406 4980186.31880185,-8231367.43081065 4979982.17443372,-8231396.69199339 4980227.59327083,-8231365.02893734 4980355.83678553)))'::geometry,
    'POLYGON((-8228215.28157305 4980103.20610541,-8228215.28157305 4970163.3879457,-8238155.09973276 4970163.3879457,-8238155.09973276 4980103.20610541,-8228215.28157305 4980103.20610541))'::geometry
));

POLYGON((-8238155.09973276 4970163.3879457,-8238155.09973276 4980103.20610541,-8231383.67109061 4980103.20610541,-8231367.43081065 4979982
.17443372,-8231381.86136498 4980103.20610541,-8228215.28157305 4980103.20610541,-8228215.28157305 4970163.3879457,-8238155.09973276 4970163
.3879457))

The result seen is not the intersection (as expected) but the result of subtracting the geometry area to the box.

In truth this series of chained calls come from St_AsMVTGeom where this functions are called inside the same C function, but seeing the individual calls makes it easier to spot the issue.

I've traced this behaviour back to GEOS (​https://github.com/Algunenano/geos/commit/1737fa1a34e6cc30fb26b2898f607d5db3b176a6).

Does anyone know if GEOSClipByRect is supposed to handle invalid geometries correctly or to just give an invalid output as the one seen above?

[Algunenano]

BTW, here is the call to St_AsMVTGeom reproducing the issue

With clipping on you get the whole extent as the geometry output:

SELECT St_AsText(ST_AsMVTGeom(
 'SRID=3857;MULTIPOLYGON(((-8231365.02893734 4980355.83678553,-8231366.94866817 4980350.2976052,-8231368.82883367 4980344.52250402,-8231370.96255986 4980333.47364302,-8231372.09200747 4980327.69798564,-8231372.84505222 4980323.68014243,-8231374.2255656 4980316.77455645,-8231375.10554947 4980309.86865218,-8231380.6311111 4980275.9676398,-8231394.82332406 4980186.31880185,-8231394.8569833 4980133.57928934,-8231367.43081065 4979982.17443372,-8231372.85765532 4979985.17751813,-8231395.25669799 4980132.93828551,-8231396.69199339 4980227.59327083,-8231377.30092207 4980340.77515211,-8231368.35061694 4980357.62467295,-8231365.02893734 4980355.83678553)))'::geometry,
 'SRID=3857;POLYGON((-8238077.16046316 4980025.2668358,-8238077.16046316 4970241.3272153,-8228293.22084265 4970241.3272153,-8228293.22084265 4980025.2668358,-8238077.16046316 4980025.2668358))'::geometry,
 256,
 1,
 true
 )) as the_geom_webmercator;
                    the_geom_webmercator                     
-------------------------------------------------------------
 MULTIPOLYGON(((175 -1,257 -1,257 257,-1 257,-1 -1,175 -1)))
(1 row)

With clipping off you get something that resembles the initial geometry (keep in mind that MVT flips the y coordinate):

SELECT St_AsText(ST_AsMVTGeom(
 'SRID=3857;MULTIPOLYGON(((-8231365.02893734 4980355.83678553,-8231366.94866817 4980350.2976052,-8231368.82883367 4980344.52250402,-8231370.96255986 4980333.47364302,-8231372.09200747 4980327.69798564,-8231372.84505222 4980323.68014243,-8231374.2255656 4980316.77455645,-8231375.10554947 4980309.86865218,-8231380.6311111 4980275.9676398,-8231394.82332406 4980186.31880185,-8231394.8569833 4980133.57928934,-8231367.43081065 4979982.17443372,-8231372.85765532 4979985.17751813,-8231395.25669799 4980132.93828551,-8231396.69199339 4980227.59327083,-8231377.30092207 4980340.77515211,-8231368.35061694 4980357.62467295,-8231365.02893734 4980355.83678553)))'::geometry,
 'SRID=3857;POLYGON((-8238077.16046316 4980025.2668358,-8238077.16046316 4970241.3272153,-8228293.22084265 4970241.3272153,-8228293.22084265 4980025.2668358,-8238077.16046316 4980025.2668358))'::geometry,
 256,
 1,
 false
 )) as the_geom_webmercator;
                                                                                     the_geom_webmercator                                  
                                                   
-------------------------------------------------------------------------------------------------------------------------------------------
---------------------------------------------------
 MULTIPOLYGON(((175.090909090909 -4.45454545454545,175 -5,176 -9,175.090909090909 -4.45454545454545)),((175.090909090909 -4.45454545454545,
176 1,175 -4,175.090909090909 -4.45454545454545)))

[komzpa]

Does anyone know if GEOSClipByRect is supposed to handle invalid geometries correctly or to just give an invalid output as the one seen above?

I would expect GEOSClipByRect to return correct part of geometry, but not resolving the ambiguities inside it. Basically, I'd expect that an assert that result's box is within the clipping box, to hold.

[Algunenano]

I would expect GEOSClipByRect to return correct part of geometry, but not resolving the ambiguities inside it. Basically, I'd expect that an assert that result's box is within the clipping box, to hold.

Following this way of thinking, I would expect the bbox of the result to be both within the geometry bbox and the clipping box.

[Algunenano]

I've been going deeper into this issue and ended up reviewing the code in Geos that does the clipping.

Here is the geometry that exemplifies this issue: 

If you take the order used by the orange arrow (right hand order) then when the clipping happens you'll get the correct part of the geometry; on the other hand, if you get the green arrow (left hand order) you will get the opposite (the whole box minus that part of the geometry). The issue comes from the fact that this geometry is invalid so the part of the geometry getting clipped is in the wrong winding order (left hand).

I've been thinking about how to fix it and the way I've come up with is to get the points where the geometries cross and check whether the intermediate point (or the corner of the box if the points are in different sides of the box) is inside the geometry or not. Doing this involves a Point in Polygon check which is not fast, and forcing this check for all geometries (even the ones that are valid) would mean a big hit in performance in a function whose main purpose is to be a fast and dirty alternative to St_Intersection when the second geometry is a box.

With this in mind, I've decided to explore an alternative solution: check if the resulting geometry's bbox is inside the bbox of the geometry (as it should normally) and, if not, try to fix the geometry and relaunch again. This should have a negligible impact in performance for the most common case (valid input) and fix this issue.

My main concern is whether to do this only for St_AsMVTGeom or in lwgeom_clip_by_rect itself If done only in St_AsMVTGeom we could, instead of fixing the geometry, just drop it (as done by Mapnik); and if done in lwgeom_clip_by_rect it would help everyone relying on that function (St_SnaptoGrid) with the caveat of hiding geometry changes to the user, which is something that was discussed during 2.5 and postponed to 3.0.

I'm open to ideas, but I think that for now (and backporting to 2.4 and 2.5) I'd take the easy way out (drop that geometry in St_AsMVTGeom) and, whenever some decision is taken as to how to handle invalid geometries inside Postgis, adapt St_AsMVTGeom (or lwgeom_clip_by_rect) do follow that behaviour.

[Algunenano]

WIP in ​https://github.com/postgis/postgis/commit/38cfcd4767075f806f32d780973be448c1e67e70

[Algunenano]

Fix in trunk in r16856.

[Raul Marin]

In 16857:

St_AsMVTGeom: Drop invalid geometries after simplification

References #4183

[Raul Marin]

In 16858:

St_AsMVTGeom: Drop invalid geometries after simplification

Closes #4183