Repository: geopandas/geopandas

Issue:
Reprojection: loss of elevation when using pygeos 
With https://github.com/geopandas/geopandas/commit/5d1181af6bcca88351d84bb63fc0fbe0e32ce055 and pygeos support enabled, to_crs() transforms all geometries to 2-dimensions.

_vectorized.transform uses pygeos.get_coordinates, and the documentation mentions: Three-dimensional data is ignored.

Discussion:
martinfleis: I think that we should try to mimic the behaviour of shapely 1.8, meaning we would include Z for geometries that have Z and don't include it when the geometry don't have it. Like this:

```py
        has_z = shapely.has_z(data)

        result = np.empty(data.shape, dtype=data.dtype)

        coords = shapely.get_coordinates(data[~has_z], include_z=False)
        new_coords_z = func(coords[:, 0], coords[:, 1])
        result[~has_z] = shapely.set_coordinates(
            data[~has_z].copy(), np.array(new_coords_z).T
        )

        coords_z = shapely.get_coordinates(data[has_z], include_z=True)
        new_coords_z = func(coords_z[:, 0], coords_z[:, 1], coords_z[:, 2])
        result[has_z] = shapely.set_coordinates(
            data[has_z].copy(), np.array(new_coords_z).T
        )
```
derinwalters: Thank you, @jorisvandenbossche . This response was incredibly fast and exactly what I was hoping for. Question: if we want to allow overriding the include_z inference, don't geodataframe.py and geoseries.py to_crs method definitions need the include_z argument exposed and subsequently drilled down as well?
jorisvandenbossche: > I think that we should try to mimic the behaviour of shapely 1.8, meaning we would include Z for geometries that have Z and don't include it when the geometry don't have it. 

OK, I was assuming to make this a on/off choice for the _whole_ geometry column, basically not allowing mixed dimensions in the output (knowing this is still a breaking change in behaviour compared to before). 
But it's true we can do the effort to keep the mixed geometries. You example code is not that complex.

Now, I have to check more, but it might also be this already works somewhat out of the box (the following is with the current PR):

```
In [8]: arr = geopandas.array.from_shapely([shapely.Point(1, 2, 3), shapely.Point(2, 3), shapely.Point(3, 4, 5)], crs="EPSG:432
   ...: 6")

In [9]: arr
Out[9]: 
<GeometryArray>
[<POINT Z (1 2 3)>, <POINT (2 3)>, <POINT Z (3 4 5)>]
Length: 3, dtype: geometry

In [10]: arr.to_crs("EPSG:3857", include_z=True)
Out[10]: 
<GeometryArray>
[<POINT Z (111319.491 222684.209 3)>,     <POINT (222638.982 334111.171)>,
  <POINT Z (333958.472 445640.11 5)>]
Length: 3, dtype: geometry
```

I _think_ this is because `get_coordinates` fills with NaNs for the geometries without z:

```
In [12]: shapely.get_coordinates(arr.data, include_z=True)
Out[12]: 
array([[ 1.,  2.,  3.],
       [ 2.,  3., nan],
       [ 3.,  4.,  5.]])
```

And then the transform function keeps the NaN, and when recreating the original geometries those NaNs are just ignored (since the original geometry was 2D).

jorisvandenbossche: Although, my test above uses a CRS transformation that is not defined for Z values, and so leaves those values alone. If you have a transformation that includes the z values, NaN elevation values might actually give problems. 
Is that what you mean with getting empty Points? (if the NaN z values gets used and propagates to x/y?)
martinfleis: Consider this snippet:

