diff --git a/autotest/ogr/ogr_gpkg.py b/autotest/ogr/ogr_gpkg.py
index 93d23426290e..67fb786af1ba 100755
--- a/autotest/ogr/ogr_gpkg.py
+++ b/autotest/ogr/ogr_gpkg.py
@@ -10147,8 +10147,13 @@ 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)
@@ -10156,7 +10161,14 @@ def test_ogr_gpkg_geom_coord_precision(tmp_vsimem, with_metadata):
     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)
 
diff --git a/doc/source/drivers/vector/gpkg.rst b/doc/source/drivers/vector/gpkg.rst
index 49e93ac4f8e7..a6b67133b513 100644
--- a/doc/source/drivers/vector/gpkg.rst
+++ b/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,6 +646,8 @@ 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
 -----------------------------
 
@@ -639,10 +655,36 @@ Geometry coordinate precision
 
 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}'
diff --git a/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h b/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h
index 4367680810f0..f8c7688201aa 100644
--- a/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h
+++ b/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)
     {
diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp
index e6c4d5a92d28..7b122f5e090a 100644
--- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp
+++ b/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)
diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp
index e68d6fa37ccc..e5e3043bad72 100644
--- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp
+++ b/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 "
diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagelayer.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagelayer.cpp
index 4bc3f00ee879..82c578364829 100644
--- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagelayer.cpp
+++ b/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();
diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp
index e9f24d9a21ef..45092bfcbe45 100644
--- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp
+++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp
@@ -1477,7 +1477,16 @@ OGRErr OGRGeoPackageTableLayer::ReadTableDefinition()
                         }
                         m_poFeatureDefn->GetGeomFieldDefn(iGeomCol)
                             ->SetCoordinatePrecision(sCoordPrec);
-                        m_sBinaryPrecision.SetFrom(sCoordPrec);
+                        if (CPLTestBool(CPLGetXMLValue(
+                                psXMLNode.get(), "discard_coord_lsb", "false")))
+                        {
+                            m_sBinaryPrecision.SetFrom(sCoordPrec);
+                            m_bUndoDiscardCoordLSBOnReading =
+                                CPLTestBool(CPLGetXMLValue(
+                                    psXMLNode.get(),
+                                    "undo_discard_coord_lsb_on_reading",
+                                    "false"));
+                        }
                     }
                 }
             }
@@ -5594,6 +5603,7 @@ void OGRGeoPackageTableLayer::SetOpeningParameters(
 void OGRGeoPackageTableLayer::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)
 {
@@ -5606,6 +5616,7 @@ void OGRGeoPackageTableLayer::SetCreationParameters(
     m_bTableCreatedInTransaction = m_poDS->IsInTransaction();
     m_bHasTriedDetectingFID64 = true;
     m_pszFidColumn = CPLStrdup(pszFIDColumnName);
+    m_bUndoDiscardCoordLSBOnReading = bUndoDiscardCoordLSBOnReading;
 
     if (eGType != wkbNone)
     {
@@ -5617,7 +5628,9 @@ void OGRGeoPackageTableLayer::SetCreationParameters(
         oGeomFieldDefn.SetSpatialRef(poSRS);
         oGeomFieldDefn.SetNullable(bGeomNullable);
         oGeomFieldDefn.SetCoordinatePrecision(oCoordPrec);
-        m_sBinaryPrecision.SetFrom(oCoordPrec);
+
+        if (bDiscardCoordLSB)
+            m_sBinaryPrecision.SetFrom(oCoordPrec);
 
         // Save coordinate precision in gpkg_metadata/gpkg_metadata_reference
         if ((oCoordPrec.dfXYResolution != OGRGeomCoordinatePrecision::UNKNOWN ||
@@ -5636,6 +5649,11 @@ void OGRGeoPackageTableLayer::SetCreationParameters(
             if (oCoordPrec.dfMResolution != OGRGeomCoordinatePrecision::UNKNOWN)
                 osCoordPrecision += CPLSPrintf(" m_resolution=\"%g\"",
                                                oCoordPrec.dfMResolution);
+            osCoordPrecision += CPLSPrintf(" discard_coord_lsb=\"%s\"",
+                                           bDiscardCoordLSB ? "true" : "false");
+            osCoordPrecision +=
+                CPLSPrintf(" undo_discard_coord_lsb_on_reading=\"%s\"",
+                           m_bUndoDiscardCoordLSBOnReading ? "true" : "false");
             osCoordPrecision += " />";
 
             char *pszSQL = sqlite3_mprintf(
@@ -7742,13 +7760,41 @@ void OGR_GPKG_FillArrowArray_Step(sqlite3_context *pContext, int /*argc*/,
             if (nBlobSize >= 8 && pabyBlob && pabyBlob[0] == 'G' &&
                 pabyBlob[1] == 'P')
             {
-                /* Read header */
-                OGRErr err = GPkgHeaderFromWKB(pabyBlob, nBlobSize, &oHeader);
-                if (err == OGRERR_NONE)
+                if (psFillArrowArray->poLayer->m_bUndoDiscardCoordLSBOnReading)
                 {
-                    /* WKB pointer */
-                    pabyWkb = pabyBlob + oHeader.nHeaderLen;
-                    nWKBSize = nBlobSize - oHeader.nHeaderLen;
+                    OGRGeometry *poGeomPtr =
+                        GPkgGeometryToOGR(pabyBlob, nBlobSize, nullptr);
+                    if (poGeomPtr)
+                    {
+                        poGeomPtr->roundCoordinates(
+                            psFillArrowArray->poFeatureDefn->GetGeomFieldDefn(0)
+                                ->GetCoordinatePrecision());
+                        nWKBSize = poGeomPtr->WkbSize();
+                        abyWkb.resize(nWKBSize);
+                        if (poGeomPtr->exportToWkb(wkbNDR, abyWkb.data(),
+                                                   wkbVariantIso) !=
+                            OGRERR_NONE)
+                        {
+                            nWKBSize = 0;
+                        }
+                        else
+                        {
+                            pabyWkb = abyWkb.data();
+                        }
+                        delete poGeomPtr;
+                    }
+                }
+                else
+                {
+                    /* Read header */
+                    OGRErr err =
+                        GPkgHeaderFromWKB(pabyBlob, nBlobSize, &oHeader);
+                    if (err == OGRERR_NONE)
+                    {
+                        /* WKB pointer */
+                        pabyWkb = pabyBlob + oHeader.nHeaderLen;
+                        nWKBSize = nBlobSize - oHeader.nHeaderLen;
+                    }
                 }
             }
             else if (nBlobSize > 0 && pabyBlob)