Simplify method provides wrong result

[kamad3]

In some random cases the gdal.Geometry.Simplify() which uses the GEOS Simplify() provides wrong results. I tested it both with geos 3.7.2 and 3.8.0, the results were the same. The details from the gdal guys can be found here: ​https://github.com/OSGeo/gdal/issues/2044

EXAMPLE

The original polygon is:

POLYGON ((21.44786 47.65323,21.38886 47.66023,21.37386 47.65123,21.36586 47.66323,21.34886 47.66423,21.33286 47.65623,21.32186 47.66123,21.30886 47.65923,21.30186 47.66623,21.30786 47.66823,21.30286 47.67523,21.26786 47.67823,21.25586 47.69623,21.23886 47.70023,21.25186 47.70123,21.25186 47.70523,21.24486 47.71023,21.23286 47.70723,21.23286 47.71123,21.25186 47.71523,21.23886 47.72223,21.24186 47.72523,21.23586 47.73023,21.24286 47.73023,21.23586 47.73123,21.23686 47.73523,21.24986 47.73423,21.24486 47.73323,21.24986 47.72423,21.26586 47.71823,21.25886 47.72923,21.26086 47.73723,21.24986 47.74023,21.25986 47.74823,21.25386 47.75123,21.27786 47.75423,21.28086 47.76323,21.26786 47.76123,21.29286 47.77723,21.28086 47.78123,21.26086 47.77623,21.28986 47.78523,21.29486 47.78123,21.29586 47.78623,21.32086 47.77123,21.32686 47.78723,21.31686 47.79023,21.31186 47.80023,21.32386 47.79023,21.34186 47.79223,21.34886 47.78723,21.34686 47.78323,21.35686 47.78223,21.33886 47.79923,21.32886 47.80523,21.32186 47.80223,21.31486 47.81023,21.32786 47.81123,21.33986 47.80223,21.33886 47.81123,21.32686 47.82023,21.32586 47.82723,21.32786 47.82323,21.33886 47.82623,21.34186 47.82123,21.36386 47.82223,21.36786 47.81823,21.36086 47.81423,21.37386 47.81423,21.37386 47.81023,21.37486 47.81423,21.38686 47.81423,21.39486 47.80723,21.39586 47.81423,21.40686 47.81723,21.42286 47.80823,21.42286 47.81323,21.43586 47.81623,21.44486 47.80723,21.44486 47.81723,21.44986 47.80523,21.46386 47.81023,21.46186 47.80323,21.47986 47.80723,21.48186 47.80223,21.47686 47.80123,21.49486 47.80123,21.50486 47.78523,21.48686 47.78023,21.49386 47.77023,21.50686 47.77223,21.49986 47.76623,21.49686 47.75123,21.50486 47.75323,21.50686 47.74623,21.50986 47.75223,21.50886 47.74123,21.51286 47.74623,21.51986 47.74223,21.51386 47.72323,21.52586 47.71723,21.52586 47.71123,21.51886 47.70823,21.51986 47.71223,21.51386 47.71223,21.50886 47.70723,21.50686 47.69223,21.49686 47.69023,21.50286 47.68823,21.50186 47.67323,21.48586 47.66623,21.48386 47.65623,21.45386 47.65423,21.45686 47.65923,21.44786 47.65323))

This polygon with Simplify(0.0036) provides the following wrong polygon:

POLYGON ((21.3280682018927 47.8232867823344,21.32786 47.82323,21.32586 47.82723,21.3280682018927 47.8232867823344))

I checked with QGIS, it generates a proper one:

