GPKG: map Null CRS to a new srs_id=99999, and add a SRID layer creati…

…on option

New layer creation option:
```
-  .. 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.
```

```
Coordinate Reference Systems
----------------------------

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``: 0
- ``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``: ``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 nature".
```

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"] == 0
+        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"] == "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)

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``: 0
+- ``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``: ``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 nature".
 
 Level of support of GeoPackage Extensions
 -----------------------------------------

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)