Merge pull request #7989 from rouault/fix_7986

 OpenFileGDB: correctly read POINT EMPTY (fixes #7986)

autotest/ogr/data/filegdb/empty_polygon.gdb/timestamps

@@ -0,0 +1 @@
+����������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������

autotest/ogr/ogr_fgdb.py

@@ -3113,3 +3113,16 @@ def test_ogr_filegdb_incompatible_geometry_types(fgdb_drv, layer_geom_type, wkt)
         shutil.rmtree(dirname)
     except OSError:
         pass
+
+
+###############################################################################
+# Test reading an empty polygon
+
+
+def test_ogr_filegdb_read_empty_polygon():
+
+    # Dataset generated by OpenFileGDB driver
+    ds = ogr.Open("data/filegdb/empty_polygon.gdb")
+    lyr = ds.GetLayer(0)
+    f = lyr.GetNextFeature()
+    assert f.GetGeometryRef().IsEmpty()

autotest/ogr/ogr_geom.py

@@ -4066,3 +4066,35 @@ def test_ogr_isclockwise():
     geom = ogr.CreateGeometryFromWkt("POLYGON((0 0,1 0,1 1,0 0))")
     ring = geom.GetGeometryRef(0)
     assert not ring.IsClockwise()
+
+
+###############################################################################
+
+
+@pytest.mark.parametrize(
+    "geom_type",
+    [
+        ogr.wkbPoint,
+        ogr.wkbPoint25D,
+        ogr.wkbPointM,
+        ogr.wkbPointZM,
+        ogr.wkbLineString,
+        ogr.wkbPolygon,
+        ogr.wkbMultiPoint,
+        ogr.wkbMultiLineString,
+        ogr.wkbMultiPolygon,
+        ogr.wkbGeometryCollection,
+        ogr.wkbCircularString,
+        ogr.wkbCompoundCurve,
+        ogr.wkbCurvePolygon,
+        ogr.wkbMultiCurve,
+        ogr.wkbMultiSurface,
+        ogr.wkbPolyhedralSurface,
+        ogr.wkbTIN,
+    ],
+)
+def test_ogr_geom_CreateGeometry(geom_type):
+
+    geom = ogr.Geometry(geom_type)
+    assert geom.GetGeometryType() == geom_type
+    assert geom.IsEmpty()

autotest/ogr/ogr_openfilegdb_write.py

@@ -4259,3 +4259,49 @@ def test_ogr_openfilegdb_write_compound_crs(write_wkid, write_vcswkid):
 
     finally:
         gdal.RmdirRecursive(dirname)
+
+
+###############################################################################
+# Test writing empty geometries
+
+
+@pytest.mark.parametrize(
+    "geom_type",
+    [
+        ogr.wkbPoint,
+        ogr.wkbPoint25D,
+        ogr.wkbPointM,
+        ogr.wkbPointZM,
+        ogr.wkbMultiLineString,
+        ogr.wkbMultiLineString25D,
+        ogr.wkbMultiLineStringM,
+        ogr.wkbMultiLineStringZM,
+        ogr.wkbMultiPolygon,
+        ogr.wkbMultiPolygon25D,
+        ogr.wkbMultiPolygonM,
+        ogr.wkbMultiPolygonZM,
+    ],
+)
+def test_ogr_openfilegdb_write_empty_geoms(geom_type):
+
+    dirname = "/vsimem/test_ogr_openfilegdb_write_empty_geoms.gdb"
+    try:
+        ds = ogr.GetDriverByName("OpenFileGDB").CreateDataSource(dirname)
+        lyr = ds.CreateLayer("test", geom_type=geom_type)
+        f = ogr.Feature(lyr.GetLayerDefn())
+        g = ogr.Geometry(geom_type)
+        f.SetGeometry(g)
+        with gdaltest.config_option("OGR_OPENFILEGDB_WRITE_EMPTY_GEOMETRY", "YES"):
+            lyr.CreateFeature(f)
+        ds = None
+
+        ds = ogr.Open(dirname)
+        lyr = ds.GetLayer(0)
+        assert lyr.GetGeomType() == geom_type
+        f = lyr.GetNextFeature()
+        g = f.GetGeometryRef()
+        assert g.GetGeometryType() == geom_type
+        assert g.IsEmpty()
+
+    finally:
+        gdal.RmdirRecursive(dirname)

ogr/ogrgeometryfactory.cpp

@@ -554,59 +554,85 @@ OGRGeometry *
 OGRGeometryFactory::createGeometry(OGRwkbGeometryType eGeometryType)
 
 {
+    OGRGeometry *poGeom = nullptr;
     switch (wkbFlatten(eGeometryType))
     {
         case wkbPoint:
-            return new (std::nothrow) OGRPoint();
+            poGeom = new (std::nothrow) OGRPoint();
+            break;
 
         case wkbLineString:
-            return new (std::nothrow) OGRLineString();
+            poGeom = new (std::nothrow) OGRLineString();
+            break;
 
         case wkbPolygon:
-            return new (std::nothrow) OGRPolygon();
+            poGeom = new (std::nothrow) OGRPolygon();
+            break;
 
         case wkbGeometryCollection:
-            return new (std::nothrow) OGRGeometryCollection();
+            poGeom = new (std::nothrow) OGRGeometryCollection();
+            break;
 
         case wkbMultiPolygon:
-            return new (std::nothrow) OGRMultiPolygon();
+            poGeom = new (std::nothrow) OGRMultiPolygon();
+            break;
 
         case wkbMultiPoint:
-            return new (std::nothrow) OGRMultiPoint();
+            poGeom = new (std::nothrow) OGRMultiPoint();
+            break;
 
         case wkbMultiLineString:
-            return new (std::nothrow) OGRMultiLineString();
+            poGeom = new (std::nothrow) OGRMultiLineString();
+            break;
 
         case wkbLinearRing:
-            return new (std::nothrow) OGRLinearRing();
+            poGeom = new (std::nothrow) OGRLinearRing();
+            break;
 
         case wkbCircularString:
-            return new (std::nothrow) OGRCircularString();
+            poGeom = new (std::nothrow) OGRCircularString();
+            break;
 
         case wkbCompoundCurve:
-            return new (std::nothrow) OGRCompoundCurve();
+            poGeom = new (std::nothrow) OGRCompoundCurve();
+            break;
 
         case wkbCurvePolygon:
-            return new (std::nothrow) OGRCurvePolygon();
+            poGeom = new (std::nothrow) OGRCurvePolygon();
+            break;
 
         case wkbMultiCurve:
-            return new (std::nothrow) OGRMultiCurve();
+            poGeom = new (std::nothrow) OGRMultiCurve();
+            break;
 
         case wkbMultiSurface:
-            return new (std::nothrow) OGRMultiSurface();
+            poGeom = new (std::nothrow) OGRMultiSurface();
+            break;
 
         case wkbTriangle:
-            return new (std::nothrow) OGRTriangle();
+            poGeom = new (std::nothrow) OGRTriangle();
+            break;
 
         case wkbPolyhedralSurface:
-            return new (std::nothrow) OGRPolyhedralSurface();
+            poGeom = new (std::nothrow) OGRPolyhedralSurface();
+            break;
 
         case wkbTIN:
-            return new (std::nothrow) OGRTriangulatedSurface();
+            poGeom = new (std::nothrow) OGRTriangulatedSurface();
+            break;
 
         default:
-            return nullptr;
+            CPLAssert(false);
+            break;
     }
+    if (poGeom)
+    {
+        if (OGR_GT_HasZ(eGeometryType))
+            poGeom->set3D(true);
+        if (OGR_GT_HasM(eGeometryType))
+            poGeom->setMeasured(true);
+    }
+    return poGeom;
 }
 
 /************************************************************************/

ogr/ogrpgeogeometry.cpp

@@ -2861,6 +2861,10 @@ OGRErr OGRCreateFromShapeBin(GByte *pabyShape, OGRGeometry **ppoGeom,
                     }
                 }
             }
