Zarr: fix incorrect DataAxisToSRSAxisMapping for EPSG geographic CRS

Resolves issue of #11380

Co-authored-by: Samuel Kogler <[email protected]>

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
 

frmts/zarr/zarr_array.cpp

@@ -2724,6 +2724,7 @@ void ZarrArray::ParseSpecialAttributes(
             if (item.IsValid())
             {
                 poSRS = std::make_shared<OGRSpatialReference>();
+                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<OGRSpatialReference>();
+                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<int>{2, 1})
+            const auto &oMapping = poSRS->GetDataAxisToSRSAxisMapping();
+            if (oMapping == std::vector<int>{2, 1} ||
+                oMapping == std::vector<int>{2, 1, 3})
                 poSRS->SetDataAxisToSRSAxisMapping({iDimY, iDimX});
-            else if (poSRS->GetDataAxisToSRSAxisMapping() ==
-                     std::vector<int>{1, 2})
+            else if (oMapping == std::vector<int>{1, 2} ||
+                     oMapping == std::vector<int>{1, 2, 3})
                 poSRS->SetDataAxisToSRSAxisMapping({iDimX, iDimY});
         }