```py
s = geopandas.GeoSeries(
    [
        Point(1, 2),
        LineString([(1, 2, 3), (4, 5, 6)]),
    ],
    crs=2056
)
data = s.values.data

transformer = Transformer.from_crs(s.crs, 4326, always_xy=True)
func = transformer.transform

coords_z = pygeos.get_coordinates(data, include_z=True)
new_coords_z = func(coords_z[:, 0], coords_z[:, 1], coords_z[:, 2])
result = pygeos.set_coordinates(data.copy(), np.array(new_coords_z).T)
```
```
array([<pygeos.Geometry POINT EMPTY>,
       <pygeos.Geometry LINESTRING Z (-19.918 32.125 3, -19.918 32.125 6)>],
      dtype=object)
```

If you include Z which is NaN, pyproj actually returns an array of NaNs resulting in an empty point. If you mask non-Z and use 2d transformation, all is fine.
jorisvandenbossche: Trying to find an example of 3D CRS, the following still seems to work (it gives the same x/y values regardless of include_z=True/False, and the z value is actually transformed for the 3D Point):

```
In [17]: arr = geopandas.array.from_shapely([shapely.Point(4, 50, 10), shapely.Point(4, 51)], crs="EPSG:4979")

In [18]: arr.to_crs("EPSG:6190", include_z=True)
Out[18]: 
<GeometryArray>
[<POINT Z (123563.583 76585.173 -33.228)>, <POINT (124114.915 187817.059)>]
Length: 2, dtype: geometry

In [19]: arr.to_crs("EPSG:6190", include_z=False)
Out[19]: 
<GeometryArray>
[<POINT (123563.583 76585.173)>, <POINT (124114.915 187817.059)>]
Length: 2, dtype: geometry
```
martinfleis: Use the one above. That fails.

```py
arr = geopandas.array.from_shapely([shapely.Point(4, 50, 10), shapely.Point(4, 51)], crs="EPSG:2056")
arr.to_crs("EPSG:4326", include_z=True)
```
```
<GeometryArray>
[<POINT Z (-19.918 32.125 10)>, <POINT EMPTY>]
Length: 2, dtype: geometry
```

jorisvandenbossche: OK, I see, it very much depends on the exact CRS you are using what PROJ exactly does with NaN z values: (I was just using the wrong example ;))

```
In [40]: transformer = Transformer.from_crs(4326, 2056, always_xy=True)

In [41]: transformer.transform(1, 2, np.nan)
Out[41]: (nan, nan, nan)

In [42]: transformer = Transformer.from_crs(4326, 3857, always_xy=True)

In [43]: transformer.transform(1, 2, np.nan)
Out[43]: (111319.49079327357, 222684.20850554405, nan)
```
martinfleis: Yeah, I was a bit lucky to catch that selecting 2056->4326. But if some transformations result in this, which I'd assume comes from PROJ, we should make an effort to mask Z geometries as I suggested in https://github.com/geopandas/geopandas/pull/2641#issuecomment-1303127373. The overhead of that doesn't seem to be significant, timings of `has_z` check over the whole array are in ms.
jorisvandenbossche: Updated to handle mixed dimensions as suggested!

