Repository: geopandas/geopandas

Issue:
TST: failures with shapely and GEOS 3.9
CI started coming with GEOS 3.9 today which caused a lot of errors if we use shapely backend. Most of them are related to overlay, so I guess it is caused by the new `OverlayNG`, which returns slightly different geometries.

https://github.com/geopandas/geopandas/runs/1854075806

Discussion:
martinfleis: Thanks a lot!

> Eventually we might also want to add it to the GeometryArray/GeoSeries, but for that I would like to wait on Shapely including it (Toblerity/Shapely#1090)

Agree.

> We can skip the geometry equality check for now with GEOS 3.9 (we already check area and bounds, and that passes), and later re-produce the expected files with GEOS 3.9

Let's do that. I have also observed some differences in how GEOS 3.9 handles specific situations (e.g. order of geometries from polygonize is different) and concluded that it is the best solution to tailor tests for that.

jorisvandenbossche: > I have also observed some differences in how GEOS 3.9 handles specific situations (e.g. order of geometries from polygonize is different)

If it's only the order of geometries, that's something that `normalize` should solve. But in this overlay case with NYBB, also the number of geometries of a MultiPolygon is different, which is not something normalize fixes.
martinfleis: It is not an order of geometries within MultiPolygon, it is an order of Polygons from `polygonize`.

> the number of geometries of a MultiPolygon is different

That is weird. Have you looked at the difference in detail?
martinfleis: It seems that our `sort` keyword in STRTree is now obsolete since the result is sorted in any case with GEOS 3.9 - https://github.com/geopandas/geopandas/pull/1826/checks?check_run_id=1901961534#step:7:2912
jorisvandenbossche: > > the number of geometries of a MultiPolygon is different
>
> That is weird. Have you looked at the difference in detail?

That's actually already the case before GEOS 3.9 as well (that's the reason we *already* were skipping many polygons when asserting the actual geometries in the nybb test)

Looking at the actual vs expected result on GEOS 3.9 in a notebook, it are all small but real differences at the border of the circles. 

New result:

```
result.plot(figsize=(10, 10))
```

