GPKG: only quantize geometries if DISCARD_COORD_LSB layer creation op…

…tion is set to YES. Also add a UNDO_DISCARD_COORD_LSB_ON_READING layer creation option

autotest/ogr/ogr_gpkg.py

@@ -10147,16 +10147,28 @@ def test_ogr_gpkg_creation_with_foreign_key_constraint_enabled(tmp_vsimem):
 
 
 @gdaltest.enable_exceptions()
-@pytest.mark.parametrize("with_metadata", [True, False])
-def test_ogr_gpkg_geom_coord_precision(tmp_vsimem, with_metadata):
+@pytest.mark.parametrize(
+    "with_metadata,DISCARD_COORD_LSB,UNDO_DISCARD_COORD_LSB_ON_READING",
+    [(True, "YES", "YES"), (False, "YES", "NO"), (False, "NO", "NO")],
+)
+def test_ogr_gpkg_geom_coord_precision(
+    tmp_vsimem, with_metadata, DISCARD_COORD_LSB, UNDO_DISCARD_COORD_LSB_ON_READING
+):
 
     filename = str(tmp_vsimem / "test.gpkg")
     ds = gdal.GetDriverByName("GPKG").Create(filename, 0, 0, 0, gdal.GDT_Unknown)
     geom_fld = ogr.GeomFieldDefn("my_geom", ogr.wkbUnknown)
     prec = ogr.CreateGeomCoordinatePrecision()
     prec.Set(1e-5, 1e-3, 1e-2)
     geom_fld.SetCoordinatePrecision(prec)
-    lyr = ds.CreateLayerFromGeomFieldDefn("test", geom_fld)
+    lyr = ds.CreateLayerFromGeomFieldDefn(
+        "test",
+        geom_fld,
+        [
+            "DISCARD_COORD_LSB=" + DISCARD_COORD_LSB,
+            "UNDO_DISCARD_COORD_LSB_ON_READING=" + UNDO_DISCARD_COORD_LSB_ON_READING,
+        ],
+    )
     geom_fld = lyr.GetLayerDefn().GetGeomFieldDefn(0)
     prec = geom_fld.GetCoordinatePrecision()
     assert prec.GetXYResolution() == 1e-5
@@ -10183,15 +10195,71 @@ def test_ogr_gpkg_geom_coord_precision(tmp_vsimem, with_metadata):
     assert ds.GetMetadata() == {}
     assert lyr.GetMetadata() == ({"FOO": "BAR"} if with_metadata else {})
     f = lyr.GetNextFeature()
+
+    g = f.GetGeometryRef()
+    assert g.GetX(0) == pytest.approx(1.23456789, abs=1e-5)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetX(0) != pytest.approx(1.23456789, abs=1e-8)
+    if UNDO_DISCARD_COORD_LSB_ON_READING == "YES":
+        assert g.GetX(0) == pytest.approx(1.23457, abs=1e-8)
+    else:
+        assert g.GetX(0) != pytest.approx(1.23457, abs=1e-8)
+
+    assert g.GetY(0) == pytest.approx(2.34567891, abs=1e-5)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetY(0) != pytest.approx(2.34567891, abs=1e-8)
+    if UNDO_DISCARD_COORD_LSB_ON_READING == "YES":
+        assert g.GetY(0) == pytest.approx(2.34568, abs=1e-8)
+
+    assert g.GetZ(0) == pytest.approx(9.87654321, abs=1e-3)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetZ(0) != pytest.approx(9.87654321, abs=1e-8)
+    if UNDO_DISCARD_COORD_LSB_ON_READING == "YES":
+        assert g.GetZ(0) == pytest.approx(9.876, abs=1e-8)
+
+    assert g.GetM(0) == pytest.approx(-1.23456789, abs=1e-2)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetM(0) != pytest.approx(-1.23456789, abs=1e-8)
+    if UNDO_DISCARD_COORD_LSB_ON_READING == "YES":
+        assert g.GetM(0) == pytest.approx(-1.23, abs=1e-8)
+
+    # Test Arrow interface
+    lyr.ResetReading()
+    mem_ds = ogr.GetDriverByName("Memory").CreateDataSource("")
+    mem_lyr = mem_ds.CreateLayer("test", geom_type=ogr.wkbNone)
+    mem_lyr.CreateGeomField(ogr.GeomFieldDefn("my_geom"))
+    mem_lyr.WriteArrow(lyr)
+    f = mem_lyr.GetNextFeature()
+
     g = f.GetGeometryRef()
     assert g.GetX(0) == pytest.approx(1.23456789, abs=1e-5)
