Add S102 raster read-only driver for S-102 bathymetric products (depe…

…nds on libhdf5)

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,93 @@
+#!/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)
+
+    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,177 @@
+#!/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")

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

doc/source/drivers/raster/s102.rst

@@ -0,0 +1,65 @@
+.. _raster.s102:
+
+================================================================================
+S102 -- S-102 Bathymetric Surface Product
+================================================================================
+
+.. shortname:: S102
+
+.. build_dependencies:: libhdf5
+
+.. versionadded:: 3.8
+
+This driver provides read-only support for bathymetry data in the S-102 format,
+which is a specific product profile in an HDF5 file
+
+S-102 files have two image bands representing depth (band 1),
+uncertainty (band 2) values for each cell in a raster grid area.
+
+Note that positive values of depth mean values *below* the reference surface
+of the vertical datum. The :oo:`DEPTH_OR_ELEVATION` open option can be set
+to ``ELEVATION`` to expose depth values as elevation values, by negating their sign
+(i.e. positive values of elevation mean values above the reference surface)
+
+Georeferencing is reported.
+
+Nodata, minimum and maximum values for each band are also reported.
+
+Driver capabilities
+-------------------
+
+.. supports_georeferencing::
+
+.. supports_virtualio::
+
+Open options
+------------
+
+- .. oo:: DEPTH_OR_ELEVATION
+     :choices: DEPTH, ELEVATION
+     :default: DEPTH
+
+     Whether to report depth or elevation.
+     Positive values of depth mean values *below* the reference surface of the
+     vertical datum.
+     Positive values of elevation mean values *above* the reference surface of the
+     vertical datum (which is the convention used by the :ref:`BAG driver <raster.bag>`)
+
+- .. oo:: NORTH_UP
+     :choices: YES, NO
+     :default: YES
+
+     Whether the top line of the dataset should be the northern-most one.
+
+     This is the default behavior of most GDAL formats, but the native
+     organization of the data in S-102 products is to have the first line of
+     the grid being the southern-most one. This native organization can be
+     exposed by the driver by setting this option to NO (in which case the
+     6th term of the geotransform matrix will be positive)
+
+See Also
+--------
+
+-  Implemented as :source_file:`frmts/hdf5/s102dataset.cpp`.
+-  `S-102 Bathymetric Surface Product Specification <https://registry.iho.int/productspec/view.do?idx=199&product_ID=S-102&statusS=5&domainS=20&category=product_ID&searchValue=>`__
+-  :ref:`BAG driver <raster.bag>`

frmts/drivers.ini

@@ -136,6 +136,7 @@ USGSDEM
 GXF
 KEA
 BAG
+S102
 HDF5
 HDF5Image
 NWT_GRD