diff --git a/autotest/gcore/multidim.py b/autotest/gcore/multidim.py
index f16a34f286b5..eb40ea755643 100644
--- a/autotest/gcore/multidim.py
+++ b/autotest/gcore/multidim.py
@@ -811,3 +811,160 @@ def test_multidim_asclassicsubdataset_band_metadata():
         "OTHER_STRING_ATTR": "string_attr_value",
         "AUX_VAR": "bar",
     }
+
+
+@gdaltest.enable_exceptions()
+def test_multidim_SubsetDimensionFromSelection():
+
+    drv = gdal.GetDriverByName("MEM")
+    mem_ds = drv.CreateMultiDimensional("myds")
+    rg = mem_ds.GetRootGroup()
+
+    numeric_attr = rg.CreateAttribute(
+        "numeric_attr", [], gdal.ExtendedDataType.Create(gdal.GDT_Float64)
+    )
+    assert numeric_attr.Write(1) == gdal.CE_None
+
+    subgroup = rg.CreateGroup("band_properties")
+
+    dimBand = rg.CreateDimension("band", None, None, 3)
+    validity = subgroup.CreateMDArray(
+        "validity", [dimBand], gdal.ExtendedDataType.Create(gdal.GDT_Int32)
+    )
+    validity.Write(array.array("i", [1, 0, 1]))
+    dimBand.SetIndexingVariable(validity)
+
+    dimX = rg.CreateDimension("X", None, None, 2)
+    X = rg.CreateMDArray("X", [dimX], gdal.ExtendedDataType.Create(gdal.GDT_Float64))
+    X.Write(array.array("d", [10, 20]))
+    dimX.SetIndexingVariable(X)
+
+    dimY = rg.CreateDimension("Y", None, None, 2)
+
+    ar = rg.CreateMDArray(
+        "ar", [dimBand, dimY, dimX], gdal.ExtendedDataType.Create(gdal.GDT_Float64)
+    )
+    assert ar.Write(array.array("d", [i for i in range(3 * 2 * 2)])) == gdal.CE_None
+    ar.SetUnit("foo")
+    ar.SetNoDataValueDouble(5)
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(4326)
+    ar.SetSpatialRef(srs)
+
+    ar_numeric_attr = ar.CreateAttribute(
+        "ar_numeric_attr", [], gdal.ExtendedDataType.Create(gdal.GDT_Float64)
+    )
+    assert ar_numeric_attr.Write(1) == gdal.CE_None
+
+    subgroup = rg.CreateGroup("subgroup")
+
+    indexed_twice_by_band = subgroup.CreateMDArray(
+        "indexed_twice_by_band",
+        [dimBand, dimX, dimBand],
+        gdal.ExtendedDataType.Create(gdal.GDT_Float64),
+    )
+    assert (
+        indexed_twice_by_band.Write(array.array("d", [i for i in range(3 * 3 * 2)]))
+        == gdal.CE_None
+    )
+
+    rg.CreateMDArray("non_numeric_ar", [dimBand], gdal.ExtendedDataType.CreateString())
+
+    too_large_dim = rg.CreateDimension(
+        "too_large_dim", None, None, 10 * 1024 * 1024 + 1
+    )
+    rg.CreateMDArray(
+        "too_large_dim_ar", [too_large_dim], gdal.ExtendedDataType.Create(gdal.GDT_Byte)
+    )
+
+    same_value = rg.CreateMDArray(
+        "same_value", [dimX], gdal.ExtendedDataType.Create(gdal.GDT_Float64)
+    )
+    same_value.Write(array.array("d", [10, 10]))
+
+    with pytest.raises(Exception, match="Invalid value for selection"):
+        rg.SubsetDimensionFromSelection("")
+    with pytest.raises(Exception, match="Invalid value for selection"):
+        rg.SubsetDimensionFromSelection("/band_properties/validity")
+    with pytest.raises(Exception, match="Non-numeric value in selection criterion"):
+        rg.SubsetDimensionFromSelection("/band_properties/validity='foo'")
+    with pytest.raises(Exception, match="Cannot find array /i/do/not/exist"):
+        rg.SubsetDimensionFromSelection("/i/do/not/exist=1")
+    with pytest.raises(Exception, match="Array /ar is not single dimensional"):
+        rg.SubsetDimensionFromSelection("/ar=1")
+    with pytest.raises(Exception, match="Array /non_numeric_ar is not of numeric type"):
+        rg.SubsetDimensionFromSelection("/non_numeric_ar=1")
+    with pytest.raises(Exception, match="Too many values in /too_large_dim_ar"):
+        rg.SubsetDimensionFromSelection("/too_large_dim_ar=1")
+    with pytest.raises(
+        Exception, match="No value in /band_properties/validity matching 2.000000"
+    ):
+        rg.SubsetDimensionFromSelection("/band_properties/validity=2")
+
+    # In C++, should return the same object as rg, but we can't easily check
+    rg.SubsetDimensionFromSelection("/same_value=10")
+
+    rg_subset = rg.SubsetDimensionFromSelection("/band_properties/validity=1")
+    assert rg_subset
+    dims = rg_subset.GetDimensions()
+    assert len(dims) == 4
+    bandFound = False
+    for dim in dims:
+        if dim.GetName() == "band":
+            bandFound = True
+            assert dim.GetFullName() == "/band"
+            assert dim.GetSize() == 2
+            var = dim.GetIndexingVariable()
+            assert var
+            assert var.GetName() == "validity"
+            assert var.GetFullName() == "/band_properties/validity"
+    assert bandFound
+
+    assert rg_subset.GetGroupNames() == rg.GetGroupNames()
+    with pytest.raises(Exception, match="Group i_do_not_exist does not exist"):
+        assert rg_subset.OpenGroup("i_do_not_exist") is None
+    assert len(rg_subset.GetAttributes()) == 1
+    assert rg_subset.GetAttribute("numeric_attr") is not None
+
+    subgroup_subset = rg_subset.OpenGroup("band_properties")
+    assert subgroup_subset.GetFullName() == "/band_properties"
+    assert subgroup_subset.GetMDArrayNames() == ["validity"]
+    validity_subset = subgroup_subset.OpenMDArray("validity")
+    assert validity_subset.GetFullName() == "/band_properties/validity"
+    assert validity_subset.Read() == array.array("i", [1, 1])
+
+    assert rg_subset.GetMDArrayNames() == rg.GetMDArrayNames()
+    with pytest.raises(Exception, match="Array i_do_not_exist does not exist"):
+        assert rg_subset.OpenMDArray("i_do_not_exist") is None
+    ar_subset = rg_subset.OpenMDArray("ar")
+    assert ar_subset.GetFullName() == "/ar"
+    assert [dim.GetFullName() for dim in ar_subset.GetDimensions()] == [
+        "/band",
+        "/Y",
+        "/X",
+    ]
+    assert ar_subset.GetDataType() == ar.GetDataType()
+    assert ar_subset.Read() == array.array(
+        "d", [0.0, 1.0, 2.0, 3.0, 8.0, 9.0, 10.0, 11.0]
+    )
+    assert ar_subset.Read(
+        array_start_idx=[1, 0, 0], array_step=[-1, 1, 1]
+    ) == array.array("d", [8.0, 9.0, 10.0, 11.0, 0.0, 1.0, 2.0, 3.0])
+    assert len(ar_subset.GetAttributes()) == 1
+    assert ar_subset.GetAttribute("ar_numeric_attr") is not None
+    assert ar_subset.GetUnit() == ar.GetUnit()
+    assert ar_subset.GetSpatialRef().IsSame(ar.GetSpatialRef())
+    assert ar_subset.GetNoDataValueAsRaw() == ar.GetNoDataValueAsRaw()
+    assert ar_subset.GetBlockSize() == [0, 1, 0]
+
+    subgroup_subset = rg_subset.OpenGroup("subgroup")
+    indexed_twice_by_band_subset = subgroup_subset.OpenMDArray("indexed_twice_by_band")
+    assert (
+        indexed_twice_by_band_subset.GetFullName() == "/subgroup/indexed_twice_by_band"
+    )
+    assert indexed_twice_by_band_subset.Read() == array.array(
+        "d", [0.0, 2.0, 3.0, 5.0, 12.0, 14.0, 15.0, 17.0]
+    )
+    assert indexed_twice_by_band_subset.Read(
+        array_start_idx=[1, 0, 0], array_step=[-1, 1, 1]
+    ) == array.array("d", [12.0, 14.0, 15.0, 17.0, 0.0, 2.0, 3.0, 5.0])
diff --git a/autotest/gdrivers/netcdf_multidim.py b/autotest/gdrivers/netcdf_multidim.py
index 49e390ea947d..634e3683a231 100755
--- a/autotest/gdrivers/netcdf_multidim.py
+++ b/autotest/gdrivers/netcdf_multidim.py
@@ -3751,7 +3751,7 @@ def test_netcdf_multidim_serialize_statistics_asclassicdataset(tmp_path):
     filename = str(
         tmp_path / "test_netcdf_multidim_serialize_statistics_asclassicdataset.nc"
     )
-    shutil.copy("data/netcdf/byte_no_cf.nc", filename)
+    shutil.copy("data/netcdf/byte.nc", filename)
 
     def test():
         ds = gdal.OpenEx(filename, gdal.OF_MULTIDIM_RASTER | gdal.OF_UPDATE)
@@ -3772,6 +3772,16 @@ def test():
         classic_ds = view.AsClassicDataset(1, 0)
         classic_ds.GetRasterBand(1).ComputeStatistics(False)
 
+        rg_subset = rg.SubsetDimensionFromSelection("/x=440750")
+        ds = rg_subset.OpenMDArray("Band1").AsClassicDataset(1, 0)
+        ds.GetRasterBand(1).ComputeStatistics(False)
+
+    def test2():
+        ds = gdal.OpenEx(filename, gdal.OF_MULTIDIM_RASTER | gdal.OF_UPDATE)
+        rg = ds.GetRootGroup()
+        rg_subset = rg.SubsetDimensionFromSelection("/x=440750")
+        rg_subset.OpenMDArray("Band1").GetStatistics(False, force=True)
+
     def reopen():
 
         aux_xml = open(filename + ".aux.xml", "rb").read().decode("UTF-8")
@@ -3783,6 +3793,11 @@ def reopen():
             '<DerivedDataset name="AsClassicDataset(1,0) view of Sliced view of /Band1 ([0:10,...])">'
             in aux_xml
         )
+        assert (
+            '<DerivedDataset name="AsClassicDataset(1,0) view of /Band1 with context Selection /x=440750">'
+            in aux_xml
+        )
+        assert '<Array name="/Band1" context="Selection /x=440750">' in aux_xml
 
         ds = gdal.OpenEx(filename, gdal.OF_MULTIDIM_RASTER)
         rg = ds.GetRootGroup()
@@ -3808,5 +3823,16 @@ def reopen():
             -1.0,
         ]
 
+        rg_subset = rg.SubsetDimensionFromSelection("/x=440750")
+
+        assert rg_subset.OpenMDArray("Band1").GetStatistics(False, False).min == 107
+        assert rg_subset.OpenMDArray("Band1").GetStatistics(False, False).max == 197
+
+        ds = rg_subset.OpenMDArray("Band1").AsClassicDataset(1, 0)
+        assert ds.GetRasterBand(1).GetStatistics(False, False) == pytest.approx(
+            [107.0, 197.0, 149.7, 26.595300336714]
+        )
+
     test()
+    test2()
     reopen()
diff --git a/frmts/hdf4/hdf4multidim.cpp b/frmts/hdf4/hdf4multidim.cpp
index 2f2686156b8d..dda53384f658 100644
--- a/frmts/hdf4/hdf4multidim.cpp
+++ b/frmts/hdf4/hdf4multidim.cpp
@@ -92,10 +92,21 @@ class HDF4Group final : public GDALGroup
     std::shared_ptr<HDF4SharedResources> m_poShared;
     std::shared_ptr<HDF4SDSGroup> m_poGDALGroup{};
 
