[gdal-dev] WGS84 -> WGS84 transform flips axes

Thu Feb 6 02:30:35 PST 2020

To give a somewhat more concise example:

from osgeo import gdal, osr, ogr

projection = osr.SpatialReference()
projection.ImportFromEPSG(4326)

temp_ds = gdal.GetDriverByName('Memory').Create('', 0, 0, 0)
layer = temp_ds.CreateLayer("foo", projection, geom_type=ogr.wkbPolygon)
projFromLayer = layer.GetSpatialRef()

xform = osr.CoordinateTransformation(projection, projFromLayer)
xform.TransformPoint(1.0, 0.0)
# Outputs (0.0, 1.0, 0.0)

xform_no_layer = osr.CoordinateTransformation(projection, projection)
xform_no_layer.TransformPoint(1.0, 0.0)
# Outputs (1.0, 0.0, 0.0)

Dr. Daniel Evans
Software Developer

Skype<sip:daniel.evans at jbarisk.com>

From: Daniel Evans
Sent: 06 February 2020 10:21
To: gdal-dev at lists.osgeo.org
Subject: WGS84 -> WGS84 transform flips axes

I am currently looking at upgrading from GDAL 2.x to GDAL 3.x. Most of my software tests pass, but I have a few failing because a transform from WGS84 to WGS84 is flipping the axes on transformed geometries.

I have a recollection of reading about GDAL 3 (or PROJ 6) resulting in axis swapping in some cases, but can't find the discussion/Github issue again.

The two SRS objects appear identical, and seem to report the same axis order, but a point at (1, 0) gets transformed to (0, 1) - Python code below:

>>> from_proj =  ds.GetLayer().GetSpatialRef()
>>> to_proj = osr.SpatialReference()
>>> to_proj.ImportFromEPSG(4326)
>>> to_proj.GetAxisName(None, 0)
'Geodetic latitude'
>>> to_proj.GetAxisName(None, 1)
'Geodetic longitude'
>>> from_proj.GetAxisName(None, 0)
'Geodetic latitude'
>>> from_proj.GetAxisName(None, 1)
'Geodetic longitude'
>>> xform = osr.CoordinateTransformation(from_proj, to_proj)
>>> xform.TransformPoint(1.0, 0.0)
(0.0, 1.0, 0.0)

If I follow from_proj back to its origin, it is also imported by importing from an EPSG code, which is again 4326:

>>> projection = osr.SpatialReference()
>>> projection.ImportFromEPSG(epsg)
>>> temp_ds = gdal.GetDriverByName('Memory').Create('', 0, 0, 0)
>>> layer = temp_ds.CreateLayer(name, projection, geom_type=geom_type)

The two SRS objects appear identical when exported to WKT:

>>> from_proj.ExportToWkt()
'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"]]'
>>> to_proj.ExportToWkt()
'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"]]'

And if I re-import the SRSes from those WKT strings, the transform starts to behave:

>>> from_proj_reimported = osr.SpatialReference(); from_proj_reimported.ImportFromWkt(from_proj.ExportToWkt())
>>> to_proj_reimported = osr.SpatialReference(); to_proj_reimported.ImportFromWkt(to_proj.ExportToWkt())
>>> new_xform = osr.CoordinateTransformation(from_proj_reimported, to_proj_reimported)
>>> new_xform.TransformPoint(1.0, 0.0)
(1.0, 0.0, 0.0)

How are these two SRS objects ending up with flipped axes, how can I actually tell that the axes are flipped, and what do I need to do to un-flip them (short of exporting and reimporting every time I use them)?

Dr. Daniel Evans
Software Developer

Skype<sip:daniel.evans at jbarisk.com>

T +44 (0) 1756 799919
www.jbarisk.com<http://www.jbarisk.com>

[Visit our website]<http://www.jbarisk.com> [http://www.jbagroup.co.uk/imgstore/JBA-Email-Sig-Icons-LINKEDIN.png]  [Follow us on Twitter] <https://twitter.com/jbarisk>
Our postal address and registered office is JBA Risk Management Limited, 1 Broughton Park, Old Lane North, Broughton, Skipton, North Yorkshire BD23 3FD.

Find out more about us here: www.jbarisk.com<http://www.jbarisk.com/> and follow us on Twitter @JBARisk<http://twitter.com/JBARisk> and LinkedIn<https://www.linkedin.com/company/2370847?trk=tyah&trkInfo=clickedVertical%3Acompany%2CclickedEntityId%3A2370847%2Cidx%3A2-1-2%2CtarId%3A1447414259786%2Ctas%3AJBA%20RISK%20MANAGEMENT>

The JBA Group supports the JBA Trust.

All JBA Risk Management's email messages contain confidential information and are intended only for the individual(s) named. If you are not the named addressee you should not disseminate, distribute or copy this e-mail.
Please notify the sender immediately by email if you have received this email by mistake and delete this email from your system.

JBA Risk Management Limited is registered in England, company number 07732946, 1 Broughton Park, Old Lane North, Broughton, Skipton, North Yorkshire, BD23 3FD, Telephone: +441756799919

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200206/3ac98533/attachment-0001.html>