diff --git a/autotest/gdrivers/gpkg.py b/autotest/gdrivers/gpkg.py
index 14baa64162c7..1a9bb542a4cf 100755
--- a/autotest/gdrivers/gpkg.py
+++ b/autotest/gdrivers/gpkg.py
@@ -1383,8 +1383,7 @@ def test_gpkg_15():
     out_ds = None
 
     out_ds = gdal.Open("/vsimem/tmp.gpkg")
-    assert out_ds.GetSpatialRef().IsLocal()
-    assert out_ds.GetProjectionRef().find("Undefined Cartesian SRS") >= 0
+    assert out_ds.GetSpatialRef() is None
     # Test setting on read-only dataset
     with gdal.quiet_errors():
         ret = out_ds.SetProjection("")
diff --git a/autotest/ogr/ogr_gpkg.py b/autotest/ogr/ogr_gpkg.py
index 63b4954320c8..394ca882730c 100755
--- a/autotest/ogr/ogr_gpkg.py
+++ b/autotest/ogr/ogr_gpkg.py
@@ -206,7 +206,9 @@ def point_no_spi_but_with_dashes(gpkg_ds):
 @pytest.fixture()
 def point_with_spi_and_dashes(gpkg_ds):
 
-    lyr = gpkg_ds.CreateLayer("point-with-spi-and-dashes", geom_type=ogr.wkbPoint)
+    lyr = gpkg_ds.CreateLayer(
+        "point-with-spi-and-dashes", geom_type=ogr.wkbPoint, options=["SRID=0"]
+    )
     assert lyr.TestCapability(ogr.OLCFastSpatialFilter) == 1
     feat = ogr.Feature(lyr.GetLayerDefn())
     feat.SetGeometry(ogr.CreateGeometryFromWkt("POINT(1000 30000000)"))
@@ -2172,6 +2174,131 @@ def test_ogr_gpkg_write_srs_undefined_Cartesian(tmp_path):
     gpkg_ds = None
 
 
