Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reprojection issue #110

Closed
sindizzy opened this issue Aug 14, 2023 · 3 comments
Closed

Reprojection issue #110

sindizzy opened this issue Aug 14, 2023 · 3 comments
Assignees
Labels
duplicate This issue or pull request already exists not-a-bug This is intended behavior and not a bug

Comments

@sindizzy
Copy link

sindizzy commented Aug 14, 2023

Describe the bug
Using TransformTo() method results in an exception: System.ApplicationException: Invalid coordinate.
Possible issue with GDAL: OSGeo/gdal#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 sindizzy added the bug Something isn't working label Aug 14, 2023
@sindizzy
Copy link
Author

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
Copy link
Author

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 MaxRev-Dev added duplicate This issue or pull request already exists not-a-bug This is intended behavior and not a bug and removed bug Something isn't working labels Aug 14, 2023
@MaxRev-Dev
Copy link
Owner

Indeed, this behaviour is related to incompatible changes between major GDAL versions. Here's another related issue - #13.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
duplicate This issue or pull request already exists not-a-bug This is intended behavior and not a bug
Projects
None yet
Development

No branches or pull requests

2 participants