BUG: fix to_parquet/to_feather to write ISO WKB flavor for 3D geometries (if shapely >= 2 is available) (#2760)

CHANGELOG.md

@@ -24,6 +24,7 @@ Bug fixes:
 - Escape special characters to avoid TemplateSyntaxError in ``explore()`` (#2657)
 - Fix `to_parquet`/`to_feather` to not write an invalid bbox (with NaNs) in the
   metadata in case of an empty GeoDataFrame (#2653)
+- Fix `to_parquet`/`to_feather` to use correct WKB flavor for 3D geometries (#2654)
 
 Notes on (optional) dependencies:
 

geopandas/geodataframe.py

@@ -951,6 +951,7 @@ individually so that features may have different properties
             The default is to return a binary bytes object.
         kwargs
             Additional keyword args will be passed to
+            :func:`shapely.to_wkb` if shapely >= 2 is installed or
             :func:`pygeos.to_wkb` if pygeos is installed.
 
         Returns

geopandas/geoseries.py

@@ -1209,6 +1209,7 @@ e": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [3
             The default is to return a binary bytes object.
         kwargs
             Additional keyword args will be passed to
+            :func:`shapely.to_wkb` if shapely >= 2 is installed or
             :func:`pygeos.to_wkb` if pygeos is installed.
 
         Returns

geopandas/io/arrow.py

@@ -5,6 +5,7 @@ import warnings
 import numpy as np
 from pandas import DataFrame, Series
 
+import geopandas._compat as compat
 from geopandas._compat import import_optional_dependency
 from geopandas.array import from_wkb
 from geopandas import GeoDataFrame
@@ -246,7 +247,23 @@ def _geopandas_to_arrow(df, index=None, schema_version=None):
     # create geo metadata before altering incoming data frame
     geo_metadata = _create_metadata(df, schema_version=schema_version)
 
-    df = df.to_wkb()
+    kwargs = {}
+    if compat.USE_SHAPELY_20:
+        kwargs = dict(flavor="iso")
+    else:
+        for col in df.columns[df.dtypes == "geometry"]:
+            series = df[col]
+            if series.has_z.any():
+                warnings.warn(
+                    "The GeoDataFrame contains 3D geometries, and when using "
+                    "shapely < 2.0, such geometries will be written not exactly "
+                    "following to the GeoParquet spec (not using ISO WKB). For "
+                    "most use cases this should not be a problem (GeoPandas can "
+                    "read such files fine).",
+                    stacklevel=2,
+                )
+                break
+    df = df.to_wkb(**kwargs)
 
     table = Table.from_pandas(df, preserve_index=index)
 

geopandas/io/tests/test_arrow.py

@@ -15,6 +15,7 @@ from shapely.geometry import box, Point, MultiPolygon
 
 
 import geopandas
+import geopandas._compat as compat
 from geopandas import GeoDataFrame, read_file, read_parquet, read_feather
 from geopandas.array import to_wkb
 from geopandas.datasets import get_path
@@ -692,6 +693,28 @@ def test_write_read_default_crs(tmpdir, format):
     assert df.crs.equals(pyproj.CRS("OGC:CRS84"))
 
 
+def test_write_iso_wkb(tmpdir):
+    gdf = geopandas.GeoDataFrame(
+        geometry=geopandas.GeoSeries.from_wkt(["POINT Z (1 2 3)"])
+    )
+    if compat.USE_SHAPELY_20:
+        gdf.to_parquet(tmpdir / "test.parquet")
+    else:
+        with pytest.warns(UserWarning, match="The GeoDataFrame contains 3D geometries"):
+            gdf.to_parquet(tmpdir / "test.parquet")
+
+    from pyarrow.parquet import read_table
+
+    table = read_table(tmpdir / "test.parquet")
+    wkb = table["geometry"][0].as_py().hex()
+
+    if compat.USE_SHAPELY_20:
+        # correct ISO flavor
+        assert wkb == "01e9030000000000000000f03f00000000000000400000000000000840"
+    else:
+        assert wkb == "0101000080000000000000f03f00000000000000400000000000000840"
+
+
 @pytest.mark.parametrize(
     "format,schema_version",
     product(["feather", "parquet"], [None] + SUPPORTED_VERSIONS),
@@ -699,7 +722,6 @@ def test_write_read_default_crs(tmpdir, format):
 def test_write_spec_version(tmpdir, format, schema_version):
     if format == "feather":
         from pyarrow.feather import read_table
-
     else:
         from pyarrow.parquet import read_table