+###############################################################################
+# Test writing a None SRS
+
+
+@pytest.mark.parametrize("crs_wkt_extension", [True, False])
+@gdaltest.enable_exceptions()
+def test_ogr_gpkg_write_no_srs(tmp_path, crs_wkt_extension):
+
+    fname = tmp_path / "test_ogr_gpkg_write_no_srs.gpkg"
+
+    options = []
+    if crs_wkt_extension:
+        options += ["CRS_WKT_EXTENSION=YES"]
+    gpkg_ds = gdaltest.gpkg_dr.CreateDataSource(fname, options=options)
+    assert gpkg_ds is not None
+
+    lyr = gpkg_ds.CreateLayer("layer1", geom_type=ogr.wkbPoint, srs=None)
+    assert lyr.GetSpatialRef() is None
+
+    lyr = gpkg_ds.CreateLayer("layer2", geom_type=ogr.wkbPoint, srs=None)
+    assert lyr.GetSpatialRef() is None
+
+    srs = osr.SpatialReference()
+    srs.SetFromUserInput(
+        'LOCAL_CS["Undefined SRS",LOCAL_DATUM["unknown",32767],UNIT["unknown",0],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
+    )
+    gpkg_ds.CreateLayer("layer3", geom_type=ogr.wkbPoint, srs=srs)
+
+    gpkg_ds = None
+    gpkg_ds = ogr.Open(fname)
+
+    # Check no unexpected SRS entries have been inserted into gpkg_spatial_ref_sys
+    with gpkg_ds.ExecuteSQL("SELECT COUNT(*) FROM gpkg_spatial_ref_sys") as sql_lyr:
+        assert sql_lyr.GetNextFeature().GetField(0) == 4
+
+    # Check no new SRS entries have been inserted into gpkg_spatial_ref_sys
+    with gpkg_ds.ExecuteSQL(
+        "SELECT * FROM gpkg_spatial_ref_sys WHERE srs_id = 99999"
+    ) as sql_lyr:
+        f = sql_lyr.GetNextFeature()
+        assert f["srs_name"] == "Undefined SRS"
+        assert f["organization"] == "GDAL"
+        assert f["organization_coordsys_id"] == 99999
+        assert (
+            f["definition"]
+            == 'LOCAL_CS["Undefined SRS",LOCAL_DATUM["unknown",32767],UNIT["unknown",0],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
+        )
+        if crs_wkt_extension:
+            assert (
+                f["definition_12_063"]
+                == 'ENGCRS["Undefined SRS",EDATUM["unknown"],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["unknown",0]],AXIS["northing",north,ORDER[2],LENGTHUNIT["unknown",0]]]'
+            )
+        assert f["description"] == "Custom undefined coordinate reference system"
+
+    for lyr in gpkg_ds:
+        assert lyr.GetSpatialRef() is None
+
+    gpkg_ds = None
+
+
+###############################################################################
+# Test SRID layer creation option
+
+
+@gdaltest.enable_exceptions()
+@pytest.mark.parametrize(
+    "srid,expected_wkt",
+    [
+        (-1, 'ENGCRS["Undefined Cartesian SRS",'),
+        (0, 'GEOGCRS["Undefined geographic SRS",'),
+        (4326, 'GEOGCRS["WGS 84",'),
+        (4258, 'GEOGCRS["ETRS89",'),
+        (1, None),
+        (99999, None),
+        (123456, None),
+    ],
+)
+def test_ogr_gpkg_SRID_creation_option(tmp_path, srid, expected_wkt):
+
+    fname = tmp_path / "test_ogr_gpkg_SRID_creation_option.gpkg"
+
+    gpkg_ds = gdaltest.gpkg_dr.CreateDataSource(fname)
+    assert gpkg_ds is not None
+
+    if srid in (1, 123456):
+        with gdal.quiet_errors():
+            lyr = gpkg_ds.CreateLayer(
+                "layer1", geom_type=ogr.wkbPoint, options=[f"SRID={srid}"]
+            )
+        assert (
+            gdal.GetLastErrorMsg()
+            == f"No entry in gpkg_spatial_ref_sys matching SRID={srid}"
+        )
+    else:
+        lyr = gpkg_ds.CreateLayer(
+            "layer1", geom_type=ogr.wkbPoint, options=[f"SRID={srid}"]
+        )
+    got_srs = lyr.GetSpatialRef()
+    if expected_wkt is None:
+        assert got_srs is None
+    else:
+        assert got_srs.ExportToWkt(["FORMAT=WKT2_2019"]).startswith(expected_wkt)
+
+    gpkg_ds = None
+    gpkg_ds = ogr.Open(fname)
+
+    if srid in (1, 123456):
+        with gdal.quiet_errors():
+            lyr = gpkg_ds.GetLayer(0)
+            got_srs = lyr.GetSpatialRef()
+        assert (
+            gdal.GetLastErrorMsg()
+            == f"unable to read srs_id '{srid}' from gpkg_spatial_ref_sys"
+        )
+    else:
+        lyr = gpkg_ds.GetLayer(0)
+        got_srs = lyr.GetSpatialRef()
+    if expected_wkt is None:
+        assert got_srs is None
+    else:
+        assert got_srs.ExportToWkt(["FORMAT=WKT2_2019"]).startswith(expected_wkt)
+
+    gpkg_ds = None
+
+
 ###############################################################################
 # Test maximum width of text fields
 
@@ -6619,7 +6746,7 @@ def test_ogr_gpkg_views(tmp_vsimem, tmp_path):
 
     filename = tmp_vsimem / "test_ogr_gpkg_views.gpkg"
     ds = gdaltest.gpkg_dr.CreateDataSource(filename)
-    lyr = ds.CreateLayer("foo", geom_type=ogr.wkbPoint)
+    lyr = ds.CreateLayer("foo", geom_type=ogr.wkbPoint, options=["SRID=0"])
     lyr.CreateField(ogr.FieldDefn("str_field"))
     f = ogr.Feature(lyr.GetLayerDefn())
     f.SetGeometry(ogr.CreateGeometryFromWkt("POINT(0 0)"))
@@ -7846,8 +7973,7 @@ def test_ogr_gpkg_alter_geom_field_defn(tmp_vsimem, tmp_path):
     ds = ogr.Open(filename, update=1)
     lyr = ds.GetLayer(0)
     srs = lyr.GetSpatialRef()
-    assert srs is not None
-    assert srs.GetName() == "Undefined geographic SRS"
+    assert srs is None
 
     new_geom_field_defn = ogr.GeomFieldDefn("", ogr.wkbNone)
     new_geom_field_defn.SetSpatialRef(srs_4326)