-    assert g.GetX(0) != pytest.approx(1.23456789, abs=1e-8)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetX(0) != pytest.approx(1.23456789, abs=1e-8)
+    if UNDO_DISCARD_COORD_LSB_ON_READING == "YES":
+        assert g.GetX(0) == pytest.approx(1.23457, abs=1e-8)
+    else:
+        assert g.GetX(0) != pytest.approx(1.23457, abs=1e-8)
+
     assert g.GetY(0) == pytest.approx(2.34567891, abs=1e-5)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetY(0) != pytest.approx(2.34567891, abs=1e-8)
+    if UNDO_DISCARD_COORD_LSB_ON_READING == "YES":
+        assert g.GetY(0) == pytest.approx(2.34568, abs=1e-8)
+
     assert g.GetZ(0) == pytest.approx(9.87654321, abs=1e-3)
-    assert g.GetZ(0) != pytest.approx(9.87654321, abs=1e-8)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetZ(0) != pytest.approx(9.87654321, abs=1e-8)
+    if UNDO_DISCARD_COORD_LSB_ON_READING == "YES":
+        assert g.GetZ(0) == pytest.approx(9.876, abs=1e-8)
+
     assert g.GetM(0) == pytest.approx(-1.23456789, abs=1e-2)
-    assert g.GetM(0) != pytest.approx(-1.23456789, abs=1e-8)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetM(0) != pytest.approx(-1.23456789, abs=1e-8)
+    if UNDO_DISCARD_COORD_LSB_ON_READING == "YES":
+        assert g.GetM(0) == pytest.approx(-1.23, abs=1e-8)
 
