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

[idanmiara]

What does this PR do?

This PR solves two issues, I wanted to separate them but apparently for writing proper tests that don't depend on numpy I needed to merge them.

I've added many tests for to make it right this time...

osr_util.py depends on numpy. this PR remove this dependency.

osgeo_utils/auxiliary/array_util.py - add ArrayLike, ScalarLike, ArrayOrScalarLike types (replaces NumpyCompatibleArray, NumpyCompatibleArrayOrReal)

osr_util.py, gdallocationinfo.py, numpy_util.py - use new types
autotest/pyscripts/test_gdal_utils.py - test new types

osr_util.py - add default_axis_order flag to allow the user to change the default
some applications use gis_order, this PR will allow them to override the default if they wish by adding the following code to their package/module __init__.py

from osgeo.osr import OAMS_TRADITIONAL_GIS_ORDER
from osgeo_utils.auxiliary.osr_util import set_default_axis_order
set_default_axis_order(OAMS_TRADITIONAL_GIS_ORDER)

Because get_srs is used in many inner functions, it would be very hard to propagate the gis_order default without this PR.

this PR also fixes a bug in get_srs() for the following case:

    pj4326_gis = osr_util.get_srs(4326, gis_order=True)
    assert pj4326_gis.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER

    # check that srs object is also affected if explicitly set
    srs = osr_util.get_srs(pj4326_gis, gis_order=False)
    assert srs.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT

Also fixing the issue that calling get_srs() with an osr.SpatialReference object modified the input object.
With this PR it creates a copy and returns it.

Related PR

python/typing#593
#3780

[idanmiara]

@rouault I can't figure out why the following 3 tests fail.
https://github.com/OSGeo/gdal/pull/3778/files#diff-9f53635932ec3d0256fdf3dba9f904de01f9bdb5aca0ddd8d8f805f199e41e40R148
I've narrowed down the issue as written in the comment

[rouault]

I've narrowed down the issue as written in the comment

Hum, this is quite subtle, so let's see at:

$ python
>>> from osgeo import osr
>>> sr = osr.SpatialReference()
>>> sr.ImportFromEPSG(4326)
>>> sr2 = osr.SpatialReference()
>>> sr2.SetFromUserInput('+proj=longlat +datum=WGS84 +no_defs')

>>> sr.ExportToWkt()
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'

>>> sr2.ExportToWkt()
'GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]'

>>> sr.GetDataAxisToSRSAxisMapping()
[1, 2]
>>> sr2.GetDataAxisToSRSAxisMapping()
[1, 2]


>>> sr.IsSame(sr2)
1

Why IsSame() returns 1 whereas the axis order in the proj4 case is long/lat and in the importFromEPSG() is lat/long ? The reason is that the IsSame() behaviour by default is CRITERION=EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS, that is it ignores axis order difference for GeographicCRS. The reason for that was probably to diminish the level of backward incompatibility when we switched to PROJ 6 during the GDAL 2 --> GDAL 3 transition.

But if now I do

>>> sr.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
>>> sr.GetDataAxisToSRSAxisMapping()
[2, 1]
>>> sr.IsSame(sr2)
0

they are no longer consider equivalent for IsSame(). The reason is that the first test in IsSame() checks that the GetDataAxisToSRSAxisMapping() of the 2 CRS is the same.

Admitedly, this is quite confusing. Perhaps in the CRITERION=EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS case when we compare 2 Geographic CRS, we should also ignore the GetDataAxisToSRSAxisMapping() setting, but there might be some unforeseen impacts in doing that

[idanmiara]

I've narrowed down the issue as written in the comment

Hum, this is quite subtle, so let's see at:

$ python
>>> from osgeo import osr
>>> sr = osr.SpatialReference()
>>> sr.ImportFromEPSG(4326)
>>> sr2 = osr.SpatialReference()
>>> sr2.SetFromUserInput('+proj=longlat +datum=WGS84 +no_defs')

>>> sr.ExportToWkt()
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'

>>> sr2.ExportToWkt()
'GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]'

>>> sr.GetDataAxisToSRSAxisMapping()
[1, 2]
>>> sr2.GetDataAxisToSRSAxisMapping()
[1, 2]


>>> sr.IsSame(sr2)
1

Why IsSame() returns 1 whereas the axis order in the proj4 case is long/lat and in the importFromEPSG() is lat/long ? The reason is that the IsSame() behaviour by default is CRITERION=EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS, that is it ignores axis order difference for GeographicCRS. The reason for that was probably to diminish the level of backward incompatibility when we switched to PROJ 6 during the GDAL 2 --> GDAL 3 transition.

But if now I do

>>> sr.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
>>> sr.GetDataAxisToSRSAxisMapping()
[2, 1]
>>> sr.IsSame(sr2)
0

