Merge pull request OSGeo#8569 from rouault/s102

Add S102 raster read-only driver for S-102 bathymetric products (depends on libhdf5)

apps/gdalmdiminfo_lib.cpp

@@ -692,10 +692,14 @@ DumpDimensions(const std::shared_ptr<GDALGroup> &rootGroup,
             }
             else
             {
+                std::set<std::string> alreadyDumpedDimensionsLocal(
+                    alreadyDumpedDimensions);
+                alreadyDumpedDimensionsLocal.insert(osFullname);
+
                 auto indexingVariableContext(serializer.MakeObjectContext());
                 serializer.AddObjKey(poIndexingVariable->GetName());
                 DumpArray(rootGroup, poIndexingVariable, serializer, psOptions,
-                          alreadyDumpedDimensions,
+                          alreadyDumpedDimensionsLocal,
                           /* bOutputObjType = */ false,
                           /* bOutputName = */ false);
             }

autotest/gdrivers/data/s102/MD_test_s102_v2.1.xml

@@ -0,0 +1 @@
+<nothing/>

autotest/gdrivers/data/s102/MD_test_s102_v2.2.xml

@@ -0,0 +1 @@
+<nothing/>

autotest/gdrivers/data/s102/generate_test.py

@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+###############################################################################
+# $Id$
+#
+# Project:  GDAL/OGR Test Suite
+# Purpose:  Generate test_s102.h5
+# Author:   Even Rouault <even dot rouault at spatialys.com>
+#
+###############################################################################
+# Copyright (c) 2023, Even Rouault <even dot 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.
+###############################################################################
+
+import os
+
+import h5py
+import numpy as np
+
+
+def generate(filename, version):
+    f = h5py.File(os.path.join(os.path.dirname(__file__), f"{filename}.h5"), "w")
+    BathymetryCoverage = f.create_group("BathymetryCoverage")
+    BathymetryCoverage_01 = BathymetryCoverage.create_group("BathymetryCoverage.01")
+    Group_001 = BathymetryCoverage_01.create_group("Group_001")
+
+    values_struct_type = np.dtype(
+        [
+            ("depth", "f4"),
+            ("uncertainty", "f4"),
+        ]
+    )
+    values = Group_001.create_dataset("values", (2, 3), dtype=values_struct_type)
+    data = np.array(
+        [(0, 100), (1, 101), (2, 102), (1e6, 103), (4, 1e6), (5, 105)],
+        dtype=values_struct_type,
+    ).reshape(values.shape)
+    values[...] = data
+
+    Group_001.attrs["minimumUncertainty"] = np.float32(100)
+    Group_001.attrs["maximumUncertainty"] = np.float32(105)
+    Group_001.attrs["minimumDepth"] = np.float32(0)
+    Group_001.attrs["maximumDepth"] = np.float32(5)
+
+    BathymetryCoverage_01.attrs["gridOriginLongitude"] = np.float64(2)
+    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")
+
+    f.attrs["issueDate"] = "2023-12-31"
+    f.attrs["geographicIdentifier"] = "Somewhere"
+    f.attrs["verticalDatum"] = np.int16(12)
+    if version == "INT.IHO.S-102.2.2":
+        f.attrs["horizontalCRS"] = np.int32(4326)
+        f.attrs["verticalCS"] = np.int32(6498)  # Depth, metres, orientation down
+        f.attrs["verticalCoordinateBase"] = np.int32(2)
+        f.attrs["verticalDatumReference"] = np.int32(1)
+    else:
+        assert version == "INT.IHO.S-102.2.1"
+        f.attrs["horizontalDatumReference"] = "EPSG"
+        f.attrs["horizontalDatumValue"] = np.int32(4326)
+    f.attrs["productSpecification"] = version
+    f.attrs[
+        "producer"
+    ] = "Generated by autotest/gdrivers/data/s102/generate_test.py (not strictly fully S102 compliant)"
+    f.attrs["metadata"] = f"MD_{filename}.xml"
+
+    open(os.path.join(os.path.dirname(__file__), f.attrs["metadata"]), "wb").write(
+        b"<nothing/>"
+    )
+
+
+generate("test_s102_v2.1", "INT.IHO.S-102.2.1")
+generate("test_s102_v2.2", "INT.IHO.S-102.2.2")

autotest/gdrivers/s102.py

