Incompatible projection reporting of EPSG:27700 with GDAL 3+

[mattenvsys]

Raster files known to be projected in EPSG:27700 with GDAL <= 2.4 are reported differently in recent versions of gdalinfo. This causes them to be interpreted as "Unknown CRS" in QGIS 3.10.3 (running against PROJ 6 and GDAL 3) and become incorrectly positioned.

Unsure of where to post this issue. Seems to affect projection of raster/vector files, but unsure if it stems from GDAL/PROJ/other. QGIS devs suggested to post it here. We're not building from source so are unable to test different combinations of GDAL and PROJ versions to pinpoint the issue further.

The following projection info on the same file (attached here 27700.zip):

1. GDAL 2.4.3 / PROJ 5.2.0

PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646,
                AUTHORITY["EPSG","7001"]],
            TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4277"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","27700"]]

2. GDAL 3.0.4 / PROJ 6.3.1
   (note the additional sections of BOUNDCRS, SOURCECRS, TARGETCRS and ABRIDGEDTRANSFORMATION)

BOUNDCRS[
    SOURCECRS[
        PROJCRS["OSGB 1936 / British National Grid",
            BASEGEOGCRS["OSGB 1936",
                DATUM["OSGB 1936",
                    ELLIPSOID["Airy 1830",6377563.396,299.3249646,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4277]],
            CONVERSION["British National Grid",
                METHOD["Transverse Mercator",
                    ID["EPSG",9807]],
                PARAMETER["Latitude of natural origin",49,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8801]],
                PARAMETER["Longitude of natural origin",-2,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["Scale factor at natural origin",0.9996012717,
                    SCALEUNIT["unity",1],
                    ID["EPSG",8805]],
                PARAMETER["False easting",400000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",-100000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8807]]],
            CS[Cartesian,2],
                AXIS["(E)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1]],
                AXIS["(N)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1]],
            USAGE[
                SCOPE["unknown"],
                AREA["UK - Britain and UKCS 49┬░46'N to 61┬░01'N, 7┬░33'W to 3┬░33'E"],
                BBOX[49.75,-9.2,61.14,2.88]],
            ID["EPSG",27700]]],
    TARGETCRS[
        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]],
            USAGE[
                SCOPE["unknown"],
                AREA["World"],
                BBOX[-90,-180,90,180]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",446.448,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",-125.157,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",542.06,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0.15,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0.247,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0.842,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",0.999979511,
            ID["EPSG",8611]]]]

The latter does not interpret as EPSG:27700 in QGIS 3.10.3, but rather as "Unknown CRS". Its position is offset from other files that have projections read correctly.

How To Reproduce

1. Download test file 27700.zip
2. Load into a new instance of QGIS 3.10.3 on Windows
3. The project and file projection should be labelled as "Unknown CRS"

Expected Behaviour

The projection in QGIS should be recognised as "EPSG:27700"

Operating System

Windows 10 64bit build 1903

GDAL Version

3.0.4 from QGIS 3.10.3 64bit on Windows

[rouault]

As both a GDAL and QGIS dev :-), I don't think this is a GDAL bug. See my comment in qgis/QGIS#34993 (comment)

[rouault]

As mentionned in qgis/QGIS#34993 (comment), "I've just written the https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/samples/gdal_remove_towgs84.py script. You can download it and run it with GDAL 2 or GDAL 3 to remove, in-place, TOWGS84[] clauses from GeoTIFF file whose CRS contains both a EPSG code and a TOWGS84[] clause. Fixing such files is the way forward."

Considering this as the fix for this issue