+            else
+            {
+                *ppoGeom = new OGRPolygon();
+            }
         }  // Polygon.
 
         /* --------------------------------------------------------------------

ogr/ogrsf_frmts/openfilegdb/filegdbindex_write.cpp

@@ -1180,7 +1180,7 @@ bool FileGDBTable::CreateSpatialIndex()
             {
                 auto poGeom = std::unique_ptr<OGRGeometry>(
                     poGeomConverter->GetAsGeometry(psField));
-                if (poGeom != nullptr)
+                if (poGeom != nullptr && !poGeom->IsEmpty())
                 {
                     aSetValues.clear();
                     const auto eGeomType =

ogr/ogrsf_frmts/openfilegdb/filegdbtable.cpp

@@ -2528,6 +2528,8 @@ int FileGDBTable::DoesGeometryIntersectsFilterEnvelope(const OGRField *psField)
         {
             GUIntBig x, y;
             ReadVarUInt64NoCheck(pabyCur, x);
+            if (x == 0)  // POINT EMPTY
+                return FALSE;
             x--;
             if (x < m_nFilterXMin || x > m_nFilterXMax)
                 return FALSE;
@@ -3403,28 +3405,29 @@ FileGDBOGRGeometryConverterImpl::GetAsGeometry(const OGRField *psField)
             ReadVarUInt64NoCheck(pabyCur, x);
             ReadVarUInt64NoCheck(pabyCur, y);
 
-            const double dfX =
-                CPLUnsanitizedAdd<GUIntBig>(x, -1) / poGeomField->GetXYScale() +
-                poGeomField->GetXOrigin();
-            const double dfY =
-                CPLUnsanitizedAdd<GUIntBig>(y, -1) / poGeomField->GetXYScale() +
-                poGeomField->GetYOrigin();
+            const double dfX = x == 0 ? std::numeric_limits<double>::quiet_NaN()
+                                      : (x - 1U) / poGeomField->GetXYScale() +
+                                            poGeomField->GetXOrigin();
+            const double dfY = y == 0 ? std::numeric_limits<double>::quiet_NaN()
+                                      : (y - 1U) / poGeomField->GetXYScale() +
+                                            poGeomField->GetYOrigin();
             if (bHasZ)
             {
                 ReadVarUInt64NoCheck(pabyCur, z);
                 const double dfZScale = SanitizeScale(poGeomField->GetZScale());
                 const double dfZ =
-                    CPLUnsanitizedAdd<GUIntBig>(z, -1) / dfZScale +
-                    poGeomField->GetZOrigin();
+                    z == 0 ? std::numeric_limits<double>::quiet_NaN()
+                           : (z - 1U) / dfZScale + poGeomField->GetZOrigin();
                 if (bHasM)
                 {
                     GUIntBig m = 0;
                     ReadVarUInt64NoCheck(pabyCur, m);
                     const double dfMScale =
                         SanitizeScale(poGeomField->GetMScale());
                     const double dfM =
-                        CPLUnsanitizedAdd<GUIntBig>(m, -1) / dfMScale +
-                        poGeomField->GetMOrigin();
+                        m == 0
+                            ? std::numeric_limits<double>::quiet_NaN()
+                            : (m - 1U) / dfMScale + poGeomField->GetMOrigin();
                     return new OGRPoint(dfX, dfY, dfZ, dfM);
                 }
                 return new OGRPoint(dfX, dfY, dfZ);
@@ -3436,8 +3439,8 @@ FileGDBOGRGeometryConverterImpl::GetAsGeometry(const OGRField *psField)
                 ReadVarUInt64NoCheck(pabyCur, m);
                 const double dfMScale = SanitizeScale(poGeomField->GetMScale());
                 const double dfM =
-                    CPLUnsanitizedAdd<GUIntBig>(m, -1) / dfMScale +
-                    poGeomField->GetMOrigin();
+                    m == 0 ? std::numeric_limits<double>::quiet_NaN()
+                           : (m - 1U) / dfMScale + poGeomField->GetMOrigin();
                 poPoint->setM(dfM);
                 return poPoint;
             }

ogr/ogrsf_frmts/openfilegdb/filegdbtable_write.cpp

@@ -578,38 +578,54 @@ bool FileGDBTable::EncodeGeometry(const FileGDBGeomField *poGeomField,
                 }
             }
             const auto poPoint = poGeom->toPoint();
-            double dfVal;
-
-            dfVal = (poPoint->getX() - poGeomField->GetXOrigin()) *
-                        poGeomField->GetXYScale() +
-                    1;
-            CHECK_CAN_BE_ENCODED_ON_VARUINT(dfVal, "Cannot encode value");
-            WriteVarUInt(m_abyGeomBuffer, static_cast<uint64_t>(dfVal + 0.5));
-
-            dfVal = (poPoint->getY() - poGeomField->GetYOrigin()) *
-                        poGeomField->GetXYScale() +
-                    1;
-            CHECK_CAN_BE_ENCODED_ON_VARUINT(dfVal, "Cannot encode Y value");
-            WriteVarUInt(m_abyGeomBuffer, static_cast<uint64_t>(dfVal + 0.5));
-
-            if (bIs3D)
+            if (poPoint->IsEmpty())
             {
-                dfVal = (poPoint->getZ() - poGeomField->GetZOrigin()) *
-                            poGeomField->GetZScale() +
+                WriteUInt8(m_abyGeomBuffer, 0);
+                WriteUInt8(m_abyGeomBuffer, 0);
+                if (bIs3D)
+                    WriteUInt8(m_abyGeomBuffer, 0);
+                if (bIsMeasured)
+                    WriteUInt8(m_abyGeomBuffer, 0);
+            }
+            else
+            {
+                double dfVal;
+
+                dfVal = (poPoint->getX() - poGeomField->GetXOrigin()) *
+                            poGeomField->GetXYScale() +
                         1;
-                CHECK_CAN_BE_ENCODED_ON_VARUINT(dfVal, "Cannot encode Z value");
+                CHECK_CAN_BE_ENCODED_ON_VARUINT(dfVal, "Cannot encode value");
                 WriteVarUInt(m_abyGeomBuffer,
                              static_cast<uint64_t>(dfVal + 0.5));
-            }
 
-            if (bIsMeasured)
-            {
-                dfVal = (poPoint->getM() - poGeomField->GetMOrigin()) *
-                            poGeomField->GetMScale() +
+                dfVal = (poPoint->getY() - poGeomField->GetYOrigin()) *
+                            poGeomField->GetXYScale() +
                         1;
-                CHECK_CAN_BE_ENCODED_ON_VARUINT(dfVal, "Cannot encode M value");
+                CHECK_CAN_BE_ENCODED_ON_VARUINT(dfVal, "Cannot encode Y value");
                 WriteVarUInt(m_abyGeomBuffer,
                              static_cast<uint64_t>(dfVal + 0.5));
+
+                if (bIs3D)
+                {
+                    dfVal = (poPoint->getZ() - poGeomField->GetZOrigin()) *
+                                poGeomField->GetZScale() +
+                            1;
+                    CHECK_CAN_BE_ENCODED_ON_VARUINT(dfVal,
+                                                    "Cannot encode Z value");
+                    WriteVarUInt(m_abyGeomBuffer,
+                                 static_cast<uint64_t>(dfVal + 0.5));
+                }
+
+                if (bIsMeasured)
+                {
+                    dfVal = (poPoint->getM() - poGeomField->GetMOrigin()) *
+                                poGeomField->GetMScale() +
+                            1;
+                    CHECK_CAN_BE_ENCODED_ON_VARUINT(dfVal,
+                                                    "Cannot encode M value");
+                    WriteVarUInt(m_abyGeomBuffer,
+                                 static_cast<uint64_t>(dfVal + 0.5));
+                }
             }
 
             return true;

ogr/ogrsf_frmts/openfilegdb/ogropenfilegdblayer_write.cpp

@@ -2049,8 +2049,12 @@ bool OGROpenFileGDBLayer::PrepareFileGDBFeature(OGRFeature *poFeature,
         }
 
         // Treat empty geometries as NULL, like the FileGDB driver
-        if (poGeom->IsEmpty())
+        if (poGeom->IsEmpty() &&
+            !CPLTestBool(CPLGetConfigOption(
+                "OGR_OPENFILEGDB_WRITE_EMPTY_GEOMETRY", "NO")))
+        {
             poGeom = nullptr;
+        }
     }
 
     if (m_iAreaField >= 0)