S102: add multidimensional support

autotest/gdrivers/data/s102/generate_test.py

@@ -63,6 +63,8 @@ def generate(filename, version):
     BathymetryCoverage_01.attrs["gridOriginLatitude"] = np.float64(48)
     BathymetryCoverage_01.attrs["gridSpacingLongitudinal"] = np.float64(0.4)
     BathymetryCoverage_01.attrs["gridSpacingLatitudinal"] = np.float64(0.5)
+    BathymetryCoverage_01.attrs["numPointsLongitudinal"] = np.uint32(values.shape[1])
+    BathymetryCoverage_01.attrs["numPointsLatitudinal"] = np.uint32(values.shape[0])
 
     f.create_group("Group_F")
 

autotest/gdrivers/s102.py

@@ -175,3 +175,30 @@ def test_s102_identify_fallback_through_HDF5_driver():
     assert ds.GetDriver().GetDescription() == "S102"
     del ds
     assert not os.path.exists("data/s102/test_s102_v2.1.h5.aux.xml")
+
+
+###############################################################################
+
+
+def test_s102_multidim():
+    ds = gdal.OpenEx("data/s102/test_s102_v2.1.h5", gdal.OF_MULTIDIM_RASTER)
+    rg = ds.GetRootGroup()
+    ar = rg.OpenMDArrayFromFullname(
+        "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values"
+    )
+    assert ar.GetSpatialRef().GetAuthorityCode(None) == "4326"
+
+    assert ar.GetDimensions()[0].GetName() == "Y"
+    y = ar.GetDimensions()[0].GetIndexingVariable()
+    y_data = struct.unpack("d" * y.GetDimensions()[0].GetSize(), y.Read())
+    assert y_data[0] == 48.0
+    assert y_data[-1] == 48.5
+
+    assert ar.GetDimensions()[1].GetName() == "X"
+    x = ar.GetDimensions()[1].GetIndexingVariable()
+    x_data = struct.unpack("d" * x.GetDimensions()[0].GetSize(), x.Read())
+    assert x_data[0] == 2.0
+    assert x_data[-1] == 2.8
+
+    # Check that it doesn't go into infinite recursion
+    gdal.MultiDimInfo(ds)

frmts/hdf5/CMakeLists.txt

@@ -10,6 +10,8 @@ set(SOURCE
     bagdataset.cpp
     hdf5multidim.cpp
     hdf5eosparser.cpp
+    s100.cpp
+    s100.h
     s102dataset.cpp
 )
 

frmts/hdf5/hdf5dataset.h

@@ -131,8 +131,7 @@ class HDF5SharedResources
         m_oMapEOSGridNameToDimensions{};
     std::map<std::string, std::vector<std::shared_ptr<GDALDimension>>>
         m_oMapEOSSwathNameToDimensions{};
-    std::map<std::string, std::shared_ptr<GDALMDArrayRegularlySpaced>>
-        m_oRefKeeper;
+    std::map<std::string, std::shared_ptr<GDALMDArray>> m_oRefKeeper;
 
     explicit HDF5SharedResources(const std::string &osFilename);
 
@@ -199,7 +198,7 @@ class HDF5SharedResources
         return oIter->second;
     }
 
-    void KeepRef(const std::shared_ptr<GDALMDArrayRegularlySpaced> &poArray)
+    void KeepRef(const std::shared_ptr<GDALMDArray> &poArray)
     {
         m_oRefKeeper[poArray->GetFullName()] = poArray;
     }

frmts/hdf5/hdf5multidim.cpp

@@ -27,6 +27,7 @@
 
 #include "hdf5dataset.h"
 #include "hdf5eosparser.h"
+#include "s100.h"
 
 #include <algorithm>
 #include <set>
@@ -1233,6 +1234,52 @@ void HDF5Array::InstantiateDimensions(const std::string &osParentName,
             }
             return;
         }
