It seems that the Migration Guide has the answer: > OSRImportFromEPSG() takes into account official axis order. > Traditional GIS-friendly axis order can be restored with > OGRSpatialReference::SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
In my previous example, if one compares the responses to GetAxisMappingStrategy(), the original SRS and that returned by the layer are different: >>> projection.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT # True >>> projFromLayer.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER # True Is the change to this value in this circumstance intentional? If so, does it mean that I have to repeatedly enforce the axis mapping strategy, because I don't know whether it will have been modified or not, depending on the provenance of the SRS object? Dr. Daniel Evans Software Developer Skype<sip:[email protected]> From: Daniel Evans Sent: 06 February 2020 10:31 To: [email protected] Subject: RE: WGS84 -> WGS84 transform flips axes 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]<mailto:[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
