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:[email protected]>
From: Daniel Evans
Sent: 06 February 2020 10:21
To: [email protected]
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:[email protected]>
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
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
