From ac0aeacf6a1de8fa608a74e754c6a9455cf80d93 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 28 Nov 2024 19:03:48 +0100 Subject: [PATCH] Zarr: fix incorrect DataAxisToSRSAxisMapping for EPSG geographic CRS Resolves issue of https://github.com/OSGeo/gdal/pull/11380 Co-authored-by: Samuel Kogler --- autotest/gdrivers/zarr_driver.py | 72 +++++++++++++++++++++++++------- frmts/zarr/zarr_array.cpp | 10 +++-- 2 files changed, 64 insertions(+), 18 deletions(-) diff --git a/autotest/gdrivers/zarr_driver.py b/autotest/gdrivers/zarr_driver.py index 79217d38542d..caed01d26f2e 100644 --- a/autotest/gdrivers/zarr_driver.py +++ b/autotest/gdrivers/zarr_driver.py @@ -818,15 +818,21 @@ def test_zarr_read_crs(crs_member): assert rg ar = rg.OpenMDArray(rg.GetMDArrayNames()[0]) srs = ar.GetSpatialRef() - if ( - not (osr.GetPROJVersionMajor() > 6 or osr.GetPROJVersionMinor() >= 2) - and crs_member == "projjson" - ): - assert srs is None - else: - assert srs is not None - assert srs.GetAuthorityCode(None) == "4326" - assert len(ar.GetAttributes()) == 0 + assert srs is not None + assert srs.GetAuthorityCode(None) == "4326" + # Mapping is 1, 2 since the slowest varying axis in multidim + # mode is the lines, which matches latitude as the first axis of the CRS. + assert srs.GetDataAxisToSRSAxisMapping() == [1, 2] + assert len(ar.GetAttributes()) == 0 + + # Open as classic CRS + ds = gdal.Open("/vsimem/test.zarr") + srs = ds.GetSpatialRef() + assert srs.GetAuthorityCode(None) == "4326" + # Inverted mapping in classic raster mode compared to multidim mode, + # because the first "axis" in our data model is columns. + assert srs.GetDataAxisToSRSAxisMapping() == [2, 1] + finally: gdal.RmdirRecursive("/vsimem/test.zarr") @@ -3469,7 +3475,8 @@ def set_crs(): assert ar crs = osr.SpatialReference() crs.ImportFromEPSG(4326) - crs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + # lat first + crs.SetDataAxisToSRSAxisMapping([1, 2]) crs.SetCoordinateEpoch(2021.2) assert ar.SetSpatialRef(crs) == gdal.CE_None @@ -3496,7 +3503,7 @@ def check_crs(): crs = ar.GetSpatialRef() assert crs is not None assert crs.GetAuthorityCode(None) == "4326" - assert crs.GetDataAxisToSRSAxisMapping() == [2, 1] + assert crs.GetDataAxisToSRSAxisMapping() == [1, 2] assert crs.GetCoordinateEpoch() == 2021.2 check_crs() @@ -3505,6 +3512,8 @@ def check_crs_classic_dataset(): ds = gdal.Open("/vsimem/test.zarr") crs = ds.GetSpatialRef() assert crs is not None + assert crs.GetAuthorityCode(None) == "4326" + assert crs.GetDataAxisToSRSAxisMapping() == [2, 1] check_crs_classic_dataset() @@ -5505,10 +5514,9 @@ def test_zarr_read_cf1(): ds = gdal.Open("data/zarr/byte_cf1.zarr") assert ds - assert ( - ds.GetSpatialRef().ExportToProj4() - == "+proj=utm +zone=11 +ellps=clrk66 +units=m +no_defs" - ) + srs = ds.GetSpatialRef() + assert srs.ExportToProj4() == "+proj=utm +zone=11 +ellps=clrk66 +units=m +no_defs" + assert srs.GetDataAxisToSRSAxisMapping() == [1, 2] ############################################################################### @@ -5525,6 +5533,40 @@ def test_zarr_read_cf1_zarrv3(): ) +############################################################################### + + +@gdaltest.enable_exceptions() +@pytest.mark.require_proj(9) +def test_zarr_write_WGS84_and_EGM96_height(tmp_vsimem): + + tmp_filename = str(tmp_vsimem / "out.zarr") + with gdal.GetDriverByName("Zarr").Create(tmp_filename, 1, 1) as ds: + srs = osr.SpatialReference() + srs.ImportFromEPSG(9707) + ds.SetSpatialRef(srs) + with gdal.Open(tmp_filename) as ds: + srs = ds.GetSpatialRef() + assert srs.GetAuthorityCode(None) == "9707" + assert srs.GetDataAxisToSRSAxisMapping() == [2, 1] + + +############################################################################### + + +@gdaltest.enable_exceptions() +def test_zarr_write_UTM31N_and_EGM96_height(tmp_vsimem): + + tmp_filename = str(tmp_vsimem / "out.zarr") + with gdal.GetDriverByName("Zarr").Create(tmp_filename, 1, 1) as ds: + srs = osr.SpatialReference() + srs.SetFromUserInput("EPSG:32631+5773") + ds.SetSpatialRef(srs) + with gdal.Open(tmp_filename) as ds: + srs = ds.GetSpatialRef() + assert srs.GetDataAxisToSRSAxisMapping() == [1, 2] + + ############################################################################### # Test bug fix for https://github.com/OSGeo/gdal/issues/11016 diff --git a/frmts/zarr/zarr_array.cpp b/frmts/zarr/zarr_array.cpp index 6d12705201f5..5f038d8945db 100644 --- a/frmts/zarr/zarr_array.cpp +++ b/frmts/zarr/zarr_array.cpp @@ -2724,6 +2724,7 @@ void ZarrArray::ParseSpecialAttributes( if (item.IsValid()) { poSRS = std::make_shared(); + poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); if (poSRS->SetFromUserInput( item.ToString().c_str(), OGRSpatialReference:: @@ -2748,6 +2749,7 @@ void ZarrArray::ParseSpecialAttributes( if (gridMappingArray) { poSRS = std::make_shared(); + poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); CPLStringList aosKeyValues; for (const auto &poAttr : gridMappingArray->GetAttributes()) { @@ -2798,10 +2800,12 @@ void ZarrArray::ParseSpecialAttributes( } if (iDimX > 0 && iDimY > 0) { - if (poSRS->GetDataAxisToSRSAxisMapping() == std::vector{2, 1}) + const auto &oMapping = poSRS->GetDataAxisToSRSAxisMapping(); + if (oMapping == std::vector{2, 1} || + oMapping == std::vector{2, 1, 3}) poSRS->SetDataAxisToSRSAxisMapping({iDimY, iDimX}); - else if (poSRS->GetDataAxisToSRSAxisMapping() == - std::vector{1, 2}) + else if (oMapping == std::vector{1, 2} || + oMapping == std::vector{1, 2, 3}) poSRS->SetDataAxisToSRSAxisMapping({iDimX, iDimY}); }