Merge pull request #3778 from talos-gis/default_gis_order

`osr_util.py` - bugfix and add `default_axis_order` flag; `arrays.py` - add types to remove dependency on `numpy`

autotest/pyscripts/test_gdal_utils.py

@@ -27,7 +27,7 @@
 # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 # DEALINGS IN THE SOFTWARE.
 ###############################################################################
-
+import array
 import os
 from pathlib import Path
 
@@ -38,7 +38,7 @@
 
 pytest.importorskip('osgeo_utils')
 
-from osgeo_utils.auxiliary import util, raster_creation, base
+from osgeo_utils.auxiliary import util, raster_creation, base, array_util
 
 temp_files = []
 
@@ -89,13 +89,13 @@ def test_utils_py_1():
     ds = util.open_ds(filename)
     compression = util.get_image_structure_metadata(filename, 'COMPRESSION')
     assert compression == 'DEFLATE'
-    ovr_count = util.get_ovr_count(ds)+1
+    ovr_count = util.get_ovr_count(ds) + 1
     assert ovr_count == 3
     pixel_size = util.get_pixel_size(ds)
     assert pixel_size == (10, -10)
 
     for i in range(-ovr_count, ovr_count):
-        assert util.get_ovr_idx(filename, ovr_idx=i) == (i if i >= 0 else ovr_count+i)
+        assert util.get_ovr_idx(filename, ovr_idx=i) == (i if i >= 0 else ovr_count + i)
 
     for res, ovr in [(5, 0), (10, 0), (11, 0), (19.99, 0), (20, 1), (20.1, 1), (39, 1), (40, 2), (41, 2), (400, 2)]:
         assert util.get_ovr_idx(filename, ovr_res=res) == ovr
@@ -110,6 +110,32 @@ def test_utils_py_1():
     ds_list = None
 
 
+def test_utils_arrays():
+    scalars = [7, 5.2]
+
+    for scalar in scalars:
+        assert isinstance(scalar, array_util.ScalarLike.__args__)
+        assert isinstance(scalar, array_util.ArrayOrScalarLike.__args__)
+
+    for vec in (scalars, tuple(scalars), array.array('d', scalars), array.array('i', [2, 3])):
+        assert isinstance(vec, array_util.ArrayLike.__args__)
+        assert isinstance(vec, array_util.ArrayOrScalarLike.__args__)
+
+    for not_vec in (None, {1: 2}):
+        assert not isinstance(not_vec, array_util.ArrayLike.__args__)
+        assert not isinstance(not_vec, array_util.ArrayOrScalarLike.__args__)
+
+
+def test_utils_np_arrays():
+    np = pytest.importorskip('numpy')
+    vec_2d = [[1, 2, 3], [4, 5, 6]]
+
+    for dtype in (np.int8, np.int32, np.float64):
+        for vec in (vec_2d[0], vec_2d):
+            arr = np.array(vec, dtype=dtype)
+            assert isinstance(arr, array_util.ArrayLike.__args__)
+
+
 def test_utils_py_cleanup():
     for filename in temp_files:
         try:

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/array_util.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# ******************************************************************************
+#
+#  Project:  GDAL utils.auxiliary
+#  Purpose:  array and scalar related types functions
+#  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 array
+import math
+from typing import Union, Sequence, TYPE_CHECKING
+
+ScalarLike = Union[
+    int,
+    float,
+]
+
+if TYPE_CHECKING:
+    # avoid "TypeError: Subscripted generics cannot be used with class and instance checks" while using isinstance()
+    ArrayLike = Union[
+        ScalarLike,
+        Sequence[ScalarLike],
+        array.array,
+    ]
+else:
+    ArrayLike = Union[
+        ScalarLike,
+        Sequence,
+        array.array,
+    ]
+
+try:
+    from numpy import ndarray
+    ArrayLike = Union[ArrayLike, ndarray]
+except ImportError:
+    pass
+
+ArrayOrScalarLike = Union[ArrayLike, ScalarLike]
+
+
+def array_dist(x: ArrayOrScalarLike, y: ArrayOrScalarLike, is_max: bool = True) -> ScalarLike:
+    if isinstance(x, ScalarLike.__args__):
+        return abs(x-y)
+    try:
+        from osgeo_utils.auxiliary.numpy_util import array_dist as np_array_dist
+        return np_array_dist(x=x, y=y, is_max=is_max)
+    except ImportError as e:
+        return max(abs(a-b) for a, b in zip(x, y)) if is_max else max(abs(a-b) for a, b in zip(x, y))

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

@@ -27,15 +27,11 @@
 #  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 #  DEALINGS IN THE SOFTWARE.
 # ******************************************************************************
-import array
-from numbers import Real
-from typing import Union, Sequence
 
 from osgeo import gdal, gdal_array
 import numpy as np
 
-NumpyCompatibleArray = Union[np.ndarray, array.array, Sequence]
-NumpyCompatibleArrayOrReal = Union[NumpyCompatibleArray, Real]
+from osgeo_utils.auxiliary.array_util import ArrayOrScalarLike, ScalarLike
 
 
 def GDALTypeCodeToNumericTypeCodeEx(buf_type, signed_byte, default=None):
@@ -53,3 +49,14 @@ def GDALTypeCodeAndNumericTypeCodeFromDataSet(ds):
     signed_byte = ds.GetRasterBand(1).GetMetadataItem('PIXELTYPE', 'IMAGE_STRUCTURE') == 'SIGNEDBYTE'
     np_typecode = GDALTypeCodeToNumericTypeCodeEx(buf_type, signed_byte=signed_byte, default=np.float32)
     return buf_type, np_typecode
+
+
+def array_dist(x: ArrayOrScalarLike, y: ArrayOrScalarLike, is_max: bool = True) -> ScalarLike:
+    if isinstance(x, ScalarLike.__args__):
+        return abs(x-y)
+    if not isinstance(x, np.ndarray):
+        x = np.array(x)
+    if not isinstance(y, np.ndarray):
+        y = np.array(y)
+    diff = np.abs(x-y)
+    return np.max(diff) if is_max else np.min(diff)