+    lyr.ResetReading()
+    f = lyr.GetNextFeature()
     # Update geometry to check existing precision settings are used
     f.SetGeometry(
         ogr.CreateGeometryFromWkt(
@@ -10225,12 +10293,15 @@ def test_ogr_gpkg_geom_coord_precision(tmp_vsimem, with_metadata):
     f = lyr.GetNextFeature()
     g = f.GetGeometryRef()
     assert g.GetX(0) == pytest.approx(-1.23456789, abs=1e-5)
-    assert g.GetX(0) != pytest.approx(-1.23456789, abs=1e-8)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetX(0) != pytest.approx(-1.23456789, abs=1e-8)
     assert g.GetY(0) == pytest.approx(-2.34567891, abs=1e-5)
     assert g.GetZ(0) == pytest.approx(-9.87654321, abs=1e-3)
-    assert g.GetZ(0) != pytest.approx(-9.87654321, abs=1e-8)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetZ(0) != pytest.approx(-9.87654321, abs=1e-8)
     assert g.GetM(0) == pytest.approx(1.23456789, abs=1e-2)
-    assert g.GetM(0) != pytest.approx(1.23456789, abs=1e-8)
+    if DISCARD_COORD_LSB == "YES":
+        assert g.GetM(0) != pytest.approx(1.23456789, abs=1e-8)
 
     ds.DeleteLayer(0)
 

doc/source/drivers/vector/gpkg.rst

@@ -351,6 +351,20 @@ The following layer creation options are available:
       geometry column can be NULL. Can be set to NO so that geometry is
       required.
 
+-  .. lco:: DISCARD_COORD_LSB
+      :choices: YES, NO
+      :default: NO
+      :since: 3.9
+
+      Whether the geometry coordinate precision should be used to set to zero non-significant least-significant bits of geometries. Helps when further compression is used. See :ref:`ogr_gpkg_geometry_coordinate_precision` for more details.
+
+-  .. lco:: UNDO_DISCARD_COORD_LSB_ON_READING
+      :choices: YES, NO
+      :default: NO
+      :since: 3.9
+
+      Whether to ask GDAL to take into coordinate precision to undo the effects of DISCARD_COORD_LSB. See :ref:`ogr_gpkg_geometry_coordinate_precision` for more details.
+
 -  .. lco:: FID
       :default: fid
 
@@ -632,17 +646,45 @@ Update of an existing file is not supported.
 Creation involves the creation of a temporary file. Sufficiently large files
 will be automatically compressed using the SOZip optimization.
 
+.. _ogr_gpkg_geometry_coordinate_precision:
+
 Geometry coordinate precision
 -----------------------------
 
 .. versionadded:: GDAL 3.9
 
 The GeoPackage driver supports reading and writing the geometry coordinate
 precision, using the :cpp:class:`OGRGeomCoordinatePrecision` settings of the
-:cpp:class:`OGRGeomFieldDefn`, as well as setting to zero the less-significant
-bits which are not relevant for the requested precision. This can improve
-further lossless compression stages, for example when putting a GeoPackage in
-an archive.
+:cpp:class:`OGRGeomFieldDefn`. By default, the geometry coordinate precision
+is only noted in metadata, and does not cause geometries that are written to
+be modified to comply with this precision.
+
+Several settings may be combined to apply further processing:
+
+* if the :config:`OGR_APPLY_GEOM_SET_PRECISION` configuration option is set to
+  ``YES``, the :cpp:func:`OGRGeometry::SetPrecision` method will be applied
+  when calling the CreateFeature() and SetFeature() method of the driver, to
+  round X and Y coordinates to the specified precision, and fix potential
+  geometry invalidities resulting from the rounding.
+
+* if the ``DISCARD_COORD_LSB`` layer creation option is set to YES, the
+  less-significant bits of the WKB geometry encoding which are not relevant for
+  the requested precision are set to zero. This can improve further lossless
+  compression stages, for example when putting a GeoPackage in an archive.
+  Note however that when reading back such geometries and displaying them
+  to the maximum precision, they will not "exactly" match the original
+  :cpp:class:`OGRGeomCoordinatePrecision` settings. However, they will round
+  back to it.
+  The value of the ``DISCARD_COORD_LSB`` layer creation option is written in
+  the dataset metadata, and will be re-used for later edition sessions.
+
+* if the ``UNDO_DISCARD_COORD_LSB_ON_READING`` layer creation option is set to
+  YES (only makes sense if the ``DISCARD_COORD_LSB`` layer creation option is
+  also set to YES), when *reading* back geometries from a dataset, the
+  :cpp:func:`OGRGeometry::roundCoordinates` method will be applied so that
+  the geometry coordinates exactly match the original specified coordinate
+  precision. That option will only be honored by GDAL 3.9 or later.
+
 
 Implementation details: the coordinate precision is stored in a record in each
 of the ``gpkg_metadata`` and ``gpkg_metadata_reference`` table, with the
@@ -651,7 +693,7 @@ specification:
 
 - gpkg_metadata.md_standard_uri = 'http://gdal.org'
 - gpkg_metadata.mime_type = 'text/xml'
-- gpkg_metadata.metadata = '<CoordinatePrecision xy_resolution="{xy_resolution}" z_resolution="{z_resolution}" m_resolution="{m_resolution}" />'
+- gpkg_metadata.metadata = '<CoordinatePrecision xy_resolution="{xy_resolution}" z_resolution="{z_resolution}" m_resolution="{m_resolution}" discard_coord_lsb={true or false} undo_discard_coord_lsb_on_reading={true or false} />'
 - gpkg_metadata_reference.reference_scope = 'column'
 - gpkg_metadata_reference.table_name = '{table_name}'
 - gpkg_metadata_reference.column_name = '{geometry_column_name}'

ogr/ogrsf_frmts/gpkg/ogr_geopackage.h

@@ -574,6 +574,9 @@ class OGRGeoPackageLayer CPL_NON_FINAL : public OGRLayer,
     int m_iGeomCol = -1;
     std::vector<int> m_anFieldOrdinals{};
 
+    //! Whether to call OGRGeometry::SetPrecision() when reading back geometries from the database
+    bool m_bUndoDiscardCoordLSBOnReading = false;
+
     void ClearStatement();
     virtual OGRErr ResetStatement() = 0;
 
@@ -884,13 +887,12 @@ class OGRGeoPackageTableLayer final : public OGRGeoPackageLayer
                               const char *pszObjectType, bool bIsInGpkgContents,
                               bool bIsSpatial, const char *pszGeomColName,
                               const char *pszGeomType, bool bHasZ, bool bHasM);
-    void SetCreationParameters(OGRwkbGeometryType eGType,
-                               const char *pszGeomColumnName, int bGeomNullable,
-                               OGRSpatialReference *poSRS,
-                               const OGRGeomCoordinatePrecision &oCoordPrec,
-                               const char *pszFIDColumnName,
-                               const char *pszIdentifier,
-                               const char *pszDescription);
+    void SetCreationParameters(
+        OGRwkbGeometryType eGType, const char *pszGeomColumnName,
+        int bGeomNullable, OGRSpatialReference *poSRS,
+        const OGRGeomCoordinatePrecision &oCoordPrec, bool bDiscardCoordLSB,
+        bool bUndoDiscardCoordLSBOnReading, const char *pszFIDColumnName,
+        const char *pszIdentifier, const char *pszDescription);
     void SetDeferredSpatialIndexCreation(bool bFlag);
     void SetASpatialVariant(GPKGASpatialVariant eASpatialVariant)
     {

ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp

@@ -6774,6 +6774,10 @@ GDALGeoPackageDataset::ICreateLayer(const char *pszLayerName,
         eGType, pszGeomColumnName, bGeomNullable, poSRS,
         poSrcGeomFieldDefn ? poSrcGeomFieldDefn->GetCoordinatePrecision()
                            : OGRGeomCoordinatePrecision(),
+        CPLTestBool(
+            CSLFetchNameValueDef(papszOptions, "DISCARD_COORD_LSB", "NO")),
+        CPLTestBool(CSLFetchNameValueDef(
+            papszOptions, "UNDO_DISCARD_COORD_LSB_ON_READING", "NO")),
         pszFIDColumnName, pszIdentifier,
         CSLFetchNameValue(papszOptions, "DESCRIPTION"));
     if (poSRS)

ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp

@@ -675,6 +675,13 @@ void RegisterOGRGeoPackage()
         "  <Option name='GEOMETRY_NULLABLE' type='boolean' "
         "description='Whether the values of the geometry column can be NULL' "
         "default='YES'/>"
+        "  <Option name='DISCARD_COORD_LSB' type='boolean' "
+        "description='Whether the geometry coordinate precision should be used "
+        "to set to zero non-significant least-significant bits of geometries. "
+        "Helps when further compression is used' default='NO'/>"
+        "  <Option name='UNDO_DISCARD_COORD_LSB_ON_READING' type='boolean' "
+        "description='Whether to ask GDAL to take into coordinate precision to "
+        "undo the effects of DISCARD_COORD_LSB' default='NO'/>"
         "  <Option name='FID' type='string' description='Name of the FID "
         "column to create' default='fid'/>"
         "  <Option name='OVERWRITE' type='boolean' description='Whether to "

ogr/ogrsf_frmts/gpkg/ogrgeopackagelayer.cpp

@@ -397,8 +397,16 @@ OGRFeature *OGRGeoPackageLayer::TranslateFeature(sqlite3_stmt *hStmt)
                              "Unable to read geometry");
                 }
             }
