Merge pull request #9595 from rouault/fix_9580

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

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("")

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)

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

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
 -----------------------------------------

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)