OGR: Add OGRLayer::GetExtent3D() (#8806)

* OGR: Add OGRLayer::GetExtent3D()
* Specialized implementation in Shapefile driver
* Add OLCFastGetExtent3D

Refs #8570

autotest/ogr/ogr_csv.py

@@ -27,6 +27,7 @@
 # Boston, MA 02111-1307, USA.
 ###############################################################################
 
+import math
 import pathlib
 import sys
 
@@ -2927,6 +2928,66 @@ def test_ogr_csv_separator_open_option(tmp_vsimem, sep, sep_opt_value, other_sep
     assert f["baz"] == "3"
 
 
+def test_ogr_csv_getextent3d(tmp_vsimem):
+
+    gdal.FileFromMemBuffer(
+        tmp_vsimem / "test.csv",
+        "id,WKT\n1,POINT Z(1 1 1)\n1,POINT Z(2 2 2)",
+    )
+
+    gdal.ErrorReset()
+    ds = gdal.OpenEx(
+        tmp_vsimem / "test.csv", gdal.OF_VECTOR, open_options=["SEPARATOR=COMMA"]
+    )
+    assert gdal.GetLastErrorMsg() == ""
+    lyr = ds.GetLayer(0)
+    dfn = lyr.GetLayerDefn()
+    assert dfn.GetGeomFieldCount() == 1
+    ext2d = lyr.GetExtent()
+    assert ext2d == (1.0, 2.0, 1.0, 2.0)
+    ext3d = lyr.GetExtent3D()
+    assert ext3d == (1.0, 2.0, 1.0, 2.0, 1.0, 2.0)
+
+    # Test 2D
+    gdal.FileFromMemBuffer(
+        tmp_vsimem / "test.csv",
+        "id,WKT\n1,POINT(1 1)\n2,POINT(2 2)",
+    )
+    gdal.ErrorReset()
+    ds = gdal.OpenEx(
+        tmp_vsimem / "test.csv", gdal.OF_VECTOR, open_options=["SEPARATOR=COMMA"]
+    )
+    assert gdal.GetLastErrorMsg() == ""
+    lyr = ds.GetLayer(0)
+    dfn = lyr.GetLayerDefn()
+    assert dfn.GetGeomFieldCount() == 1
+    ext2d = lyr.GetExtent()
+    assert ext2d == (1.0, 2.0, 1.0, 2.0)
+    ext3d = lyr.GetExtent3D()
+    assert ext3d[:4] == (1.0, 2.0, 1.0, 2.0)
+    assert math.isnan(ext3d[4])
+    assert math.isnan(ext3d[5])
+
+    # Test mixed 2D
+    gdal.FileFromMemBuffer(
+        tmp_vsimem / "test.csv",
+        "id,WKT\n1,POINT Z(1 1 1)\n2,POINT(2 2)",
+    )
+    gdal.ErrorReset()
+    ds = gdal.OpenEx(
+        tmp_vsimem / "test.csv", gdal.OF_VECTOR, open_options=["SEPARATOR=COMMA"]
+    )
+    assert gdal.GetLastErrorMsg() == ""
+    lyr = ds.GetLayer(0)
+    assert not lyr.TestCapability(ogr.OLCFastGetExtent3D)
+    dfn = lyr.GetLayerDefn()
+    assert dfn.GetGeomFieldCount() == 1
+    ext2d = lyr.GetExtent()
+    assert ext2d == (1.0, 2.0, 1.0, 2.0)
+    ext3d = lyr.GetExtent3D()
+    assert ext3d == (1.0, 2.0, 1.0, 2.0, 1.0, 1.0)
+
+
 ###############################################################################
 
 

autotest/ogr/ogr_geojson.py

@@ -4743,3 +4743,93 @@ def test_ogr_geojson_foreign_members_feature(tmp_vsimem):
         """{\n"type": "FeatureCollection",\n"name": "test",\n"features": [\n{ "type": "Feature", "properties": { }, "geometry": null, "foo":"bar"}\n]\n}"""
         in data
     )
+
+
+###############################################################################
+
+
+def test_ogr_json_getextent3d(tmp_vsimem):
+
+    jdata = r"""{
+            "type": "FeatureCollection",
+            "features": [
+                {
+                    "type": "Feature",
+                    "properties": {},
+                    "geometry": {
+                        "type": "%s",
+                        "coordinates": %s
+                    }
+                },
+                {
+                    "type": "Feature",
+                    "properties": {},
+                    "geometry": {
+                        "type": "%s",
+                        "coordinates": %s
+                    }
+                }
+            ]
+        }"""
+
+    gdal.FileFromMemBuffer(
+        tmp_vsimem / "test.json",
+        jdata % ("Point", "[1, 1, 1]", "Point", "[2, 2, 2]"),
+    )
+
+    gdal.ErrorReset()
+    ds = gdal.OpenEx(
+        tmp_vsimem / "test.json",
+        gdal.OF_VECTOR,
+    )
+    assert gdal.GetLastErrorMsg() == ""
+    lyr = ds.GetLayer(0)
+    dfn = lyr.GetLayerDefn()
+    assert dfn.GetGeomFieldCount() == 1
+    ext2d = lyr.GetExtent()
+    assert ext2d == (1.0, 2.0, 1.0, 2.0)
+    ext3d = lyr.GetExtent3D()
+    assert ext3d == (1.0, 2.0, 1.0, 2.0, 1.0, 2.0)
+
+    # Test 2D
+    gdal.FileFromMemBuffer(
+        tmp_vsimem / "test.json",
+        jdata % ("Point", "[1, 1]", "Point", "[2, 2]"),
+    )
+
+    ds = gdal.OpenEx(
+        tmp_vsimem / "test.json",
+        gdal.OF_VECTOR,
+    )
+
+    assert gdal.GetLastErrorMsg() == ""
+    lyr = ds.GetLayer(0)
+    assert not lyr.TestCapability(ogr.OLCFastGetExtent3D)
+    dfn = lyr.GetLayerDefn()
+    assert dfn.GetGeomFieldCount() == 1
+    ext2d = lyr.GetExtent()
+    assert ext2d == (1.0, 2.0, 1.0, 2.0)
+    ext3d = lyr.GetExtent3D()
+    assert ext3d[:4] == (1.0, 2.0, 1.0, 2.0)
+    assert math.isnan(ext3d[4])
+    assert math.isnan(ext3d[5])
+
+    # Test mixed 2D
+    gdal.FileFromMemBuffer(
+        tmp_vsimem / "test.json",
+        jdata % ("Point", "[1, 1, 1]", "Point", "[2, 2]"),
+    )
+
+    ds = gdal.OpenEx(
+        tmp_vsimem / "test.json",
+        gdal.OF_VECTOR,
+    )
+
+    assert gdal.GetLastErrorMsg() == ""
+    lyr = ds.GetLayer(0)
+    dfn = lyr.GetLayerDefn()
+    assert dfn.GetGeomFieldCount() == 1
+    ext2d = lyr.GetExtent()
+    assert ext2d == (1.0, 2.0, 1.0, 2.0)
+    ext3d = lyr.GetExtent3D()
+    assert ext3d == (1.0, 2.0, 1.0, 2.0, 1.0, 1.0)

autotest/ogr/ogr_shape.py

@@ -27,6 +27,7 @@
 # Boston, MA 02111-1307, USA.
 ###############################################################################
 
+import math
 import os
 import shutil
 import struct
@@ -1853,6 +1854,11 @@ def test_ogr_shape_48(tmp_vsimem):
     lyr.SetFeature(feat)
     extent = lyr.GetExtent()
     assert extent == (1, 3, 2, 4), "did not get expected extent (1)"
+    extent3D = lyr.GetExtent3D()
+    assert lyr.TestCapability(ogr.OLCFastGetExtent3D)
+    assert extent3D[:4] == (1, 3, 2, 4), "did not get expected extent 3D"
+    assert math.isnan(extent3D[4])
+    assert math.isnan(extent3D[5])
 
     ds.ExecuteSQL("RECOMPUTE EXTENT ON ogr_shape_48")
     extent = lyr.GetExtent()
@@ -1898,13 +1904,15 @@ def test_ogr_shape_48(tmp_vsimem):
     )
     lyr.CreateFeature(feat)
     feat.SetGeometry(
-        ogr.CreateGeometryFromWkt("POLYGON((0 0 2,0 1 2,1 1 2,1 0 2,0 0 2))")
+        ogr.CreateGeometryFromWkt("POLYGON((0 0 2,0 1 1,1 1 2,1 0 2,0 0 3))")
     )
     lyr.SetFeature(feat)
     ds.ExecuteSQL("RECOMPUTE EXTENT ON ogr_shape_48")
