From 507a09d7a77086777277221a22042a41bd1294fa Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Sat, 16 Mar 2024 00:42:18 +0100
Subject: [PATCH] OGRGeometry::exportToGEOS(): retain Z values for
 POLYHEDRALSURFACE Z / TIN Z / GEOMETRYCOLLECTION Z of surfaces (master only)

---
 autotest/ogr/ogr_geos.py               | 22 ++++++++++++++++++++++
 autotest/utilities/test_ogr2ogr_lib.py |  4 ++--
 ogr/ogrgeometry.cpp                    | 20 +++++++++++++++-----
 3 files changed, 39 insertions(+), 7 deletions(-)

diff --git a/autotest/ogr/ogr_geos.py b/autotest/ogr/ogr_geos.py
index 2d010efdc94f..343d28bf7967 100755
--- a/autotest/ogr/ogr_geos.py
+++ b/autotest/ogr/ogr_geos.py
@@ -634,3 +634,25 @@ def test_ogr_geos_set_precision():
     g = ogr.CreateGeometryFromWkt("LINESTRING (1 1,9 9)")
     g = g.SetPrecision(10)
     assert g.ExportToWkt() == "LINESTRING (0 0,10 10)"
+
+
+###############################################################################
+
+
+def test_ogr_geos_set_unary_union_TINZ():
+
+    g = ogr.CreateGeometryFromWkt("TIN Z (((0 0 10,0 1 10,1 1 10,0 0 10)))")
+    g = g.UnaryUnion()
+    assert g.ExportToIsoWkt() == "POLYGON Z ((0 0 10,0 1 10,1 1 10,0 0 10))"
+
+
+###############################################################################
+
+
+def test_ogr_geos_set_unary_union_GEOMETRYCOLLECTIONZ_POLYGONZ():
+
+    g = ogr.CreateGeometryFromWkt(
+        "GEOMETRYCOLLECTION Z (POLYGON Z ((0 0 10,0 1 10,1 1 10,0 0 10)))"
+    )
+    g = g.UnaryUnion()
+    assert g.ExportToIsoWkt() == "POLYGON Z ((0 0 10,0 1 10,1 1 10,0 0 10))"
diff --git a/autotest/utilities/test_ogr2ogr_lib.py b/autotest/utilities/test_ogr2ogr_lib.py
index 1f1a6999e1e4..27ae450d0fc0 100755
--- a/autotest/utilities/test_ogr2ogr_lib.py
+++ b/autotest/utilities/test_ogr2ogr_lib.py
@@ -1248,10 +1248,10 @@ def test_ogr2ogr_lib_clipsrc_3d_polygon(tmp_vsimem):
     assert lyr.GetFeatureCount() == 2
 
     feat = lyr.GetNextFeature()
-    ogrtest.check_feature_geometry(feat, "LINESTRING (0 0, 5 5)")
+    ogrtest.check_feature_geometry(feat, "LINESTRING Z (0 0 0, 5 5 5)")
 
     feat = lyr.GetNextFeature()
-    ogrtest.check_feature_geometry(feat, "LINESTRING (5 5, 10 0)")
+    ogrtest.check_feature_geometry(feat, "LINESTRING Z (5 5 5, 10 0 10)")
 
     ds = None
 
diff --git a/ogr/ogrgeometry.cpp b/ogr/ogrgeometry.cpp
index 476772b8a793..f08f81312f74 100644
--- a/ogr/ogrgeometry.cpp
+++ b/ogr/ogrgeometry.cpp
@@ -3318,7 +3318,10 @@ OGRGeometry::exportToGEOS(UNUSED_IF_NO_GEOS GEOSContextHandle_t hGEOSCtxt) const
     else if (eType == wkbPolyhedralSurface || eType == wkbTIN)
     {
         OGRGeometry *poGC = OGRGeometryFactory::forceTo(
-            poLinearGeom->clone(), wkbGeometryCollection, nullptr);
+            poLinearGeom->clone(),
+            OGR_GT_SetModifier(wkbGeometryCollection, poLinearGeom->Is3D(),
+                               poLinearGeom->IsMeasured()),
+            nullptr);
         hGeom = convertToGEOSGeom(hGEOSCtxt, poGC);
         delete poGC;
     }
@@ -3326,10 +3329,11 @@ OGRGeometry::exportToGEOS(UNUSED_IF_NO_GEOS GEOSContextHandle_t hGEOSCtxt) const
     {
         bool bCanConvertToMultiPoly = true;
         // bool bMustConvertToMultiPoly = true;
-        OGRGeometryCollection *poGC = poLinearGeom->toGeometryCollection();
+        const OGRGeometryCollection *poGC =
+            poLinearGeom->toGeometryCollection();
         for (int iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++)
         {
-            OGRwkbGeometryType eSubGeomType =
+            const OGRwkbGeometryType eSubGeomType =
                 wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType());
             if (eSubGeomType == wkbPolyhedralSurface || eSubGeomType == wkbTIN)
             {
@@ -3345,9 +3349,15 @@ OGRGeometry::exportToGEOS(UNUSED_IF_NO_GEOS GEOSContextHandle_t hGEOSCtxt) const
         if (bCanConvertToMultiPoly /* && bMustConvertToMultiPoly */)
         {
             OGRGeometry *poMultiPolygon = OGRGeometryFactory::forceTo(
-                poLinearGeom->clone(), wkbMultiPolygon, nullptr);
+                poLinearGeom->clone(),
+                OGR_GT_SetModifier(wkbMultiPolygon, poLinearGeom->Is3D(),
+                                   poLinearGeom->IsMeasured()),
+                nullptr);
             OGRGeometry *poGCDest = OGRGeometryFactory::forceTo(
-                poMultiPolygon, wkbGeometryCollection, nullptr);
+                poMultiPolygon,
+                OGR_GT_SetModifier(wkbGeometryCollection, poLinearGeom->Is3D(),
+                                   poLinearGeom->IsMeasured()),
+                nullptr);
             hGeom = convertToGEOSGeom(hGEOSCtxt, poGCDest);
             delete poGCDest;
         }