GDAL 3.0 Coordinate transformation (backwards compatibility (?))

[RutgerK]

Expected behavior and actual behavior.

When using the Python bindings, the order of the input arguments for the TransformPoint method of a CoordinateTransformation appears to be flipped compared to GDAL 2.4.

I would expect the format to be:
x,y = ct.TransformPoint(x,y)

But GDAL 3.0 behaves like:
x,y = ct.TransformPoint(y,x)

Especially the mismatch between the input and output order is confusing.

I'm not completely sure this is a bug, since the migration guide mentions some changes with 'axis', but the new behavior seems very strange to me. I encountered this when debugging geom.TransformTo which caused my geometry to have flipped lat/lon.

Steps to reproduce the problem.

import gdal

srs_4326 = gdal.osr.SpatialReference()
srs_4326.ImportFromEPSG(4326)

srs_3857 = gdal.osr.SpatialReference()
srs_3857.ImportFromEPSG(3857)

ct_4326_to_3857 = osr.CoordinateTransformation(srs_4326, srs_3857)

lat = 50.0
lon = -120.0

mapx, mapy, z = ct_4326_to_3857.TransformPoint(lon, lat)

With GDAL 3.0.0, the above code returns:
(inf, inf)

With GDAL 2.4.1, the above code returns:
(-13358338.89519283, 6446275.8410171615)

Using GDAL 3.0.0, but changing the input order to lat/lon returns:
mapx, mapy, z = ct_4326_to_3857.TransformPoint(lat, lon)
(-13358338.895192828, 6446275.841017158)

Both 3.0 and 2.4 show the expected behavior when using gdaltransform.exe from the command line, this is the output of 3.0:

The docstring of TransformPoint changed a little bit from 2.4 to 3.0, but both docstrings suggest that the order of inputs should be (x,y,z).

Some more comparisons (with varying x/y orders can be found in these notebooks:

GDAL 2.4:
https://nbviewer.jupyter.org/gist/RutgerK/75ee3016e58999b60469fd479655cf18

GDAL 3.0:
https://nbviewer.jupyter.org/gist/RutgerK/9c8d5255f566bafcb4bb543cc0f4eb7c

Since I've only used builds from conda, I can't rule out that it's related to the package/build process.

Operating system

Windows 10 (64bit)

GDAL version and provenance

GDAL 3.0.0 from conda-forge

Specifically:
https://conda.anaconda.org/conda-forge/win-64/libgdal-3.0.0-h47faea2_1.tar.bz2
https://conda.anaconda.org/conda-forge/win-64/gdal-3.0.0-py37hdf5ee75_1.tar.bz2

[rouault]

The change is intended.
The migration guide suggests a fix/workaround by calling srs_4326.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)

[naught101]

Just to be clear, the if code needs to run with an arbitrary version of GDAL, you can do this:

import osgeo

# ... so something to create and SRS:
srs = SpatialReference()

if int(osgeo.__version__[0]) >= 3:
    # GDAL 3 changes axis order: https://github.com/OSGeo/gdal/issues/1546
    srs.SetAxisMappingStrategy(osgeo.osr.OAMS_TRADITIONAL_GIS_ORDER)

it would probably better to invert it the other way, eventually, and depend on gdal>=3.0, but I guess there's no way to force forward compatibility with older versions?

[alexgrayeccc]

The weird thing is if you use Transform on a geometry you get the same problem. I now have to reverse the X,Y in the Point geometry I create for the transformation to work. I can't imagine having to do this on polygon geometries, parsing out each point and reversing the X,Y to project it back correctly...

 inSpatialRef = osr.SpatialReference()
    inSpatialRef.ImportFromEPSG(WSG84)
    esri_spatial_ref = osr.SpatialReference(wkt='PROJCS["PS100",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic_North_Pole"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-100.0],PARAMETER["Standard_Parallel_1",65.0],UNIT["Meter",1.0]]')
    esri_spatial_ref.MorphFromESRI()  
    NewPoint = ogr.Geometry(ogr.wkbPoint)
    Ylat = 63.04444613939555
    Xlon = -65.64445686919069
    NewPoint.AddPoint(Ylat, Xlon)

    coordTransform = osr.CoordinateTransformation(inSpatialRef, esri_spatial_ref)

    NewPoint.Transform(coordTransform)

    print("X:")
    print (NewPoint.GetX())
    print("Y:")
    print (NewPoint.GetY())

[janosrusiczki]

If you're wondering why phyghtmap 2.23 reverses coordinates, using latitude for longitude and viceversa, it's this issue. Until there's release with a fix for phyghtmap I have provided steps on how to install a patched version. Thanks @naught101 for the copy/paste ready fix. 😄