-            if (poGeom != nullptr)
+            if (poGeom)
+            {
+                if (m_bUndoDiscardCoordLSBOnReading)
+                {
+                    poGeom->roundCoordinates(
+                        poGeomFieldDefn->GetCoordinatePrecision());
+                }
                 poGeom->assignSpatialReference(poSrs);
+            }
+
             poFeature->SetGeometryDirectly(poGeom);
         }
     }
@@ -633,7 +641,8 @@ int OGRGeoPackageLayer::GetNextArrowArray(struct ArrowArrayStream *stream,
                 const GByte *pabyGpkg = static_cast<const GByte *>(
                     sqlite3_column_blob(hStmt, m_iGeomCol));
                 if (m_poFilterGeom == nullptr && iGpkgSize >= 8 && pabyGpkg &&
-                    pabyGpkg[0] == 'G' && pabyGpkg[1] == 'P')
+                    pabyGpkg[0] == 'G' && pabyGpkg[1] == 'P' &&
+                    !m_bUndoDiscardCoordLSBOnReading)
                 {
                     GPkgHeader oHeader;
 
@@ -663,6 +672,12 @@ int OGRGeoPackageLayer::GetNextArrowArray(struct ArrowArrayStream *stream,
                         }
                         poGeom.reset(poGeomPtr);
                     }
+                    else if (m_bUndoDiscardCoordLSBOnReading)
+                    {
+                        poGeom->roundCoordinates(
+                            m_poFeatureDefn->GetGeomFieldDefn(0)
+                                ->GetCoordinatePrecision());
+                    }
                     if (poGeom != nullptr)
                     {
                         nWKBSize = poGeom->WkbSize();