Polygon ((21.44785999999999859 47.65323000000000064, 21.38886000000000109 47.66022999999999854, 21.37386000000000053 47.6512299999999982, 21.36586000000000141 47.66322999999999865, 21.34885999999999839 47.66423000000000343, 21.33286000000000016 47.65623000000000076, 21.30885999999999925 47.65923000000000087, 21.30186000000000135 47.66622999999999877, 21.30786000000000158 47.66823000000000121, 21.30285999999999902 47.67522999999999911, 21.26785999999999888 47.67822999999999922, 21.25585999999999842 47.6962299999999999, 21.23885999999999896 47.70022999999999769, 21.25186000000000064 47.70123000000000246, 21.25186000000000064 47.70523000000000025, 21.24485999999999919 47.7102300000000028, 21.23285999999999873 47.70723000000000269, 21.23285999999999873 47.71123000000000047, 21.25186000000000064 47.71522999999999826, 21.23885999999999896 47.72223000000000326, 21.24185999999999908 47.72523000000000337, 21.23585999999999885 47.73022999999999882, 21.2428600000000003 47.73022999999999882, 21.23585999999999885 47.73122999999999649, 21.23686000000000007 47.73523000000000138, 21.24986000000000175 47.73422999999999661, 21.24485999999999919 47.73322999999999894, 21.24986000000000175 47.7242299999999986, 21.26585999999999999 47.71822999999999837, 21.25885999999999854 47.72923000000000116, 21.26086000000000098 47.73722999999999672, 21.24986000000000175 47.74022999999999683, 21.25985999999999976 47.74822999999999951, 21.25385999999999953 47.75122999999999962, 21.27786000000000044 47.75422999999999973, 21.28086000000000055 47.76323000000000008, 21.26785999999999888 47.76122999999999763, 21.29286000000000101 47.77723000000000297, 21.28086000000000055 47.78123000000000076, 21.26086000000000098 47.7762299999999982, 21.28986000000000089 47.78522999999999854, 21.2948599999999999 47.78123000000000076, 21.29586000000000112 47.78623000000000332, 21.3208599999999997 47.77123000000000275, 21.32685999999999993 47.78723000000000098, 21.31685999999999837 47.7902300000000011, 21.31185999999999936 47.80022999999999911, 21.32385999999999981 47.7902300000000011, 21.3418600000000005 47.79223000000000354, 21.34885999999999839 47.78723000000000098, 21.3468599999999995 47.7832300000000032, 21.35686000000000107 47.78222999999999843, 21.33886000000000038 47.79923000000000144, 21.32885999999999882 47.80523000000000167, 21.32186000000000092 47.80223000000000155, 21.31485999999999947 47.81022999999999712, 21.32786000000000115 47.81123000000000189, 21.33986000000000161 47.80223000000000155, 21.33886000000000038 47.81123000000000189, 21.32685999999999993 47.82023000000000224, 21.32585999999999871 47.82723000000000013, 21.32786000000000115 47.82323000000000235, 21.33886000000000038 47.82623000000000246, 21.3418600000000005 47.8212299999999999, 21.36385999999999896 47.82222999999999757, 21.3678600000000003 47.81822999999999979, 21.36085999999999885 47.81423000000000201, 21.37386000000000053 47.81423000000000201, 21.37386000000000053 47.81022999999999712, 21.38685999999999865 47.81423000000000201, 21.39486000000000132 47.80722999999999701, 21.39585999999999899 47.81423000000000201, 21.40686000000000178 47.81723000000000212, 21.42286000000000001 47.80823000000000178, 21.42286000000000001 47.81322999999999723, 21.43586000000000169 47.81622999999999735, 21.44485999999999848 47.80722999999999701, 21.44485999999999848 47.81723000000000212, 21.44986000000000104 47.80523000000000167, 21.46386000000000038 47.81022999999999712, 21.46186000000000149 47.80322999999999922, 21.47985999999999862 47.80722999999999701, 21.48186000000000107 47.80223000000000155, 21.47685999999999851 47.80122999999999678, 21.49485999999999919 47.80122999999999678, 21.50486000000000075 47.78522999999999854, 21.48686000000000007 47.78023000000000309, 21.49386000000000152 47.77022999999999797, 21.50685999999999964 47.77223000000000042, 21.49986000000000175 47.76623000000000019, 21.49686000000000163 47.75122999999999962, 21.50486000000000075 47.75323000000000206, 21.50685999999999964 47.74622999999999706, 21.50985999999999976 47.75222999999999729, 21.50885999999999854 47.74123000000000161, 21.51285999999999987 47.74622999999999706, 21.51986000000000132 47.74222999999999928, 21.51386000000000109 47.72323000000000093, 21.52586000000000155 47.7172300000000007, 21.52586000000000155 47.71123000000000047, 21.51386000000000109 47.71222999999999814, 21.50885999999999854 47.70723000000000269, 21.50685999999999964 47.69223000000000212, 21.49686000000000163 47.69022999999999968, 21.50285999999999831 47.68822999999999723, 21.50186000000000064 47.67322999999999666, 21.48585999999999885 47.66622999999999877, 21.48385999999999996 47.65623000000000076, 21.45385999999999882 47.65422999999999831, 21.45685999999999893 47.65923000000000087, 21.44785999999999859 47.65323000000000064))

[mdavis]

JTS has the same bug. Ticket ​here.

[mdavis]

Copying the comments on the GDAL ticket for posterity:

There's a few things going on here.

• The bad output is caused by GEOS trying to rectify an invalid output polygon from simplification by using buffer(0). Obviously that trick doesn't always work... ( buffer sometimes keeps the "wrong" lobe of a self-crossing polygon)
• QGIS appears to be using the GEOS TopologyPreservingSimplifier, which preserves topology by not simplifying vertices which would cause invalid output.
• The PostGIS ST_Simplify uses an internal implementation which can produce invalid output. But it also offers ST_SimplifyPreserveTopology.
• In JTS there is a flag on DouglasPeuckerSimplifier which allows skipping the area-topology-fixing step. This will give the same output as PostGIS ST_Simplify. Unfortunately the GEOS code doesn't have this flag.

My suggestion for dealing with this is:

• The method of forcing valid polygonal output for Simplify is problematic, and no simple fix is obvious. So for now this hack should probably be disabled. Either the flag should be added, or (more simply) Simplify should be changed to change its default behaviour.
• If valid output is required then either use TopologyPreservingSimplifier or MakeValid

[dbaston]

Agreed. Given that we expose a GEOSTopologyPreserveSimplify, I think we should just let GEOSSimplify produce an invalid output.

[mdavis]

Fixed in JTS by ​https://github.com/locationtech/jts/pull/655

[mdavis]

Fixed by efe09d6/git

See also #732