diff --git a/autotest/ogr/ogr_mssqlspatial.py b/autotest/ogr/ogr_mssqlspatial.py
index 76595301570c..9a42741cd9f4 100755
--- a/autotest/ogr/ogr_mssqlspatial.py
+++ b/autotest/ogr/ogr_mssqlspatial.py
@@ -678,7 +678,14 @@ def test_binary_field_bcp(mssql_ds):
             res_ds = gdal.VectorTranslate(
                 mssql_ds.GetDescription(),
                 filename,
-                options=["-nln", "test_binary_field", "-lco", "OVERWRITE=YES"],
+                options=[
+                    "-nln",
+                    "test_binary_field",
+                    "-lco",
+                    "OVERWRITE=YES",
+                    "-a_srs",
+                    "EPSG:4326",
+                ],
             )
 
             assert res_ds is not None
diff --git a/doc/source/drivers/vector/gpkg.rst b/doc/source/drivers/vector/gpkg.rst
index 8f08c714f9f3..cf33fb484284 100644
--- a/doc/source/drivers/vector/gpkg.rst
+++ b/doc/source/drivers/vector/gpkg.rst
@@ -357,6 +357,16 @@ The following layer creation options are available:
       geometry column can be NULL. Can be set to NO so that geometry is
       required.
 
+-  .. lco:: SRID
+      :choices: <integer>
+      :since: 3.9
+
+      Forced ``srs_id`` of the entry in the ``gpkg_spatial_ref_sys`` table to point to.
+      This may be -1 ("Undefined Cartesian SRS"), 0 ("Undefined geographic SRS"),
+      99999 ("Undefined SRS"), a valid EPSG CRS code or an existing entry of the
+      ``gpkg_spatial_ref_sys`` table. If pointing to a non-existing entry, only a warning
+      will be emitted.
+
 -  .. lco:: DISCARD_COORD_LSB
       :choices: YES, NO
       :default: NO
@@ -594,12 +604,27 @@ also supported by GeoPackage and stored in the ``gpkg_spatial_ref_sys`` table.
 
 Two special hard-coded CRS are reserved per the GeoPackage specification:
 
-- SRID 0, for a Undefined Geographic CRS. This one is selected by default if
-  creating a spatial layer without any explicit CRS
+- srs_id=0, for a Undefined Geographic CRS. For GDAL 3.8 or earlier, this one is
+  selected by default if creating a spatial layer without any explicit CRS
 
-- SRID -1, for a Undefined Projected CRS. It might be selected by creating a
+- srs_id=-1, for a Undefined Projected CRS. It might be selected by creating a
   layer with a CRS instantiated from the following WKT string:
-  ``LOCAL_CS["Undefined cartesian SRS"]``. (GDAL >= 3.3)
+  ``LOCAL_CS["Undefined Cartesian SRS"]``. (GDAL >= 3.3)
+
+Starting with GDAL 3.9, a layer without any explicit CRS is mapped from/to a
+custom entry of srs_id=99999 with the following properties:
+
+- ``srs_name``: ``Undefined SRS``
+- ``organization``: ``GDAL``
+- ``organization_coordsys_id``: 99999
+- ``definition``: ``LOCAL_CS["Undefined SRS",LOCAL_DATUM["unknown",32767],UNIT["unknown",0],AXIS["Easting",EAST],AXIS["Northing",NORTH]]``
+- ``definition_12_063`` (when the CRS WKT extension is used): ``ENGCRS["Undefined SRS",EDATUM["unknown"],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["unknown",0]],AXIS["northing",north,ORDER[2],LENGTHUNIT["unknown",0]]]``
+- ``description``: ``Custom undefined coordinate reference system``
+
+Note that the use of a LOCAL_CS / EngineeringCRS is mostly to provide a valid
+CRS definition to comply with the requirements of the GeoPackage specification
+and to be compatible of other applications (or GDAL 3.8 or earlier), but the
+semantics of that entry is intented to be "undefined SRS of any kind".
 
 Level of support of GeoPackage Extensions
 -----------------------------------------
diff --git a/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h b/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h
index a25517eb0446..617bd019b92e 100644
--- a/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h
+++ b/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h
@@ -374,10 +374,13 @@ class GDALGeoPackageDataset final : public OGRSQLiteBaseDataSource,
         return nSoftTransactionLevel > 0;
     }
 
-    int GetSrsId(const OGRSpatialReference &oSRS);
+    // At least 100000 to avoid conflicting with EPSG codes
+    static constexpr int FIRST_CUSTOM_SRSID = 100000;
+
+    int GetSrsId(const OGRSpatialReference *poSRS);
     const char *GetSrsName(const OGRSpatialReference &oSRS);
-    OGRSpatialReference *GetSpatialRef(int iSrsId,
-                                       bool bFallbackToEPSG = false);
+    OGRSpatialReference *GetSpatialRef(int iSrsId, bool bFallbackToEPSG = false,
+                                       bool bEmitErrorIfNotFound = true);
     OGRErr CreateExtensionsTableIfNecessary();
     bool HasExtensionsTable();
 