-  public:
+  protected:
     HDF4Group(const std::string &osParentName, const std::string &osName,
               const std::shared_ptr<HDF4SharedResources> &poShared);
 
+  public:
+    static std::shared_ptr<HDF4Group>
+    Create(const std::string &osParentName, const std::string &osName,
+           const std::shared_ptr<HDF4SharedResources> &poShared)
+    {
+        auto poGroup = std::shared_ptr<HDF4Group>(
+            new HDF4Group(osParentName, osName, poShared));
+        poGroup->SetSelf(poGroup);
+        return poGroup;
+    }
+
     std::vector<std::shared_ptr<GDALAttribute>>
     GetAttributes(CSLConstList papszOptions = nullptr) const override;
 
@@ -3190,7 +3201,7 @@ void HDF4Dataset::OpenMultiDim(const char *pszFilename,
 
     hSD = -1;
 
-    m_poRootGroup = std::make_shared<HDF4Group>(std::string(), "/", poShared);
+    m_poRootGroup = HDF4Group::Create(std::string(), "/", poShared);
 
     SetDescription(pszFilename);
 
diff --git a/frmts/hdf5/hdf5multidim.cpp b/frmts/hdf5/hdf5multidim.cpp
index 275187ce7144..e7115be82a58 100644
--- a/frmts/hdf5/hdf5multidim.cpp
+++ b/frmts/hdf5/hdf5multidim.cpp
@@ -64,7 +64,7 @@ class HDF5Group final : public GDALGroup
     static herr_t GetAttributesCallback(hid_t hGroup, const char *pszObjName,
                                         void *);
 
-  public:
+  protected:
     HDF5Group(
         const std::string &osParentName, const std::string &osName,
         const std::shared_ptr<HDF5SharedResources> &poShared,
@@ -85,6 +85,19 @@ class HDF5Group final : public GDALGroup
         }
     }
 
+  public:
+    static std::shared_ptr<HDF5Group> Create(
+        const std::string &osParentName, const std::string &osName,
+        const std::shared_ptr<HDF5SharedResources> &poShared,
+        const std::set<std::pair<unsigned long, unsigned long>> &oSetParentIds,
+        hid_t hGroup, unsigned long objIds[2])
+    {
+        auto poGroup = std::shared_ptr<HDF5Group>(new HDF5Group(
+            osParentName, osName, poShared, oSetParentIds, hGroup, objIds));
+        poGroup->SetSelf(poGroup);
+        return poGroup;
+    }
+
     ~HDF5Group()
     {
         H5Gclose(m_hGroup);
@@ -537,7 +550,7 @@ std::shared_ptr<HDF5Group> HDF5SharedResources::GetRootGroup()
 
     auto poSharedResources = m_poSelf.lock();
     CPLAssert(poSharedResources != nullptr);
-    return std::make_shared<HDF5Group>(
+    return HDF5Group::Create(
         std::string(), "/", poSharedResources,
         std::set<std::pair<unsigned long, unsigned long>>(), hGroup,
         oStatbuf.objno);
@@ -778,9 +791,8 @@ std::shared_ptr<GDALGroup> HDF5Group::OpenGroup(const std::string &osName,
     {
         return nullptr;
     }
-    return std::make_shared<HDF5Group>(GetFullName(), osName, m_poShared,
-                                       m_oSetParentIds, hSubGroup,
-                                       oStatbuf.objno);
+    return HDF5Group::Create(GetFullName(), osName, m_poShared, m_oSetParentIds,
+                             hSubGroup, oStatbuf.objno);
 }
 
 /************************************************************************/
diff --git a/frmts/mem/memdataset.cpp b/frmts/mem/memdataset.cpp
index 2083b24989dd..5c46c01201b6 100644
--- a/frmts/mem/memdataset.cpp
+++ b/frmts/mem/memdataset.cpp
@@ -1452,8 +1452,9 @@ std::shared_ptr<GDALGroup> MEMGroup::OpenGroup(const std::string &osName,
 std::shared_ptr<MEMGroup> MEMGroup::Create(const std::string &osParentName,
                                            const char *pszName)
 {
-    auto newGroup(std::make_shared<MEMGroup>(osParentName, pszName));
-    newGroup->m_pSelf = newGroup;
+    auto newGroup(
+        std::shared_ptr<MEMGroup>(new MEMGroup(osParentName, pszName)));
+    newGroup->SetSelf(newGroup);
     return newGroup;
 }
 
@@ -1479,7 +1480,7 @@ std::shared_ptr<GDALGroup> MEMGroup::CreateGroup(const std::string &osName,
         return nullptr;
     }
     auto newGroup = MEMGroup::Create(GetFullName(), osName.c_str());
-    newGroup->m_pParent = m_pSelf;
+    newGroup->m_pParent = std::dynamic_pointer_cast<MEMGroup>(m_pSelf.lock());
     m_oMapGroups[osName] = newGroup;
     return newGroup;
 }
@@ -1717,8 +1718,9 @@ MEMGroup::CreateAttribute(const std::string &osName,
                  "An attribute with same name already exists");
         return nullptr;
     }
-    auto newAttr(
-        MEMAttribute::Create(m_pSelf.lock(), osName, anDimensions, oDataType));
+    auto newAttr(MEMAttribute::Create(
+        std::dynamic_pointer_cast<MEMGroup>(m_pSelf.lock()), osName,
+        anDimensions, oDataType));
     if (!newAttr)
         return nullptr;
     m_oMapAttributes[osName] = newAttr;
@@ -3049,8 +3051,9 @@ MEMGroup::CreateDimension(const std::string &osName, const std::string &osType,
                  "A dimension with same name already exists");
         return nullptr;
     }
-    auto newDim(MEMDimension::Create(m_pSelf.lock(), osName, osType,
-                                     osDirection, nSize));
+    auto newDim(MEMDimension::Create(
+        std::dynamic_pointer_cast<MEMGroup>(m_pSelf.lock()), osName, osType,
+        osDirection, nSize));
     m_oMapDimensions[osName] = newDim;
     return newDim;
 }
diff --git a/frmts/mem/memmultidim.h b/frmts/mem/memmultidim.h
index 64e27a2b3baa..d5e346c13a48 100644
--- a/frmts/mem/memmultidim.h
+++ b/frmts/mem/memmultidim.h
@@ -66,7 +66,6 @@ class CPL_DLL MEMGroup CPL_NON_FINAL : public GDALGroup,
     std::map<CPLString, std::shared_ptr<GDALGroup>> m_oMapGroups{};
     std::map<CPLString, std::shared_ptr<GDALMDArray>> m_oMapMDArrays{};
     std::map<CPLString, std::shared_ptr<GDALDimension>> m_oMapDimensions{};
-    std::weak_ptr<MEMGroup> m_pSelf{};
     std::weak_ptr<MEMGroup> m_pParent{};
 
   protected:
@@ -82,7 +81,6 @@ class CPL_DLL MEMGroup CPL_NON_FINAL : public GDALGroup,
 
     void NotifyChildrenOfDeletion() override;
 
-  public:
     MEMGroup(const std::string &osParentName, const char *pszName)
         : GDALGroup(osParentName, pszName ? pszName : "")
     {
@@ -90,6 +88,7 @@ class CPL_DLL MEMGroup CPL_NON_FINAL : public GDALGroup,
             m_osFullName = osParentName;
     }
 
+  public:
     static std::shared_ptr<MEMGroup> Create(const std::string &osParentName,
                                             const char *pszName);
 
diff --git a/frmts/netcdf/netcdfmultidim.cpp b/frmts/netcdf/netcdfmultidim.cpp
index 1e84dcbf6d49..ab5571c29576 100644
--- a/frmts/netcdf/netcdfmultidim.cpp
+++ b/frmts/netcdf/netcdfmultidim.cpp
@@ -228,7 +228,6 @@ class netCDFGroup final : public GDALGroup, public netCDFAttributeHolder
     std::shared_ptr<netCDFSharedResources> m_poShared;
     int m_gid = 0;
     CPLStringList m_aosStructuralInfo{};
-    std::weak_ptr<netCDFGroup> m_poSelf{};
     std::weak_ptr<netCDFGroup> m_poParent{};
     std::set<GDALGroup *> m_oSetGroups{};
     std::set<GDALDimension *> m_oSetDimensions{};
@@ -274,11 +273,15 @@ class netCDFGroup final : public GDALGroup, public netCDFAttributeHolder
 
     void NotifyChildrenOfRenaming() override;
 
-  public:
     netCDFGroup(const std::shared_ptr<netCDFSharedResources> &poShared,
                 int gid);
+
+  public:
     ~netCDFGroup();
 
+    static std::shared_ptr<netCDFGroup>
+    Create(const std::shared_ptr<netCDFSharedResources> &poShared, int cdfid);
+
     static std::shared_ptr<netCDFGroup>
     Create(const std::shared_ptr<netCDFSharedResources> &poShared,
            const std::shared_ptr<netCDFGroup> &poParent, int nSubGroupId);
@@ -344,11 +347,16 @@ class netCDFVirtualGroupBySameDimension final : public GDALGroup
     std::shared_ptr<netCDFGroup> m_poGroup;
     std::string m_osDimName{};
 
-  public:
+  protected:
     netCDFVirtualGroupBySameDimension(
         const std::shared_ptr<netCDFGroup> &poGroup,
         const std::string &osDimName);
 
+  public:
+    static std::shared_ptr<netCDFVirtualGroupBySameDimension>
+    Create(const std::shared_ptr<netCDFGroup> &poGroup,
+           const std::string &osDimName);
+
     std::vector<std::string>
     GetMDArrayNames(CSLConstList papszOptions) const override;
     std::shared_ptr<GDALMDArray>
@@ -792,14 +800,28 @@ netCDFGroup::~netCDFGroup()
 /*                              Create()                                */
 /************************************************************************/
 
+/* static */
+std::shared_ptr<netCDFGroup>
+netCDFGroup::Create(const std::shared_ptr<netCDFSharedResources> &poShared,
+                    int cdfid)
+{
+    auto poGroup =
+        std::shared_ptr<netCDFGroup>(new netCDFGroup(poShared, cdfid));
+    poGroup->SetSelf(poGroup);
+    return poGroup;
+}
+
+/************************************************************************/
+/*                              Create()                                */
+/************************************************************************/
+
 /* static */
 std::shared_ptr<netCDFGroup>
 netCDFGroup::Create(const std::shared_ptr<netCDFSharedResources> &poShared,
                     const std::shared_ptr<netCDFGroup> &poParent,
                     int nSubGroupId)
 {
-    auto poSubGroup = std::make_shared<netCDFGroup>(poShared, nSubGroupId);
-    poSubGroup->m_poSelf = poSubGroup;
+    auto poSubGroup = netCDFGroup::Create(poShared, nSubGroupId);
     poSubGroup->m_poParent = poParent;
     if (poParent)
         poParent->RegisterSubGroup(poSubGroup.get());
@@ -827,7 +849,9 @@ netCDFGroup::CreateGroup(const std::string &osName,
     NCDF_ERR(ret);
     if (ret != NC_NOERR)
         return nullptr;
-    return netCDFGroup::Create(m_poShared, m_poSelf.lock(), nSubGroupId);
+    return netCDFGroup::Create(
+        m_poShared, std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()),
+        nSubGroupId);
 }
 
 /************************************************************************/
@@ -853,8 +877,9 @@ netCDFGroup::CreateDimension(const std::string &osName,
                         static_cast<size_t>(bUnlimited ? 0 : nSize), &nDimId));
     if (nDimId < 0)
         return nullptr;