-    # FIXME: when we have a GetExtent3D
     extent = lyr.GetExtent()
     assert extent == (0, 1, 0, 1), "did not get expected extent (4)"
+    extent3D = lyr.GetExtent3D()
+    assert lyr.TestCapability(ogr.OLCFastGetExtent3D)
+    assert extent3D == (0, 1, 0, 1, 1, 3), "did not get expected extent 3D"
     ds = None
     ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource(
         tmp_vsimem / "ogr_shape_48.shp"

ogr/ogr_api.h

@@ -765,6 +765,17 @@ GIntBig CPL_DLL OGR_L_GetFeatureCount(OGRLayerH, int);
 OGRErr CPL_DLL OGR_L_GetExtent(OGRLayerH, OGREnvelope *, int);
 OGRErr CPL_DLL OGR_L_GetExtentEx(OGRLayerH, int iGeomField,
                                  OGREnvelope *psExtent, int bForce);
+/**
+ * @brief OGR_L_GetExtent3D computes the 3D extent of the layer.
+ * @param hLayer the layer to consider.
+ * @param iGeomField 0-based index of the geometry field to consider.
+ * @param psExtent3D the computed 3D extent of the layer.
+ * @param bForce if TRUE, the extent will be computed even if all the
+ *        layer features have to be fetched.
+ * @return OGRERR_NONE on success or an error code in case of failure.
+ */
+OGRErr CPL_DLL OGR_L_GetExtent3D(OGRLayerH hLayer, int iGeomField,
+                                 OGREnvelope3D *psExtent3D, int bForce);
 int CPL_DLL OGR_L_TestCapability(OGRLayerH, const char *);
 OGRErr CPL_DLL OGR_L_CreateField(OGRLayerH, OGRFieldDefnH, int);
 OGRErr CPL_DLL OGR_L_CreateGeomField(OGRLayerH hLayer,

ogr/ogr_core.h

@@ -48,9 +48,12 @@
 extern "C++"
 {
 #if !defined(DOXYGEN_SKIP)
+#include <cmath>
 #include <limits>
 #endif
 
+    class OGREnvelope3D;
+
     /**
      * Simple container for a bounding region (rectangle)
      */
@@ -228,6 +231,12 @@ extern "C++"
         /** Assignment operator */
         OGREnvelope3D &operator=(const OGREnvelope3D &) = default;
 
+        /** Returns TRUE if MinZ and MaxZ are both valid numbers. */
+        bool Is3D() const
+        {
+            return !std::isnan(MinZ) && !std::isnan(MaxZ);
+        }
+
         /** Minimum Z value */
         double MinZ;
 
@@ -256,8 +265,19 @@ extern "C++"
             MaxX = MAX(MaxX, sOther.MaxX);
             MinY = MIN(MinY, sOther.MinY);
             MaxY = MAX(MaxY, sOther.MaxY);
-            MinZ = MIN(MinZ, sOther.MinZ);
-            MaxZ = MAX(MaxZ, sOther.MaxZ);
+            if (!std::isnan(sOther.MinZ) && !std::isnan(sOther.MaxZ))
+            {
+                if (std::isnan(MinZ) || std::isnan(MaxZ))
+                {
+                    MinZ = sOther.MinZ;
+                    MaxZ = sOther.MaxZ;
+                }
+                else
+                {
+                    MinZ = MIN(MinZ, sOther.MinZ);
+                    MaxZ = MAX(MaxZ, sOther.MaxZ);
+                }
+            }
         }
 
         /** Update the current object by computing its union with the other
@@ -993,6 +1013,8 @@ int CPL_DLL OGRParseDate(const char *pszInput, OGRField *psOutput,
                         */
 #define OLCFastGetExtent                                                       \
     "FastGetExtent" /**< Layer capability for fast extent retrieval */
+#define OLCFastGetExtent3D                                                     \
+    "FastGetExtent3D" /**< Layer capability for fast 3D extent retrieval */
 #define OLCCreateField                                                         \
     "CreateField" /**< Layer capability for field creation                     \
                    */

ogr/ograpispy.cpp

@@ -790,6 +790,15 @@ void OGRAPISpy_L_GetExtentEx(OGRLayerH hLayer, int iGeomField, int bForce)
     OGRAPISpyFileClose();
 }
 
+void OGRAPISpy_L_GetExtent3D(OGRLayerH hLayer, int iGeomField, int bForce)
+{
+    CPLMutexHolderD(&hMutex);
+    OGRAPISpyFlushDefered();
+    fprintf(fpSpyFile, "%s.GetExtent3D(geom_field=%d, force=%d)\n",
+            OGRAPISpyGetLayerVar(hLayer).c_str(), iGeomField, bForce);
+    OGRAPISpyFileClose();
+}
+
 void OGRAPISpy_L_SetAttributeFilter(OGRLayerH hLayer, const char *pszFilter)
 {
     CPLMutexHolderD(&hMutex);

ogr/ograpispy.h

@@ -104,6 +104,7 @@ void OGRAPISpy_Dataset_RollbackTransaction(GDALDatasetH hDS);
 void OGRAPISpy_L_GetFeatureCount(OGRLayerH hLayer, int bForce);
 void OGRAPISpy_L_GetExtent(OGRLayerH hLayer, int bForce);
 void OGRAPISpy_L_GetExtentEx(OGRLayerH hLayer, int iGeomField, int bForce);
+void OGRAPISpy_L_GetExtent3D(OGRLayerH hLayer, int iGeomField, int bForce);
 void OGRAPISpy_L_SetAttributeFilter(OGRLayerH hLayer, const char *pszFilter);
 void OGRAPISpy_L_GetFeature(OGRLayerH hLayer, GIntBig nFeatureId);
 void OGRAPISpy_L_SetNextByIndex(OGRLayerH hLayer, GIntBig nIndex);