they are no longer consider equivalent for IsSame(). The reason is that the first test in IsSame() checks that the GetDataAxisToSRSAxisMapping() of the 2 CRS is the same.

Admitedly, this is quite confusing. Perhaps in the CRITERION=EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS case when we compare 2 Geographic CRS, we should also ignore the GetDataAxisToSRSAxisMapping() setting, but there might be some unforeseen impacts in doing that

Thanks for looking into it!
I think I'll better comment out these tests as it's out of the scope of this PR.
We can try to tackle it in another PR...

[idanmiara]

Just thinking about how to solve this issue (outside of this PR)...

So if I got it strait you mean that sr.ImportFromEPSG(4326) and sr2.SetFromUserInput('+proj=longlat +datum=WGS84 +no_defs') don't have the same axis order but IsSame ignores that thus outputs equality.
After applying sr.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) they do have the same axis order, but IsSame doesn't pick that up correctly thus outputs inequality although they do equal.

I think that we could take the following approach to make this right, but is this an issue of PROJ or GDAL?
I see that this flag is mentioned in PROJ

Anyway, if it's a GDAL issue then maybe we can solve it as follows:
EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS : as you said, it should ignore the axis mapping altogether for gogcrs
(for non-gogcrs this flag should have no effect - as if the criteria is EQUIVALENT for the sake of checking the axis order).
EQUIVALENT: we need to check if the axis order after applying the mapping is the same (which it is in the case for the case above)
STRICT: should say that they are not equal as the original order is different.
I that we should keep the EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS default for keep backwards compatibility.

[rouault]

I think that we could take the following approach to make this right, but is this an issue of PROJ or GDAL?

in GDAL, in the implementation of OGRSpatialReference::IsSame()

Anyway, if it's a GDAL issue then maybe we can solve it as follows:

That sounds like a good plan.

[idanmiara]

Any other comments about this PR?
Do you think I should leave both input options for def get_srs(srs: AnySRS, gis_order: Optional[bool] = None, axis_order: Optional[OAMS_AXIS_ORDER] = None) ? or should I remove gis_order?

[rouault]

or should I remove gis_order?

I'd go for that. I see it is currently used in the 3.3 branch but we didn't make any clear statements about what is considered to be in the API or implementation details

[idanmiara]

or should I remove gis_order?

I'd go for that. I see it is currently used in the 3.3 branch but we didn't make any clear statements about what is considered to be in the API or implementation details

Right, I was thinking how far should I go to keep backwards compatibility inside the auxiliary submodule, as it does adds clutter...

[rouault]

Right, I was thinking how far should I go to keep backwards compatibility inside the auxiliary submodule

It might be useful to have a clear statement about what is in the API and what is not at some point. And what is in the API should ideally come with some not-only-contained-in-py-files documentation.

[idanmiara]

Right, I was thinking how far should I go to keep backwards compatibility inside the auxiliary submodule

It might be useful to have a clear statement about what is in the API and what is not at some point. And what is in the API should ideally come with some not-only-contained-in-py-files documentation.

Right, I'll try to formalize something like that.

[idanmiara]

or should I remove gis_order?

I'd go for that. I see it is currently used in the 3.3 branch but we didn't make any clear statements about what is considered to be in the API or implementation details

Ok, done, I think this PR is ready now.

[idanmiara]

Thanks @rouault!
I've labeld it as backport 3.3, How does it work now? I can't find that the bot has done anything related.

[rouault]

I can't find that the bot has done anything related.

the backport bot is a regular github action, and it follows the queuing logic of github actions jobs. So a bit of patience: it is queued at https://github.com/OSGeo/gdal/actions/workflows/backport.yml

[github-actions]

The backport to release/3.3 failed:

The process '/usr/bin/git' failed with exit code 1

To backport manually, run these commands in your terminal:

# Fetch latest updates from GitHub
git fetch
# Create a new working tree
git worktree add .worktrees/backport-release/3.3 release/3.3
# Navigate to the new working tree
cd .worktrees/backport-release/3.3
# Create a new branch
git switch --create backport-3778-to-release/3.3
# Cherry-pick the merged commit of this pull request and resolve the conflicts
git cherry-pick --mainline 1 7746b3c4af0007509af58f25c9404109e4351c87
# Push it to GitHub
git push --set-upstream origin backport-3778-to-release/3.3
# Go back to the original working tree
cd ../..
# Delete the working tree
git worktree remove .worktrees/backport-release/3.3

Then, create a pull request where the base branch is release/3.3 and the compare/head branch is backport-3778-to-release/3.3.

[rouault]

ok, so the automatic backport failed. @idanmiara I let you create a backport PR that cherry-picks the relevant commits

[idanmiara]

ok, so the automatic backport failed. @idanmiara I let you create a backport PR that cherry-picks the relevant commits

Thanks.
#3793