-    return netCDFDimension::Create(m_poShared, m_poSelf.lock(), m_gid, nDimId,
-                                   static_cast<size_t>(nSize), osType);
+    return netCDFDimension::Create(
+        m_poShared, std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()),
+        m_gid, nDimId, static_cast<size_t>(nSize), osType);
 }
 
 /************************************************************************/
@@ -1032,9 +1057,10 @@ std::shared_ptr<GDALMDArray> netCDFGroup::CreateMDArray(
             if (nc_inq_dimid(m_gid, dim->GetName().c_str(), &nDimId) ==
                 NC_NOERR)
             {
-                netCDFDim =
-                    netCDFDimension::Create(m_poShared, m_poSelf.lock(), m_gid,
-                                            nDimId, 0, dim->GetType());
+                netCDFDim = netCDFDimension::Create(
+                    m_poShared,
+                    std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()),
+                    m_gid, nDimId, 0, dim->GetType());
                 if (netCDFDim->GetSize() != dim->GetSize())
                 {
                     CPLError(CE_Warning, CPLE_AppDefined,
@@ -1180,8 +1206,9 @@ std::shared_ptr<GDALMDArray> netCDFGroup::CreateMDArray(
             return nullptr;
     }
 
-    return netCDFVariable::Create(m_poShared, m_poSelf.lock(), m_gid, nVarId,
-                                  dims, papszOptions, true);
+    return netCDFVariable::Create(
+        m_poShared, std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()),
+        m_gid, nVarId, dims, papszOptions, true);
 }
 
 /************************************************************************/
@@ -1192,9 +1219,9 @@ std::shared_ptr<GDALAttribute> netCDFGroup::CreateAttribute(
     const std::string &osName, const std::vector<GUInt64> &anDimensions,
     const GDALExtendedDataType &oDataType, CSLConstList papszOptions)
 {
-    return netCDFAttribute::Create(m_poShared, m_poSelf.lock(), m_gid,
-                                   NC_GLOBAL, osName, anDimensions, oDataType,
-                                   papszOptions);
+    return netCDFAttribute::Create(
+        m_poShared, std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()),
+        m_gid, NC_GLOBAL, osName, anDimensions, oDataType, papszOptions);
 }
 
 /************************************************************************/
@@ -1298,9 +1325,8 @@ netCDFGroup::OpenGroup(const std::string &osName,
             {
                 if (osCandidateGroupName == osName)
                 {
-                    auto poThisGroup =
-                        std::make_shared<netCDFGroup>(m_poShared, m_gid);
-                    return std::make_shared<netCDFVirtualGroupBySameDimension>(
+                    auto poThisGroup = netCDFGroup::Create(m_poShared, m_gid);
+                    return netCDFVirtualGroupBySameDimension::Create(
                         poThisGroup, osName);
                 }
             }
@@ -1311,7 +1337,9 @@ netCDFGroup::OpenGroup(const std::string &osName,
     if (nc_inq_grp_ncid(m_gid, osName.c_str(), &nSubGroupId) != NC_NOERR ||
         nSubGroupId <= 0)
         return nullptr;
-    return netCDFGroup::Create(m_poShared, m_poSelf.lock(), nSubGroupId);
+    return netCDFGroup::Create(
+        m_poShared, std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()),
+        nSubGroupId);
 }
 
 /************************************************************************/
@@ -1495,8 +1523,9 @@ netCDFGroup::OpenMDArray(const std::string &osName,
     if (nc_inq_varid(m_gid, osName.c_str(), &nVarId) != NC_NOERR)
         return nullptr;
     auto poVar = netCDFVariable::Create(
-        m_poShared, m_poSelf.lock(), m_gid, nVarId,
-        std::vector<std::shared_ptr<GDALDimension>>(), nullptr, false);
+        m_poShared, std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()),
+        m_gid, nVarId, std::vector<std::shared_ptr<GDALDimension>>(), nullptr,
+        false);
     if (poVar)
     {
         poVar->SetUseDefaultFillAsNoData(CPLTestBool(CSLFetchNameValueDef(
@@ -1525,9 +1554,10 @@ netCDFGroup::GetDimensions(CSLConstList) const
         auto poCachedDim = m_poShared->GetCachedDimension(dimids[i]);
         if (poCachedDim == nullptr)
         {
-            poCachedDim =
-                netCDFDimension::Create(m_poShared, m_poSelf.lock(), m_gid,
-                                        dimids[i], 0, std::string());
+            poCachedDim = netCDFDimension::Create(
+                m_poShared,
+                std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()), m_gid,
+                dimids[i], 0, std::string());
             m_poShared->CacheDimension(dimids[i], poCachedDim);
         }
         res.emplace_back(poCachedDim);
@@ -1576,8 +1606,9 @@ netCDFGroup::GetAttribute(const std::string &osName) const
         }
         return nullptr;
     }
-    return netCDFAttribute::Create(m_poShared, m_poSelf.lock(), m_gid,
-                                   NC_GLOBAL, osName);
+    return netCDFAttribute::Create(
+        m_poShared, std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()),
+        m_gid, NC_GLOBAL, osName);
 }
 
 /************************************************************************/
@@ -1600,7 +1631,9 @@ netCDFGroup::GetAttributes(CSLConstList) const
         if (!EQUAL(szAttrName, "_NCProperties"))
         {
             res.emplace_back(netCDFAttribute::Create(
-                m_poShared, m_poSelf.lock(), m_gid, NC_GLOBAL, szAttrName));
+                m_poShared,
+                std::dynamic_pointer_cast<netCDFGroup>(m_pSelf.lock()), m_gid,
+                NC_GLOBAL, szAttrName));
         }
     }
 
@@ -1656,6 +1689,20 @@ netCDFVirtualGroupBySameDimension::netCDFVirtualGroupBySameDimension(
 {
 }
 
+/************************************************************************/
+/*                              Create()                                */
+/************************************************************************/
+
+/* static */ std::shared_ptr<netCDFVirtualGroupBySameDimension>
+netCDFVirtualGroupBySameDimension::Create(
+    const std::shared_ptr<netCDFGroup> &poGroup, const std::string &osDimName)
+{
+    auto poNewGroup = std::shared_ptr<netCDFVirtualGroupBySameDimension>(
+        new netCDFVirtualGroupBySameDimension(poGroup, osDimName));
+    poNewGroup->SetSelf(poNewGroup);
+    return poNewGroup;
+}
+
 /************************************************************************/
 /*                         GetMDArrayNames()                            */
 /************************************************************************/
@@ -5413,7 +5460,7 @@ GDALDataset *netCDFDataset::OpenMultiDim(GDALOpenInfo *poOpenInfo)
         return nullptr;
     }
 
-    poDS->m_poRootGroup.reset(new netCDFGroup(poSharedResources, cdfid));
+    poDS->m_poRootGroup = netCDFGroup::Create(poSharedResources, cdfid);
 
     poDS->TryLoadXML();
 
diff --git a/frmts/vrt/vrtdataset.cpp b/frmts/vrt/vrtdataset.cpp
index 5d0f22c5aae2..de5b20cbf785 100644
--- a/frmts/vrt/vrtdataset.cpp
+++ b/frmts/vrt/vrtdataset.cpp
@@ -587,7 +587,7 @@ CPLErr VRTDataset::XMLInit(CPLXMLNode *psTree, const char *pszVRTPathIn)
             return CE_Failure;
         }
 
-        m_poRootGroup = std::make_shared<VRTGroup>(std::string(), "/");
+        m_poRootGroup = VRTGroup::Create(std::string(), "/");
         m_poRootGroup->SetIsRootGroup();
         if (!m_poRootGroup->XMLInit(m_poRootGroup, m_poRootGroup, psGroup,
                                     pszVRTPathIn))
@@ -1779,7 +1779,7 @@ VRTDataset::CreateMultiDimensional(const char *pszFilename,
     VRTDataset *poDS = new VRTDataset(0, 0);
     poDS->eAccess = GA_Update;
     poDS->SetDescription(pszFilename);
-    poDS->m_poRootGroup = std::make_shared<VRTGroup>(std::string(), "/");
+    poDS->m_poRootGroup = VRTGroup::Create(std::string(), "/");
     poDS->m_poRootGroup->SetIsRootGroup();
     poDS->m_poRootGroup->SetFilename(pszFilename);
     poDS->m_poRootGroup->SetDirty();
diff --git a/frmts/vrt/vrtdataset.h b/frmts/vrt/vrtdataset.h
index 16e37f44bd58..0d36eaea232b 100644
--- a/frmts/vrt/vrtdataset.h
+++ b/frmts/vrt/vrtdataset.h
@@ -1416,9 +1416,24 @@ class VRTGroup final : public GDALGroup
     void SetRootGroupRef(const std::weak_ptr<Ref> &rgRef);
     std::weak_ptr<Ref> GetRootGroupRef() const;
 
-  public:
+  protected:
+    friend class VRTMDArray;
+    friend std::shared_ptr<GDALMDArray>
+    VRTDerivedArrayCreate(const char *pszVRTPath, const CPLXMLNode *psTree);
+
     explicit VRTGroup(const char *pszVRTPath);
     VRTGroup(const std::string &osParentName, const std::string &osName);
+
+  public:
+    static std::shared_ptr<VRTGroup> Create(const std::string &osParentName,
+                                            const std::string &osName)
+    {
+        auto poGroup =
+            std::shared_ptr<VRTGroup>(new VRTGroup(osParentName, osName));
+        poGroup->SetSelf(poGroup);
+        return poGroup;
+    }
+
     ~VRTGroup();
 
     bool XMLInit(const std::shared_ptr<VRTGroup> &poRoot,
diff --git a/frmts/vrt/vrtmultidim.cpp b/frmts/vrt/vrtmultidim.cpp
index fbe068734cea..3dc94c3fc9e0 100644
--- a/frmts/vrt/vrtmultidim.cpp
+++ b/frmts/vrt/vrtmultidim.cpp
@@ -480,7 +480,7 @@ std::shared_ptr<GDALGroup> VRTGroup::CreateGroup(const std::string &osName,
         return nullptr;
     }
     SetDirty();
-    auto newGroup(std::make_shared<VRTGroup>(GetFullName(), osName.c_str()));
+    auto newGroup(VRTGroup::Create(GetFullName(), osName.c_str()));
     newGroup->SetRootGroupRef(GetRootGroupRef());
     m_oMapGroups[osName] = newGroup;
     return newGroup;
@@ -1106,7 +1106,7 @@ std::shared_ptr<VRTMDArray> VRTMDArray::Create(const char *pszVRTPath,
                                                const CPLXMLNode *psNode)
 {
     auto poDummyGroup =
-        std::make_shared<VRTGroup>(pszVRTPath ? pszVRTPath : "");
+        std::shared_ptr<VRTGroup>(new VRTGroup(pszVRTPath ? pszVRTPath : ""));
     auto poArray = Create(poDummyGroup, std::string(), psNode);
     if (poArray)
         poArray->m_poDummyOwningGroup = poDummyGroup;
@@ -2753,8 +2753,8 @@ CPLXMLNode *VRTArraySource::SerializeToXML(const char * /*pszVRTPath*/)
 /*                     VRTDerivedArrayCreate()                          */
 /************************************************************************/
 
-static std::shared_ptr<GDALMDArray>
-VRTDerivedArrayCreate(const char *pszVRTPath, const CPLXMLNode *psTree)
+std::shared_ptr<GDALMDArray> VRTDerivedArrayCreate(const char *pszVRTPath,
+                                                   const CPLXMLNode *psTree)
 {
     auto poArray = ParseArray(psTree, pszVRTPath, "DerivedArray");
 
@@ -2822,8 +2822,8 @@ VRTDerivedArrayCreate(const char *pszVRTPath, const CPLXMLNode *psTree)
                      CPLGetXMLNode(psStep, "Resample"))
         {
             std::vector<std::shared_ptr<GDALDimension>> apoNewDims;
-            auto poDummyGroup =
-                std::make_shared<VRTGroup>(pszVRTPath ? pszVRTPath : "");
+            auto poDummyGroup = std::shared_ptr<VRTGroup>(
+                new VRTGroup(pszVRTPath ? pszVRTPath : ""));
             for (const CPLXMLNode *psDimension =
                      CPLGetXMLNode(psResample, "Dimension");
                  psDimension; psDimension = psDimension->psNext)
diff --git a/frmts/zarr/zarr.h b/frmts/zarr/zarr.h
index 505b28e810cb..d2c395327243 100644
--- a/frmts/zarr/zarr.h
+++ b/frmts/zarr/zarr.h
@@ -341,7 +341,6 @@ class ZarrGroupBase CPL_NON_FINAL : public GDALGroup
     mutable bool m_bDimensionsInstantiated = false;
     bool m_bUpdatable = false;
     bool m_bDimSizeInUpdate = false;
-    std::weak_ptr<ZarrGroupBase> m_pSelf{};
 
     virtual void ExploreDirectory() const = 0;
     virtual void LoadAttributes() const = 0;
@@ -365,11 +364,6 @@ class ZarrGroupBase CPL_NON_FINAL : public GDALGroup
   public:
     ~ZarrGroupBase() override;
 
-    void SetSelf(const std::weak_ptr<ZarrGroupBase> &self)
-    {
-        m_pSelf = self;
-    }
-
     std::shared_ptr<GDALAttribute>
     GetAttribute(const std::string &osName) const override
     {
diff --git a/frmts/zarr/zarr_group.cpp b/frmts/zarr/zarr_group.cpp
index 2c50542d0fef..34db203d769a 100644
--- a/frmts/zarr/zarr_group.cpp
+++ b/frmts/zarr/zarr_group.cpp
@@ -76,7 +76,8 @@ void ZarrGroupBase::RegisterArray(const std::shared_ptr<ZarrArray> &array) const
     {
         m_aosArrays.emplace_back(array->GetName());
     }
-    array->RegisterGroup(m_pSelf);
+    array->RegisterGroup(
+        std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock()));
 }
 
 /************************************************************************/
@@ -314,9 +315,10 @@ std::shared_ptr<GDALDimension> ZarrGroupBase::CreateDimension(
                  "A dimension with same name already exists");
         return nullptr;
     }
-    auto newDim(std::make_shared<ZarrDimension>(m_poSharedResource, m_pSelf,
-                                                GetFullName(), osName, osType,
-                                                osDirection, nSize));
+    auto newDim(std::make_shared<ZarrDimension>(
+        m_poSharedResource,
+        std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock()), GetFullName(),
+        osName, osType, osDirection, nSize));
     newDim->SetXArrayDimension();
     m_oMapDimensions[osName] = newDim;
     return newDim;
@@ -507,7 +509,8 @@ bool ZarrGroupBase::Rename(const std::string &osNewName)
         {
             pParent->m_oMapGroups.erase(oIter);
             CPLAssert(m_pSelf.lock());
-            pParent->m_oMapGroups[osNewName] = m_pSelf.lock();
+            pParent->m_oMapGroups[osNewName] =
+                std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock());
         }
 
         for (auto &osName : pParent->m_aosGroups)
