Multidim GDALGroup::CopyFrom(): modify logic to copy first indexing v…

…ariables. Fixes issue with gdalmdimtranslate from Zarr to netCDF

autotest/gdrivers/netcdf_multidim.py

@@ -3321,3 +3321,65 @@ def reopen():
         reopen()
     finally:
         gdal.Unlink(filename)
+
+
+###############################################################################
+
+
+def test_netcdf_multidim_copy_group_with_indexing_variable_after_regular_var():
+
+    outfilename = "tmp/out.nc"
+
+    try:
+
+        def create():
+            ds = gdal.GetDriverByName("MEM").CreateMultiDimensional("")
+            rg = ds.GetRootGroup()
+            dim_y = rg.CreateDimension("y", gdal.DIM_TYPE_HORIZONTAL_Y, None, 2)
+            dim_x = rg.CreateDimension("x", gdal.DIM_TYPE_HORIZONTAL_X, None, 2)
+            var = rg.CreateMDArray(
+                "var", [dim_y, dim_x], gdal.ExtendedDataType.Create(gdal.GDT_Int16)
+            )
+
+            srs = osr.SpatialReference()
+            srs.ImportFromEPSG(32631)
+            assert var.SetSpatialRef(srs) == gdal.CE_None
+
+            assert var.Write(struct.pack("h" * 4, 1, 2, 3, 4)) == gdal.CE_None
+
+            x = rg.CreateMDArray(
+                "x", [dim_x], gdal.ExtendedDataType.Create(gdal.GDT_Float64)
+            )
+            dim_x.SetIndexingVariable(x)
+            assert x.Write(struct.pack("d" * 2, 1, 2)) == gdal.CE_None
+
+            y = rg.CreateMDArray(
+                "y", [dim_y], gdal.ExtendedDataType.Create(gdal.GDT_Float64)
+            )
+            dim_y.SetIndexingVariable(y)
+            assert y.Write(struct.pack("d" * 2, 1, 2)) == gdal.CE_None
+
+            return ds
+
+        def copy(src_ds):
+            rg = src_ds.GetRootGroup()
+            assert rg.GetMDArrayNames() == ["var", "x", "y"]
+
+            gdal.ErrorReset()
+            assert gdal.MultiDimTranslate(outfilename, src_ds)
+            assert gdal.GetLastErrorMsg() == ""
+            src_ds = None
+
+            out_ds = gdal.OpenEx(outfilename, gdal.OF_MULTIDIM_RASTER)
+            rg = out_ds.GetRootGroup()
+            ar = rg.OpenMDArray("var")
+            dim_y = ar.GetDimensions()[0]
+            assert dim_y.GetType() == gdal.DIM_TYPE_HORIZONTAL_Y
+            dim_x = ar.GetDimensions()[1]
+            assert dim_x.GetType() == gdal.DIM_TYPE_HORIZONTAL_X
+
+        copy(create())
+
+    finally:
+        if os.path.exists(outfilename):
+            gdal.Unlink(outfilename)

gcore/gdalmultidim.cpp

@@ -763,12 +763,12 @@ bool GDALGroup::CopyFrom(const std::shared_ptr<GDALGroup> &poDstRootGroup,
                 return false;
         }
 
-        auto arrayNames = poSrcGroup->GetMDArrayNames();
-        for (const auto &name : arrayNames)
+        const auto CopyArray =
+            [this, &poSrcDS, &poDstRootGroup, &mapExistingDstDims,
+             &mapSrcVariableNameToIndexedDimName, pfnProgress, pProgressData,
+             papszOptions, bStrict, &nCurCost,
+             nTotalCost](const std::shared_ptr<GDALMDArray> &srcArray)
         {
-            auto srcArray = poSrcGroup->OpenMDArray(name);
-            EXIT_OR_CONTINUE_IF_NULL(srcArray);
-
             // Map source dimensions to target dimensions
             std::vector<std::shared_ptr<GDALDimension>> dstArrayDims;
             const auto &srcArrayDims(srcArray->GetDimensions());
@@ -797,8 +797,8 @@ bool GDALGroup::CopyFrom(const std::shared_ptr<GDALGroup> &poDstRootGroup,
                         }
                         else
                         {
-                            std::string newDimNamePrefix(name + '_' +
-                                                         dim->GetName());
+                            std::string newDimNamePrefix(srcArray->GetName() +
+                                                         '_' + dim->GetName());
                             newDimName = newDimNamePrefix;
                             int nIterCount = 2;
                             while (mapExistingDstDims.find(newDimName) !=
@@ -975,7 +975,8 @@ bool GDALGroup::CopyFrom(const std::shared_ptr<GDALGroup> &poDstRootGroup,
                     CreateMDArray(srcArray->GetName(), dstArrayDims,
                                   GDALExtendedDataType::Create(eAutoScaleType),
                                   aosArrayCO.List());
-                EXIT_OR_CONTINUE_IF_NULL(dstArray);
+                if (!dstArray)
+                    return !bStrict;
 
                 if (srcArray->GetRawNoDataValue() != nullptr)
                 {
@@ -1020,7 +1021,8 @@ bool GDALGroup::CopyFrom(const std::shared_ptr<GDALGroup> &poDstRootGroup,
             {
                 dstArray = CreateMDArray(srcArray->GetName(), dstArrayDims,
                                          srcArrayType, aosArrayCO.List());
-                EXIT_OR_CONTINUE_IF_NULL(dstArray);
+                if (!dstArray)
+                    return !bStrict;
 
                 if (!dstArray->CopyFrom(poSrcDS, srcArray.get(), bStrict,
                                         nCurCost, nTotalCost, pfnProgress,
@@ -1044,9 +1046,43 @@ bool GDALGroup::CopyFrom(const std::shared_ptr<GDALGroup> &poDstRootGroup,
                         dstArray);
                 }
             }
+
+            return true;
+        };
+
+        const auto arrayNames = poSrcGroup->GetMDArrayNames();
+
+        // Start by copying arrays that are indexing variables of dimensions
+        for (const auto &name : arrayNames)
+        {
+            auto srcArray = poSrcGroup->OpenMDArray(name);
+            EXIT_OR_CONTINUE_IF_NULL(srcArray);
+
+            const auto oIterDimName =
+                mapSrcVariableNameToIndexedDimName.find(srcArray->GetName());
+            if (oIterDimName != mapSrcVariableNameToIndexedDimName.end())
+            {
+                if (!CopyArray(srcArray))
+                    return false;
+            }
+        }
+
+        // Then copy regular arrays
+        for (const auto &name : arrayNames)
+        {
+            auto srcArray = poSrcGroup->OpenMDArray(name);
+            EXIT_OR_CONTINUE_IF_NULL(srcArray);
+
+            const auto oIterDimName =
+                mapSrcVariableNameToIndexedDimName.find(srcArray->GetName());
+            if (oIterDimName == mapSrcVariableNameToIndexedDimName.end())
+            {
+                if (!CopyArray(srcArray))
+                    return false;
+            }
         }
 
-        auto groupNames = poSrcGroup->GetGroupNames();
+        const auto groupNames = poSrcGroup->GetGroupNames();
         for (const auto &name : groupNames)
         {
             auto srcSubGroup = poSrcGroup->OpenGroup(name);