Reprojection issue

[sindizzy]

Describe the bug Using TransformTo() method results in an exception: System.ApplicationException: Invalid coordinate. Possible issue with GDAL: https://github.com/OSGeo/gdal/issues/1546

Steps to reproduce I am trying to reproject a point from one CRS to another CRS and am using this simple example. From what I can deduce most of these functions act like a typical Windows API where if the result is successful then it returns a zero.

I have translated to C# as follows:

var inSr = new SpatialReference(null);
var inResI = inSr.ImportFromEPSG(4326); // WGS84/Geographic
var inResE = inSr.ExportToWkt(out string inWkt, null);
Debug.Print("inResI={0} inResE={1}", inResI, inResE);
Debug.Print("inWkt={0}", inWkt);

var outSr = new SpatialReference(null);
var outResI = outSr.ImportFromEPSG(32756); // WGS84 UTM Zone 56 South
var outResE = outSr.ExportToWkt(out string outWkt, null);
Debug.Print("outResI={0} outResE={1}", outResI, outResE);
Debug.Print("outWkt={0}", outWkt);

var pt = new Geometry(wkbGeometryType.wkbPoint);
Debug.Print("4326 coords: {0},{1}", 150.700532, -35.145618);
pt.AddPoint(150.700532, -35.145618, 0);
pt.AssignSpatialReference(inSr);
var res = pt.TransformTo(outSr);
Debug.Print("32756 coords: res={0} {1},{2}", res, pt.GetX, pt.GetY);

and the results in Debug

inResI=0 inResE=0
inWkt=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"]]
outResI=0 outResE=0
outWkt=PROJCS["WGS 84 / UTM zone 56S",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"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
4326 coords: 150.700532,-35.145618
Exception thrown: 'System.ApplicationException' in MaxRev.Gdal.Core.dll
System.ApplicationException: Invalid coordinate
   at OSGeo.OGR.Geometry.TransformTo(SpatialReference reference)

Environment information:

• OS (version): Windows 10 Enterprise 21H2
• Package version (core): 3.7.0.217
• Package version (runtime): Windows 3.7.0.100

[sindizzy]

The library returns

IsConfigured = True
GDAL Version = 3.7.0
GDAl Info = GDAL 3.7.0, released 2023/05/11
GDAl Num = 3070000

So it does seem like its an intended change in how certain functions work and there is a migration guide.

[sindizzy]

OK so I answered my own question but I will leave this here for those that are encountering the same issue. here is how I adjusted my code for this change in GDAL:

var inSr = new SpatialReference(null);
var inResI = inSr.ImportFromEPSG(4326); // WGS84/Geographic
var inResE = inSr.ExportToWkt(out string inWkt, null);
Debug.Print("inResI={0} inResE={1}", inResI, inResE);
Debug.Print("inWkt={0}", inWkt);

var outSr = new SpatialReference(null);
var outResI = outSr.ImportFromEPSG(32756); // WGS84 UTM Zone 56 South
var outResE = outSr.ExportToWkt(out string outWkt, null);
Debug.Print("outResI={0} outResE={1}", outResI, outResE);
Debug.Print("outWkt={0}", outWkt);

// GDAL changed the way that coordinates are provided to certain geoprocessing functions
// "Starting with GDAL 3.0, the axis order mandated by the authority defining a CRS is by default honoured by the OGRCoordinateTransformation class,
// and always exported in WKT1. Consequently CRS created with the "EPSG:4326" or "WGS84" strings use the latitude first, longitude second axis order.
var verCur = Version.Parse(Gdal.VersionInfo("RELEASE_NAME"));
var verChk = Version.Parse("3.0.0.0");
if (verCur >= verChk)
{
    Debug.Print("Current GDAL is >= 3.0.0.0 so set Axis Mapping Strategy");
    inSr.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);
    outSr.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);
}

var pt = new Geometry(wkbGeometryType.wkbPoint);
Debug.Print("4326 coords: {0},{1}", 150.700532, -35.145618);
pt.AddPoint(150.700532, -35.145618, 0);
pt.AssignSpatialReference(inSr);
var res = pt.TransformTo(outSr);
Debug.Print("32756 coords: res={0} {1},{2}", res, pt.GetX(0), pt.GetY(0));

and the debug output:

inResI=0 inResE=0
inWkt=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"]]
outResI=0 outResE=0
outWkt=PROJCS["WGS 84 / UTM zone 56S",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"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
Current GDAL is >= 3.0.0.0 so set Axis Mapping Strategy
4326 coords: 150.700532,-35.145618
32756 coords: res=0 290523.0254286239,6108387.720559284

[MaxRev-Dev]

Indeed, this behaviour is related to incompatible changes between major GDAL versions. Here's another related issue - https://github.com/MaxRev-Dev/gdal.netcore/issues/13.