Wrong interpretation of EPSG code from a PRJ file in QGIS-dev (proj v7.1.0)

[tobwen]

Describe the bug QGIS-dev with proj 7.1.0 currently interprets a PRJ-definition the wrong way. Instead of EPSG:3035 the string is detected as IGNF:ETRS89LAEA - ETRS89 Lambert Azimutal Equal Area - Projected

Normally, proj 7.1.0 shouldn't behave like this, since the confidence is >= 70% and there's only one matching SRS. Please check the additional context for this.

How to Reproduce

1. Download this demo dataset: test.zip
2. Load the dataset in QGIS dev
3. Check the coordinate system

QGIS and OS versions

QGIS version
3.13.0-Master
QGIS code revision
8ba73c6ff2
Compiled against Qt
5.11.2
Running against Qt
5.11.2
Compiled against GDAL/OGR
3.2.0dev
Running against GDAL/OGR
3.2.0dev
Compiled against GEOS
3.8.1-CAPI-1.13.3
Running against GEOS
3.8.1-CAPI-1.13.3
Compiled against SQLite
3.29.0
Running against SQLite
3.29.0
PostgreSQL Client Version
11.5
SpatiaLite Version
4.3.0
QWT Version
6.1.3
QScintilla2 Version
2.10.8
Compiled against PROJ
7.1.0
Running against PROJ
Rel. 7.1.0, August 1st, 2020
OS Version
Windows 10 (10.0)
This copy of QGIS writes debugging output.

Additional context According to the code, >= 70% should be enough: https://github.com/qgis/QGIS/blob/4ab9f17ee7b2b906300043c0ef9c3e1975bc4748/src/core/qgsprojutils.cpp#L235

projinfo --identify 'PROJCS["ETRS89_LAEA_Europe",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["false_easting",4321000.0],PARAMETER["false_northing",3210000.0],PARAMETER["central_meridian",10.0],PARAMETER["latitude_of_origin",52.0],UNIT["Meter",1.0]]'

Gives a confidence of 70% that this string is EPSG:3035 (and this is the only matching CRS).

[rouault]

I think I've the explanation. QGIS doesn't directly process the PJ* object built from the .prj file, but likely the one built from the WKT representation of the SRS exported by GDAL. When GDAL exports the WKT, it adds axis specification (since that's compulsory), but as the axis order infered from ESRI WKT is in easting, northing order and EPSG:3035 is northing, easting, later identification fails to find EPSG:3035 but finds instead IGNF:ETRS89LAEA that has easting, northing axis order

[tobwen]

Some background: This was the report of the alternatives used in the European Union (contains definitions by formulas etc on pages 25f.): https://esdac.jrc.ec.europa.eu/ESDB_Archive/ESDB_data_1k_raster_intro/annoni2005eurgrids.pdf

And that's the final INSPIRE Data Specification on Coordinate Reference Systems – Technical Guidelines from the working group, which links to EPSG:3035: https://inspire.ec.europa.eu/id/document/tg/rs https://inspire.ec.europa.eu/file/1726/download?token=3OGur2Ln http://www.opengis.net/def/crs/EPSG/0/3035

Comparisons

After ringing some doors, I found some guys with different software suites, which exported shapefiles in EPSG:3035 for me.

The current PRJ-string from the original dataset

PROJCS["ETRS89_LAEA_Europe",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["false_easting",4321000.0],PARAMETER["false_northing",3210000.0],PARAMETER["central_meridian",10.0],PARAMETER["latitude_of_origin",52.0],UNIT["Meter",1.0]] QGIS loads this as IGFN:ETRS89LAEA ❌

Export from ArcGIS Pro

PROJCS["ETRS89_LAEA_Europe",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["false_easting",4321000.0],PARAMETER["false_northing",3210000.0],PARAMETER["central_meridian",10.0],PARAMETER["latitude_of_origin",52.0],UNIT["Meter",1.0]] QGIS loads this as IGFN:ETRS89LAEA ❌

Export from GlobalMapper

PROJCS["ETRS89_LAEA_Europe",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["central_meridian",10],PARAMETER["latitude_of_origin",52],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["Meter",1]] QGIS loads this as IGFN:ETRS89LAEA ❌

Export from ogr2ogr 3.0.4/3.2.0-dev

PROJCS["ETRS_1989_LAEA",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",4321000.0],PARAMETER["False_Northing",3210000.0],PARAMETER["Central_Meridian",10.0],PARAMETER["Latitude_Of_Origin",52.0],UNIT["Meter",1.0]] QGIS loads this als EPSG:3035 ✔️

Export from QGIS 3.13-Master

PROJCS["ETRS_1989_LAEA",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",4321000.0],PARAMETER["False_Northing",3210000.0],PARAMETER["Central_Meridian",10.0],PARAMETER["Latitude_Of_Origin",52.0],UNIT["Meter",1.0]] QGIS loads this als EPSG:3035 ✔️

Conclusion

The ESRI definitions seem to make trouble here. While shapefile is a product by ESRI, I wonder if QGIS shouldn't take care of the characteristics of the PRJ-files. Since ArcGIS is scriptable, one could create a lookup table to match between numeric EPSG code and the output in shapefile to get a perfect match.

[rouault]

My above analysis was wrong. This was a PROJ pure issue. Fixed per https://github.com/OSGeo/PROJ/pull/2240