fix lat/lon axis order when using gdal 3 (#158)

Author: letmaik

gis4wrf/core/crs.py

@@ -180,6 +180,7 @@ def lonlat_srs(self) -> osr.SpatialReference:
         srs = self.srs
         datum = self.srs.GetAttrValue('datum')
         srs_out = osr.SpatialReference()
+        fix_axis_order(srs_out)
         srs_out.SetGeogCS('', datum, '', srs.GetSemiMajor(), srs.GetInvFlattening())
         assert not srs_out.EPSGTreatsAsLatLong(), 'expected lon/lat axis order'
         return srs_out
@@ -190,4 +191,9 @@ def transform_point(srs_in: osr.SpatialReference, srs_out: osr.SpatialReference,
         point_geom.AddPoint(point.x, point.y)
         transform = osr.CoordinateTransformation(srs_in, srs_out)
         point_geom.Transform(transform)
-        return Coordinate2D(x=point_geom.GetX(), y=point_geom.GetY())
+        return Coordinate2D(x=point_geom.GetX(), y=point_geom.GetY())
+
+def fix_axis_order(srs):
+    # https://github.com/OSGeo/gdal/blob/release/3.0/gdal/MIGRATION_GUIDE.TXT
+    if hasattr(osr, 'OAMS_TRADITIONAL_GIS_ORDER'):
+        srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

tests/core/test_crs.py

@@ -0,0 +1,42 @@
+# GIS4WRF (https://doi.org/10.5281/zenodo.1288569)
+# Copyright (c) 2019 D. Meyer and M. Riechert. Licensed under MIT.
+
+import pytest
+
+from gis4wrf.core import (
+    CRS, LonLat, Coordinate2D
+)
+
+@pytest.mark.parametrize('crs_name', ['lonlat', 'lambert', 'mercator', 'polar', 'albers_nad83'])
+def test_geo_roundtrip(crs_name: str):
+    if crs_name == 'lonlat':
+        crs = CRS.create_lonlat()
+    elif crs_name == 'lambert':
+        crs = CRS.create_lambert(truelat1=3.5, truelat2=7, origin=LonLat(lon=4, lat=0))
+    elif crs_name == 'mercator':
+        crs = CRS.create_mercator(truelat1=3.5, origin_lon=4)
+    elif crs_name == 'polar':
+        crs = CRS.create_polar(truelat1=3.5, origin_lon=4)
+    elif crs_name == 'albers_nad83':
+        crs = CRS.create_albers_nad83(truelat1=3.5, truelat2=7, origin=LonLat(lon=4, lat=0))
+    lonlat = LonLat(lon=10, lat=30)
+    xy = crs.to_xy(lonlat)
+    lonlat2 = crs.to_lonlat(xy)
+    assert lonlat.lon == pytest.approx(lonlat2.lon)
+    assert lonlat.lat == pytest.approx(lonlat2.lat)
+
+@pytest.mark.parametrize('crs_name', ['lambert', 'mercator', 'polar', 'albers_nad83'])
+def test_projection_origin(crs_name: str):
+    origin_lonlat = LonLat(lon=4, lat=0)
+    if crs_name == 'lambert':
+        crs = CRS.create_lambert(truelat1=3.5, truelat2=7, origin=origin_lonlat)
+    elif crs_name == 'mercator':
+        crs = CRS.create_mercator(truelat1=3.5, origin_lon=origin_lonlat.lon)
+    elif crs_name == 'polar':
+        origin_lonlat = LonLat(lon=origin_lonlat.lon, lat=90)
+        crs = CRS.create_polar(truelat1=3.5, origin_lon=origin_lonlat.lon)
+    elif crs_name == 'albers_nad83':
+        crs = CRS.create_albers_nad83(truelat1=3.5, truelat2=7, origin=origin_lonlat)
+    origin_xy = crs.to_xy(origin_lonlat)
+    assert origin_xy.x == pytest.approx(0)
+    assert origin_xy.y == pytest.approx(0, abs=1e-9)