GDA2020 identification failure

[nyalldawson]

A shapefile identified by ArcGIS as having GDA2020 CRS (EPSG:7844) has the .prj contents:

GEOGCS["GDA2020",DATUM["D_GDA2020",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]

Plugging this into projinfo --identify does not identify this as EPSG:7844:

PROJ.4 string:
+proj=longlat +ellps=GRS80 +no_defs +type=crs

WKT2:2019 string:
GEOGCRS["GDA2020",
    DATUM["D_GDA2020",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1],
            ID["EPSG",7019]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["Degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["Degree",0.0174532925199433]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["Degree",0.0174532925199433]]]

Identification match count: 30
IGNF:ETRS89G: 60 %
IGNF:RGWF96GDD: 60 %
IGNF:RGWF96G: 60 %
EPSG:7086: 60 %
IGNF:RGAF09GDD: 60 %
IGNF:RGAF09G: 60 %
EPSG:7084: 60 %
IGNF:RGF93GDD: 60 %
IGNF:RGF93G: 60 %
EPSG:7041: 60 %
IGNF:RGFG95GDD: 60 %
IGNF:RGFG95G: 60 %
EPSG:7039: 60 %
IGNF:RGM04GDD: 60 %
IGNF:RGM04G: 60 %
IGNF:RGPFGDD: 60 %
IGNF:RGPFG: 60 %
EPSG:7037: 60 %
IGNF:RGR92GDD: 60 %
IGNF:RGR92G: 60 %
EPSG:7035: 60 %
IGNF:RGSPM06GDD: 60 %
IGNF:RGSPM06G: 60 %
EPSG:7133: 60 %
IGNF:RGTAAF07GDD: 60 %
IGNF:RGTAAF07G: 60 %
EPSG:8902: 60 %
EPSG:7844: 25 %
EPSG:7842: 25 %
EPSG:7843: 25 %

Seems the different vs projinfo EPSG:7844 -o WKT1:ESRI is:

GDA2020 (prj) vs GCS_GDA2020 (projinfo)
D_GDA2020 (prj) vs GDA2020 ) (projinfo)

Indeed, changing the prj WKT string to

GEOGCS["GDA2020",DATUM["GDA2020",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]

correctly gives a EPSG:7844 100% match, so the GCS_GDA2020 difference seems irrelevant.

This looks similar to #1829, in that the spheroid name is prefixed with "D_". Some tool out there is creating shapefiles with WKT using these prefixes, so I think proj should be tolerant to this and add D_ variants to its list of known aliases.

[nyalldawson]

And now I find https://gdal.org/tutorials/wktproblems.html#datum-names -- so it's evidently a known difference with ESRI WKT handling. Perhaps the bug here is that proj's identify isn't correctly handling ESRI WKT 1 variants?

[rouault]

so it's evidently a known difference with ESRI WKT handling

WKT handling has been rewritten from scratch in PROJ compared to GDAL 2, so no wonder there are differences :-)

Perhaps the bug here is that proj's identify isn't correctly handling ESRI WKT 1 variants?

We could argue the bug is in ESRI database itself. Looking at their official database of geographic CRS:

$ grep -r GDA ~/projection-engine-db-doc/csv/pe_list_geogcs.csv 
4283,4283,"PE_GCS_GDA_1994","GCS_GDA_1994","GEOGCS[""GCS_GDA_1994"",DATUM[""D_GDA_1994"",SPHEROID[""GRS_1980"",6378137.0,298.257222101]],PRIMEM[""Greenwich"",0.0],UNIT[""Degree"",0.0174532925199433]]","Geocentric Datum of Australia 1994","EPSG","9.8(3.0.1)","no","Australia - GDA",-60.56,-8.47,93.41,173.35
7844,7844,"PE_GCS_GDA2020","GDA2020","GEOGCS[""GDA2020"",DATUM[""GDA2020"",SPHEROID[""GRS_1980"",6378137.0,298.257222101]],PRIMEM[""Greenwich"",0.0],UNIT[""Degree"",0.0174532925199433]]","GDA2020","EPSG","9.8(10.5.0)","no","Australia - GDA",-60.56,-8.47,93.41,173.35

it seems the official name for datum GDA94 is D_GDA_1994, but for GDAL2020 it is GDA2020 without leading D_

Confirmed also by

$ grep -r GDA ~/projection-engine-db-doc/csv/pe_list_datum.csv 
1168,1168,"PE_D_GDA2020","GDA2020","DATUM[""GDA2020"",SPHEROID[""GRS_1980"",6378137.0,298.257222101]]","GDA2020","EPSG","8.9.4(10.5.0)","no"
6283,6283,"PE_D_GDA_1994","D_GDA_1994","DATUM[""D_GDA_1994"",SPHEROID[""GRS_1980"",6378137.0,298.257222101]]","Geocentric Datum of Australia 1994","EPSG","8.3.4(3.0.1)","no"

Or perhaps they have changing datum names in different ArcGIS versions, and the DB fails to capture old names...
Workaround to follow

[nyalldawson]

WKT handling has been rewritten from scratch in PROJ compared to GDAL 2, so no wonder there are differences :-)

Right -- the emphasis there was on "a known difference with ESRI WKT handling" - I was initially unsure if these files were being created by some other tool. I wasn't pointing the blame at PROJ! 😉

Looking at https://github.com/OSGeo/gdal/blob/d5f6bb89b0e0db0a06489419050d432e8f58ff47/gdal/data/gdal_datum.csv - there was quite a comprehensive db of known ESRI datum names. Perhaps those should be included in the proj alias tables too?

[rouault]

Perhaps those should be included in the proj alias tables too?

I'd expect most of those to be hopefully already included in the alias table, but there are perhaps holes. Good exercise for a willing contributor :-)

[nyalldawson]

I'm willing to give that a shot -- is there anyway from public API to query the alias table or obtain a datum from name?

[rouault]

obtain a datum from name?

==> proj_create_from_name()
If I remember well it may use the alias table in master. But I'd suggest to play with the database directly. You may use ogr2ogr on the gdal_datum.csv file to ingest it in proj.db, and use sqlite to do joins. Or something like that

[nyalldawson]

Ok, there's only three entries in the gdal 2 file which don't have corresponding alias entries:

DATUM_CODE	DATUM_NAME	DATUM_TYPE	ORIGIN_DESCRIPTION	REALIZATION_EPOCH	ELLIPSOID_CODE	PRIME_MERIDIAN_CODE	AREA_OF_USE_CODE	DATUM_SCOPE	REMARKS	INFORMATION_SOURCE	DATA_SOURCE	REVISION_DATE	CHANGE_ID	DEPRECATED	ESRI_DATUM_NAME
5101	Ordnance Datum Newlyn	vertical	Mean Sea Level at Newlyn between 1915 and 1921. Initially realised through 1921 and then 1956 levelling network adjustments. From 2002 redefined to be realised through OSGM geoid models.				2792	Topographic mapping, geodetic survey.	Orthometric heights.	Ordnance Survey of Great Britain.	IOGP	2016/07/13	2004.100 2015.069	0	D_Ordnance_Datum_Newlyn
5100	Mean Sea Level	vertical	The average height of the surface of the sea at a tide station for all stages of the tide over a 19-year period, usually determined from hourly height readings measured from a fixed predetermined reference level.				1262	Hydrography.	Approximates geoid. Users are advised to not use this generic vertical datum but to define explicit realizations of MSL by specifying location and epoch, for instance "MSL at xxx during yyyy-yyyy".	IHO Dictionary, S-32, 5th Edition, 3156.	IOGP	2017/07/14	2011.047 2017.029	0	D_Mean_Sea_Level
6801	CH1903 (Bern)	geodetic	Fundamental point: Old Bern observatory. Latitude: 46°57'08.660"N, longitude: 0°E (of Bern).	1903-01-01	7004	8907	1286	Topographic mapping.		Bundesamt für Landestopographie	OGP	2008/06/24	2003.361 2008.045	0	D_Bern_1898

The first two don't have any corresponding entries in geodetic_datum, so that leaves just the missing "D_Bern_1898" alias for EPSG:6801.

[lgouzin]

is this a good place to say that in QGIS, 3.32, shp layer in 7844 gets loaded weirdly, even ones saved by QGIS itself. It projects the data at the right place, but the epsg is not given when hovering a layer in the browser. In Source, the CRS shows only has GCS_GDA2020. It does bring problems in pyqgis when using layer.crs().authid() which is returned has empty, only layer.crs().description() returns GCS_GDA2020

[jjimenezshaw]

@lgouzin is this the CRS defintion you have?

projinfo EPSG:7844 -o wkt1_esri -q
GEOGCS["GCS_GDA2020",DATUM["GDA2020",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

The WKT string in the unit test added in #1915 for GDA2020 is different from the string produced directly by PROJ in ESRI format. I don't know if it was in purpose or not.

--identify on the string shown above (produced by PROJ) returns

Identification match count: 1
EPSG:7844: 70 %