OGR reports incorrect SRS for GeoJSON

Question

OGR (GDAL 3.2.2) appears to incorrectly read the SRS of GeoJSON features and I would like to know if there is anything I can do to prevent this.

I have two GeoJSON files and both represent the same feature.

File 1 includes an explicit CRS property:

"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }

File 2 does not include an explicit CRS property, which I believe means CRS84 is implied (RFC-7946).

Both files are opened by OGR's GeoJSON driver, loaded into a layer, and have their SRS inspected in the following manner.

driver = GetDriverByName("GeoJSON")
datasource = driver.Open("...geojson")
layer = datasource.GetLayerByIndex(0)
print(layer.GetSpatialRef().ExportToWkt())

Both files report the following WKT (it ends with AUTHORITY["EPSG","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"]]'

However, if I explicitly create a CRS84 SRS via OGR (using exactly the same text as File 1's CRS property), it reports a different WKT.

crs84_srs: Final = SpatialReference()
crs84_srs.SetFromUserInput("urn:ogc:def:crs:OGC:1.3:CRS84")
print(crs84_srs.ExportToWkt())

This reports the following (does not end with EPSG reference)

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

Question: Why does OGR treat GeoJSON as EPSG:4326 instead of CRS84, and is there anything I can do to force it to respect the explicit or implicit presence of CRS84 in the input file? I have reviewed the driver options but I see no likely candidates. The difference between EPSG:4326 and CRS84 is subtle but important - for example is affects axis ordering when calling layer.GetExtent().

Answer 1

CRS84 is the same as EPSG:4326 but with reversed axis order, so lon-lat vs. lat-lon, respectively. But GeoJSON by definition has coordinates always in lon-lat order and therefore CRS84 = EPSG:4326. Also the default coordinate system in GeoJSON is WGS84 lot-lat coordinates. All three alternatives CRS84, EPSG:4326, and no crs at all mean the same thing.

Probably the best way to write interoperable GeoJSON it to use creation option RFC7946=YES https://gdal.org/drivers/vector/geojson.html. Then the result will have no crs and any software that reads GeoJSON knows what it means.

Ogrinfo reports the CRS of GeoJSON that is created according to RFC7946 this way:

ogrinfo crs.json -al
INFO: Open of `crs.json'
      using driver `GeoJSON' successful.

Layer name: SELECT
Geometry: Polygon
Feature Count: 1
Extent: (-34.786396, 0.417717) - (-31.168008, 2.388746)
Layer SRS WKT:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1

So GDAL seems to know how to handle the situation. It reports that the CRS is EPSG:4326 but in the data the coordinates are in longitude-latitude order. The result is the same than having CRS84 without re-mapped axis.

I have not checked from the source code what function ogrinfo is using when it gets the result Data axis to CRS axis mapping: 2,1 and if the info comes from the GeoJSON driver or some other place.