From 7205e50a961ecef62b4562247af1683e6e07888f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 2 Apr 2024 15:42:58 +0200 Subject: [PATCH 1/2] sqlite_rtree_bulk_load.c: generate a CREATE VIRTUAL TABLE ... USING rtree statement that exactly matches the OGC GeoPackage conformance test suite --- .../sqlite_rtree_bulk_load.c | 32 ++++++++++++++++--- 1 file changed, 28 insertions(+), 4 deletions(-) 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) { From c104b3146e6cf0abe90854608d7c1839080f5cf5 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 2 Apr 2024 15:51:47 +0200 Subject: [PATCH 2/2] GPKG: map Null CRS to a new srs_id=99999, and add a SRID layer creation option New layer creation option: ``` - .. lco:: SRID :choices: :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``: 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 nature". ``` --- autotest/gdrivers/gpkg.py | 3 +- autotest/ogr/ogr_gpkg.py | 134 +++++++++++++++- autotest/ogr/ogr_mssqlspatial.py | 9 +- doc/source/drivers/vector/gpkg.rst | 33 +++- ogr/ogrsf_frmts/gpkg/ogr_geopackage.h | 18 ++- .../gpkg/ogrgeopackagedatasource.cpp | 143 +++++++++++++----- ogr/ogrsf_frmts/gpkg/ogrgeopackagedriver.cpp | 3 + .../gpkg/ogrgeopackagetablelayer.cpp | 68 +++++++-- 8 files changed, 342 insertions(+), 69 deletions(-) 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: + :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::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 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 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() "