@@ -918,10 +921,11 @@ class OGRGeoPackageTableLayer final : public OGRGeoPackageLayer
                               const char *pszGeomType, bool bHasZ, bool bHasM);
     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);
+        int bGeomNullable, const OGRSpatialReference *poSRS,
+        const char *pszSRID, 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 91312d5e6e70..9701755e14de 100644
--- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp
+++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp
@@ -251,8 +251,9 @@ static OGRErr GDALGPKGImportFromEPSG(OGRSpatialReference *poSpatialRef,
     return eErr;
 }
 
-OGRSpatialReference *GDALGeoPackageDataset::GetSpatialRef(int iSrsId,
-                                                          bool bFallbackToEPSG)
+OGRSpatialReference *
+GDALGeoPackageDataset::GetSpatialRef(int iSrsId, bool bFallbackToEPSG,
+                                     bool bEmitErrorIfNotFound)
 {
     std::map<int, OGRSpatialReference *>::const_iterator oIter =
         m_oMapSrsIdToSrs.find(iSrsId);
@@ -288,7 +289,8 @@ OGRSpatialReference *GDALGeoPackageDataset::GetSpatialRef(int iSrsId,
     }
 
     CPLString oSQL;
-    oSQL.Printf("SELECT definition, organization, organization_coordsys_id%s%s "
+    oSQL.Printf("SELECT srs_name, definition, organization, "
+                "organization_coordsys_id%s%s "
                 "FROM gpkg_spatial_ref_sys WHERE "
                 "srs_id = %d LIMIT 2",
                 m_bHasDefinition12_063 ? ", definition_12_063" : "",
@@ -311,7 +313,7 @@ OGRSpatialReference *GDALGeoPackageDataset::GetSpatialRef(int iSrsId,
             }
             poSRS->Release();
         }
-        else
+        else if (bEmitErrorIfNotFound)
         {
             CPLError(CE_Warning, CPLE_AppDefined,
                      "unable to read srs_id '%d' from gpkg_spatial_ref_sys",
@@ -321,17 +323,23 @@ OGRSpatialReference *GDALGeoPackageDataset::GetSpatialRef(int iSrsId,
         return nullptr;
     }
 
-    const char *pszWkt = oResult->GetValue(0, 0);
+    const char *pszName = oResult->GetValue(0, 0);
+    if (pszName && EQUAL(pszName, "Undefined SRS"))
+    {
+        m_oMapSrsIdToSrs[iSrsId] = nullptr;
+        return nullptr;
+    }
+    const char *pszWkt = oResult->GetValue(1, 0);
     if (pszWkt == nullptr)
         return nullptr;
-    const char *pszOrganization = oResult->GetValue(1, 0);
-    const char *pszOrganizationCoordsysID = oResult->GetValue(2, 0);
+    const char *pszOrganization = oResult->GetValue(2, 0);
+    const char *pszOrganizationCoordsysID = oResult->GetValue(3, 0);
     const char *pszWkt2 =
-        m_bHasDefinition12_063 ? oResult->GetValue(3, 0) : nullptr;
+        m_bHasDefinition12_063 ? oResult->GetValue(4, 0) : nullptr;
     if (pszWkt2 && !EQUAL(pszWkt2, "undefined"))
         pszWkt = pszWkt2;
     const char *pszCoordinateEpoch =
-        m_bHasEpochColumn ? oResult->GetValue(4, 0) : nullptr;
+        m_bHasEpochColumn ? oResult->GetValue(5, 0) : nullptr;
     const double dfCoordinateEpoch =
         pszCoordinateEpoch ? CPLAtof(pszCoordinateEpoch) : 0.0;
 
@@ -537,14 +545,75 @@ bool GDALGeoPackageDataset::ConvertGpkgSpatialRefSysToExtensionWkt2(
     return bRet;
 }
 
-int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference &oSRS)
-{
-    std::unique_ptr<OGRSpatialReference> poSRS(oSRS.Clone());
+int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference *poSRSIn)
+{
+    const char *pszName = poSRSIn ? poSRSIn->GetName() : nullptr;
+    if (!poSRSIn || poSRSIn->IsEmpty() ||
+        (pszName && EQUAL(pszName, "Undefined SRS")))
+    {
+        OGRErr err = OGRERR_NONE;
+        const int nSRSId = SQLGetInteger(
+            hDB,
+            "SELECT srs_id FROM gpkg_spatial_ref_sys WHERE srs_name = "
+            "'Undefined SRS' AND organization = 'GDAL'",
+            &err);
+        if (err == OGRERR_NONE)
+            return nSRSId;
+
+        // The below WKT definitions are somehow questionable (using a unknown
+        // unit). For GDAL >= 3.9, they won't be used. They will only be used
+        // for earlier versions.
+        const char *pszSQL;
+#define UNDEFINED_CRS_SRS_ID 99999
+        static_assert(UNDEFINED_CRS_SRS_ID == FIRST_CUSTOM_SRSID - 1);
+#define STRINGIFY(x) #x
+#define XSTRINGIFY(x) STRINGIFY(x)
+        if (m_bHasDefinition12_063)
+        {
+            /* clang-format off */
+            pszSQL =
+                "INSERT INTO gpkg_spatial_ref_sys "
+                "(srs_name,srs_id,organization,organization_coordsys_id,"
+                "definition, definition_12_063, description) VALUES "
+                "('Undefined SRS'," XSTRINGIFY(UNDEFINED_CRS_SRS_ID) ",'GDAL',"
+                XSTRINGIFY(UNDEFINED_CRS_SRS_ID) ","
+                "'LOCAL_CS[\"Undefined SRS\",LOCAL_DATUM[\"unknown\",32767],"
+                "UNIT[\"unknown\",0],AXIS[\"Easting\",EAST],"
+                "AXIS[\"Northing\",NORTH]]',"
+                "'ENGCRS[\"Undefined SRS\",EDATUM[\"unknown\"],CS[Cartesian,2],"
+                "AXIS[\"easting\",east,ORDER[1],LENGTHUNIT[\"unknown\",0]],"
+                "AXIS[\"northing\",north,ORDER[2],LENGTHUNIT[\"unknown\",0]]]',"
+                "'Custom undefined coordinate reference system')";
+            /* clang-format on */
+        }
+        else
+        {
+            /* clang-format off */
+            pszSQL =
+                "INSERT INTO gpkg_spatial_ref_sys "
+                "(srs_name,srs_id,organization,organization_coordsys_id,"
+                "definition, description) VALUES "
+                "('Undefined SRS'," XSTRINGIFY(UNDEFINED_CRS_SRS_ID) ",'GDAL',"
+                XSTRINGIFY(UNDEFINED_CRS_SRS_ID) ","
+                "'LOCAL_CS[\"Undefined SRS\",LOCAL_DATUM[\"unknown\",32767],"
+                "UNIT[\"unknown\",0],AXIS[\"Easting\",EAST],"
+                "AXIS[\"Northing\",NORTH]]',"
+                "'Custom undefined coordinate reference system')";
+            /* clang-format on */
+        }
+        if (SQLCommand(hDB, pszSQL) == OGRERR_NONE)
+            return UNDEFINED_CRS_SRS_ID;
+#undef UNDEFINED_CRS_SRS_ID
+#undef XSTRINGIFY
+#undef STRINGIFY
+        return -1;
+    }
+
+    std::unique_ptr<OGRSpatialReference> poSRS(poSRSIn->Clone());
 
     if (poSRS->IsGeographic() || poSRS->IsLocal())
     {
         // See corresponding tests in GDALGeoPackageDataset::GetSpatialRef
-        const char *pszName = poSRS->GetName();
         if (pszName != nullptr && strlen(pszName) > 0)
         {
             if (EQUAL(pszName, "Undefined geographic SRS"))
@@ -575,7 +644,7 @@ int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference &oSRS)
             }
         }
 
-        poSRS->SetCoordinateEpoch(oSRS.GetCoordinateEpoch());
+        poSRS->SetCoordinateEpoch(poSRSIn->GetCoordinateEpoch());
     }
 
     // Check whether the EPSG authority code is already mapped to a
@@ -610,7 +679,7 @@ int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference &oSRS)
     }
 
     if (pszAuthorityName != nullptr && strlen(pszAuthorityName) > 0 &&
-        oSRS.GetCoordinateEpoch() == 0)
+        poSRSIn->GetCoordinateEpoch() == 0)
     {
         pszSQL =
             sqlite3_mprintf("SELECT srs_id FROM gpkg_spatial_ref_sys WHERE "
@@ -658,10 +727,10 @@ int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference &oSRS)
     const char *const apszOptionsWkt2_2019[] = {"FORMAT=WKT2_2019", nullptr};
 
     std::string osEpochTest;
-    if (oSRS.GetCoordinateEpoch() > 0 && m_bHasEpochColumn)
+    if (poSRSIn->GetCoordinateEpoch() > 0 && m_bHasEpochColumn)
     {
         osEpochTest =
-            CPLSPrintf(" AND epoch = %.18g", oSRS.GetCoordinateEpoch());
+            CPLSPrintf(" AND epoch = %.18g", poSRSIn->GetCoordinateEpoch());
     }
 
     if (!(poSRS->IsGeographic() && poSRS->GetAxesCount() == 3))
@@ -698,7 +767,7 @@ int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference &oSRS)
         return DEFAULT_SRID;
     }
 
