From a20194b3ded9edc5a23e05d30fbe3c5b4b070dbe Mon Sep 17 00:00:00 2001 From: Idan Miara Date: Sun, 2 May 2021 21:32:27 +0300 Subject: [PATCH 1/2] gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/array_util.py - add ArrayLike, ScalarLike, ArrayOrScalarLike types (to be used instead of NumpyCompatibleArray, NumpyCompatibleArrayOrReal), removes numpy dependency osr_util.py, gdallocationinfo.py, numpy_util.py - use new types autotest/pyscripts/test_gdal_utils.py - test new types --- autotest/pyscripts/test_gdal_utils.py | 34 +++++++-- .../osgeo_utils/auxiliary/array_util.py | 69 +++++++++++++++++++ .../osgeo_utils/auxiliary/numpy_util.py | 17 +++-- .../osgeo_utils/auxiliary/osr_util.py | 6 +- .../osgeo_utils/samples/gdallocationinfo.py | 12 ++-- 5 files changed, 120 insertions(+), 18 deletions(-) create mode 100644 gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/array_util.py diff --git a/autotest/pyscripts/test_gdal_utils.py b/autotest/pyscripts/test_gdal_utils.py index 431fd44ab752..60bcc245d498 100644 --- a/autotest/pyscripts/test_gdal_utils.py +++ b/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: diff --git a/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/array_util.py b/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/array_util.py new file mode 100644 index 000000000000..ed8ba8be5608 --- /dev/null +++ b/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 +# +# ****************************************************************************** +# Copyright (c) 2021, Idan Miara +# +# 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)) diff --git a/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/numpy_util.py b/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/numpy_util.py index 765025c43a48..e2b6097ec0a0 100644 --- a/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/numpy_util.py +++ b/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) diff --git a/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py b/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py index 4b311d4250a3..82f373acddae 100644 --- a/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py +++ b/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py @@ -32,7 +32,7 @@ import osgeo from osgeo import osr, ogr, gdal -from osgeo_utils.auxiliary.numpy_util import NumpyCompatibleArray +from osgeo_utils.auxiliary.array_util import ArrayLike AnySRS = Union[str, int, osr.SpatialReference, gdal.Dataset] @@ -78,8 +78,8 @@ def get_transform(src_srs: AnySRS, tgt_srs: AnySRS) -> Optional[osr.CoordinateTr def transform_points(ct: Optional[osr.CoordinateTransformation], - x: NumpyCompatibleArray, y: NumpyCompatibleArray, z: Optional[NumpyCompatibleArray] = None) -> \ - Tuple[NumpyCompatibleArray, NumpyCompatibleArray, Optional[NumpyCompatibleArray]]: + x: ArrayLike, y: ArrayLike, z: Optional[ArrayLike] = None) -> \ + Tuple[ArrayLike, ArrayLike, Optional[ArrayLike]]: if ct is not None: if z is None: for idx, (x0, y0) in enumerate(zip(x, y)): diff --git a/gdal/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py b/gdal/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py index 103f61d7dc60..3e3eafbd4ad1 100644 --- a/gdal/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py +++ b/gdal/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py @@ -40,8 +40,8 @@ from osgeo import gdalconst, osr, gdal from osgeo_utils.auxiliary.base import is_path_like -from osgeo_utils.auxiliary.numpy_util import GDALTypeCodeAndNumericTypeCodeFromDataSet, NumpyCompatibleArrayOrReal, \ - NumpyCompatibleArray +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.util import PathOrDS, open_ds, get_bands, get_scales_and_offsets, get_band_nums from osgeo_utils.auxiliary.gdal_argparse import GDALArgumentParser, GDALScript @@ -67,7 +67,7 @@ class LocationInfoOutput(Enum): def gdallocationinfo(filename_or_ds: PathOrDS, - x: NumpyCompatibleArrayOrReal, y: NumpyCompatibleArrayOrReal, + x: ArrayOrScalarLike, y: ArrayOrScalarLike, gis_order: bool = False, open_options: Optional[dict] = None, ovr_idx: Optional[int] = None, @@ -82,9 +82,9 @@ def gdallocationinfo(filename_or_ds: PathOrDS, filename = filename_or_ds if is_path_like(filename_or_ds) else '' if ds is None: raise Exception(f'Could not open {filename}.') - if not isinstance(x, NumpyCompatibleArray.__args__): + if not isinstance(x, ArrayLike.__args__): x = [x] - if not isinstance(y, NumpyCompatibleArray.__args__): + if not isinstance(y, ArrayLike.__args__): y = [y] if len(x) != len(y): raise Exception(f'len(x)={len(x)} should be the same as len(y)={len(y)}') @@ -161,7 +161,7 @@ def gdallocationinfo(filename_or_ds: PathOrDS, def gdallocationinfo_util(filename_or_ds: PathOrDS, - x: NumpyCompatibleArrayOrReal, y: NumpyCompatibleArrayOrReal, + x: ArrayOrScalarLike, y: ArrayOrScalarLike, open_options: Optional[dict] = None, band_nums: Optional[Sequence[int]] = None, resample_alg=gdalconst.GRIORA_NearestNeighbour, From d07f99bb2d27dd9fb82d10f076c7a70d3274b887 Mon Sep 17 00:00:00 2001 From: Idan Miara Date: Sun, 2 May 2021 14:27:23 +0300 Subject: [PATCH 2/2] gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py - get_srs() - 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 | 164 ++++++++++++++++++ .../osgeo_utils/auxiliary/osr_util.py | 93 +++++++--- .../osgeo_utils/samples/gdallocationinfo.py | 26 ++- 3 files changed, 253 insertions(+), 30 deletions(-) create mode 100644 autotest/pyscripts/test_osr_util.py diff --git a/autotest/pyscripts/test_osr_util.py b/autotest/pyscripts/test_osr_util.py new file mode 100644 index 000000000000..2d4a785ec85a --- /dev/null +++ b/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 +# +############################################################################### +# Copyright (c) 2021, Idan Miara +# +# 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 + diff --git a/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py b/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py index 82f373acddae..59d3678905d2 100644 --- a/gdal/swig/python/gdal-utils/osgeo_utils/auxiliary/osr_util.py +++ b/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]: diff --git a/gdal/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py b/gdal/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py index 3e3eafbd4ad1..aef403f23808 100644 --- a/gdal/swig/python/gdal-utils/osgeo_utils/samples/gdallocationinfo.py +++ b/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