Question: if we handle mixed dimensions automatically (and so preserve the dimensionality of each individual input geometry), then I don't think there is any need to expose a `ìnclude_z` keyword in the public `to_crs` methods? (cc @derinwalters) 
The only thing I could think of is that one could set `include_z=False` to also convert your XYZ geometries to XY while transforming. But then there are APIs to do that specifically (shapely has `force_2d()` / `force_3d()` functions, and so we should better expose those, xref https://github.com/geopandas/geopandas/issues/2010), which I think is the preferred way. 
derinwalters: > Question: if we handle mixed dimensions automatically (and so preserve the dimensionality of each individual input geometry), then I don't think there is any need to expose a `ìnclude_z` keyword in the public `to_crs` methods? 

Handling mixed dimensions automatically as proposed, rather than having to specify explicitly via include_z, is my preferred mechanism as well. I agree Shapely provides sufficient API for controlling other things like dropping dimensions and that it would be convenient to expose many methods in Geopandas (ex. geometry 3d->2d for polyline generation).
jorisvandenbossche: @martinfleis I think this was ready

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index bd6ede1d..f56072dd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,11 @@
 ## Development version
 
 
+Bug fixes:
+
+- Correctly handle geometries with Z dimension in ``to_crs()`` when using PyGEOS or
+  Shapely >= 2.0 (previously the z coordinates were lost) (#1345).
+
 ## Version 0.12.1 (October 29, 2022)
 
 Small bug-fix release removing the shapely<2 pin in the installation requirements.
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 45d383dd..b1700b49 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -1096,15 +1096,24 @@ def bounds(data):
 
 
 def transform(data, func):
-    if compat.USE_SHAPELY_20:
-        coords = shapely.get_coordinates(data)
-        new_coords = func(coords[:, 0], coords[:, 1])
-        result = shapely.set_coordinates(data.copy(), np.array(new_coords).T)
-        return result
-    if compat.USE_PYGEOS:
-        coords = pygeos.get_coordinates(data)
-        new_coords = func(coords[:, 0], coords[:, 1])
-        result = pygeos.set_coordinates(data.copy(), np.array(new_coords).T)
+    if compat.USE_SHAPELY_20 or compat.USE_PYGEOS:
+        if compat.USE_SHAPELY_20:
+            has_z = shapely.has_z(data)
+            from shapely import get_coordinates, set_coordinates
+        else:
+            has_z = pygeos.has_z(data)
+            from pygeos import get_coordinates, set_coordinates
+
+        result = np.empty_like(data)
+
+        coords = get_coordinates(data[~has_z], include_z=False)
+        new_coords_z = func(coords[:, 0], coords[:, 1])
+        result[~has_z] = set_coordinates(data[~has_z].copy(), np.array(new_coords_z).T)
+
+        coords_z = get_coordinates(data[has_z], include_z=True)
+        new_coords_z = func(coords_z[:, 0], coords_z[:, 1], coords_z[:, 2])
+        result[has_z] = set_coordinates(data[has_z].copy(), np.array(new_coords_z).T)
+
         return result
     else:
         from shapely.ops import transform


Test patch:
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index 1d9200a2..b8983a62 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -58,6 +58,13 @@ def test_to_crs_transform__missing_data():
     assert_geodataframe_equal(df, utm, check_less_precise=True)
 
 
+def test_to_crs_transform__empty_data():
+    df = df_epsg26918().iloc[:0]
+    lonlat = df.to_crs(epsg=4326)
+    utm = lonlat.to_crs(epsg=26918)
+    assert_geodataframe_equal(df, utm, check_less_precise=True)
+
+
 def test_to_crs_inplace():
     df = df_epsg26918()
     lonlat = df.to_crs(epsg=4326)
@@ -77,6 +84,26 @@ def test_to_crs_geo_column_name():
     assert_geodataframe_equal(df, utm, check_less_precise=True)
 
 
+def test_to_crs_dimension_z():
+    # preserve z dimension
+    arr = points_from_xy([1, 2], [2, 3], [3, 4], crs=4326)
+    assert arr.has_z.all()
+    result = arr.to_crs(epsg=3857)
+    assert result.has_z.all()
+
+
+def test_to_crs_dimension_mixed():
+    s = GeoSeries([Point(1, 2), LineString([(1, 2, 3), (4, 5, 6)])], crs=2056)
+    result = s.to_crs(epsg=4326)
+    assert not result[0].is_empty
+    assert result.has_z.tolist() == [False, True]
+    roundtrip = result.to_crs(epsg=2056)
+    # TODO replace with assert_geoseries_equal once we expose tolerance keyword
+    # assert_geoseries_equal(roundtrip, s, check_less_precise=True)
+    for a, b in zip(roundtrip, s):
+        np.testing.assert_allclose(a.coords[:], b.coords[:], atol=0.01)
+
+
 # -----------------------------------------------------------------------------
 # Test different supported formats for CRS specification