![image](https://user-images.githubusercontent.com/1020496/108121485-2b197180-70a3-11eb-9c41-bd23b2e76a8c.png)

Difference with result/expected:

```
diffs = result.difference(expected)
diffs[~diffs.is_empty].plot(edgecolor='k', figsize=(10, 10))
```

![image](https://user-images.githubusercontent.com/1020496/108121580-4dab8a80-70a3-11eb-9a63-9ead4227b760.png)

Seeing this image, with mostly small change at the border (even in areas that are otherwise fully contained in a polygon of the boroughs), I am wondering: could this be caused because the *buffer* is creating a slightly different circular polygon? Because the expected results were generated with buffers created with an older GEOS version, so if with the latest GEOS we get slightly different buffers, it's normal we can't expect the same result now. 

I don't know why the output of buffer would have changed, but so that's actually the same issue as with the `test_polygon_patch` where the number of vertices of the buffer changed (now I think of it)

jorisvandenbossche: (I should test my hypothesis by creating and saving the buffers to a file with an older GEOS version, and reading those in instead of creating on the fly in the test)
martinfleis: Good call with buffers, you're most likely right that it creates them differently now. Let's skip for now and change these tests later to be a bit more forgiving.
jorisvandenbossche: Turns out we already have the buffers included as a shapefile ;) So easy enough to use those instead, and that also solves the additional differences with GEOS 3.9.

That only leaves the pre-existing differences with the QGIS output. I looked at the difference for "union", and it's a really tiny difference:

![image](https://user-images.githubusercontent.com/1020496/108173450-6eed9480-70fe-11eb-9088-5f3845527aaa.png)

![image](https://user-images.githubusercontent.com/1020496/108173502-7ad95680-70fe-11eb-8e82-0d457a14e5f1.png)

So the right figures show the result from geopandas, and the expected QGIS result doesn't have this one sliver that sticks out.
jorisvandenbossche: > It seems that our sort keyword in STRTree is now obsolete since the result is sorted in any case with GEOS 3.9 

I don't think it's actually sorted, but just a different order (our tests currently relied on testing the order we got we previous GEOS versions for `sort=False`, and that now changed apparently):

```
>           assert_array_equal(res, expected)
E           AssertionError: 
E           Arrays are not equal
E           
E           Mismatched elements: 2 / 3 (66.7%)
E           Max absolute difference: 1
E           Max relative difference: 1.
E            x: array([1, 0, 2])
E            y: array([0, 1, 2])
```

(it's the expected that was "sorted", and the result is no longer)

Updated the tests to first check the sorted elements are the same, and then let the exact equality also fail for pygeos.
jorisvandenbossche: Aha, all green now! (for some reason this morning the dev build hang for a couple of hours...)

Patch:
diff --git a/geopandas/_compat.py b/geopandas/_compat.py
index 138cf8ab..ac42f766 100644
--- a/geopandas/_compat.py
+++ b/geopandas/_compat.py
@@ -6,6 +6,8 @@ import warnings
 
 import pandas as pd
 import shapely
+import shapely.geos
+
 
 # -----------------------------------------------------------------------------
 # pandas compat
@@ -26,6 +28,8 @@ SHAPELY_GE_17 = str(shapely.__version__) >= LooseVersion("1.7.0")
 SHAPELY_GE_18 = str(shapely.__version__) >= LooseVersion("1.8")
 SHAPELY_GE_20 = str(shapely.__version__) >= LooseVersion("2.0")
 
+GEOS_GE_390 = shapely.geos.geos_version >= (3, 9, 0)
+
 
 HAS_PYGEOS = None
 USE_PYGEOS = None
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index da3bde19..9868d255 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -828,6 +828,34 @@ def simplify(data, tolerance, preserve_topology=True):
         return out
 
 
+def _shapely_normalize(geom):
+    """
+    Small helper function for now because it is not yet available in Shapely.
+    """
+    from shapely.geos import lgeos
+    from shapely.geometry.base import geom_factory
+    from ctypes import c_void_p, c_int
+
+    lgeos._lgeos.GEOSNormalize_r.restype = c_int
+    lgeos._lgeos.GEOSNormalize_r.argtypes = [c_void_p, c_void_p]
+
+    geom_cloned = lgeos.GEOSGeom_clone(geom._geom)
+    lgeos._lgeos.GEOSNormalize_r(lgeos.geos_handle, geom_cloned)
+    return geom_factory(geom_cloned)
+
+
+def normalize(data):
+    if compat.USE_PYGEOS:
+        return pygeos.normalize(data)
+    else:
+        out = np.empty(len(data), dtype=object)
+        with compat.ignore_shapely2_warnings():
+            out[:] = [
+                _shapely_normalize(geom) if geom is not None else None for geom in data
+            ]
+        return out
+
+
 def project(data, other, normalized=False):
     if compat.USE_PYGEOS:
         return pygeos.line_locate_point(data, other, normalize=normalized)


Test patch:
diff --git a/geopandas/testing.py b/geopandas/testing.py
index cc27a147..437af0d3 100644
--- a/geopandas/testing.py
+++ b/geopandas/testing.py
@@ -7,6 +7,7 @@ import pandas as pd
 
 from geopandas import GeoDataFrame, GeoSeries
 from geopandas.array import GeometryDtype
+from geopandas import _vectorized
 
 
 def _isna(this):
@@ -67,6 +68,7 @@ def assert_geoseries_equal(
     check_less_precise=False,
     check_geom_type=False,
     check_crs=True,
+    normalize=False,
 ):
     """
     Test util for checking that two GeoSeries are equal.
@@ -89,6 +91,10 @@ def assert_geoseries_equal(
     check_crs: bool, default True
         If `check_series_type` is True, then also check that the
         crs matches.
+    normalize: bool, default False
+        If True, normalize the geometries before comparing equality.
+        Typically useful with ``check_less_precise=True``, which uses
+        ``geom_almost_equals`` and requires exact coordinate order.
     """
     assert len(left) == len(right), "%d != %d" % (len(left), len(right))
 
@@ -120,6 +126,10 @@ def assert_geoseries_equal(
             right.type,
         )
 
+    if normalize:
+        left = GeoSeries(_vectorized.normalize(left.array.data))
+        right = GeoSeries(_vectorized.normalize(right.array.data))
+
     if not check_crs:
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", "CRS mismatch", UserWarning)
@@ -145,6 +155,7 @@ def assert_geodataframe_equal(
     check_less_precise=False,
     check_geom_type=False,
     check_crs=True,
+    normalize=False,
 ):
     """
     Check that two GeoDataFrames are equal/
@@ -168,6 +179,10 @@ def assert_geodataframe_equal(
     check_crs: bool, default True
         If `check_frame_type` is True, then also check that the
         crs matches.
+    normalize: bool, default False
+        If True, normalize the geometries before comparing equality.
+        Typically useful with ``check_less_precise=True``, which uses
+        ``geom_almost_equals`` and requires exact coordinate order.
     """
     try:
         # added from pandas 0.20
@@ -217,6 +232,7 @@ def assert_geodataframe_equal(
             assert_geoseries_equal(
                 left[col],
                 right[col],
+                normalize=normalize,
                 check_dtype=check_dtype,
                 check_less_precise=check_less_precise,
                 check_geom_type=check_geom_type,
diff --git a/geopandas/tests/test_overlay.py b/geopandas/tests/test_overlay.py
index 149e0957..6ba0c398 100644
--- a/geopandas/tests/test_overlay.py
+++ b/geopandas/tests/test_overlay.py
@@ -107,19 +107,23 @@ def test_overlay(dfs_index, how):
 def test_overlay_nybb(how):
     polydf = read_file(geopandas.datasets.get_path("nybb"))
 
-    # construct circles dataframe
-    N = 10
-    b = [int(x) for x in polydf.total_bounds]
-    polydf2 = GeoDataFrame(
-        [
-            {"geometry": Point(x, y).buffer(10000), "value1": x + y, "value2": x - y}
-            for x, y in zip(
-                range(b[0], b[2], int((b[2] - b[0]) / N)),
-                range(b[1], b[3], int((b[3] - b[1]) / N)),
-            )
-        ],
-        crs=polydf.crs,
-    )
+    # The circles have been constructed and saved at the time the expected
+    # results were created (exact output of buffer algorithm can slightly
+    # change over time -> use saved ones)
+    # # construct circles dataframe
+    # N = 10
+    # b = [int(x) for x in polydf.total_bounds]
+    # polydf2 = GeoDataFrame(
+    #     [
+    #         {"geometry": Point(x, y).buffer(10000), "value1": x + y, "value2": x - y}
+    #         for x, y in zip(
+    #             range(b[0], b[2], int((b[2] - b[0]) / N)),
+    #             range(b[1], b[3], int((b[3] - b[1]) / N)),
+    #         )
+    #     ],
+    #     crs=polydf.crs,
+    # )
+    polydf2 = read_file(os.path.join(DATA, "nybb_qgis", "polydf2.shp"))
 
     result = overlay(polydf, polydf2, how=how)
 
@@ -186,12 +190,23 @@ def test_overlay_nybb(how):
         result.geometry.bounds, expected.geometry.bounds, check_less_precise=True
     )
 
-    # now drop multipolygons
-    result.geometry[result.geometry.geom_type == "MultiPolygon"] = None
-    expected.geometry[expected.geometry.geom_type == "MultiPolygon"] = None
+    # There are two cases where the multipolygon have a different number
+    # of sub-geometries -> not solved by normalize (and thus drop for now)
+    if how == "symmetric_difference":
+        expected.loc[9, "geometry"] = None
+        result.loc[9, "geometry"] = None
+
+    if how == "union":
+        expected.loc[24, "geometry"] = None
+        result.loc[24, "geometry"] = None
 
     assert_geodataframe_equal(
-        result, expected, check_crs=False, check_column_type=False
+        result,
+        expected,
+        normalize=True,
+        check_crs=False,
+        check_column_type=False,
+        check_less_precise=True,
     )
 
 
@@ -244,7 +259,11 @@ def test_overlay_overlap(how):
         result = result.sort_values(["col1", "col2"]).reset_index(drop=True)
 
     assert_geodataframe_equal(
-        result, expected, check_column_type=False, check_less_precise=True
+        result,
+        expected,
+        normalize=True,
+        check_column_type=False,
+        check_less_precise=True,
     )
 
 
@@ -469,6 +488,7 @@ def test_overlay_strict(how, keep_geom_type, geom_types):
         assert_geodataframe_equal(
             result,
             expected,
+            normalize=True,
             check_column_type=False,
             check_less_precise=True,
             check_crs=False,
diff --git a/geopandas/tests/test_plotting.py b/geopandas/tests/test_plotting.py
index e64474fa..cb518550 100644
--- a/geopandas/tests/test_plotting.py
+++ b/geopandas/tests/test_plotting.py
@@ -1490,7 +1490,10 @@ def test_polygon_patch():
     patch = _PolygonPatch(polygon)
     assert isinstance(patch, PathPatch)
     path = patch.get_path()
-    assert len(path.vertices) == len(path.codes) == 198
+    if compat.GEOS_GE_390:
+        assert len(path.vertices) == len(path.codes) == 195
+    else:
+        assert len(path.vertices) == len(path.codes) == 198
 
 
 def _check_colors(N, actual_colors, expected_colors, alpha=None):
diff --git a/geopandas/tests/test_sindex.py b/geopandas/tests/test_sindex.py
index a7391967..9ef56961 100644
--- a/geopandas/tests/test_sindex.py
+++ b/geopandas/tests/test_sindex.py
@@ -412,10 +412,14 @@ class TestPygeosInterface:
 
         test_geo = test_df.geometry.values.data[0]
         res = tree_df.sindex.query(test_geo, sort=sort)
+
+        # asserting the same elements
+        assert sorted(res) == sorted(expected)
+        # asserting the exact array can fail if sort=False
         try:
             assert_array_equal(res, expected)
         except AssertionError as e:
-            if not compat.USE_PYGEOS and sort is False:
+            if sort is False:
                 pytest.xfail(
                     "rtree results are known to be unordered, see "
                     "https://github.com/geopandas/geopandas/issues/1337\n"
@@ -649,10 +653,15 @@ class TestPygeosInterface:
         test_df = geopandas.GeoDataFrame(geometry=test_polys)
 
         res = tree_df.sindex.query_bulk(test_df.geometry, sort=sort)
+
+        # asserting the same elements
+        assert sorted(res[0]) == sorted(expected[0])
+        assert sorted(res[1]) == sorted(expected[1])
+        # asserting the exact array can fail if sort=False
         try:
             assert_array_equal(res, expected)
         except AssertionError as e:
-            if not compat.USE_PYGEOS and sort is False:
+            if sort is False:
                 pytest.xfail(
                     "rtree results are known to be unordered, see "
                     "https://github.com/geopandas/geopandas/issues/1337\n"