@@ -0,0 +1,204 @@
+#!/usr/bin/env pytest
+# -*- coding: utf-8 -*-
+###############################################################################
+# $Id$
+#
+# Project:  GDAL/OGR Test Suite
+# Purpose:  Test read functionality for S102 driver.
+# Author:   Even Rouault <even dot rouault at spatialys.com>
+#
+###############################################################################
+# Copyright (c) 2023, Even Rouault <even dot 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.
+###############################################################################
+
+import os
+import struct
+
+import gdaltest
+import pytest
+
+from osgeo import gdal
+
+pytestmark = pytest.mark.require_driver("S102")
+
+
+###############################################################################
+
+
+@pytest.mark.parametrize(
+    "filename", ["data/s102/test_s102_v2.1.h5", "data/s102/test_s102_v2.2.h5"]
+)
+def test_s102_basic(filename):
+    ds = gdal.Open(filename)
+    assert ds.RasterCount == 2
+    assert ds.RasterXSize == 3
+    assert ds.RasterYSize == 2
+    assert ds.GetSpatialRef().GetAuthorityCode(None) == "4326"
+    assert ds.GetGeoTransform() == pytest.approx((1.8, 0.4, 0.0, 48.75, 0.0, -0.5))
+    assert ds.GetMetadata_Dict() == {
+        "AREA_OR_POINT": "Point",
+        "VERTICAL_DATUM_ABBREV": "MLLW",
+        "VERTICAL_DATUM_MEANING": "meanLowerLowWater",
+        "geographicIdentifier": "Somewhere",
+        "issueDate": "2023-12-31",
+        "producer": "Generated by autotest/gdrivers/data/s102/generate_test.py (not strictly fully S102 compliant)",
+    }
+
+    depth = ds.GetRasterBand(1)
+    assert depth.GetDescription() == "depth"
+    assert depth.GetNoDataValue() == 1e6
+    assert depth.GetMinimum() == 0
+    assert depth.GetMaximum() == 5
+    assert depth.GetUnitType() == "metre"
+    assert struct.unpack("f" * 6, depth.ReadRaster()) == (1e6, 4, 5, 0, 1, 2)
+
+    uncertainty = ds.GetRasterBand(2)
+    assert uncertainty.GetDescription() == "uncertainty"
+    assert uncertainty.GetNoDataValue() == 1e6
+    assert uncertainty.GetMinimum() == 100
+    assert uncertainty.GetMaximum() == 105
+    assert uncertainty.GetUnitType() == "metre"
+    assert struct.unpack("f" * 6, uncertainty.ReadRaster()) == (
+        103,
+        1e6,
+        105,
+        100,
+        101,
+        102,
+    )
+
+    assert "MD_" in ds.GetFileList()[1]
+
+    del ds
+    assert not os.path.exists(f"{filename}.aux.xml")
+
+
+###############################################################################
+
+
+def test_s102_elevation():
+    ds = gdal.OpenEx(
+        "data/s102/test_s102_v2.1.h5", open_options=["DEPTH_OR_ELEVATION=ELEVATION"]
+    )
+    assert ds.RasterCount == 2
+    assert ds.RasterXSize == 3
+    assert ds.RasterYSize == 2
+    assert ds.GetSpatialRef().GetAuthorityCode(None) == "4326"
+    assert ds.GetGeoTransform() == pytest.approx((1.8, 0.4, 0.0, 48.75, 0.0, -0.5))
+
+    elevation = ds.GetRasterBand(1)
+    assert elevation.GetDescription() == "elevation"
+    assert elevation.GetNoDataValue() == 1e6
+    assert elevation.GetMinimum() == -5
+    assert elevation.GetMaximum() == 0
+    assert struct.unpack("f" * 6, elevation.ReadRaster()) == (1e6, -4, -5, 0, -1, -2)
+
+    uncertainty = ds.GetRasterBand(2)
+    assert uncertainty.GetDescription() == "uncertainty"
+    assert uncertainty.GetNoDataValue() == 1e6
+    assert uncertainty.GetMinimum() == 100
+    assert uncertainty.GetMaximum() == 105
+    assert struct.unpack("f" * 6, uncertainty.ReadRaster()) == (
+        103,
+        1e6,
+        105,
+        100,
+        101,
+        102,
+    )
+    del ds
+    assert not os.path.exists("data/s102/test_s102_v2.1.h5.aux.xml")
+
+
+###############################################################################
+
+
+def test_s102_north_up_no():
+    ds = gdal.OpenEx("data/s102/test_s102_v2.1.h5", open_options=["NORTH_UP=NO"])
+    assert ds.RasterCount == 2
+    assert ds.RasterXSize == 3
+    assert ds.RasterYSize == 2
+    assert ds.GetSpatialRef().GetAuthorityCode(None) == "4326"
+    assert ds.GetGeoTransform() == pytest.approx((1.8, 0.4, 0.0, 47.75, 0.0, 0.5))
+
+    depth = ds.GetRasterBand(1)
+    assert depth.GetDescription() == "depth"
+    assert depth.GetNoDataValue() == 1e6
+    assert depth.GetMinimum() == 0
+    assert depth.GetMaximum() == 5
+    assert struct.unpack("f" * 6, depth.ReadRaster()) == (0, 1, 2, 1e6, 4, 5)
+
+    uncertainty = ds.GetRasterBand(2)
+    assert uncertainty.GetDescription() == "uncertainty"
+    assert uncertainty.GetNoDataValue() == 1e6
+    assert uncertainty.GetMinimum() == 100
+    assert uncertainty.GetMaximum() == 105
+    assert struct.unpack("f" * 6, uncertainty.ReadRaster()) == (
+        100,
+        101,
+        102,
+        103,
+        1e6,
+        105,
+    )
+    del ds
+    assert not os.path.exists("data/s102/test_s102_v2.1.h5.aux.xml")
+
+
+###############################################################################
+
+
+@pytest.mark.require_driver("HDF5")
+def test_s102_identify_fallback_through_HDF5_driver():
+
+    with gdaltest.config_option("GDAL_S102_IDENTIFY", "NO"):
+        ds = gdal.Open("data/s102/test_s102_v2.1.h5")
+    assert ds
+    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)

doc/source/drivers/raster/bag.rst

@@ -474,3 +474,4 @@ See Also
    1.6 <https://github.com/OpenNavigationSurface/BAG/raw/master/docs/BAG_FSD_Release_1.6.3.doc>`__
 -  `Variable resolution grid extension for BAG
    files <https://github.com/OpenNavigationSurface/BAG/raw/master/docs/VariableResolution/2017-08-10_VariableResolution.docx>`__
+-  :ref:`S-102 driver <raster.s102>`

doc/source/drivers/raster/index.rst

@@ -157,6 +157,7 @@ Raster drivers
    rpftoc
    rraster
    rs2
+   s102
    safe
    sar_ceos
    sdat