search

OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.89k stars 2.54k forks source link

Geometry.Simplify provides wrong results #2044

Closed kamad3 closed 4 years ago

kamad3 commented 4 years ago

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))

rouault commented 4 years ago

GDAL mostly relies on GEOS Simplify(). Seems to work for me with both geos 3.5.0 and 3.8. What is yours ?

>>> from osgeo import ogr
>>> g = ogr.CreateGeometryFromWkt('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.3786 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))')
>>> print g.Simplify(0.00036)
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.3786 47.66823,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))
kamad3 commented 4 years ago

I used gdal 2.4.3 with conda package manager, my geos version is 3.7.2

Could you please try it with Simplify(0.0036) instead of Simplify(0.00036)?

Anyway I tried to upgrade to gdal 3.0.2, but with this version I receive the following error messages at another part of my code.

Code is: src = osr.SpatialReference() tgt = osr.SpatialReference() src.ImportFromEPSG(4326) tgt.ImportFromEPSG(23700) transform = osr.CoordinateTransformation(src, tgt)

Errors are: ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db ERROR 1: PROJ: proj_create: Cannot find proj.db ERROR 1: PROJ: proj_create: unrecognized format / unknown name

In spite of the error messages simplify(0.0036) is running (gdal 3.0.2, geos 3.8.0), but I receive the same output: POLYGON ((21.3280682018927 47.8232867823344,21.32786 47.82323,21.32586 47.82723,21.3280682018927 47.8232867823344))

rouault commented 4 years ago

With simplify(0.0036), I also receive the same output as you, which is different from the one from PostGIS. Anyway GEOS does the processing here, so nothing that GDAL can do about directly. Interestingly, I found that PostGIS ST_Simplify() does not use GEOSSimplify_r() but its own lwgeom_simplify_in_place() implementation in https://github.com/postgis/postgis/blob/570fbdb10de728e56957f3b7b449ac9532a53a47/liblwgeom/lwgeom.c#L1715 Perhaps @pramsey has some background about this (if it is known that GEOSSimplify() and ST_Simplify() have different behaviour, and if this was the reason for PostGIS to not use it. Or perhaps this was just for more performance)

jratike80 commented 4 years ago

JTS 1.14, used with OpenJUMP creates this small polygon

POLYGON ((
        21.328068201892744 47.823286782334385, 
        21.32786 47.82323, 
        21.32586 47.82723, 
        21.328068201892744 47.823286782334385
    ))

The use case is not very realistic because the original geometry has long spikes and rather few vertices which makes the geometry not so suitable for simplification with DP. Visvalingam-Whyatt method might work better. Naturally the current result is very wrong.

rouault commented 4 years ago

Thanks @jratike80 for investigating. @kamad3 : I'd suggest you to file a ticket at https://trac.osgeo.org/geos/newticket (GEOS) and/or https://github.com/locationtech/jts/issues (JTS). Nothing we can do on GDAL side. Closing as not a GDAL bug.

kamad3 commented 4 years ago

Thank you guys for the investigation. I try to implement the recommeneded VW algorithm and raise a ticket on geos.

dr-jts commented 4 years ago

There's a few things going on here.

The lack of the "disable valid topology hack" flag in GEOS is a bit of a showstopper. Probably the GEOS DouglsPeuckerSimplifier code should be changed to skip this by default (thus matching the PostGIS behaviour. If valid polygons are required then either TopologyPreservingSimplifier can be used. Or, post facto processing with GEOS MakeValid.