GeoJSON does not write urn:ogc:def:crs:OGC::CRS84

[snowman2]

Expected behavior and actual behavior.

With the GeoJSON driver and GDAL 3.0.2, when writing EPSG:4326 to the file, it writes urn:ogc:def:crs:OGC::CRS84 and returns urn:ogc:def:crs:OGC::CRS84. However, when you write urn:ogc:def:crs:OGC::CRS84 to the file, it does not write the CRS and returns EPSG:4326.

It would be nice if there was consistency such that the roundtrip writing/reading of the CRS preserved the original CRS. But, if there is a reason for the behavior, that would be good to know as well. Thanks!

Steps to reproduce the problem.

Write EPSG:4326:

>>> outSpatialRef = osr.SpatialReference()
>>> outSpatialRef.SetFromUserInput("EPSG:4326")
0
>>> outDataSource = outDriver.CreateDataSource('test.geojson')
>>> lyr = outDataSource.CreateLayer('foo', srs=outSpatialRef)
>>> lyr = None
>>> outDataSource = None
>>> with open("test.geojson") as tj:
...     data = json.load(tj)
... 
>>> list(data)
['type', 'name', 'crs', 'features']
>>> data['crs']
{'type': 'name', 'properties': {'name': 'urn:ogc:def:crs:OGC:1.3:CRS84'}}
>>> ds = ogr.Open("test.geojson")
>>> lyr = ds.GetLayer(0)
>>> sr_got = lyr.GetSpatialRef()
>>> sr_got.ExportToWkt()
'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["Longitude",EAST],AXIS["Latitude",NORTH]]'
>>> outSpatialRef.ExportToWkt()
'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"]]'

Write urn:ogc:def:crs:OGC::CRS84:

>>> from osgeo import ogr, osr
>>> outSpatialRef = osr.SpatialReference()
>>> outSpatialRef.SetFromUserInput("urn:ogc:def:crs:OGC::CRS84")
0
>>> outDriver = ogr.GetDriverByName('GeoJSON')
>>> outDataSource = outDriver.CreateDataSource('test.geojson')
>>> lyr = outDataSource.CreateLayer('foo', srs=outSpatialRef)
>>> lyr = None
>>> outDataSource = None
>>> ds = ogr.Open("test.geojson")
>>> lyr = ds.GetLayer(0)
>>> sr_got = lyr.GetSpatialRef()
>>> sr_got.ExportToWkt()
'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"]]'
>>> outSpatialRef.ExportToWkt()
'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["Longitude",EAST],AXIS["Latitude",NORTH]]'
>>> import json
>>> with open("test.geojson") as tj:
...     data = json.load(tj)
... 
>>> list(data)
['type', 'name', 'features']

Operating system

Ubuntu 18.04.3 LTS

Linux snowal-lx2 4.15.0-66-generic #75-Ubuntu SMP Tue Oct 1 05:24:09 UTC 2019 x86_64 x86_64 x86_64 GNU/Linux

GDAL version and provenance

# Name                    Version                   Build  Channel
gdal                      3.0.2            py38hbb6b9fb_3    conda-forge
libgdal                   3.0.2                h5534617_3    conda-forge

[rouault]

I see your point(s), but I'm not sure we want to change things in that area. Format conversion is not intended to be always round-trippable, and the last GeoJSON version RFC7946 (available as a creation option) removes any explicit support for CRS in GeoJSON.
What could possibly be improved is that if you set urn:ogc:def:crs:OGC::CRS84 as the input CRS, urn:ogc:def:crs:OGC::CRS84 is also written. But using urn:ogc:def:crs:OGC::CRS84 with OGRSpatialReference a bit unusual on the GDAL side. On reading of GeoJSON files, returning the equivalent EPSG code is the intended behaviour

[jorisvandenbossche]

On reading of GeoJSON files, returning the equivalent EPSG code is the intended behaviour

But that only happens (returning the EPSG:4326) if there is no CRS information in the geojson file.
If the file actually contains a 'crs' field with {'type': 'name', 'properties': {'name': 'urn:ogc:def:crs:OGC:1.3:CRS84'}} content, GDAL does not return the "equivalent EPSG code", but it actually returns CRS84.
Now, this feels correct to return what is written in the file, but then GDAL should maybe not write "urn:ogc:def:crs:OGC:1.3:CRS84" crs in the geojson file if the input CRS was EPSG:4326? (and only do that if the input crs was actually CRS84)