-    if (oSRS.GetCoordinateEpoch() == 0 || m_bHasEpochColumn)
+    if (poSRSIn->GetCoordinateEpoch() == 0 || m_bHasEpochColumn)
     {
         // Search if there is already an existing entry with this WKT
         if (m_bHasDefinition12_063 && (pszWKT2_2015 || pszWKT2_2019))
@@ -746,7 +815,7 @@ int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference &oSRS)
     }
 
     if (pszAuthorityName != nullptr && strlen(pszAuthorityName) > 0 &&
-        oSRS.GetCoordinateEpoch() == 0)
+        poSRSIn->GetCoordinateEpoch() == 0)
     {
         bool bTryToReuseSRSId = true;
         if (EQUAL(pszAuthorityName, "EPSG"))
@@ -796,7 +865,7 @@ int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference &oSRS)
     }
 
     // Add epoch column if needed
-    if (oSRS.GetCoordinateEpoch() > 0 && !m_bHasEpochColumn)
+    if (poSRSIn->GetCoordinateEpoch() > 0 && !m_bHasEpochColumn)
     {
         if (m_bHasDefinition12_063)
         {
@@ -851,27 +920,27 @@ int GDALGeoPackageDataset::GetSrsId(const OGRSpatialReference &oSRS)
         // Get the current maximum srid in the srs table.
         const int nMaxSRSId = SQLGetInteger(
             hDB, "SELECT MAX(srs_id) FROM gpkg_spatial_ref_sys", nullptr);
-        // At least 100000 to avoid conflicting with EPSG codes
-        nSRSId = std::max(100000, nMaxSRSId + 1);
+        nSRSId = std::max(FIRST_CUSTOM_SRSID, nMaxSRSId + 1);
     }
 
     std::string osEpochColumn;
     std::string osEpochVal;
-    if (oSRS.GetCoordinateEpoch() > 0)
+    if (poSRSIn->GetCoordinateEpoch() > 0)
     {
         osEpochColumn = ", epoch";
-        osEpochVal = CPLSPrintf(", %.18g", oSRS.GetCoordinateEpoch());
+        osEpochVal = CPLSPrintf(", %.18g", poSRSIn->GetCoordinateEpoch());
     }
 
     // Add new SRS row to gpkg_spatial_ref_sys.
     if (m_bHasDefinition12_063)
     {
         // Force WKT2_2019 when we have a dynamic CRS and coordinate epoch
-        const char *pszWKT2 =
-            oSRS.IsDynamic() && oSRS.GetCoordinateEpoch() > 0 && pszWKT2_2019
-                ? pszWKT2_2019.get()
-            : pszWKT2_2015 ? pszWKT2_2015.get()
-                           : pszWKT2_2019.get();
+        const char *pszWKT2 = poSRSIn->IsDynamic() &&
+                                      poSRSIn->GetCoordinateEpoch() > 0 &&
+                                      pszWKT2_2019
+                                  ? pszWKT2_2019.get()
+                              : pszWKT2_2015 ? pszWKT2_2015.get()
+                                             : pszWKT2_2019.get();
 
         if (pszAuthorityName != nullptr && nAuthorityCode > 0)
         {
@@ -3106,16 +3175,7 @@ CPLErr GDALGeoPackageDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
         return CE_Failure;
     }
 
-    int nSRID = -1;
-    if (poSRS == nullptr || poSRS->IsEmpty())
-    {
-        // nSRID = -1;
-    }
-    else
-    {
-        nSRID = GetSrsId(*poSRS);
-    }
-
+    const int nSRID = GetSrsId(poSRS);
     const auto poTS = GetTilingScheme(m_osTilingScheme);
     if (poTS && nSRID != poTS->nEPSGCode)
     {
@@ -6774,6 +6834,7 @@ GDALGeoPackageDataset::ICreateLayer(const char *pszLayerName,
     }
     poLayer->SetCreationParameters(
         eGType, pszGeomColumnName, bGeomNullable, poSRS,
+        CSLFetchNameValue(papszOptions, "SRID"),
         poSrcGeomFieldDefn ? poSrcGeomFieldDefn->GetCoordinatePrecision()
                            : OGRGeomCoordinatePrecision(),
         CPLTestBool(
@@ -8561,7 +8622,7 @@ static void OGRGeoPackageImportFromEPSG(sqlite3_context *pContext, int /*argc*/,
         return;
     }
 
-    sqlite3_result_int(pContext, poDS->GetSrsId(oSRS));
+    sqlite3_result_int(pContext, poDS->GetSrsId(&oSRS));
 }
 
 /************************************************************************/
@@ -9052,7 +9113,7 @@ static void GPKG_ogr_layer_Extent(sqlite3_context *pContext, int /*argc*/,
     poRing->addPoint(sExtent.MinX, sExtent.MinY);
 
     const auto poSRS = poLayer->GetSpatialRef();
-    const int nSRID = poSRS ? poDS->GetSrsId(*poSRS) : 0;
+    const int nSRID = poDS->GetSrsId(poSRS);
     size_t nBLOBDestLen = 0;
     GByte *pabyDestBLOB =
         GPkgGeometryFromOGR(&oPoly, nSRID, nullptr, &nBLOBDestLen);
diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp
index e5e3043bad72..5ff549b492ac 100644
--- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp
+++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp
@@ -675,6 +675,9 @@ void RegisterOGRGeoPackage()
         "  <Option name='GEOMETRY_NULLABLE' type='boolean' "
         "description='Whether the values of the geometry column can be NULL' "
         "default='YES'/>"
+        "  <Option name='SRID' type='integer' description='Forced srs_id of "
+        "the "
+        "entry in the gpkg_spatial_ref_sys table to point to'/>"
         "  <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. "
diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp
index b8508415a663..5fc36cc05699 100644
--- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp
+++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp
@@ -1973,8 +1973,7 @@ OGRGeoPackageTableLayer::CreateGeomField(const OGRGeomFieldDefn *poGeomFieldIn,
     }
 
     const OGRSpatialReference *poSRS = oGeomField.GetSpatialRef();
-    if (poSRS != nullptr)
-        m_iSrs = m_poDS->GetSrsId(*poSRS);
+    m_iSrs = m_poDS->GetSrsId(poSRS);
 
     /* -------------------------------------------------------------------- */
     /*      Create the new field.                                           */
@@ -5605,10 +5604,10 @@ 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)
+    const OGRSpatialReference *poSRS, const char *pszSRID,
+    const OGRGeomCoordinatePrecision &oCoordPrec, bool bDiscardCoordLSB,
+    bool bUndoDiscardCoordLSBOnReading, const char *pszFIDColumnName,
+    const char *pszIdentifier, const char *pszDescription)
 {
     m_bIsSpatial = eGType != wkbNone;
     m_bIsInGpkgContents =
@@ -5626,9 +5625,59 @@ void OGRGeoPackageTableLayer::SetCreationParameters(
         m_nZFlag = wkbHasZ(eGType) ? 1 : 0;
         m_nMFlag = wkbHasM(eGType) ? 1 : 0;
         OGRGeomFieldDefn oGeomFieldDefn(pszGeomColumnName, eGType);
-        if (poSRS)
-            m_iSrs = m_poDS->GetSrsId(*poSRS);
+
         oGeomFieldDefn.SetSpatialRef(poSRS);
+        if (pszSRID)
+        {
+            m_iSrs = atoi(pszSRID);
+            if (m_iSrs == GDALGeoPackageDataset::FIRST_CUSTOM_SRSID - 1)
+            {
+                m_iSrs = m_poDS->GetSrsId(nullptr);
+                oGeomFieldDefn.SetSpatialRef(nullptr);
+            }
+            else
+            {
+                auto poGotSRS =
+                    m_poDS->GetSpatialRef(m_iSrs, /* bFallbackToEPSG = */ false,
+                                          /* bEmitErrorIfNotFound = */ false);
+                if (poGotSRS)
+                {
+                    oGeomFieldDefn.SetSpatialRef(poGotSRS);
+                    poGotSRS->Release();
+                }
+                else
+                {
+                    bool bOK = false;
+                    OGRSpatialReference *poSRSTmp = new OGRSpatialReference();
+                    if (m_iSrs < 32767)
+                    {
+                        CPLErrorHandlerPusher oErrorHandler(
+                            CPLQuietErrorHandler);
+                        CPLErrorStateBackuper oBackuper;
+                        if (poSRSTmp->importFromEPSG(m_iSrs) == OGRERR_NONE)
+                        {
+                            bOK = true;
+                            poSRSTmp->SetAxisMappingStrategy(
+                                OAMS_TRADITIONAL_GIS_ORDER);
+                            m_iSrs = m_poDS->GetSrsId(poSRSTmp);
+                            oGeomFieldDefn.SetSpatialRef(poSRSTmp);
+                        }
+                    }
+                    if (!bOK)
+                    {
+                        CPLError(
+                            CE_Warning, CPLE_AppDefined,
+                            "No entry in gpkg_spatial_ref_sys matching SRID=%s",
+                            pszSRID);
+                    }
+                    poSRSTmp->Release();
+                }
+            }
+        }
+        else
+        {
+            m_iSrs = m_poDS->GetSrsId(poSRS);
+        }
         oGeomFieldDefn.SetNullable(bGeomNullable);
         oGeomFieldDefn.SetCoordinatePrecision(oCoordPrec);
 
@@ -7271,8 +7320,7 @@ OGRErr OGRGeoPackageTableLayer::AlterGeomFieldDefn(
             if (m_poDS->SoftStartTransaction() != OGRERR_NONE)
                 return OGRERR_FAILURE;
 
-            const int nNewSRID =
-                poNewSRS ? m_poDS->GetSrsId(*(poNewSRS.get())) : 0;
+            const int nNewSRID = m_poDS->GetSrsId(poNewSRS.get());
 
             // Replace the old SRID by the new ones in geometry blobs
             uint32_t nNewSRID_LSB = nNewSRID;
diff --git a/ogr/ogrsf_frmts/sqlite/sqlite_rtree_bulk_load/sqlite_rtree_bulk_load.c b/ogr/ogrsf_frmts/sqlite/sqlite_rtree_bulk_load/sqlite_rtree_bulk_load.c
index 598d5047bd81..910399a304d1 100644
--- a/ogr/ogrsf_frmts/sqlite/sqlite_rtree_bulk_load/sqlite_rtree_bulk_load.c
+++ b/ogr/ogrsf_frmts/sqlite/sqlite_rtree_bulk_load/sqlite_rtree_bulk_load.c
@@ -813,6 +813,16 @@ static bool insert_into_db(const struct rtree_insert_context* ctxt,
     return true;
 }
 
+static bool IsLowercaseAlpha(const char* s)
+{
+    for (; *s != 0; ++s) {
+        if (!(*s >= 'a' && *s <= 'z')) {
+            return false;
+        }
+    }
+    return true;
+}
+
 bool SQLITE_RTREE_BL_SYMBOL(sqlite_rtree_bl_serialize)(
                                const struct sqlite_rtree_bl *tr,
                                sqlite3* hDB,
@@ -827,10 +837,24 @@ bool SQLITE_RTREE_BL_SYMBOL(sqlite_rtree_bl_serialize)(
         *p_error_msg = NULL;
     }
 
-    char* sql = sqlite3_mprintf(
-        "CREATE VIRTUAL TABLE \"%w\" USING rtree(\"%w\", \"%w\", \"%w\", \"%w\", \"%w\")",
-        rtree_name, rowid_colname, minx_colname, maxx_colname, miny_colname,
-        maxy_colname);
+    char* sql;
+    if (IsLowercaseAlpha(rowid_colname) &&
+        IsLowercaseAlpha(minx_colname) &&
+        IsLowercaseAlpha(maxx_colname) &&
+        IsLowercaseAlpha(miny_colname) &&
+        IsLowercaseAlpha(maxy_colname)) {
+        /* To make OGC GeoPackage compliance test happy... */
+        sql = sqlite3_mprintf(
+            "CREATE VIRTUAL TABLE \"%w\" USING rtree(%s, %s, %s, %s, %s)",
+            rtree_name, rowid_colname, minx_colname, maxx_colname, miny_colname,
+            maxy_colname);
+    }
+    else {
+        sql = sqlite3_mprintf(
+            "CREATE VIRTUAL TABLE \"%w\" USING rtree(\"%w\", \"%w\", \"%w\", \"%w\", \"%w\")",
+            rtree_name, rowid_colname, minx_colname, maxx_colname, miny_colname,
+            maxy_colname);
+    }
     int ret = sqlite3_exec(hDB, sql, NULL, NULL, p_error_msg);
     sqlite3_free(sql);
     if (ret != SQLITE_OK) {