diff --git a/frmts/zarr/zarr_v2_array.cpp b/frmts/zarr/zarr_v2_array.cpp
index d9677c771f7c..cb5e6ac5c545 100644
--- a/frmts/zarr/zarr_v2_array.cpp
+++ b/frmts/zarr/zarr_v2_array.cpp
@@ -1402,8 +1402,10 @@ ZarrV2Group::LoadArray(const std::string &osArrayName,
             return nullptr;
         }
         aoDims.emplace_back(std::make_shared<ZarrDimension>(
-            m_poSharedResource, m_pSelf, std::string(), CPLSPrintf("dim%d", i),
-            std::string(), std::string(), nSize));
+            m_poSharedResource,
+            std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock()),
+            std::string(), CPLSPrintf("dim%d", i), std::string(), std::string(),
+            nSize));
     }
 
     // XArray extension
@@ -1511,8 +1513,9 @@ ZarrV2Group::LoadArray(const std::string &osArrayName,
         }
 
         auto poDimLocal = std::make_shared<ZarrDimension>(
-            m_poSharedResource, m_pSelf, GetFullName(), osDimName, osType,
-            osDirection, poDim->GetSize());
+            m_poSharedResource,
+            std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock()),
+            GetFullName(), osDimName, osType, osDirection, poDim->GetSize());
         poDimLocal->SetXArrayDimension();
         m_oMapDimensions[osDimName] = poDimLocal;
         poDim = poDimLocal;
@@ -1553,7 +1556,8 @@ ZarrV2Group::LoadArray(const std::string &osArrayName,
         const auto arrayDims = nczarrArrayDimrefs.ToArray();
         if (arrayDims.Size() == oShape.Size())
         {
-            auto poRG = m_pSelf.lock();
+            auto poRG =
+                std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock());
             CPLAssert(poRG != nullptr);
             while (true)
             {
@@ -1617,7 +1621,10 @@ ZarrV2Group::LoadArray(const std::string &osArrayName,
                                     auto poDimLocal =
                                         std::make_shared<ZarrDimension>(
                                             m_poSharedResource,
-                                            poDimParentGroup->m_pSelf,
+                                            std::dynamic_pointer_cast<
+                                                ZarrGroupBase>(
+                                                poDimParentGroup->m_pSelf
+                                                    .lock()),
                                             poDimParentGroup->GetFullName(),
                                             poDim->GetName(), osType,
                                             osDirection, poDim->GetSize());
diff --git a/frmts/zarr/zarr_v2_group.cpp b/frmts/zarr/zarr_v2_group.cpp
index f16b2f5cf0a3..eb70ade14516 100644
--- a/frmts/zarr/zarr_v2_group.cpp
+++ b/frmts/zarr/zarr_v2_group.cpp
@@ -189,7 +189,8 @@ ZarrV2Group::OpenZarrGroup(const std::string &osName, CSLConstList) const
 
             auto poSubGroup =
                 ZarrV2Group::Create(m_poSharedResource, GetFullName(), osName);
-            poSubGroup->m_poParent = m_pSelf;
+            poSubGroup->m_poParent =
+                std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock());
             poSubGroup->SetUpdatable(m_bUpdatable);
             poSubGroup->SetDirectoryName(osSubDir);
             m_oMapGroups[osName] = poSubGroup;
@@ -251,7 +252,8 @@ ZarrV2Group::GetOrCreateSubGroup(const std::string &osSubGroupFullname)
     poSubGroup =
         ZarrV2Group::Create(m_poSharedResource, poBelongingGroup->GetFullName(),
                             osSubGroupFullname.substr(nLastSlashPos + 1));
-    poSubGroup->m_poParent = poBelongingGroup->m_pSelf;
+    poSubGroup->m_poParent = std::dynamic_pointer_cast<ZarrGroupBase>(
+        poBelongingGroup->m_pSelf.lock());
     poSubGroup->SetDirectoryName(
         CPLFormFilename(poBelongingGroup->m_osDirectoryName.c_str(),
                         poSubGroup->GetName().c_str(), nullptr));
@@ -629,7 +631,8 @@ ZarrV2Group::CreateGroup(const std::string &osName,
                                 osDirectoryName);
     if (!poGroup)
         return nullptr;
-    poGroup->m_poParent = m_pSelf;
+    poGroup->m_poParent =
+        std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock());
     m_oMapGroups[osName] = poGroup;
     m_aosGroups.emplace_back(osName);
     return poGroup;
diff --git a/frmts/zarr/zarr_v3_array.cpp b/frmts/zarr/zarr_v3_array.cpp
index 7b7fdf713770..23bc60c2a591 100644
--- a/frmts/zarr/zarr_v3_array.cpp
+++ b/frmts/zarr/zarr_v3_array.cpp
@@ -1243,8 +1243,10 @@ ZarrV3Group::LoadArray(const std::string &osArrayName,
             return nullptr;
         }
         aoDims.emplace_back(std::make_shared<ZarrDimension>(
-            m_poSharedResource, m_pSelf, std::string(), CPLSPrintf("dim%d", i),
-            std::string(), std::string(), nSize));
+            m_poSharedResource,
+            std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock()),
+            std::string(), CPLSPrintf("dim%d", i), std::string(), std::string(),
+            nSize));
     }
 
     // Deal with dimension_names
@@ -1331,8 +1333,9 @@ ZarrV3Group::LoadArray(const std::string &osArrayName,
         }
 
         auto poDimLocal = std::make_shared<ZarrDimension>(
-            m_poSharedResource, m_pSelf, GetFullName(), osDimName, osType,
-            osDirection, poDim->GetSize());
+            m_poSharedResource,
+            std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock()),
+            GetFullName(), osDimName, osType, osDirection, poDim->GetSize());
         poDimLocal->SetXArrayDimension();
         m_oMapDimensions[osDimName] = poDimLocal;
         poDim = poDimLocal;
diff --git a/frmts/zarr/zarr_v3_group.cpp b/frmts/zarr/zarr_v3_group.cpp
index 41f262602808..d5a42a6f933c 100644
--- a/frmts/zarr/zarr_v3_group.cpp
+++ b/frmts/zarr/zarr_v3_group.cpp
@@ -254,7 +254,8 @@ ZarrV3Group::OpenZarrGroup(const std::string &osName, CSLConstList) const
             }
             auto poSubGroup = ZarrV3Group::Create(
                 m_poSharedResource, GetFullName(), osName, osSubDir);
-            poSubGroup->m_poParent = m_pSelf;
+            poSubGroup->m_poParent =
+                std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock());
             poSubGroup->SetUpdatable(m_bUpdatable);
             m_oMapGroups[osName] = poSubGroup;
             return poSubGroup;
@@ -267,7 +268,8 @@ ZarrV3Group::OpenZarrGroup(const std::string &osName, CSLConstList) const
     {
         auto poSubGroup = ZarrV3Group::Create(m_poSharedResource, GetFullName(),
                                               osName, osSubDir);
-        poSubGroup->m_poParent = m_pSelf;
+        poSubGroup->m_poParent =
+            std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock());
         poSubGroup->SetUpdatable(m_bUpdatable);
         m_oMapGroups[osName] = poSubGroup;
         return poSubGroup;
@@ -363,7 +365,8 @@ ZarrV3Group::CreateGroup(const std::string &osName,
                                 osDirectoryName);
     if (!poGroup)
         return nullptr;
-    poGroup->m_poParent = m_pSelf;
+    poGroup->m_poParent =
+        std::dynamic_pointer_cast<ZarrGroupBase>(m_pSelf.lock());
     m_oMapGroups[osName] = poGroup;
     m_aosGroups.emplace_back(osName);
     return poGroup;