[rouault]

But that only happens (returning the EPSG:4326) if there is no CRS information in the geojson file.
If the file actually contains a 'crs' field with {'type': 'name', 'properties': {'name': 'urn:ogc:def:crs:OGC:1.3:CRS84'}} content, GDAL does not return the "equivalent EPSG code", but it actually returns CRS84.

True indeed. Missed that. Comes from the fact that GDAL has a historic hardcoded WKT representation for urn:ogc:def:crs:OGC:1.3:CRS84 that is the same as EPSG:4326 but doesn't include the EPSG:4326 authority and use longitude, latitude axis order (on purpose). That said, this is historic behaviour, and with PROJ 6 we have now a OGC:CRS84 entry in proj.db, so if we removed that hardcoded string from SetFromUserInput() implementation, we'd got a WKT representation with a OGC:CRS84 authority. Which wouldn't probably be ideal for the sake of doing later coordinate transformations, as mixing CRS from different authorities will lead to suboptimal datum transformation propositions.

then GDAL should maybe not write "urn:ogc:def:crs:OGC:1.3:CRS84" crs in the geojson file if the input CRS was EPSG:4326?

Mostly comes from the fact that the old GeoJSON 2008 spec recommended using urn:ogc:def:crs:OGC:1.3:CRS84 over EPSG:4326: https://geojson.org/geojson-spec#named-crs

Writing EPSG:4326 in the GeoJSN file could lead to ambiguities about the ordering of coordinates in the GeoJSON file.
On the reading side, when there's no crs info, OGR returns though EPSG:4326 (with the Data axis to CRS axis mapping: 2,1 info now since GDAL 3 to indicate that the data axis order is the reverse of the one mandated by the authority !), since EPSG codes are more easily interoperable for reasons given above.

All in all, a mix of format specificities, historic behaviour, data axis order mess... As there's no major functional disorder due to current behaviour (like wrong axis order), the benefit vs risk ratio of changing current behaviour is a bit low :-)

What could perhaps be done is on the writing side to write "urn:ogc:def:crs:OGC::CRS84" in the GeoJSON file if the input CRS was "urn:ogc:def:crs:OGC::CRS84", but I'd suggest rather not to use "urn:ogc:def:crs:OGC::CRS84" for CRS passed to CreateLayer(). With other formats, it could cause issues that you wouldn't have if you would use plain EPSG:4326

[jorisvandenbossche]

Thanks for the explanation! I fully agree that given the current (historic) "mess", changing something might only increase the mess ;)

But, from an application standpoint, like GeoPandas, I was thinking if we should automatically convert CRS84 into EPSG:4326 on reading the file (if the file source is geojson).
One reason is to preserve historic behaviour in GeoPandas (or actually in Fiona), where with GDAL < 3 / PROJ < 6, you get {'init': 'epsg:4326'} as the CRS information when reading a geojson file. So if we don't do any automatic conversion, users will start getting CRS84 instead of EPSG:4326 when reading a geojson file with the upcoming geopandas version (which will start properly using the new features of PROJ 6 / pyproj > 2).
Another reason is that I think getting back EPSG:4326 is also less confusing, as for better or worse, that's what users will expect to see to know they have geographic (lon/lat) data.
But is this something you would advise against?

[rouault]

I was thinking if we should automatically convert CRS84 into EPSG:4326 on reading the file (if the file source is geojson).

That does make sense. Actually I missed there was a difference in behaviour between GDAL 3 and GDAL 2 when the GeoJSON file had a crs: urn:ogc:def:crs:OGC:1.3:CRS84. GDAL 2 used to report EPSG:4326, whereas currently GDAL does not. GDAL 2 behaviour seems preferable, so I've just fixed that per cd9e06c and backported it to 3.0 branch. So with that fix, a workaround on GeoPandas/Fiona wouldn't be needed

Note: my previous explanation about GDAL 2 not including EPSG:4326 authority in its urn:ogc:def:crs:OGC:1.3:CRS84 hardcoded WKT was wrong ! It actually did included it (but omitted explicit axis order, which meant implicitly long, lat). This is something I've changed in GDAL 3 (not to lie on the authority and adding explicit long, lat order), but this had this side effect on the GeoJSON driver. Pfff :-)

[snowman2]

Thanks @rouault!

[jorisvandenbossche]

Thanks a lot on the quick follow-up @rouault !