Merge pull request OSGeo#5339 from rouault/OSR_DEFAULT_AXIS_MAPPING_STRATEGY

Add a OSR_DEFAULT_AXIS_MAPPING_STRATEGY configuration option

Author: rouault

autotest/cpp/test_osr_ct.cpp

@@ -180,9 +180,11 @@ namespace tut
     void object::test<4>()
     {
         OGRSpatialReference oSRSSource;
+        oSRSSource.SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT);
         oSRSSource.importFromEPSG(4267);
 
         OGRSpatialReference oSRSTarget;
+        oSRSTarget.SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT);
         oSRSTarget.importFromEPSG(4269);
 
         auto poCT = std::unique_ptr<OGRCoordinateTransformation>(

autotest/gcore/vrtmisc.py

@@ -509,6 +509,7 @@ def test_vrtmisc_write_srs():
     tmpfile = '/vsimem/test_vrtmisc_write_srs.vrt'
     ds = gdal.Translate(tmpfile, 'data/byte.tif', format='VRT')
     sr = osr.SpatialReference()
+    sr.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     sr.ImportFromEPSG(4326)
     ds.SetSpatialRef(sr)
     ds = None

autotest/ogr/ogr_geojson.py

@@ -3126,6 +3126,7 @@ def test_ogr_geojson_export_geometry_axis_order():
     # EPSG:4326 and lat,long data order
     sr = osr.SpatialReference()
     sr.ImportFromEPSG(4326)
+    sr.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     g = ogr.CreateGeometryFromWkt('POINT (49 2)')
     g.AssignSpatialReference(sr)
     before_wkt = g.ExportToWkt()
@@ -3157,6 +3158,7 @@ def test_ogr_geojson_export_geometry_axis_order():
     # Projected CRS with northing, easting order
     sr = osr.SpatialReference()
     sr.ImportFromEPSG(2393)
+    sr.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     g = ogr.CreateGeometryFromWkt('POINT (49 2)')
     g.AssignSpatialReference(sr)
     assert json.loads(g.ExportToJson()) == { "type": "Point", "coordinates": [ 2.0, 49.0 ] }

autotest/ogr/ogr_gml_geom.py

@@ -1270,6 +1270,7 @@ def test_gml_write_gml3_srs():
 
     srlatlong = osr.SpatialReference()
     srlatlong.SetFromUserInput("EPSG:4326")
+    srlatlong.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
 
     geom = ogr.CreateGeometryFromWkt('POINT(500000 4500000)')
     geom.AssignSpatialReference(sr32631)

autotest/osr/osr_ct.py

@@ -79,6 +79,7 @@ def test_osr_ct_2():
     utm_srs.SetWellKnownGeogCS('WGS84')
 
     ll_srs = osr.SpatialReference()
+    ll_srs.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     ll_srs.SetWellKnownGeogCS('WGS84')
 
     ct = osr.CoordinateTransformation(ll_srs, utm_srs)
@@ -353,8 +354,10 @@ def test_osr_ct_options_area_of_interest():
 
     srs_nad27 = osr.SpatialReference()
     srs_nad27.SetFromUserInput("NAD27")
+    srs_nad27.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     srs_wgs84 = osr.SpatialReference()
     srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     options = osr.CoordinateTransformationOptions()
     assert not options.SetAreaOfInterest(-200,40,-99,41)
     assert not options.SetAreaOfInterest(-100,-100,-99,41)
@@ -543,8 +546,10 @@ def test_osr_ct_take_into_account_srs_coordinate_epoch():
 
     s = osr.SpatialReference()
     s.SetFromUserInput("EPSG:7844") # GDA2020
+    s.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
 
     t_2020 = osr.SpatialReference()
+    t_2020.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     t_2020.SetFromUserInput("EPSG:9000") # ITRF2014
     t_2020.SetCoordinateEpoch(2020)
 
@@ -556,6 +561,7 @@ def test_osr_ct_take_into_account_srs_coordinate_epoch():
     assert y == pytest.approx(150, abs=1e-10)
 
     t_2030 = osr.SpatialReference()
+    t_2030.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     t_2030.SetFromUserInput("EPSG:9000") # ITRF2014
     t_2030.SetCoordinateEpoch(2030)
 
@@ -595,6 +601,7 @@ def test_osr_ct_only_axis_order_different():
 
     t = osr.SpatialReference()
     t.ImportFromEPSG(4326)
+    t.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
 
     ct = osr.CoordinateTransformation(s_wrong_axis_order, t)
     x, y, _ = ct.TransformPoint(2, 49, 0)
@@ -624,6 +631,7 @@ def test_osr_ct_wkt_non_consistent_with_epsg_definition():
 
     t = osr.SpatialReference()
     t.ImportFromEPSG(4326)
+    t.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
 
     ct = osr.CoordinateTransformation(s_wrong_axis_order, t)
     x, y, _ = ct.TransformPoint(2, 49, 0)

autotest/osr/osr_ct_proj.py

@@ -194,10 +194,12 @@ def test_proj(src_srs, src_xyz, src_error,
             pytest.skip(f'PROJ version < {proj_version_req}')
 
     src = osr.SpatialReference()
+    src.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert src.SetFromUserInput(src_srs) == 0, \
         ('SetFromUserInput(%s) failed.' % src_srs)
 
     dst = osr.SpatialReference()
+    dst.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert dst.SetFromUserInput(dst_srs) == 0, \
         ('SetFromUserInput(%s) failed.' % dst_srs)
 
@@ -261,6 +263,7 @@ def test_proj(src_srs, src_xyz, src_error,
 )
 def test_transform_bounds_densify(density, expected):
     src = osr.SpatialReference()
+    src.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert src.ImportFromEPSG(4326) == 0
     dst = osr.SpatialReference()
     assert dst.ImportFromProj4(
@@ -316,6 +319,7 @@ def test_transform_bounds_densify_out_of_bounds__geographic_output():
        "+a=6370997 +b=6370997 +units=m +no_defs"
     ) == 0
     dst = osr.SpatialReference()
+    dst.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert dst.ImportFromEPSG(4326) == 0
     ctr = osr.CoordinateTransformation(src, dst)
     assert ctr.TransformBounds(
@@ -325,8 +329,10 @@ def test_transform_bounds_densify_out_of_bounds__geographic_output():
 
 def test_transform_bounds_antimeridian():
     src = osr.SpatialReference()
+    src.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert src.ImportFromEPSG(4167) == 0
     dst = osr.SpatialReference()
+    dst.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert dst.ImportFromEPSG(3851) == 0
     ctr = osr.CoordinateTransformation(src, dst)
     assert ctr.TransformBounds(
@@ -423,8 +429,10 @@ def test_transform_bounds__noop_geographic():
 
 def test_transform_bounds__north_pole():
     src = osr.SpatialReference()
+    src.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert src.ImportFromEPSG(32661) == 0
     dst = osr.SpatialReference()
+    dst.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert dst.ImportFromEPSG(4326) == 0
     ctr = osr.CoordinateTransformation(src, dst)
     assert ctr.TransformBounds(
@@ -459,8 +467,10 @@ def test_transform_bounds__north_pole__xy():
 
 def test_transform_bounds__south_pole():
     src = osr.SpatialReference()
+    src.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert src.ImportFromEPSG(32761) == 0
     dst = osr.SpatialReference()
+    dst.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     assert dst.ImportFromEPSG(4326) == 0
     ctr = osr.CoordinateTransformation(src, dst)
     assert ctr.TransformBounds(

autotest/osr/osr_esri.py

@@ -917,6 +917,7 @@ def test_osr_esri_30():
            '']
 
     srs_prj = osr.SpatialReference()
+    srs_prj.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     srs_prj.ImportFromESRI(prj)
 
     wkt = """GEOGCS["unknown",
@@ -926,6 +927,7 @@ def test_osr_esri_30():
     UNIT["degree",0.0174532925199433]]"""
 
     srs_wkt = osr.SpatialReference(wkt=wkt)
+    srs_wkt.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
 
     assert srs_prj.IsSame(srs_wkt)
 

autotest/osr/osr_ozi.py

@@ -79,6 +79,7 @@ def test_osr_ozi_2():
 def test_osr_ozi_3():
 
     srs = osr.SpatialReference()
+    srs.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     srs.ImportFromOzi(["OziExplorer Map Data File Version 2.2",
                        "Test_Map",
                        "Test_Map.png",

autotest/pymod/gdaltest.py

@@ -858,9 +858,11 @@ def user_srs_to_wkt(user_text):
 
 def equal_srs_from_wkt(expected_wkt, got_wkt, verbose=True):
     expected_srs = osr.SpatialReference()
+    expected_srs.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     expected_srs.ImportFromWkt(expected_wkt)
 
     got_srs = osr.SpatialReference()
+    got_srs.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
     got_srs.ImportFromWkt(got_wkt)
 
     if got_srs.IsSame(expected_srs):

autotest/pyscripts/test_osr_util.py

@@ -30,7 +30,7 @@
 
 import pytest
 
-from osgeo import osr
+from osgeo import gdal, osr
 
 # test that osgeo_utils and numpy are available, otherwise skip all tests
 pytest.importorskip('osgeo_utils')
@@ -40,7 +40,8 @@
 
 def test_gis_order():
     pj4326_gis2 = osr_util.get_srs(4326)  # tests the correct default
-    assert pj4326_gis2.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+    if gdal.GetConfigOption('OSR_DEFAULT_AXIS_MAPPING_STRATEGY') is None:
+        assert pj4326_gis2.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
 
     pj4326_auth = osr_util.get_srs(4326, axis_order=osr.OAMS_AUTHORITY_COMPLIANT)
     assert pj4326_auth.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
@@ -49,23 +50,26 @@ def test_gis_order():
     assert pj4326_gis.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
 
     assert not osr_util.are_srs_equivalent(pj4326_gis, pj4326_auth)
-    assert osr_util.are_srs_equivalent(pj4326_auth, 4326)
-    assert not osr_util.are_srs_equivalent(pj4326_gis, 4326)
+    if gdal.GetConfigOption('OSR_DEFAULT_AXIS_MAPPING_STRATEGY') is None:
+        assert osr_util.are_srs_equivalent(pj4326_auth, 4326)
+        assert not osr_util.are_srs_equivalent(pj4326_gis, 4326)
 
     pj4326_str = osr_util.get_srs_pj(pj4326_auth)
     pj4326_str2 = osr_util.get_srs_pj(pj4326_gis)
 
     # axis order is not reflected in proj strings
     assert isinstance(pj4326_str, str) and pj4326_str == pj4326_str2
 
-    assert osr_util.are_srs_equivalent(pj4326_str, 4326)
+    if gdal.GetConfigOption('OSR_DEFAULT_AXIS_MAPPING_STRATEGY') is None:
+        assert osr_util.are_srs_equivalent(pj4326_str, 4326)
     assert osr_util.are_srs_equivalent(pj4326_auth, pj4326_str)
     assert not osr_util.are_srs_equivalent(pj4326_gis, pj4326_str)
 
     osr_util.set_default_axis_order(osr.OAMS_TRADITIONAL_GIS_ORDER)  # sets gis order
 
     srs = osr_util.get_srs(4326)  # check the the default was changed
-    assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
+    if gdal.GetConfigOption('OSR_DEFAULT_AXIS_MAPPING_STRATEGY') is None:
+        assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
 
     # check that srs object is not affected by default
     srs = osr_util.get_srs(pj4326_auth)
@@ -94,10 +98,12 @@ def test_gis_order():
     osr_util.set_default_axis_order()
 
     srs = osr_util.get_srs(4326)  # check the the default was restored
-    assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+    if gdal.GetConfigOption('OSR_DEFAULT_AXIS_MAPPING_STRATEGY') is None:
+        assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
 
     srs = osr_util.get_srs(pj4326_str)
-    assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+    if gdal.GetConfigOption('OSR_DEFAULT_AXIS_MAPPING_STRATEGY') is None:
+        assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
 
     # check that srs object is not affected by default
     srs = osr_util.get_srs(pj4326_gis)  # check that srs object is also affected
@@ -118,11 +124,13 @@ def test_gis_order2():
     srs_from_epsg = osr.SpatialReference()
     srs_from_epsg.ImportFromEPSG(4326)
 
-    assert srs_from_epsg.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+    if gdal.GetConfigOption('OSR_DEFAULT_AXIS_MAPPING_STRATEGY') is None:
+        assert srs_from_epsg.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
     srs_from_str = osr.SpatialReference()
     srs_from_str.ImportFromProj4(pj4326_str)
-    assert srs_from_str.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
-    assert srs_from_epsg.IsSame(srs_from_str)
+    if gdal.GetConfigOption('OSR_DEFAULT_AXIS_MAPPING_STRATEGY') is None:
+        assert srs_from_str.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+        assert srs_from_epsg.IsSame(srs_from_str)
 
     # testing that explicitly setting OAMS_AUTHORITY_COMPLIANT does not effect equivalence
     srs_from_epsg.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)

ogr/ogrspatialreference.cpp

@@ -166,9 +166,34 @@ struct OGRSpatialReference::Private
     void                refreshAxisMapping();
 };
 
+static OSRAxisMappingStrategy GetDefaultAxisMappingStrategy()
+{
+    const char* pszDefaultAMS = CPLGetConfigOption(
+        "OSR_DEFAULT_AXIS_MAPPING_STRATEGY", nullptr);
+    if( pszDefaultAMS )
+    {
+        if( EQUAL(pszDefaultAMS, "AUTHORITY_COMPLIANT") )
+            return OAMS_AUTHORITY_COMPLIANT;
+        else if( EQUAL(pszDefaultAMS, "TRADITIONAL_GIS_ORDER") )
+            return OAMS_TRADITIONAL_GIS_ORDER;
+        else
+        {
+            CPLError(CE_Failure, CPLE_AppDefined,
+                     "Illegal value for OSR_DEFAULT_AXIS_MAPPING_STRATEGY = %s",
+                     pszDefaultAMS);
+        }
+    }
+    return OAMS_AUTHORITY_COMPLIANT;
+}
+
 OGRSpatialReference::Private::Private():
     m_poListener(std::shared_ptr<Listener>(new Listener(this)))
 {
+    // Get the default value for m_axisMappingStrategy from the
+    // OSR_DEFAULT_AXIS_MAPPING_STRATEGY configuration option, if set.
+    static const OSRAxisMappingStrategy defaultAxisMappingStrategy = GetDefaultAxisMappingStrategy();
+
+    m_axisMappingStrategy = defaultAxisMappingStrategy;
 }
 
 OGRSpatialReference::Private::~Private()
@@ -708,6 +733,18 @@ void OGRsnPrintDouble( char * pszStrBuf, size_t size, double dfValue )
  *
  * Note that newly created objects are given a reference count of one.
  *
+ * Starting with GDAL 3.0, coordinates associated with a OGRSpatialReference
+ * object are assumed to be in the order of the axis of the CRS definition (which
+ * for example means latitude first, longitude second for geographic CRS belonging
+ * to the EPSG authority). It is possible to define a data axis to CRS axis
+ * mapping strategy with the SetAxisMappingStrategy() method.
+ *
+ * Starting with GDAL 3.5, the OSR_DEFAULT_AXIS_MAPPING_STRATEGY configuration
+ * option can be set to "TRADITIONAL_GIS_ORDER" / "AUTHORITY_COMPLIANT" (the later
+ * being the default value when the option is not set) to control the value of the
+ * data axis to CRS axis mapping strategy when a OSRSpatialReference object is
+ * created. Calling SetAxisMappingStrategy() will override this default value.
+
  * The C function OSRNewSpatialReference() does the same thing as this
  * constructor.
  *
@@ -729,6 +766,18 @@ OGRSpatialReference::OGRSpatialReference( const char * pszWKT ) :
 /**
  * \brief Constructor.
  *
+ * Starting with GDAL 3.0, coordinates associated with a OGRSpatialReference
+ * object are assumed to be in the order of the axis of the CRS definition (which
+ * for example means latitude first, longitude second for geographic CRS belonging
+ * to the EPSG authority). It is possible to define a data axis to CRS axis
+ * mapping strategy with the SetAxisMappingStrategy() method.
+ *
+ * Starting with GDAL 3.5, the OSR_DEFAULT_AXIS_MAPPING_STRATEGY configuration
+ * option can be set to "TRADITIONAL_GIS_ORDER" / "AUTHORITY_COMPLIANT" (the later
+ * being the default value when the option is not set) to control the value of the
+ * data axis to CRS axis mapping strategy when a OSRSpatialReference object is
+ * created. Calling SetAxisMappingStrategy() will override this default value.
+ *
  * This function is the same as OGRSpatialReference::OGRSpatialReference()
  */
 OGRSpatialReferenceH CPL_STDCALL OSRNewSpatialReference( const char *pszWKT )
@@ -11628,6 +11677,12 @@ OSRAxisMappingStrategy OSRGetAxisMappingStrategy( OGRSpatialReferenceH hSRS )
 /************************************************************************/
 
 /** \brief Set the data axis to CRS axis mapping strategy.
+ *
+ * Starting with GDAL 3.5, the OSR_DEFAULT_AXIS_MAPPING_STRATEGY configuration
+ * option can be set to "TRADITIONAL_GIS_ORDER" / "AUTHORITY_COMPLIANT" (the later
+ * being the default value when the option is not set) to control the value of the
+ * data axis to CRS axis mapping strategy when a OSRSpatialReference object is
+ * created. Calling SetAxisMappingStrategy() will override this default value.
  *
  * See OGRSpatialReference::GetAxisMappingStrategy()
  * @since GDAL 3.0