diff --git a/gcore/CMakeLists.txt b/gcore/CMakeLists.txt
index 227768a80a2f..304b94f2268a 100644
--- a/gcore/CMakeLists.txt
+++ b/gcore/CMakeLists.txt
@@ -45,6 +45,7 @@ add_library(
   rawdataset.cpp
   gdalmultidim.cpp
   gdalmultidim_gridded.cpp
+  gdalmultidim_subsetdimension.cpp
   gdalpython.cpp
   gdalpythondriverloader.cpp
   tilematrixset.cpp
diff --git a/gcore/gdal.h b/gcore/gdal.h
index cee606c2570c..49cc651d11ab 100644
--- a/gcore/gdal.h
+++ b/gcore/gdal.h
@@ -2219,6 +2219,8 @@ GDALAttributeH CPL_DLL GDALGroupCreateAttribute(
 bool CPL_DLL GDALGroupDeleteAttribute(GDALGroupH hGroup, const char *pszName,
                                       CSLConstList papszOptions);
 bool CPL_DLL GDALGroupRename(GDALGroupH hGroup, const char *pszNewName);
+GDALGroupH CPL_DLL GDALGroupSubsetDimensionFromSelection(
+    GDALGroupH hGroup, const char *pszSelection, CSLConstList papszOptions);
 
 void CPL_DLL GDALMDArrayRelease(GDALMDArrayH hMDArray);
 const char CPL_DLL *GDALMDArrayGetName(GDALMDArrayH hArray);
diff --git a/gcore/gdal_pam.h b/gcore/gdal_pam.h
index 5b23e3d0e2f1..b40ae21826a7 100644
--- a/gcore/gdal_pam.h
+++ b/gcore/gdal_pam.h
@@ -380,22 +380,30 @@ class CPL_DLL GDALPamMultiDim
     virtual ~GDALPamMultiDim();
 
     std::shared_ptr<OGRSpatialReference>
-    GetSpatialRef(const std::string &osArrayFullName);
+    GetSpatialRef(const std::string &osArrayFullName,
+                  const std::string &osContext);
 
     void SetSpatialRef(const std::string &osArrayFullName,
+                       const std::string &osContext,
                        const OGRSpatialReference *poSRS);
 
-    CPLErr GetStatistics(const std::string &osArrayFullName, bool bApproxOK,
+    CPLErr GetStatistics(const std::string &osArrayFullName,
+                         const std::string &osContext, bool bApproxOK,
                          double *pdfMin, double *pdfMax, double *pdfMean,
                          double *pdfStdDev, GUInt64 *pnValidCount);
 
-    void SetStatistics(const std::string &osArrayFullName, bool bApproxStats,
+    void SetStatistics(const std::string &osArrayFullName,
+                       const std::string &osContext, bool bApproxStats,
                        double dfMin, double dfMax, double dfMean,
                        double dfStdDev, GUInt64 nValidCount);
 
     void ClearStatistics();
 
-    void ClearStatistics(const std::string &osArrayFullName);
+    void ClearStatistics(const std::string &osArrayFullName,
+                         const std::string &osContext);
+
+    static std::shared_ptr<GDALPamMultiDim>
+    GetPAM(const std::shared_ptr<GDALMDArray> &poParent);
 };
 
 /* ******************************************************************** */
@@ -409,7 +417,8 @@ class CPL_DLL GDALPamMDArray : public GDALMDArray
 
   protected:
     GDALPamMDArray(const std::string &osParentName, const std::string &osName,
-                   const std::shared_ptr<GDALPamMultiDim> &poPam);
+                   const std::shared_ptr<GDALPamMultiDim> &poPam,
+                   const std::string &osContext = std::string());
 
     bool SetStatistics(bool bApproxStats, double dfMin, double dfMax,
                        double dfMean, double dfStdDev, GUInt64 nValidCount,
diff --git a/gcore/gdal_priv.h b/gcore/gdal_priv.h
index 06e077feae9f..8434f692ab90 100644
--- a/gcore/gdal_priv.h
+++ b/gcore/gdal_priv.h
@@ -2328,12 +2328,21 @@ class CPL_DLL GDALGroup : public GDALIHasAttribute
   protected:
     //! @cond Doxygen_Suppress
     std::string m_osName{};
+
+    // This is actually a path of the form "/parent_path/{m_osName}"
     std::string m_osFullName{};
 
+    // Used for example by GDALSubsetGroup to distinguish a derived group
+    //from its original, without altering its name
+    const std::string m_osContext{};
+
+    std::weak_ptr<GDALGroup> m_pSelf{};
+
     //! Can be set to false by the owing group, when deleting this object
     bool m_bValid = true;
 
-    GDALGroup(const std::string &osParentName, const std::string &osName);
+    GDALGroup(const std::string &osParentName, const std::string &osName,
+              const std::string &osContext = std::string());
 
     const GDALGroup *
     GetInnerMostGroup(const std::string &osPathOrArrayOrDim,
@@ -2344,6 +2353,11 @@ class CPL_DLL GDALGroup : public GDALIHasAttribute
 
     bool CheckValidAndErrorOutIfNot() const;
 
+    void SetSelf(const std::shared_ptr<GDALGroup> &self)
+    {
+        m_pSelf = self;
+    }
+
     virtual void NotifyChildrenOfRenaming()
     {
     }
@@ -2446,12 +2460,20 @@ class CPL_DLL GDALGroup : public GDALIHasAttribute
 
     virtual bool Rename(const std::string &osNewName);
 
+    std::shared_ptr<GDALGroup>
+    SubsetDimensionFromSelection(const std::string &osSelection) const;
+
     //! @cond Doxygen_Suppress
     virtual void ParentRenamed(const std::string &osNewParentFullName);
 
     virtual void Deleted();
 
     virtual void ParentDeleted();
+
+    const std::string &GetContext() const
+    {
+        return m_osContext;
+    }
     //! @endcond
 
     //! @cond Doxygen_Suppress
@@ -2473,6 +2495,8 @@ class CPL_DLL GDALAbstractMDArray
   protected:
     //! @cond Doxygen_Suppress
     std::string m_osName{};
+
+    // This is actually a path of the form "/parent_path/{m_osName}"
     std::string m_osFullName{};
     std::weak_ptr<GDALAbstractMDArray> m_pSelf{};
 
@@ -2828,12 +2852,17 @@ class CPL_DLL GDALMDArray : virtual public GDALAbstractMDArray,
         return atInternal(indices, tail...);
     }
 
+    // Used for example by GDALSubsetGroup to distinguish a derived group
+    //from its original, without altering its name
+    const std::string m_osContext{};
+
     mutable bool m_bHasTriedCachedArray = false;
     mutable std::shared_ptr<GDALMDArray> m_poCachedArray{};
 
   protected:
     //! @cond Doxygen_Suppress
-    GDALMDArray(const std::string &osParentName, const std::string &osName);
+    GDALMDArray(const std::string &osParentName, const std::string &osName,
+                const std::string &osContext = std::string());
 
     virtual bool IAdviseRead(const GUInt64 *arrayStartIdx, const size_t *count,
                              CSLConstList papszOptions) const;
@@ -3065,6 +3094,11 @@ class CPL_DLL GDALMDArray : virtual public GDALAbstractMDArray,
     virtual std::shared_ptr<GDALMDArray>
     GetView(const std::string &viewExpr, bool bRenameDimensions,
             std::vector<ViewSpec> &viewSpecs) const;
+
+    const std::string &GetContext() const
+    {
+        return m_osContext;
+    }
     //! @endcond
 };
 
diff --git a/gcore/gdalmultidim.cpp b/gcore/gdalmultidim.cpp
index c206f0dce208..d52706bd9118 100644
--- a/gcore/gdalmultidim.cpp
+++ b/gcore/gdalmultidim.cpp
@@ -33,6 +33,7 @@
 #include <limits>
 #include <queue>
 #include <set>
+#include <utility>
 #include <time.h>
 
 #include <ctype.h>  // isalnum
@@ -49,19 +50,6 @@
 #define COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT
 #endif
 
-/************************************************************************/
-/*                          GetPAM()                                    */
-/************************************************************************/
-
-static std::shared_ptr<GDALPamMultiDim>
-GetPAM(const std::shared_ptr<GDALMDArray> &poParent)
-{
-    auto poPamArray = dynamic_cast<GDALPamMDArray *>(poParent.get());
-    if (poPamArray)
-        return poPamArray->GetPAM();
-    return nullptr;
-}
-
 /************************************************************************/
 /*                       GDALMDArrayUnscaled                            */
 /************************************************************************/
@@ -79,9 +67,9 @@ class GDALMDArrayUnscaled final : public GDALPamMDArray
     explicit GDALMDArrayUnscaled(const std::shared_ptr<GDALMDArray> &poParent)
         : GDALAbstractMDArray(std::string(),
                               "Unscaled view of " + poParent->GetFullName()),
-          GDALPamMDArray(std::string(),
-                         "Unscaled view of " + poParent->GetFullName(),
-                         ::GetPAM(poParent)),
+          GDALPamMDArray(
+              std::string(), "Unscaled view of " + poParent->GetFullName(),
+              GDALPamMultiDim::GetPAM(poParent), poParent->GetContext()),
           m_poParent(std::move(poParent)),
           m_dt(GDALExtendedDataType::Create(
               GDALDataTypeIsComplex(
@@ -348,12 +336,14 @@ bool GDALIHasAttribute::DeleteAttribute(CPL_UNUSED const std::string &osName,
 /************************************************************************/
 
 //! @cond Doxygen_Suppress
-GDALGroup::GDALGroup(const std::string &osParentName, const std::string &osName)
+GDALGroup::GDALGroup(const std::string &osParentName, const std::string &osName,
+                     const std::string &osContext)
     : m_osName(osParentName.empty() ? "/" : osName),
       m_osFullName(
           !osParentName.empty()
               ? ((osParentName == "/" ? "/" : osParentName + "/") + osName)
-              : "/")
+              : "/"),
+      m_osContext(osContext)
 {
 }
 //! @endcond
@@ -3636,10 +3626,13 @@ bool GDALAttribute::Write(const double *vals, size_t nVals)
 
 //! @cond Doxygen_Suppress
 GDALMDArray::GDALMDArray(CPL_UNUSED const std::string &osParentName,
-                         CPL_UNUSED const std::string &osName)
+                         CPL_UNUSED const std::string &osName,
+                         const std::string &osContext)
+    :
 #if !defined(COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT)
-    : GDALAbstractMDArray(osParentName, osName)
+      GDALAbstractMDArray(osParentName, osName),
 #endif
+      m_osContext(osContext)
 {
 }
 //! @endcond
@@ -4821,7 +4814,8 @@ class GDALSlicedMDArray final : public GDALPamMDArray
           GDALPamMDArray(std::string(),
                          "Sliced view of " + poParent->GetFullName() + " (" +
                              viewExpr + ")",
-                         ::GetPAM(poParent)),
+                         GDALPamMultiDim::GetPAM(poParent),
+                         poParent->GetContext()),
           m_poParent(std::move(poParent)), m_dims(std::move(dims)),
           m_mapDimIdxToParentDimIdx(std::move(mapDimIdxToParentDimIdx)),
           m_parentRanges(std::move(parentRanges)),
@@ -5287,10 +5281,10 @@ class GDALExtractFieldMDArray final : public GDALPamMDArray
         : GDALAbstractMDArray(std::string(), "Extract field " + fieldName +
                                                  " of " +
                                                  poParent->GetFullName()),
-          GDALPamMDArray(std::string(),
-                         "Extract field " + fieldName + " of " +
-                             poParent->GetFullName(),
-                         ::GetPAM(poParent)),
+          GDALPamMDArray(
+              std::string(),
+              "Extract field " + fieldName + " of " + poParent->GetFullName(),
+              GDALPamMultiDim::GetPAM(poParent), poParent->GetContext()),
           m_poParent(poParent), m_dt(srcComp->GetType()),
           m_srcCompName(srcComp->GetName())
     {
@@ -5706,7 +5700,8 @@ class GDALMDArrayTransposed final : public GDALPamMDArray
           GDALPamMDArray(std::string(),
                          "Transposed view of " + poParent->GetFullName() +
                              " along " + MappingToStr(anMapNewAxisToOldAxis),
-                         ::GetPAM(poParent)),
+                         GDALPamMultiDim::GetPAM(poParent),
+                         poParent->GetContext()),
           m_poParent(std::move(poParent)),
           m_anMapNewAxisToOldAxis(anMapNewAxisToOldAxis),
           m_dims(std::move(dims)),
@@ -6440,7 +6435,8 @@ class GDALMDArrayMask final : public GDALPamMDArray
         : GDALAbstractMDArray(std::string(),
                               "Mask of " + poParent->GetFullName()),
           GDALPamMDArray(std::string(), "Mask of " + poParent->GetFullName(),
-                         ::GetPAM(poParent)),
+                         GDALPamMultiDim::GetPAM(poParent),
+                         poParent->GetContext()),
           m_poParent(std::move(poParent))
     {
     }
@@ -7586,9 +7582,9 @@ class GDALMDArrayResampled final : public GDALPamMDArray
         const std::vector<GUInt64> &anBlockSize)
         : GDALAbstractMDArray(std::string(),
                               "Resampled view of " + poParent->GetFullName()),
-          GDALPamMDArray(std::string(),
-                         "Resampled view of " + poParent->GetFullName(),
-                         ::GetPAM(poParent)),
+          GDALPamMDArray(
+              std::string(), "Resampled view of " + poParent->GetFullName(),
+              GDALPamMultiDim::GetPAM(poParent), poParent->GetContext()),
           m_poParent(std::move(poParent)), m_apoDims(apoDims),
           m_anBlockSize(anBlockSize), m_dt(m_poParent->GetDataType())
     {
@@ -9158,9 +9154,15 @@ GDALDatasetFromArray *GDALDatasetFromArray::Create(
     if (!array->GetFilename().empty())
     {
         poDS->SetPhysicalFilename(array->GetFilename().c_str());
-        poDS->SetDerivedDatasetName(
+        std::string osDerivedDatasetName(
             CPLSPrintf("AsClassicDataset(%d,%d) view of %s", int(iXDim),
                        int(iYDim), array->GetFullName().c_str()));
+        if (!array->GetContext().empty())
+        {
+            osDerivedDatasetName += " with context ";
+            osDerivedDatasetName += array->GetContext();
+        }
+        poDS->SetDerivedDatasetName(osDerivedDatasetName.c_str());
         poDS->TryLoadXML();
     }
 
@@ -11139,6 +11141,31 @@ bool GDALGroupRename(GDALGroupH hGroup, const char *pszNewName)
     return hGroup->m_poImpl->Rename(pszNewName);
 }
 
+/************************************************************************/
+/*                 GDALGroupSubsetDimensionFromSelection()              */
+/************************************************************************/
+
+/** Return a virtual group whose one dimension has been subset according to a
+ * selection.
+ *
+ * This is the same as the C++ method GDALGroup::SubsetDimensionFromSelection().
+ *
+ * @return a virtual group, to be freed with GDALGroupRelease(), or nullptr.
+ */
+GDALGroupH
+GDALGroupSubsetDimensionFromSelection(GDALGroupH hGroup,
+                                      const char *pszSelection,
+                                      CPL_UNUSED CSLConstList papszOptions)
+{
+    VALIDATE_POINTER1(hGroup, __func__, nullptr);
+    VALIDATE_POINTER1(pszSelection, __func__, nullptr);
+    auto hNewGroup = hGroup->m_poImpl->SubsetDimensionFromSelection(
+        std::string(pszSelection));
+    if (!hNewGroup)
+        return nullptr;
+    return new GDALGroupHS(hNewGroup);
+}
+
 /************************************************************************/
 /*                        GDALMDArrayRelease()                          */
 /************************************************************************/
@@ -13336,7 +13363,8 @@ struct GDALPamMultiDim::Private
         Statistics stats{};
     };
 
-    std::map<std::string, ArrayInfo> m_oMapArray{};
+    typedef std::pair<std::string, std::string> NameContext;
+    std::map<NameContext, ArrayInfo> m_oMapArray{};
     std::vector<CPLXMLTreeCloser> m_apoOtherNodes{};
     bool m_bDirty = false;
     bool m_bLoaded = false;
@@ -13397,6 +13425,9 @@ void GDALPamMultiDim::Load()
             const char *pszName = CPLGetXMLValue(psIter, "name", nullptr);
             if (!pszName)
                 continue;
+            const char *pszContext = CPLGetXMLValue(psIter, "context", "");
+            const auto oKey =
+                std::pair<std::string, std::string>(pszName, pszContext);
 
             /* --------------------------------------------------------------------
              */
@@ -13435,7 +13466,7 @@ void GDALPamMultiDim::Load()
                 if (pszCoordinateEpoch)
                     poSRS->SetCoordinateEpoch(CPLAtof(pszCoordinateEpoch));
 
-                d->m_oMapArray[pszName].poSRS = poSRS;
+                d->m_oMapArray[oKey].poSRS = poSRS;
             }
 
             const CPLXMLNode *psStatistics =
@@ -13456,7 +13487,7 @@ void GDALPamMultiDim::Load()
                     CPLAtofM(CPLGetXMLValue(psStatistics, "StdDev", "0"));
                 sStats.nValidCount = static_cast<GUInt64>(CPLAtoGIntBig(
                     CPLGetXMLValue(psStatistics, "ValidSampleCount", "0")));
-                d->m_oMapArray[pszName].stats = sStats;
+                d->m_oMapArray[oKey].stats = sStats;
             }
         }
         else
