vector::ops::transformations::tests::test_simplify_preserve_topology macOS test regression

[metasim]

The test vector::ops::transformations::tests::test_simplify_preserve_topology is failing on macOS under GDAL 3.7.1.

Output

---- vector::ops::transformations::tests::test_simplify_preserve_topology stdout ----
thread 'vector::ops::transformations::tests::test_simplify_preserve_topology' panicked at 'assertion failed: `(left == right)`
  left: `POLYGON ((45 20,10 10,30 5,45 20),(30 20,20 15,20 25,30 20))`,
 right: `POLYGON ((20 35,10 10,30 5,45 20,20 35),(30 20,20 15,20 25,30 20))`', src/vector/ops/transformations.rs:258:9

Build Details

GDALVersionInfo {
    RELEASE_NAME: "3.7.1"
    RELEASE_DATE: "20230706"
    VERSION_NUM: "3070100"
    BUILD_INFO {
        CURL_ENABLED: "YES"
        PAM_ENABLED: "YES"
        GEOS_ENABLED: "YES"
        OGR_ENABLED: "YES"
        GEOS_VERSION: "3.12.0-CAPI-1.18.0"
        CURL_VERSION: "7.79.1"
        PROJ_RUNTIME_VERSION: "9.2.1"
        PROJ_BUILD_VERSION: "9.2.1"
        COMPILER: "clang 14.0.0 (clang-1400.0.29.202)"
    }
}

$ uname -a
Darwin  21.6.0 Darwin Kernel Version 21.6.0: Thu Jun  8 23:56:13 PDT 2023; root:xnu-8020.240.18.701.6~1/RELEASE_ARM64_T6000 arm64

[lnicola]

I can only conclude that your MacOS GDAL is cursed 🥲.

[metasim]

I'll try rolling back to an earlier version.

[lnicola]

So if I'm not mistaken, this is the original polygon:

This is the simplified version:

And this is what your GDAL outputs:

I've.. never considered what happens when you have an interior ring located like that.

[lnicola]

CC @rouault does this look like an expected output for OGR_G_SimplifyPreserveTopology? GEOS 3.10 "works", but 3.12 appears (unconfirmed, I don't have it on hand) to produce the output at the bottom. I did notice some changes related to this in the GEOS 3.12 release notes.

Possible repro:

>>> from osgeo import gdal, ogr
>>> print(gdal.VersionInfo("BUILD_INFO")
... )
PAM_ENABLED=YES
OGR_ENABLED=YES
CURL_ENABLED=YES
CURL_VERSION=7.81.0
GEOS_ENABLED=YES
GEOS_VERSION=3.10.2-CAPI-1.16.0
PROJ_BUILD_VERSION=9.2.1
PROJ_RUNTIME_VERSION=9.2.1
COMPILER=GCC 11.3.0

>>> poly = ogr.CreateGeometryFromWkt("POLYGON ((20 35,10 30,10 10,30 5,45 20,20 35),(30 20,20 15,20 25,30 20))")

# "good" output
>>> poly.SimplifyPreserveTopology(100).ExportToWkt()
'POLYGON ((20 35,10 10,30 5,45 20,20 35),(30 20,20 15,20 25,30 20))'

# "bad" output is POLYGON ((45 20,10 10,30 5,45 20),(30 20,20 15,20 25,30 20))

[lnicola]

According to @dr-jts, this is a know but in GEOS. We should drop the hole or use another polygon (from JTS?). But keeping the outer ring is probably enough.

We might still get a different answer though.

[metasim]

@lnicola Thanks for the fantastic research on this.

This is the only test I can find around this functionality in OSGeo/gdal. Given that validating GEOS or GDAL is a non-goal here, should we just replicate this simpler test, and consider that enough to verify the binding works?

[lnicola]

Yup, that sounds good.