Unexpected consequences of AutoIdentifyEPSG()

[AndrewAtAvenza]

I read the note at the top of the submission blurb and while this is somewhat axis-adjacent, but not really an axis issue even though the consequences present like an axis bug.

We made the switch to 3.3 but had the usual axis-related issues. Using OAMS_TRADITIONAL_GIS_ORDER solved most of them, but we had some odd outliers. After a great deal of digging, I managed to boil it down to a fairly simple code block that demonstrates the issue. I've included comments that indicate what's happening under the hood along the way

#include <ogr_spatialref.h>
#include <string>

int main()
{
    const char* kGda94Wkt = R"(PROJCS["GDA_1994_MGA_Zone_55",GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",10000000.0],PARAMETER["Central_Meridian",147.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]])";

    OGRSpatialReference gda94;
    gda94.importFromWkt(kGda94Wkt);
    gda94.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    const auto autoIdentifyResult = gda94.AutoIdentifyEPSG();

    char* wkt = nullptr;
    gda94.exportToWkt(&wkt);
    printf("Notice that it has added the EPSG code for the GEOSCS, even though it says it failed (%d):\n'%s'\n", autoIdentifyResult, wkt);
    CPLFree(wkt);

    // when it clones the GEOCS, it fills in the AXIS tags, but does not notice that they
    // contradict the adjacent AUTHORITY tag for EPSG:4283. This discrepency causes all sorts
    // of headaches to follow
    auto* gda94GeodeticBase = gda94.CloneGeogCS();
    assert(gda94GeodeticBase);
    gda94GeodeticBase->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

    gda94GeodeticBase->exportToWkt(&wkt);
    printf("Notice that the WKT includes EPSG 4283, which assumes lat/long, but the AXIS indicates that it's long/lat:\n'%s'", wkt);
    CPLFree(wkt);

    auto* wgs84 = OGRSpatialReference::GetWGS84SRS();
    // I recognize this is redundant, but for consistency
    wgs84->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

    // when this is constructed, it attempts to see if it can feed an EPSG code to PROJ instead
    // of the WKT; when it checks the GDA94 geodetic system, it will use EPSG 4283, and when
    // it compares them, it will think it's fine because of this line:
    //
    // oTmpSRS.SetDataAxisToSRSAxisMapping(poSRS->GetDataAxisToSRSAxisMapping());
    //
    // This effectively nullifies the axis test portion of IsSame(). I suspect it should be:
    //
    // oTmpSRS.SetAxisMappingStrategy(poSRS->SetAxisMappingStrategy());
    //
    // This would cause it to refresh the axis mapping internally and when it compares them,
    // oTmpSRS will be { 1, 2 } while poSRS will be { 2, 1 } causing it to fall back on the WKT
    // in which PROJ ignores the EPSG tag and properly handles the transformation.
    auto* wgs84ToGda94GeodeticBase = OGRCreateCoordinateTransformation(wgs84, gda94GeodeticBase);

    // rounded to whole number for testing clarity
    double x = -79;
    double y = 43;

    // Because it used the EPSG code instead of the WKT, PROJ makes this a no-op, but
    // GDAL thinks there's an axis swap because WGS84 is { 2, 1 } and the GDA94 geodetic is
    // { 1, 2 }. If PROJ reads the WKT instead, it will add an axis swap, which is then
    // properly negated by the mapping swap
    if (TRUE == wgs84ToGda94GeodeticBase->Transform(1, &x, &y)) {
        printf("Expected value: (-79, 43); actual value (%f, %f)\n", x, y);
    } else {
        printf("This shouldn't happen.\n");
    }
}

Basically, the axis-related code is fine, but there are two unrelated issues that work together to cause the ultimate problem.

1. AutoIdentifyEPSG() adds an EPSG authority based on implied AXIS that aren't present, and when the GEOCS is pulled out via CloneGeogCS(), it inserts AXIS tags that contradict the EPSG AUTHORITY tag that was added.
2. When setting up a coordinate transform, the attempt to use an ESPG code instead of WKT flaw where it doesn't properly consider that the EPSG code specified may not match the rest of the WKT, specifically the AXIS tags (see comment in code above).

Also, not for nothing, but it's a bit strange that AutoIdentifyEPSG() returns OGRERR_UNSUPPORTED_SRS, implying that it didn't do anything, when in fact it did inject some AUTHORITY tags.

I don't know what the best solution to the two problems above are (though maybe my suggestion on using SetDatAxisToSRSAxisMapping() would solve point 2) -- I suspect this can only happen with CloneGeogCS() and not other aspects of the WKT so maybe just check for contradictory AUTHORITY tags when AXIS tags are added as a conseuence of 'CloneGeogCS()`?

[AndrewAtAvenza]

The solution for us is basically not to call AutoIdentifyEPSG() but I feel like this is going to burn other people, and given the hoopla around axes, it has the potential to muddy the waters for people (like me!) who think they have a handle on the 3.x changes and then get thrown a loop like this.

[AndrewAtAvenza]

Oh wow, you're fast. Many thanks!