@@ -13486,7 +13517,12 @@ void GDALPamMultiDim::Save()
     {
         CPLXMLNode *psArrayNode =
             CPLCreateXMLNode(oTree.get(), CXT_Element, "Array");
-        CPLAddXMLAttributeAndValue(psArrayNode, "name", kv.first.c_str());
+        CPLAddXMLAttributeAndValue(psArrayNode, "name", kv.first.first.c_str());
+        if (!kv.first.second.empty())
+        {
+            CPLAddXMLAttributeAndValue(psArrayNode, "context",
+                                       kv.first.second.c_str());
+        }
         if (kv.second.poSRS)
         {
             char *pszWKT = nullptr;
@@ -13578,10 +13614,12 @@ void GDALPamMultiDim::Save()
 /************************************************************************/
 
 std::shared_ptr<OGRSpatialReference>
-GDALPamMultiDim::GetSpatialRef(const std::string &osArrayFullName)
+GDALPamMultiDim::GetSpatialRef(const std::string &osArrayFullName,
+                               const std::string &osContext)
 {
     Load();
-    auto oIter = d->m_oMapArray.find(osArrayFullName);
+    auto oIter =
+        d->m_oMapArray.find(std::make_pair(osArrayFullName, osContext));
     if (oIter != d->m_oMapArray.end())
         return oIter->second.poSRS;
     return nullptr;
@@ -13592,14 +13630,17 @@ GDALPamMultiDim::GetSpatialRef(const std::string &osArrayFullName)
 /************************************************************************/
 
 void GDALPamMultiDim::SetSpatialRef(const std::string &osArrayFullName,
+                                    const std::string &osContext,
                                     const OGRSpatialReference *poSRS)
 {
     Load();
     d->m_bDirty = true;
     if (poSRS && !poSRS->IsEmpty())
-        d->m_oMapArray[osArrayFullName].poSRS.reset(poSRS->Clone());
+        d->m_oMapArray[std::make_pair(osArrayFullName, osContext)].poSRS.reset(
+            poSRS->Clone());
     else
-        d->m_oMapArray[osArrayFullName].poSRS.reset();
+        d->m_oMapArray[std::make_pair(osArrayFullName, osContext)]
+            .poSRS.reset();
 }
 
 /************************************************************************/
@@ -13607,12 +13648,14 @@ void GDALPamMultiDim::SetSpatialRef(const std::string &osArrayFullName,
 /************************************************************************/
 
 CPLErr GDALPamMultiDim::GetStatistics(const std::string &osArrayFullName,
+                                      const std::string &osContext,
                                       bool bApproxOK, double *pdfMin,
                                       double *pdfMax, double *pdfMean,
                                       double *pdfStdDev, GUInt64 *pnValidCount)
 {
     Load();
-    auto oIter = d->m_oMapArray.find(osArrayFullName);
+    auto oIter =
+        d->m_oMapArray.find(std::make_pair(osArrayFullName, osContext));
     if (oIter == d->m_oMapArray.end())
         return CE_Failure;
     const auto &stats = oIter->second.stats;
@@ -13638,13 +13681,15 @@ CPLErr GDALPamMultiDim::GetStatistics(const std::string &osArrayFullName,
 /************************************************************************/
 
 void GDALPamMultiDim::SetStatistics(const std::string &osArrayFullName,
+                                    const std::string &osContext,
                                     bool bApproxStats, double dfMin,
                                     double dfMax, double dfMean,
                                     double dfStdDev, GUInt64 nValidCount)
 {
     Load();
     d->m_bDirty = true;
-    auto &stats = d->m_oMapArray[osArrayFullName].stats;
+    auto &stats =
+        d->m_oMapArray[std::make_pair(osArrayFullName, osContext)].stats;
     stats.bHasStats = true;
     stats.bApproxStats = bApproxStats;
     stats.dfMin = dfMin;
@@ -13658,11 +13703,13 @@ void GDALPamMultiDim::SetStatistics(const std::string &osArrayFullName,
 /*                           ClearStatistics()                          */
 /************************************************************************/
 
-void GDALPamMultiDim::ClearStatistics(const std::string &osArrayFullName)
+void GDALPamMultiDim::ClearStatistics(const std::string &osArrayFullName,
+                                      const std::string &osContext)
 {
     Load();
     d->m_bDirty = true;
-    d->m_oMapArray[osArrayFullName].stats.bHasStats = false;
+    d->m_oMapArray[std::make_pair(osArrayFullName, osContext)].stats.bHasStats =
+        false;
 }
 
 /************************************************************************/
@@ -13677,18 +13724,32 @@ void GDALPamMultiDim::ClearStatistics()
         kv.second.stats.bHasStats = false;
 }
 
+/************************************************************************/
+/*                             GetPAM()                                 */
+/************************************************************************/
+
+/*static*/ std::shared_ptr<GDALPamMultiDim>
+GDALPamMultiDim::GetPAM(const std::shared_ptr<GDALMDArray> &poParent)
+{
+    auto poPamArray = dynamic_cast<GDALPamMDArray *>(poParent.get());
+    if (poPamArray)
+        return poPamArray->GetPAM();
+    return nullptr;
+}
+
 /************************************************************************/
 /*                           GDALPamMDArray                             */
 /************************************************************************/
 
 GDALPamMDArray::GDALPamMDArray(const std::string &osParentName,
                                const std::string &osName,
-                               const std::shared_ptr<GDALPamMultiDim> &poPam)
+                               const std::shared_ptr<GDALPamMultiDim> &poPam,
+                               const std::string &osContext)
     :
 #if !defined(COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT)
       GDALAbstractMDArray(osParentName, osName),
 #endif
-      GDALMDArray(osParentName, osName), m_poPam(poPam)
+      GDALMDArray(osParentName, osName, osContext), m_poPam(poPam)
 {
 }
 
@@ -13700,7 +13761,7 @@ bool GDALPamMDArray::SetSpatialRef(const OGRSpatialReference *poSRS)
 {
     if (!m_poPam)
         return false;
-    m_poPam->SetSpatialRef(GetFullName(), poSRS);
+    m_poPam->SetSpatialRef(GetFullName(), GetContext(), poSRS);
     return true;
 }
 
@@ -13712,7 +13773,7 @@ std::shared_ptr<OGRSpatialReference> GDALPamMDArray::GetSpatialRef() const
 {
     if (!m_poPam)
         return nullptr;
-    return m_poPam->GetSpatialRef(GetFullName());
+    return m_poPam->GetSpatialRef(GetFullName(), GetContext());
 }
 
 /************************************************************************/
@@ -13726,9 +13787,9 @@ CPLErr GDALPamMDArray::GetStatistics(bool bApproxOK, bool bForce,
                                      GDALProgressFunc pfnProgress,
                                      void *pProgressData)
 {
-    if (m_poPam &&
-        m_poPam->GetStatistics(GetFullName(), bApproxOK, pdfMin, pdfMax,
-                               pdfMean, pdfStdDev, pnValidCount) == CE_None)
+    if (m_poPam && m_poPam->GetStatistics(GetFullName(), GetContext(),
+                                          bApproxOK, pdfMin, pdfMax, pdfMean,
+                                          pdfStdDev, pnValidCount) == CE_None)
     {
         return CE_None;
     }
@@ -13751,8 +13812,8 @@ bool GDALPamMDArray::SetStatistics(bool bApproxStats, double dfMin,
 {
     if (!m_poPam)
         return false;
-    m_poPam->SetStatistics(GetFullName(), bApproxStats, dfMin, dfMax, dfMean,
-                           dfStdDev, nValidCount);
+    m_poPam->SetStatistics(GetFullName(), GetContext(), bApproxStats, dfMin,
+                           dfMax, dfMean, dfStdDev, nValidCount);
     return true;
 }
 
@@ -13764,7 +13825,7 @@ void GDALPamMDArray::ClearStatistics()
 {
     if (!m_poPam)
         return;
-    m_poPam->ClearStatistics(GetFullName());
+    m_poPam->ClearStatistics(GetFullName(), GetContext());
 }
 
 //! @endcond
diff --git a/gcore/gdalmultidim_gridded.cpp b/gcore/gdalmultidim_gridded.cpp
index a65c936a4699..30cf3b1cd487 100644
--- a/gcore/gdalmultidim_gridded.cpp
+++ b/gcore/gdalmultidim_gridded.cpp
@@ -39,19 +39,6 @@
 #include <limits>
 #include <new>
 
-/************************************************************************/
-/*                          GetPAM()                                    */
-/************************************************************************/
-
-static std::shared_ptr<GDALPamMultiDim>
-GetPAM(const std::shared_ptr<GDALMDArray> &poParent)
-{
-    auto poPamArray = dynamic_cast<GDALPamMDArray *>(poParent.get());
-    if (poPamArray)
-        return poPamArray->GetPAM();
-    return nullptr;
-}
-
 /************************************************************************/
 /*                         GDALMDArrayGridded                           */
 /************************************************************************/
@@ -89,9 +76,9 @@ class GDALMDArrayGridded final : public GDALPamMDArray
         double dfResY, double dfRadius)
         : GDALAbstractMDArray(std::string(),
                               "Gridded view of " + poParent->GetFullName()),
-          GDALPamMDArray(std::string(),
-                         "Gridded view of " + poParent->GetFullName(),
-                         ::GetPAM(poParent)),
+          GDALPamMDArray(
+              std::string(), "Gridded view of " + poParent->GetFullName(),
+              GDALPamMultiDim::GetPAM(poParent), poParent->GetContext()),
           m_poParent(std::move(poParent)), m_apoDims(apoDims), m_poVarX(poVarX),
           m_poVarY(poVarY), m_poVectorDS(std::move(poVectorDS)), m_eAlg(eAlg),
           m_poGridOptions(std::move(poGridOptions)),
diff --git a/gcore/gdalmultidim_subsetdimension.cpp b/gcore/gdalmultidim_subsetdimension.cpp
new file mode 100644
index 000000000000..1440eb64ed5f
--- /dev/null
+++ b/gcore/gdalmultidim_subsetdimension.cpp
@@ -0,0 +1,563 @@
+/******************************************************************************
+ * Name:     gdalmultidim_subsetdimension.cpp
+ * Project:  GDAL Core
+ * Purpose:  GDALGroup::SubsetDimensionFromSelection() implementation
+ * Author:   Even Rouault <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2023, Even Rouault <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include "gdal_priv.h"
+#include "gdal_pam.h"
+
+#include <algorithm>
+
+/************************************************************************/
+/*                           GetParentName()                            */
+/************************************************************************/
+
+static std::string GetParentName(const std::string &osPath)
+{
+    if (osPath == "/" || osPath.rfind('/') == 0)
+        return "/";
+    return osPath.substr(0, osPath.rfind('/'));
+}
+
+/************************************************************************/
+/*                   GDALSubsetGroupSharedResources                     */
+/************************************************************************/
+
+struct GDALSubsetGroupSharedResources
+{
+    std::string m_osDimFullName{};
+    std::vector<int> m_anMapNewDimToOldDim{};
+    std::string m_osSelection{};
+    std::shared_ptr<GDALDimension> m_poNewDim{};
+    std::shared_ptr<GDALMDArray> m_poNewIndexingVar{};
+};
+
+/************************************************************************/
+/*                          CreateContext()                             */
+/************************************************************************/
+
+static std::string
+CreateContext(const std::string &osParentContext,
+              const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared)
+{
+    std::string osRet(osParentContext);
+    if (!osRet.empty())
+        osRet += ". ";
+    osRet += "Selection ";
+    osRet += poShared->m_osSelection;
+    return osRet;
+}
+
+/************************************************************************/
+/*                           GDALSubsetGroup                            */
+/************************************************************************/
+
+class GDALSubsetGroup final : public GDALGroup
+{
+    std::shared_ptr<GDALGroup> m_poParent{};
+    std::shared_ptr<GDALSubsetGroupSharedResources> m_poShared{};
+
+    GDALSubsetGroup(
+        const std::shared_ptr<GDALGroup> &poParent,
+        const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared)
+        : GDALGroup(GetParentName(poParent->GetFullName()), poParent->GetName(),
+                    CreateContext(poParent->GetContext(), poShared)),
+          m_poParent(std::move(poParent)), m_poShared(std::move(poShared))
+    {
+    }
+
+  public:
+    static std::shared_ptr<GDALGroup>
+    Create(const std::shared_ptr<GDALGroup> &poParent,
+           const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared)
+    {
+        auto poGroup = std::shared_ptr<GDALSubsetGroup>(
+            new GDALSubsetGroup(poParent, poShared));
+        poGroup->SetSelf(poGroup);
+        return poGroup;
+    }
+
+    std::vector<std::string>
+    GetMDArrayNames(CSLConstList papszOptions = nullptr) const override
+    {
+        return m_poParent->GetMDArrayNames(papszOptions);
+    }
+
+    std::shared_ptr<GDALMDArray>
+    OpenMDArray(const std::string &osName,
+                CSLConstList papszOptions = nullptr) const override;
+
+    std::vector<std::string>
+    GetGroupNames(CSLConstList papszOptions = nullptr) const override
+    {
+        return m_poParent->GetGroupNames(papszOptions);
+    }
+
+    std::shared_ptr<GDALGroup>
+    OpenGroup(const std::string &osName,
+              CSLConstList papszOptions = nullptr) const override;
+
+    std::vector<std::shared_ptr<GDALDimension>>
+    GetDimensions(CSLConstList papszOptions = nullptr) const override;
+
+    std::shared_ptr<GDALAttribute>
+    GetAttribute(const std::string &osName) const override
+    {
+        return m_poParent->GetAttribute(osName);
+    }
+
+    std::vector<std::shared_ptr<GDALAttribute>>
+    GetAttributes(CSLConstList papszOptions = nullptr) const override
+    {
+        return m_poParent->GetAttributes(papszOptions);
+    }
+};
+
+/************************************************************************/
+/*                           GDALSubsetArray                            */
+/************************************************************************/
+
+class GDALSubsetArray final : public GDALPamMDArray
+{
+  private:
+    std::shared_ptr<GDALMDArray> m_poParent{};
+    std::shared_ptr<GDALSubsetGroupSharedResources> m_poShared{};
+    std::vector<std::shared_ptr<GDALDimension>> m_apoDims{};
+    std::vector<bool> m_abPatchedDim{};
+    bool m_bPatchedDimIsFirst = false;
+
+  protected:
+    GDALSubsetArray(
+        const std::shared_ptr<GDALMDArray> &poParent,
+        const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared,
+        const std::string &osContext)
+        : GDALAbstractMDArray(GetParentName(poParent->GetFullName()),
+                              poParent->GetName()),
+          GDALPamMDArray(GetParentName(poParent->GetFullName()),
+                         poParent->GetName(), GDALPamMultiDim::GetPAM(poParent),
+                         osContext),
+          m_poParent(std::move(poParent)), m_poShared(std::move(poShared))
+    {
+        m_apoDims = m_poParent->GetDimensions();
+        for (size_t i = 0; i < m_apoDims.size(); ++i)
+        {
+            auto &poDim = m_apoDims[i];
+            if (poDim->GetFullName() == m_poShared->m_osDimFullName)
+            {
+                m_bPatchedDimIsFirst = (i == 0);
+                poDim = m_poShared->m_poNewDim;
+                m_abPatchedDim.push_back(true);
+            }
+            else
+            {
+                m_abPatchedDim.push_back(false);
+            }
+        }
+    }
+
+    bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
+               const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
+               const GDALExtendedDataType &bufferDataType,
+               void *pDstBuffer) const override;
+
+  public:
+    static std::shared_ptr<GDALSubsetArray>
+    Create(const std::shared_ptr<GDALMDArray> &poParent,
+           const std::shared_ptr<GDALSubsetGroupSharedResources> &poShared,
+           const std::string &osContext)
+    {
+        auto newAr(std::shared_ptr<GDALSubsetArray>(
+            new GDALSubsetArray(poParent, poShared, osContext)));
+        newAr->SetSelf(newAr);
+        return newAr;
+    }
+
+    bool IsWritable() const override
+    {
+        return false;
+    }
+
+    const std::string &GetFilename() const override
+    {
+        return m_poParent->GetFilename();
+    }
+
+    const std::vector<std::shared_ptr<GDALDimension>> &
+    GetDimensions() const override
+    {
+        return m_apoDims;
+    }
+
+    const GDALExtendedDataType &GetDataType() const override
+    {
+        return m_poParent->GetDataType();
+    }
+
+    const std::string &GetUnit() const override
+    {
+        return m_poParent->GetUnit();
+    }
+
+    std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
+    {
+        return m_poParent->GetSpatialRef();
+    }
+
+    const void *GetRawNoDataValue() const override
+    {
+        return m_poParent->GetRawNoDataValue();
+    }
+
+    std::vector<GUInt64> GetBlockSize() const override
+    {
+        std::vector<GUInt64> ret(m_poParent->GetBlockSize());
+        for (size_t i = 0; i < m_apoDims.size(); ++i)
+        {
+            if (m_abPatchedDim[i])
+                ret[1] = 1;
+        }
+        return ret;
+    }
+
+    std::shared_ptr<GDALAttribute>
+    GetAttribute(const std::string &osName) const override
+    {
+        return m_poParent->GetAttribute(osName);
+    }
+
+    std::vector<std::shared_ptr<GDALAttribute>>
+    GetAttributes(CSLConstList papszOptions = nullptr) const override
+    {
+        return m_poParent->GetAttributes(papszOptions);
+    }
+};
+
+/************************************************************************/
+/*                            OpenMDArray()                             */
+/************************************************************************/
+
+std::shared_ptr<GDALMDArray>
+GDALSubsetGroup::OpenMDArray(const std::string &osName,
+                             CSLConstList papszOptions) const
+{
+    auto poArray = m_poParent->OpenMDArray(osName, papszOptions);
+    if (poArray)
+    {
+        for (const auto &poDim : poArray->GetDimensions())
+        {
+            if (poDim->GetFullName() == m_poShared->m_osDimFullName)
+            {
+                return GDALSubsetArray::Create(poArray, m_poShared,
+                                               GetContext());
+            }
+        }
+    }
+    return poArray;
+}
+
+/************************************************************************/
+/*                             OpenGroup()                              */
+/************************************************************************/
+
+std::shared_ptr<GDALGroup>
+GDALSubsetGroup::OpenGroup(const std::string &osName,
+                           CSLConstList papszOptions) const
+{
+    auto poSubGroup = m_poParent->OpenGroup(osName, papszOptions);
+    if (poSubGroup)
+    {
+        poSubGroup = GDALSubsetGroup::Create(poSubGroup, m_poShared);
+    }
+    return poSubGroup;
+}
+
+/************************************************************************/
+/*                             GetDimensions()                          */
+/************************************************************************/
+
+std::vector<std::shared_ptr<GDALDimension>>
+GDALSubsetGroup::GetDimensions(CSLConstList papszOptions) const
+{
+    auto apoDims = m_poParent->GetDimensions(papszOptions);
+    for (auto &poDim : apoDims)
+    {
+        if (poDim->GetFullName() == m_poShared->m_osDimFullName)
+        {
+            poDim = m_poShared->m_poNewDim;
+        }
+    }
+    return apoDims;
+}
+
+/************************************************************************/
+/*                             IRead()                                  */
+/************************************************************************/
+
+bool GDALSubsetArray::IRead(const GUInt64 *arrayStartIdx, const size_t *count,
+                            const GInt64 *arrayStep,
+                            const GPtrDiff_t *bufferStride,
+                            const GDALExtendedDataType &bufferDataType,
+                            void *pDstBuffer) const
+{
+    const auto nDims = m_apoDims.size();
+    std::vector<GUInt64> newArrayStartIdx(nDims);
+    // the +1 in nDims + 1 is to make happy -Werror=null-dereference when
+    // doing newCount[0] = 1 and newArrayStep[0] = 1
+    std::vector<size_t> newCount(nDims + 1, 1);
+    std::vector<GInt64> newArrayStep(nDims + 1, 1);
+    const size_t nBufferDTSize = bufferDataType.GetSize();
+
+    if (m_bPatchedDimIsFirst)
+    {
+        // Optimized case when the only patched dimension is the first one.
+        std::copy_n(arrayStartIdx, nDims, newArrayStartIdx.data());
+        std::copy_n(count, nDims, newCount.data());
+        std::copy_n(arrayStep, nDims, newArrayStep.data());
+        GUInt64 arrayIdx = arrayStartIdx[0];
+        GByte *pabyDstBuffer = static_cast<GByte *>(pDstBuffer);
+        newCount[0] = 1;
+        newArrayStep[0] = 1;
+        for (size_t i = 0; i < count[0]; ++i)
+        {
+            newArrayStartIdx[0] =
+                m_poShared->m_anMapNewDimToOldDim[static_cast<int>(arrayIdx)];
+            if (!m_poParent->Read(newArrayStartIdx.data(), newCount.data(),
+                                  newArrayStep.data(), bufferStride,
+                                  bufferDataType, pabyDstBuffer))
+            {
+                return false;
+            }
+            if (arrayStep[0] > 0)
+                arrayIdx += arrayStep[0];
+            else
+                arrayIdx -= -arrayStep[0];
+            pabyDstBuffer += bufferStride[0] * nBufferDTSize;
+        }
+        return true;
+    }
+
+    // Slow/unoptimized case
+    std::vector<size_t> anStackIter(nDims);
+    std::vector<GUInt64> anStackArrayIdx(nDims);
+    std::vector<GByte *> pabyDstBufferStack(nDims + 1);
+    pabyDstBufferStack[0] = static_cast<GByte *>(pDstBuffer);
+    size_t iDim = 0;
+lbl_next_depth:
+    if (iDim == nDims)
+    {
+        if (!m_poParent->Read(newArrayStartIdx.data(), newCount.data(),
+                              newArrayStep.data(), bufferStride, bufferDataType,
+                              pabyDstBufferStack[iDim]))
+        {
+            return false;
+        }
+    }
+    else
+    {
+        anStackIter[iDim] = 0;
+        anStackArrayIdx[iDim] = arrayStartIdx[iDim];
+        while (true)
+        {
+            if (m_abPatchedDim[iDim])
+            {
+                newArrayStartIdx[iDim] =
+                    m_poShared->m_anMapNewDimToOldDim[static_cast<int>(
+                        anStackArrayIdx[iDim])];
+            }
+            else
+            {
+                newArrayStartIdx[iDim] = anStackArrayIdx[iDim];
+            }
+            ++iDim;
+            pabyDstBufferStack[iDim] = pabyDstBufferStack[iDim - 1];
+            goto lbl_next_depth;
+        lbl_return_to_caller_in_loop:
+            --iDim;
+            ++anStackIter[iDim];
+            if (anStackIter[iDim] == count[iDim])
+                break;
+            if (arrayStep[iDim] > 0)
+                anStackArrayIdx[iDim] += arrayStep[iDim];
+            else
+                anStackArrayIdx[iDim] -= -arrayStep[iDim];
+            pabyDstBufferStack[iDim] += bufferStride[iDim] * nBufferDTSize;
+        }
+    }
+    if (iDim > 0)
+        goto lbl_return_to_caller_in_loop;
+
+    return true;
+}
+
+/************************************************************************/
+/*                   SubsetDimensionFromSelection()                     */
+/************************************************************************/
+
+/** Return a virtual group whose one dimension has been subset according to a
+ * selection.
+ *
+ * The selection criterion is currently restricted to the form
+ * "/path/to/array=numeric_value" (no spaces around equal)
+ *
+ * This is similar to XArray indexing by name and label on a XArray Dataset
+ * using the sel() method.
+ * Cf https://docs.xarray.dev/en/latest/user-guide/indexing.html#quick-overview
+ *
+ * For example on a EMIT L2A product
+ * (https://github.com/nasa/EMIT-Data-Resources/blob/main/python/tutorials/Exploring_EMIT_L2A_Reflectance.ipynb),
+ * this can be used to keep only valid bands with
+ * SubsetDimensionFromSelection("/sensor_band_parameters/good_wavelengths=1")
+ *
+ * This is the same as the C function GDALGroupSubsetDimensionFromSelection().
+ *
+ * @param osSelection Selection criterion.
+ * @return a virtual group, or nullptr in case of error
+ * @since 3.8
+ */
+std::shared_ptr<GDALGroup>
+GDALGroup::SubsetDimensionFromSelection(const std::string &osSelection) const
+{
+    auto self = std::dynamic_pointer_cast<GDALGroup>(m_pSelf.lock());
+    if (!self)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "Driver implementation issue: m_pSelf not set !");
+        return nullptr;
+    }
+
+    const auto nEqualPos = osSelection.find('=');
+    if (nEqualPos == std::string::npos)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for selection");
+        return nullptr;
+    }
+    const auto osArrayName = osSelection.substr(0, nEqualPos);
+    const auto osValue = osSelection.substr(nEqualPos + 1);
+    if (CPLGetValueType(osValue.c_str()) != CPL_VALUE_INTEGER &&
+        CPLGetValueType(osValue.c_str()) != CPL_VALUE_REAL)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "Non-numeric value in selection criterion");
+        return nullptr;
+    }
+    auto poArray = OpenMDArrayFromFullname(osArrayName);
+    if (!poArray)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined, "Cannot find array %s",
+                 osArrayName.c_str());
+        return nullptr;
+    }
+    if (poArray->GetDimensionCount() != 1)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "Array %s is not single dimensional", osArrayName.c_str());
+        return nullptr;
+    }
+    if (poArray->GetDataType().GetClass() != GEDTC_NUMERIC)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined, "Array %s is not of numeric type",
+                 osArrayName.c_str());
+        return nullptr;
+    }
+
+    const auto nElts = poArray->GetTotalElementsCount();
+    if (nElts > 10 * 1024 * 1024)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined, "Too many values in %s",
+                 osArrayName.c_str());
+        return nullptr;
+    }
+    std::vector<double> values;
+    try
+    {
+        values.resize(static_cast<size_t>(nElts));
+    }
+    catch (const std::bad_alloc &e)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined, "Out of memory: %s", e.what());
+        return nullptr;
+    }
+    const GUInt64 startIdx[1] = {0};
+    const size_t count[1] = {values.size()};
+    if (!poArray->Read(startIdx, count, nullptr, nullptr,
+                       GDALExtendedDataType::Create(GDT_Float64), &values[0],
+                       values.data(), values.size() * sizeof(values[0])))
+    {
+        return nullptr;
+    }
+    const double dfSelectionValue = CPLAtof(osValue.c_str());
+    std::vector<int> anMapNewDimToOldDim;
+    for (int i = 0; i < static_cast<int>(nElts); ++i)
+    {
+        if (values[i] == dfSelectionValue)
+            anMapNewDimToOldDim.push_back(i);
+    }
+    if (anMapNewDimToOldDim.empty())
+    {
+        CPLError(CE_Failure, CPLE_AppDefined, "No value in %s matching %f",
+                 osArrayName.c_str(), dfSelectionValue);
+        return nullptr;
+    }
+    if (anMapNewDimToOldDim.size() == nElts)
+    {
+        return self;
+    }
+
+    auto poDim = poArray->GetDimensions()[0];
+    auto poShared = std::make_shared<GDALSubsetGroupSharedResources>();
+    poShared->m_osSelection = osSelection;
+    poShared->m_osDimFullName = poArray->GetDimensions()[0]->GetFullName();
+    poShared->m_anMapNewDimToOldDim = std::move(anMapNewDimToOldDim);
+
+    // Create a modified dimension of reduced size
+    auto poNewDim = std::make_shared<GDALDimensionWeakIndexingVar>(
+        GetParentName(poDim->GetFullName()), poDim->GetName(), poDim->GetType(),
+        poDim->GetDirection(), poShared->m_anMapNewDimToOldDim.size());
+    poShared->m_poNewDim = poNewDim;
+
+    auto poIndexingVar = poDim->GetIndexingVariable();
+    if (poIndexingVar)
+    {
+        // poNewIndexingVar must be created with a different GDALSubsetGroupSharedResources
+        // instance than poShared, to avoid cross reference, that would result in
+        // objects not being freed !
+        auto poSpecificShared =
+            std::make_shared<GDALSubsetGroupSharedResources>();
+        poSpecificShared->m_osSelection = osSelection;
+        poSpecificShared->m_osDimFullName =
+            poArray->GetDimensions()[0]->GetFullName();
+        poSpecificShared->m_anMapNewDimToOldDim =
+            std::move(anMapNewDimToOldDim);
+        poSpecificShared->m_poNewDim = poNewDim;
+        auto poNewIndexingVar =
+            GDALSubsetArray::Create(poIndexingVar, poSpecificShared,
+                                    CreateContext(GetContext(), poShared));
+        poNewDim->SetIndexingVariable(poNewIndexingVar);
+        poShared->m_poNewIndexingVar = poNewIndexingVar;
+    }
+
+    return GDALSubsetGroup::Create(self, poShared);
+}
diff --git a/swig/include/MultiDimensional.i b/swig/include/MultiDimensional.i
index fa5268c0aee7..753c2722ae11 100644
--- a/swig/include/MultiDimensional.i
+++ b/swig/include/MultiDimensional.i
@@ -250,6 +250,12 @@ public:
     return GDALGroupRename( self, newName ) ? CE_None : CE_Failure;
   }
 
+%newobject SubsetDimensionFromSelection;
+  GDALGroupHS *SubsetDimensionFromSelection( const char *selection,
+                                             char **options = 0 ) {
+    return GDALGroupSubsetDimensionFromSelection(self, selection, options);
+  }
+
 } /* extend */
 }; /* GDALGroupH */