+
+        // Special case for S102
+        if (nDims == 2 &&
+            GetFullName() ==
+                "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values")
+        {
+            auto poRootGroup = m_poShared->GetRootGroup();
+            if (poRootGroup)
+            {
+                m_poSRS = std::make_shared<OGRSpatialReference>();
+                if (S100ReadSRS(poRootGroup.get(), *(m_poSRS.get())))
+                {
+                    if (m_poSRS->GetDataAxisToSRSAxisMapping() ==
+                        std::vector<int>{2, 1})
+                        m_poSRS->SetDataAxisToSRSAxisMapping({1, 2});
+                    else
+                        m_poSRS->SetDataAxisToSRSAxisMapping({2, 1});
+                }
+                else
+                {
+                    m_poSRS.reset();
+                }
+
+                auto poBathymetryCoverage01 =
+                    poRootGroup->OpenGroupFromFullname(
+                        "/BathymetryCoverage/BathymetryCoverage.01");
+                if (poBathymetryCoverage01)
+                {
+                    std::vector<std::shared_ptr<GDALMDArray>> apoIndexingVars;
+                    if (S100GetDimensions(poBathymetryCoverage01.get(), m_dims,
+                                          apoIndexingVars) &&
+                        m_dims.size() == 2 &&
+                        m_dims[0]->GetSize() == anDimSizes[0] &&
+                        m_dims[1]->GetSize() == anDimSizes[1])
+                    {
+                        for (const auto &poIndexingVar : apoIndexingVars)
+                            m_poShared->KeepRef(poIndexingVar);
+                        return;
+                    }
+                    else
+                    {
+                        m_dims.clear();
+                    }
+                }
+            }
+        }
     }
 
     std::map<std::string, std::shared_ptr<GDALDimension>> oMapFullNameToDim;

frmts/hdf5/s100.cpp

