GeoPackage: Fix handling of invalid SRS ID when writing

[mloskot]

This issue is dedicated to discuss further improvements to PR #3286 suggested by @rouault in #3286 (comment), quoting:

It occurs to me that we would probably want to do something on the writing side of the driver to be able to round-trip nicely those 2 special CRS.

---

@rouault I admit I'm having hard time trying to find an anchor for the 'writing side' in the GeoPackage driver - I have very little experience with the format.

Do you mean the driver is lacking as not accepting any of the two possible undefined SRS when new dataset/layer is created?

Putting it differently, what new test cases are expected to be covered (and any required implementation on the drivier added/changed, of course). Here is what I have just been scratching out:

def test_ogr_gpkg_write_srs_undefined():

    ds = gdaltest.gpkg_dr.CreateDataSource('/vsimem/ogr_gpkg_srs_undefined.gpkg')
    assert ds is not None

    srs= osr.SpatialReference()
    srs.SetFromUserInput('GEOGCS["Undefined geographic SRS",DATUM["unknown",SPHEROID["unknown",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]')
    lyr = ds.CreateLayer('srs_test_layer', geom_type=ogr.wkbPoint, srs=srs)

    assert lyr.GetSpatialRef().ExportToWkt().find('Undefined geographic SRS') >= 0

def test_gpkg_write_srs_undefined():

    ds = gdaltest.gpkg_dr.Create('/vsimem/tmp.gpkg', 1, 1)
    ds.SetGeoTransform([0, 1, 0, 0, 0, -1])

    srs = osr.SpatialReference()
    srs.ImportFromWkt('GEOGCS["Undefined geographic SRS",DATUM["unknown",SPHEROID["unknown",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]')
    ret = ds.SetSpatialRef(srs)

    srs_got = ds.GetSpatialRef()
    assert srs_got.IsSame(srs), srs_got.Export

I'll appreciate any hints or bearings where the driver is lacking.

[rouault]

I'll appreciate any hints or bearings where the driver is lacking.

I've not checked but I suspect that if you run those cases, instead of using SRID = 0 or -1, it will add a new entry in the gpkg_spatial_ref_sys table. This is not a huge deal, but it would be more in the spirit of the spec to just use the special SRID values.

[mloskot]

Yes, new SRS is inserted, because GDALGeoPackageDataset::GetSrsId does not match the undefined definition of the srid_id = -1 to the equivalent SRS constructed from the WKT:

srs = osr.SpatialReference()
srs.ImportFromWkt('GEOGCS["Undefined geographic SRS",DATUM["unknown",SPHEROID["unknown",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]')

No authority is found, so it falls-back to definition = <WKT> match which fails and the GetSrsId ends up adding new record with srid_id = 10000 (see Python test sample below).

Undefined geographic SRS | 100000 | 100000 | GEOGCS["Undefined geographic SRS",DATUM["unknown",SPHEROID["unknown",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]] | None

Problem

I don't really know how to change the GetSrsId logic to match the 'GEOGCS["Undefined geographic SRS", ...] or the LOCALCS["Undefined cartesian SRS", ...] to their 'degenerate' definition-less equivalents from the gpkg_spatial_ref_sys table:

Python Test

Here is my current Python sample that exercises the situation:

def dump_gpkg_spatial_ref_sys(ds, tag):
    print('--- {0} ------------'.format(tag))
    lyr = ds.ExecuteSQL('SELECT * FROM gpkg_spatial_ref_sys ORDER BY srs_id')
    f = lyr.GetNextFeature()
    i = 1
    while f:
        print('{0}:'.format(i), f.GetField('srs_name'), '|', f.GetField('srs_id'), '|', f.GetField('organization_coordsys_id'), '|',  f.GetField('definition'), '|', f.GetField('description'))
        f = lyr.GetNextFeature()
        i = i + 1

def test_gpkg_write_srs_undefined():

    if gdaltest.gpkg_dr is None:
        pytest.skip()

    gdal.Unlink('/vsimem/tmp.gpkg')

    ds = gdaltest.gpkg_dr.Create('/vsimem/tmp.gpkg', 1, 1)
    ds.SetGeoTransform([0, 1, 0, 0, 0, -1])

    dump_gpkg_spatial_ref_sys(ds, 'BEFORE')

    srs = osr.SpatialReference()
    srs.ImportFromWkt('GEOGCS["Undefined geographic SRS",DATUM["unknown",SPHEROID["unknown",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]')
    ret = ds.SetSpatialRef(srs)

    dump_gpkg_spatial_ref_sys(ds, 'AFTER ')

Output:

--- BEFORE ------------
1: Undefined cartesian SRS | -1 | -1 | undefined | undefined cartesian coordinate reference system
2: Undefined geographic SRS | 0 | 0 | undefined | undefined geographic coordinate reference system
3: WGS 84 geodetic | 4326 | 4326 | GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]] | longitude/latitude coordinates in decimal degrees on the WGS 84 spheroid
--- AFTER  ------------
1: Undefined cartesian SRS | -1 | -1 | undefined | undefined cartesian coordinate reference system
2: Undefined geographic SRS | 0 | 0 | undefined | undefined geographic coordinate reference system
3: WGS 84 geodetic | 4326 | 4326 | GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]] | longitude/latitude coordinates in decimal degrees on the WGS 84 spheroid
4: Undefined geographic SRS | 100000 | 100000 | GEOGCS["Undefined geographic SRS",DATUM["unknown",SPHEROID["unknown",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]] | None

[rouault]

I don't really know how to change the GetSrsId logic to match the 'GEOGCS["Undefined geographic SRS", ...] or the LOCALCS["Undefined cartesian SRS", ...]

I would likely add an early test in the method for the nature of the CRS with OGRSpatialReference::IsGeographic() / OGRSpatialReference::IsLocal() and on its name with OGRSpatialReference::GetName()

[mloskot]

Okey.

Seeing the any in the srs_name field of the table, I wasn't certain about relying on the OGRSpatialReference::GetName.

This does seem to do the job:

    if (poSRS->IsGeographic() || poSRS->IsLocal())
    {
        const char* pszName = poSRS->GetName();
        if ( pszName != nullptr && strlen(pszName) > 0 )
        {
            if (EQUAL(pszName, "Undefined cartesian SRS"))
                return -1;

            if (EQUAL(pszName, "Undefined geographic SRS"))
                return 0;
        }
    }

I'll prepare PR.

[rouault]

Seeing the any in the srs_name field of the table, I wasn't certain about relying on the OGRSpatialReference::GetName.

that's a good point, but without testing the name, we have no way of deciding since we have no concept of SRID in OGRSpatialReference. My idea was that a ogr2ogr of a source GeoPackage with SRID=0 or -1 to a target GeoPackage would result in 0 and -1 being used in the target dataset.