netCDF multidim: report correct axis order when reading SRS with EPSG…

… geographic CRS + geoid model

autotest/gdrivers/netcdf_multidim.py

@@ -4156,3 +4156,30 @@ def check_metadata():
 
     assert os.path.exists(pam_filename)
     os.unlink(pam_filename)
+
+
+###############################################################################
+
+
+@gdaltest.enable_exceptions()
+@pytest.mark.require_proj(9)
+def test_netcdf_multidim_WGS84_and_EGM96_height(tmp_path):
+
+    tmp_filename = str(tmp_path / "out.nc")
+    with gdal.GetDriverByName("netCDF").Create(tmp_filename, 3, 3) as ds:
+        srs = osr.SpatialReference()
+        srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+        srs.ImportFromEPSG(9707)
+        ds.SetSpatialRef(srs)
+        ds.SetGeoTransform([2, 1, 0, 49, 0, -1])
+    with gdal.Open(tmp_filename) as ds:
+        srs = ds.GetSpatialRef()
+        assert srs.GetAuthorityCode(None) == "9707"
+        # We currently report a 3rd axis, which is dubious
+        assert srs.GetDataAxisToSRSAxisMapping()[0:2] == [2, 1]
+    with gdal.OpenEx(tmp_filename, gdal.OF_MULTIDIM_RASTER) as ds:
+        rg = ds.GetRootGroup()
+        ar = rg.OpenMDArray("Band1")
+        srs = ar.GetSpatialRef()
+        assert srs.GetAuthorityCode(None) == "9707"
+        assert srs.GetDataAxisToSRSAxisMapping() == [1, 2]

frmts/netcdf/netcdfmultidim.cpp

@@ -2646,8 +2646,9 @@ std::shared_ptr<OGRSpatialReference> netCDFVariable::GetSpatialRef() const
         m_poSRS.reset(poSRS->Clone());
         if (iDimX > 0 && iDimY > 0)
         {
-            if (m_poSRS->GetDataAxisToSRSAxisMapping() ==
-                std::vector<int>{2, 1})
+            const auto &oMapping = m_poSRS->GetDataAxisToSRSAxisMapping();
+            if (oMapping == std::vector<int>{2, 1} ||
+                oMapping == std::vector<int>{2, 1, 3})
                 m_poSRS->SetDataAxisToSRSAxisMapping({iDimY, iDimX});
             else
                 m_poSRS->SetDataAxisToSRSAxisMapping({iDimX, iDimY});