gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py - get_s…

…rs() - add `default_axis_order` flag (with get and set) to allow the user to change the default; add get_srs_pj(), are_srs_equivalent() support functions

autotest/pyscripts/test_osr_util.py.py - add test get_srs() and get_transform() with different axis orders

autotest/pyscripts/test_osr_util.py

@@ -0,0 +1,164 @@
+#!/usr/bin/env pytest
+# -*- coding: utf-8 -*-
+###############################################################################
+#
+# Project:  GDAL/OGR Test Suite
+# Purpose:  osgeo_utils.auxiliary.osr_utils (gdal-utils) testing
+# Author:   Idan Miara <[email protected]>
+#
+###############################################################################
+# Copyright (c) 2021, Idan Miara <[email protected]>
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included
+# in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
+###############################################################################
+
+import pytest
+
+from osgeo import osr
+
+# test that osgeo_utils and numpy are available, otherwise skip all tests
+pytest.importorskip('osgeo_utils')
+
+from osgeo_utils.auxiliary import osr_util, array_util
+
+
+def test_gis_order():
+    pj4326_gis2 = osr_util.get_srs(4326)  # tests the correct default
+    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
+
+    pj4326_gis = osr_util.get_srs(4326, axis_order=osr.OAMS_TRADITIONAL_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)
+
+    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)
+    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
+
+    # check that srs object is not affected by default
+    srs = osr_util.get_srs(pj4326_auth)
+    assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+    assert osr_util.are_srs_equivalent(srs, pj4326_auth)
+
+    # check that srs object is also affected if explicitly set
+    srs = osr_util.get_srs(pj4326_auth, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
+    assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
+
+    # check that the default does not effect explicit order
+    srs = osr_util.get_srs(4326, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
+    assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
+
+    srs = osr_util.get_srs(pj4326_str)
+    assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
+
+    srs = osr_util.get_srs(4326, axis_order=osr.OAMS_AUTHORITY_COMPLIANT)
+    assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+
+    assert not osr_util.are_srs_equivalent(pj4326_gis, pj4326_auth)
+    assert not osr_util.are_srs_equivalent(pj4326_auth, 4326)
+    assert osr_util.are_srs_equivalent(pj4326_gis, 4326)
+
+    # restore the default and repeat some tests
+    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
+
+    srs = osr_util.get_srs(pj4326_str)
+    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
+    assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
+
+    # check that srs object is also affected if explicitly set
+    srs = osr_util.get_srs(pj4326_gis, axis_order=osr.OAMS_AUTHORITY_COMPLIANT)
+    assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+
+    srs = osr_util.get_srs(4326, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
+    assert srs.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
+
+
+def test_gis_order2():
+    pj4326_gis = osr_util.get_srs(4326, axis_order=osr.OAMS_TRADITIONAL_GIS_ORDER)
+    pj4326_str = osr_util.get_srs_pj(4326)
+
+    srs_from_epsg = osr.SpatialReference()
+    srs_from_epsg.ImportFromEPSG(4326)
+
+    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)
+
+    # testing that explicitly setting OAMS_AUTHORITY_COMPLIANT does not effect equivalence
+    srs_from_epsg.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
+    srs_from_str.SetAxisMappingStrategy(osr.OAMS_AUTHORITY_COMPLIANT)
+    assert srs_from_epsg.IsSame(srs_from_str)
+
+    srs_from_epsg.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    srs_from_str.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+
+    srs_from_epsg2 = osr.SpatialReference()
+    srs_from_epsg2.ImportFromEPSG(4326)
+    srs_from_epsg2.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    assert srs_from_epsg.IsSame(srs_from_epsg2)
+
+    test_gis_order_proj_str_vs_epsg = False
+    # explicitly setting OAMS_TRADITIONAL_GIS_ORDER triggers inequality between srs from proj-string and srs from epsg
+    # if this issue is resolved these tests can be enabled
+    if test_gis_order_proj_str_vs_epsg:
+        assert srs_from_epsg.IsSame(srs_from_str)
+        assert osr_util.are_srs_equivalent(pj4326_str, 4326)
+        assert osr_util.are_srs_equivalent(pj4326_gis, pj4326_str)
+
+
+def test_transform():
+    pj_utm = osr_util.get_srs(32636)
+    utm_x = [690950.4640, 688927.6381]
+    utm_y = [3431318.8435, 3542183.4911]
+
+    for gis_order in (False, True):
+        axis_order = osr_util.get_axis_order_from_gis_order(gis_order)
+        pj4326 = osr_util.get_srs(4326, axis_order=axis_order)
+        ct = osr_util.get_transform(pj4326, pj_utm)
+        lon = [35, 35]
+        lat = [31, 32]
+        x, y = (lon, lat) if gis_order else (lat, lon)
+        osr_util.transform_points(ct, x, y)
+        d = array_util.array_dist(x, utm_x), array_util.array_dist(y, utm_y)
+        assert max(d) < 0.01
+

gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py

@@ -35,37 +35,84 @@
 from osgeo_utils.auxiliary.array_util import ArrayLike
 
 AnySRS = Union[str, int, osr.SpatialReference, gdal.Dataset]
+OAMS_AXIS_ORDER = int
+_default_axis_order: Optional[OAMS_AXIS_ORDER] = None
 
 
-def get_srs(srs: AnySRS, gis_order: bool = False) -> osr.SpatialReference:
-    """ returns an SRS object from epsg, pj_string or DataSet or SRS object """
-    if isinstance(srs, osr.SpatialReference):
-        pass
-    elif isinstance(srs, gdal.Dataset):
-        srs = get_srs_from_ds(srs)
-    elif isinstance(srs, int):
-        srs_ = osr.SpatialReference()
-        if srs_.ImportFromEPSG(srs) != ogr.OGRERR_NONE:
-            raise Exception(f'ogr error when parsing srs from epsg:{srs}')
-        srs = srs_
-    elif isinstance(srs, str):
-        srs_ = osr.SpatialReference()
-        if srs_.SetFromUserInput(srs) != ogr.OGRERR_NONE:  # accept PROJ string, WKT, PROJJSON, etc.
-            raise Exception(f'ogr error when parsing srs from proj string: {srs}')
-        srs = srs_
+def get_srs(srs_like: AnySRS,
+            axis_order: Optional[OAMS_AXIS_ORDER] = None) -> osr.SpatialReference:
+    """
+    returns an SRS object from epsg, pj_string or DataSet or SRS object
+
+    gis_order and axis_order are mutually exclusive
+    if gis_order is not None then axis_order is set according to gis_order
+    if axis_order == None -> set axis order to default_axis_order iff the input is not an osr.SpatialRefernce reference,
+    otherwise set requested axis_order
+    """
+    if isinstance(srs_like, osr.SpatialReference):
+        srs = srs_like
+    elif isinstance(srs_like, gdal.Dataset):
+        srs = osr.SpatialReference()
+        srs.ImportFromWkt(srs_like.GetProjection())
+    elif isinstance(srs_like, int):
+        srs = osr.SpatialReference()
+        if srs.ImportFromEPSG(srs_like) != ogr.OGRERR_NONE:
+            raise Exception(f'ogr error when parsing srs from epsg:{srs_like}')
+    elif isinstance(srs_like, str):
+        srs = osr.SpatialReference()
+        if srs.SetFromUserInput(srs_like) != ogr.OGRERR_NONE:  # accept PROJ string, WKT, PROJJSON, etc.
+            raise Exception(f'ogr error when parsing srs from user input: {srs_like}')
     else:
-        raise Exception(f'Unknown SRS: {srs}')
+        raise Exception(f'Unknown SRS: {srs_like}')
 
-    if gis_order and int(osgeo.__version__[0]) >= 3:
+    if axis_order is None and srs != srs_like:
+        axis_order = _default_axis_order
+    if axis_order is not None and int(osgeo.__version__[0]) >= 3:
         # GDAL 3 changes axis order: https://github.com/OSGeo/gdal/issues/1546
-        srs.SetAxisMappingStrategy(osgeo.osr.OAMS_TRADITIONAL_GIS_ORDER)
+        if not isinstance(axis_order, OAMS_AXIS_ORDER):
+            raise Exception(f'Unexpected axis_order: {axis_order}')
+        srs_axis_order = srs.GetAxisMappingStrategy()
+        if axis_order != srs_axis_order:
+            if srs == srs_like:
+                # we don't want to change the input srs, thus create a copy
+                srs = srs.Clone()
+            srs.SetAxisMappingStrategy(axis_order)
     return srs
 
 
-def get_srs_from_ds(ds: gdal.Dataset) -> osr.SpatialReference:
-    srs = osr.SpatialReference()
-    srs.ImportFromWkt(ds.GetProjection())
-    return srs
+def get_axis_order_from_gis_order(gis_order: Optional[bool]):
+    return None if gis_order is None \
+        else osr.OAMS_TRADITIONAL_GIS_ORDER if gis_order \
+        else osr.OAMS_AUTHORITY_COMPLIANT
+
+
+def get_gis_order_from_axis_order(axis_order: Optional[OAMS_AXIS_ORDER]):
+    return None if axis_order is None else axis_order == osr.OAMS_TRADITIONAL_GIS_ORDER
+
+
+def set_default_axis_order(axis_order: Optional[OAMS_AXIS_ORDER] = None) -> Optional[OAMS_AXIS_ORDER]:
+    global _default_axis_order
+    _default_axis_order = axis_order
+    return _default_axis_order
+
+
+def get_default_axis_order() -> Optional[OAMS_AXIS_ORDER]:
+    global _default_axis_order
+    return _default_axis_order
+
+
+def get_srs_pj(srs: AnySRS) -> str:
+    srs = get_srs(srs)
+    srs_pj4 = srs.ExportToProj4()
+    return srs_pj4
+
+
+def are_srs_equivalent(srs1: AnySRS, srs2: AnySRS) -> bool:
+    if srs1 == srs2:
+        return True
+    srs1 = get_srs(srs1)
+    srs2 = get_srs(srs2)
+    return srs1.IsSame(srs2)
 
 
 def get_transform(src_srs: AnySRS, tgt_srs: AnySRS) -> Optional[osr.CoordinateTransformation]:

gdal/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py

@@ -42,7 +42,8 @@
 from osgeo_utils.auxiliary.base import is_path_like
 from osgeo_utils.auxiliary.array_util import ArrayLike, ArrayOrScalarLike
 from osgeo_utils.auxiliary.numpy_util import GDALTypeCodeAndNumericTypeCodeFromDataSet
-from osgeo_utils.auxiliary.osr_util import transform_points, AnySRS, get_transform, get_srs
+from osgeo_utils.auxiliary.osr_util import transform_points, AnySRS, get_transform, get_srs, \
+    get_axis_order_from_gis_order, OAMS_AXIS_ORDER
 from osgeo_utils.auxiliary.util import PathOrDS, open_ds, get_bands, get_scales_and_offsets, get_band_nums
 from osgeo_utils.auxiliary.gdal_argparse import GDALArgumentParser, GDALScript
 
@@ -68,7 +69,7 @@ class LocationInfoOutput(Enum):
 
 def gdallocationinfo(filename_or_ds: PathOrDS,
                      x: ArrayOrScalarLike, y: ArrayOrScalarLike,
-                     gis_order: bool = False,
+                     axis_order: Optional[OAMS_AXIS_ORDER] = None,
                      open_options: Optional[dict] = None,
                      ovr_idx: Optional[int] = None,
                      band_nums: Optional[Sequence[int]] = None,
@@ -109,7 +110,7 @@ def gdallocationinfo(filename_or_ds: PathOrDS,
                 if srs == LocationInfoSRS.SameAsDS_SRS_GeogCS:
                     points_srs = ds_srs.CloneGeogCS()
                 else:
-                    points_srs = get_srs(srs, gis_order=gis_order)
+                    points_srs = get_srs(srs, axis_order=axis_order)
                 ct = get_transform(points_srs, ds_srs)
             x, y, _z = transform_points(ct, x, y)
 
@@ -308,7 +309,7 @@ def get_parser(self, argv) -> GDALArgumentParser:
                                  "instead of the base band. Note that the x,y location (if the coordinate system is "
                                  "pixel/line) must still be given with respect to the base band.")
 
-        parser.add_argument("-axis_order", dest="gis_order", choices=['gis', 'authority'], type=str,
+        parser.add_argument("-axis_order", dest="axis_order", choices=['gis', 'authority'], type=str, default=None,
                             help="X, Y Axis order: Traditional GIS, Authority complaint or otherwise utility default.")
 
         parser.add_argument("-oo", dest="open_options", metavar="NAME=VALUE",
@@ -338,9 +339,20 @@ def augment_kwargs(self, kwargs) -> dict:
         elif kwargs['wgs84']:
             kwargs['srs'] = 4326
 
-        kwargs['gis_order'] = \
-            kwargs['srs'] and (kwargs['srs'] != LocationInfoSRS.SameAsDS_SRS) if kwargs['gis_order'] is None \
-                else str(kwargs['gis_order']).lower() == 'gis'
+        axis_order = kwargs['axis_order']
+        if isinstance(axis_order, OAMS_AXIS_ORDER):
+            pass
+        else:
+            if isinstance(axis_order, str):
+                gis_order = axis_order.lower() == 'gis'
+            elif axis_order is None:
+                gis_order = kwargs['srs'] is not None and (kwargs['srs'] != LocationInfoSRS.SameAsDS_SRS)
+            elif isinstance(axis_order, bool):
+                gis_order = axis_order
+            else:
+                raise Exception(f'Unknown axis order {axis_order}')
+            axis_order = get_axis_order_from_gis_order(gis_order)
+        kwargs['axis_order'] = axis_order
 
         if kwargs['xml']:
             kwargs['output_mode'] = LocationInfoOutput.XML