@@ -0,0 +1,183 @@
+/******************************************************************************
+ *
+ * Project:  Hierarchical Data Format Release 5 (HDF5)
+ * Purpose:  Read S100 bathymetric datasets.
+ * Author:   Even Rouault <even dot rouault at spatialys dot com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys dot 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 "s100.h"
+
+/************************************************************************/
+/*                            S100ReadSRS()                             */
+/************************************************************************/
+
+bool S100ReadSRS(const GDALGroup *poRootGroup, OGRSpatialReference &oSRS)
+{
+    // Get SRS
+    auto poHorizontalCRS = poRootGroup->GetAttribute("horizontalCRS");
+    if (poHorizontalCRS &&
+        poHorizontalCRS->GetDataType().GetClass() == GEDTC_NUMERIC)
+    {
+        // horizontalCRS is v2.2
+        oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+        if (oSRS.importFromEPSG(poHorizontalCRS->ReadAsInt()) != OGRERR_NONE)
+        {
+            oSRS.Clear();
+        }
+    }
+    else
+    {
+        auto poHorizontalDatumReference =
+            poRootGroup->GetAttribute("horizontalDatumReference");
+        auto poHorizontalDatumValue =
+            poRootGroup->GetAttribute("horizontalDatumValue");
+        if (poHorizontalDatumReference && poHorizontalDatumValue)
+        {
+            const char *pszAuthName =
+                poHorizontalDatumReference->ReadAsString();
+            const char *pszAuthCode = poHorizontalDatumValue->ReadAsString();
+            if (pszAuthName && pszAuthCode)
+            {
+                oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+                if (oSRS.SetFromUserInput(
+                        (std::string(pszAuthName) + ':' + pszAuthCode).c_str(),
+                        OGRSpatialReference::
+                            SET_FROM_USER_INPUT_LIMITATIONS_get()) !=
+                    OGRERR_NONE)
+                {
+                    oSRS.Clear();
+                }
+            }
+        }
+    }
+    return !oSRS.IsEmpty();
+}
+
+/************************************************************************/
+/*                        S100GetGeoTransform()                         */
+/************************************************************************/
+
+bool S100GetGeoTransform(const GDALGroup *poGroup, double adfGeoTransform[6],
+                         bool bNorthUp)
+{
+    auto poOriginX = poGroup->GetAttribute("gridOriginLongitude");
+    auto poOriginY = poGroup->GetAttribute("gridOriginLatitude");
+    auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
+    auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
+    auto poNumPointsLongitudinal =
+        poGroup->GetAttribute("numPointsLongitudinal");
+    auto poNumPointsLatitudinal = poGroup->GetAttribute("numPointsLatitudinal");
+    if (poOriginX &&
+        poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 &&
+        poOriginY &&
+        poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 &&
+        poSpacingX &&
+        poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
+        poSpacingY &&
+        poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64 &&
+        poNumPointsLongitudinal &&
+        GDALDataTypeIsInteger(
+            poNumPointsLongitudinal->GetDataType().GetNumericDataType()) &&
+        poNumPointsLatitudinal &&
+        GDALDataTypeIsInteger(
+            poNumPointsLatitudinal->GetDataType().GetNumericDataType()))
+    {
+        adfGeoTransform[0] = poOriginX->ReadAsDouble();
+        adfGeoTransform[3] =
+            poOriginY->ReadAsDouble() +
+            (bNorthUp ? poSpacingY->ReadAsDouble() *
+                            (poNumPointsLatitudinal->ReadAsInt() - 1)
+                      : 0);
+        adfGeoTransform[1] = poSpacingX->ReadAsDouble();
+        adfGeoTransform[5] =
+            bNorthUp ? -poSpacingY->ReadAsDouble() : poSpacingY->ReadAsDouble();
+
+        // From pixel-center convention to pixel-corner convention
+        adfGeoTransform[0] -= adfGeoTransform[1] / 2;
+        adfGeoTransform[3] -= adfGeoTransform[5] / 2;
+
+        return true;
+    }
+    return false;
+}
+
+/************************************************************************/
+/*                        S100GetDimensions()                           */
+/************************************************************************/
+
+bool S100GetDimensions(
+    const GDALGroup *poGroup,
+    std::vector<std::shared_ptr<GDALDimension>> &apoDims,
+    std::vector<std::shared_ptr<GDALMDArray>> &apoIndexingVars)
+{
+    auto poOriginX = poGroup->GetAttribute("gridOriginLongitude");
+    auto poOriginY = poGroup->GetAttribute("gridOriginLatitude");
+    auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
+    auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
+    auto poNumPointsLongitudinal =
+        poGroup->GetAttribute("numPointsLongitudinal");
+    auto poNumPointsLatitudinal = poGroup->GetAttribute("numPointsLatitudinal");
+    if (poOriginX &&
+        poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 &&
+        poOriginY &&
+        poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 &&
+        poSpacingX &&
+        poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
+        poSpacingY &&
+        poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64 &&
+        poNumPointsLongitudinal &&
+        GDALDataTypeIsInteger(
+            poNumPointsLongitudinal->GetDataType().GetNumericDataType()) &&
+        poNumPointsLatitudinal &&
+        GDALDataTypeIsInteger(
+            poNumPointsLatitudinal->GetDataType().GetNumericDataType()))
+    {
+        {
+            auto poDim = std::make_shared<GDALDimensionWeakIndexingVar>(
+                std::string(), "Y", GDAL_DIM_TYPE_HORIZONTAL_Y, std::string(),
+                poNumPointsLatitudinal->ReadAsInt());
+            auto poIndexingVar = GDALMDArrayRegularlySpaced::Create(
+                std::string(), poDim->GetName(), poDim,
+                poOriginY->ReadAsDouble(), poSpacingY->ReadAsDouble(), 0);
+            poDim->SetIndexingVariable(poIndexingVar);
+            apoDims.emplace_back(poDim);
+            apoIndexingVars.emplace_back(poIndexingVar);
+        }
+
+        {
+            auto poDim = std::make_shared<GDALDimensionWeakIndexingVar>(
+                std::string(), "X", GDAL_DIM_TYPE_HORIZONTAL_X, std::string(),
+                poNumPointsLongitudinal->ReadAsInt());
+            auto poIndexingVar = GDALMDArrayRegularlySpaced::Create(
+                std::string(), poDim->GetName(), poDim,
+                poOriginX->ReadAsDouble(), poSpacingX->ReadAsDouble(), 0);
+            poDim->SetIndexingVariable(poIndexingVar);
+            apoDims.emplace_back(poDim);
+            apoIndexingVars.emplace_back(poIndexingVar);
+        }
+
+        return true;
+    }
+    return false;
+}

frmts/hdf5/s100.h

@@ -0,0 +1,47 @@
+/******************************************************************************
+ *
+ * Project:  Hierarchical Data Format Release 5 (HDF5)
+ * Purpose:  Read S100 bathymetric datasets.
+ * Author:   Even Rouault <even dot rouault at spatialys dot com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys dot 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.
+ ****************************************************************************/
+
+#ifndef S100_H
+#define S100_H
+
+#include "cpl_port.h"
+
+#include "gdal_priv.h"
+#include "ogr_spatialref.h"
+
+bool S100ReadSRS(const GDALGroup *poRootGroup, OGRSpatialReference &oSRS);
+
+bool S100GetDimensions(
+    const GDALGroup *poGroup,
+    std::vector<std::shared_ptr<GDALDimension>> &apoDims,
+    std::vector<std::shared_ptr<GDALMDArray>> &apoIndexingVars);
+
+bool S100GetGeoTransform(const GDALGroup *poGroup, double adfGeoTransform[6],
+                         bool bNorthUp);
+
+#endif  // S100_H