diff --git a/autotest/gdrivers/data/s102/MD_test_s102_v2.1.xml b/autotest/gdrivers/data/s102/MD_test_s102_v2.1.xml new file mode 100644 index 000000000000..7d5c8daa45e2 --- /dev/null +++ b/autotest/gdrivers/data/s102/MD_test_s102_v2.1.xml @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/autotest/gdrivers/data/s102/MD_test_s102_v2.2.xml b/autotest/gdrivers/data/s102/MD_test_s102_v2.2.xml new file mode 100644 index 000000000000..7d5c8daa45e2 --- /dev/null +++ b/autotest/gdrivers/data/s102/MD_test_s102_v2.2.xml @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/autotest/gdrivers/data/s102/generate_test.py b/autotest/gdrivers/data/s102/generate_test.py new file mode 100755 index 000000000000..b17b429f98a4 --- /dev/null +++ b/autotest/gdrivers/data/s102/generate_test.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +############################################################################### +# $Id$ +# +# Project: GDAL/OGR Test Suite +# Purpose: Generate test_s102.h5 +# Author: Even Rouault +# +############################################################################### +# Copyright (c) 2023, Even Rouault +# +# 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 "2.2" in version: + 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: + 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"" + ) + + +generate("test_s102_v2.1", "INT.IHO.S-102.2.1") +generate("test_s102_v2.2", "INT.IHO.S-102.2.2") diff --git a/autotest/gdrivers/data/s102/test_s102_v2.1.h5 b/autotest/gdrivers/data/s102/test_s102_v2.1.h5 new file mode 100644 index 000000000000..6bad93ae6fdf Binary files /dev/null and b/autotest/gdrivers/data/s102/test_s102_v2.1.h5 differ diff --git a/autotest/gdrivers/data/s102/test_s102_v2.2.h5 b/autotest/gdrivers/data/s102/test_s102_v2.2.h5 new file mode 100644 index 000000000000..859c80637ef9 Binary files /dev/null and b/autotest/gdrivers/data/s102/test_s102_v2.2.h5 differ diff --git a/autotest/gdrivers/s102.py b/autotest/gdrivers/s102.py new file mode 100755 index 000000000000..fe8ca7fe7a05 --- /dev/null +++ b/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 +# +############################################################################### +# Copyright (c) 2023, Even Rouault +# +# 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") diff --git a/doc/source/drivers/raster/bag.rst b/doc/source/drivers/raster/bag.rst index 8e08b4ef930a..25adc8f79c55 100644 --- a/doc/source/drivers/raster/bag.rst +++ b/doc/source/drivers/raster/bag.rst @@ -474,3 +474,4 @@ See Also 1.6 `__ - `Variable resolution grid extension for BAG files `__ +- :ref:`S-102 driver ` diff --git a/doc/source/drivers/raster/index.rst b/doc/source/drivers/raster/index.rst index 0313f17c6d76..c48f6603e207 100644 --- a/doc/source/drivers/raster/index.rst +++ b/doc/source/drivers/raster/index.rst @@ -157,6 +157,7 @@ Raster drivers rpftoc rraster rs2 + s102 safe sar_ceos sdat diff --git a/doc/source/drivers/raster/s102.rst b/doc/source/drivers/raster/s102.rst new file mode 100644 index 000000000000..25692b7d6f69 --- /dev/null +++ b/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 ``YES`` 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 `) + +- .. 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 `__ +- :ref:`BAG driver ` diff --git a/frmts/drivers.ini b/frmts/drivers.ini index 9672bedd5e53..ec782d9b076a 100644 --- a/frmts/drivers.ini +++ b/frmts/drivers.ini @@ -136,6 +136,7 @@ USGSDEM GXF KEA BAG +S102 HDF5 HDF5Image NWT_GRD diff --git a/frmts/gdalallregister.cpp b/frmts/gdalallregister.cpp index 32d17a192b19..145e6b260897 100644 --- a/frmts/gdalallregister.cpp +++ b/frmts/gdalallregister.cpp @@ -47,14 +47,14 @@ static char *szConfiguredFormats = "GDAL_FORMATS"; /************************************************************************/ /** - * \brief Register a plugin by name, returning an error if not found - * + * \brief Register a plugin by name, returning an error if not found + * * This function will call GDALDriverManager::LoadPlugin() to register a * specific plugin by name. - * + * * This method is intended to be called instead of GDALAllRegister() or * GDALRegisterPlugins() when fine tuning which drivers are needed at runtime. - * + * * @see GDALDriverManager::LoadPlugin() * @see GDALDriverManager::AutoLoadDrivers() * @since GDAL 3.8 @@ -74,15 +74,15 @@ CPLErr GDALRegisterPlugin(const char *name) /** * \brief Register drivers and support code available as a plugin. - * + * * This function will call GDALDriverManager::AutoLoadDrivers() to * register all drivers or supporting code (for example VRT pixelfunctions * or VSI adapters) that have not been compiled into the GDAL core but instead * are made available through GDAL's plugin mechanism. - * + * * This method is intended to be called instead of GDALAllRegister() when * fine tuning which drivers are needed at runtime. - * + * * @see GDALDriverManager::AutoLoadDrivers() * @since GDAL 3.8 */ @@ -513,6 +513,7 @@ void CPL_STDCALL GDALAllRegister() #ifdef FRMT_hdf5 GDALRegister_BAG(); + GDALRegister_S102(); GDALRegister_HDF5(); GDALRegister_HDF5Image(); #endif diff --git a/frmts/hdf5/CMakeLists.txt b/frmts/hdf5/CMakeLists.txt index 5a7115133d1b..100d13567202 100644 --- a/frmts/hdf5/CMakeLists.txt +++ b/frmts/hdf5/CMakeLists.txt @@ -10,6 +10,7 @@ set(SOURCE bagdataset.cpp hdf5multidim.cpp hdf5eosparser.cpp + s102dataset.cpp ) add_gdal_driver(TARGET gdal_HDF5 SOURCES ${SOURCE} PLUGIN_CAPABLE) diff --git a/frmts/hdf5/hdf5dataset.cpp b/frmts/hdf5/hdf5dataset.cpp index cffdbe852f91..ea9230623ca0 100644 --- a/frmts/hdf5/hdf5dataset.cpp +++ b/frmts/hdf5/hdf5dataset.cpp @@ -198,6 +198,7 @@ void GDALRegister_HDF5() #ifdef HDF5_PLUGIN GDALRegister_HDF5Image(); GDALRegister_BAG(); + GDALRegister_S102(); #endif } @@ -644,8 +645,6 @@ GDALDataset *HDF5Dataset::Open(GDALOpenInfo *poOpenInfo) poDS->ReadGlobalAttributes(true); - poDS->SetMetadata(poDS->m_aosMetadata.List()); - if (STARTS_WITH(poDS->m_aosMetadata.FetchNameValueDef("mission_name", ""), "Sentinel 3") && EQUAL( @@ -660,6 +659,22 @@ GDALDataset *HDF5Dataset::Open(GDALOpenInfo *poOpenInfo) return nullptr; } + // Safety belt if S102Dataset::Identify() failed + if (STARTS_WITH( + poDS->m_aosMetadata.FetchNameValueDef("productSpecification", ""), + "INT.IHO.S-102.") && + GDALGetDriverByName("S102") != nullptr) + { + delete poDS; + std::string osS102Filename("S102:\""); + osS102Filename += + CPLString(poOpenInfo->pszFilename).replaceAll("\"", "\\\""); + osS102Filename += '"'; + return GDALDataset::Open(osS102Filename.c_str(), GDAL_OF_RASTER); + } + + poDS->SetMetadata(poDS->m_aosMetadata.List()); + if (CSLCount(poDS->papszSubDatasets) / 2 >= 1) poDS->SetMetadata(poDS->papszSubDatasets, "SUBDATASETS"); diff --git a/frmts/hdf5/hdf5dataset.h b/frmts/hdf5/hdf5dataset.h index d10ef4771f19..6c255ed6577c 100644 --- a/frmts/hdf5/hdf5dataset.h +++ b/frmts/hdf5/hdf5dataset.h @@ -105,6 +105,7 @@ hid_t GDAL_HDF5Open(const std::string &osFilename); class HDF5Dataset; class HDF5EOSParser; class BAGDataset; +class S102Dataset; namespace GDAL { @@ -118,6 +119,7 @@ class HDF5SharedResources { friend class ::HDF5Dataset; friend class ::BAGDataset; + friend class ::S102Dataset; std::weak_ptr m_poSelf{}; bool m_bReadOnly = true; diff --git a/frmts/hdf5/hdf5multidim.cpp b/frmts/hdf5/hdf5multidim.cpp index f8b199326371..c27aa74aa503 100644 --- a/frmts/hdf5/hdf5multidim.cpp +++ b/frmts/hdf5/hdf5multidim.cpp @@ -972,6 +972,21 @@ HDF5Array::HDF5Array(const std::string &osParentName, const std::string &osName, HDF5Array::GetAttributes(); + // Special case for S102 nodata value that is at 1e6 + if (GetFullName() == + "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values" && + m_dt.GetClass() == GEDTC_COMPOUND && + m_dt.GetSize() == 2 * sizeof(float) && + m_dt.GetComponents().size() == 2 && + m_dt.GetComponents()[0]->GetType().GetNumericDataType() == + GDT_Float32 && + m_dt.GetComponents()[1]->GetType().GetNumericDataType() == GDT_Float32) + { + m_abyNoData.resize(m_dt.GetSize()); + float afNoData[2] = {1e6f, 1e6f}; + memcpy(m_abyNoData.data(), afNoData, m_abyNoData.size()); + } + if (bSkipFullDimensionInstantiation) { const int nDims = H5Sget_simple_extent_ndims(m_hDataSpace); diff --git a/frmts/hdf5/s102dataset.cpp b/frmts/hdf5/s102dataset.cpp new file mode 100644 index 000000000000..aa7db5b5b669 --- /dev/null +++ b/frmts/hdf5/s102dataset.cpp @@ -0,0 +1,642 @@ +/****************************************************************************** + * + * Project: Hierarchical Data Format Release 5 (HDF5) + * Purpose: Read S102 bathymetric datasets. + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2023, Even Rouault + * + * 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 "cpl_port.h" +#include "hdf5dataset.h" +#include "gh5_convenience.h" + +#include "gdal_priv.h" +#include "gdal_proxy.h" + +#include + +/************************************************************************/ +/* S102Dataset */ +/************************************************************************/ + +class S102Dataset final : public GDALPamDataset +{ + OGRSpatialReference m_oSRS{}; + bool m_bHasGT = false; + double m_adfGeoTransform[6] = {0, 1, 0, 0, 0, 1}; + std::unique_ptr m_poDepthDS{}; + std::unique_ptr m_poUncertaintyDS{}; + std::string m_osMetadataFile{}; + + public: + S102Dataset() = default; + + CPLErr GetGeoTransform(double *) override; + const OGRSpatialReference *GetSpatialRef() const override; + + char **GetFileList() override; + + static GDALDataset *Open(GDALOpenInfo *); + static int Identify(GDALOpenInfo *); +}; + +/************************************************************************/ +/* S102RasterBand */ +/************************************************************************/ + +class S102RasterBand : public GDALProxyRasterBand +{ + friend class S102Dataset; + GDALRasterBand *m_poUnderlyingBand = nullptr; + double m_dfMinimum = std::numeric_limits::quiet_NaN(); + double m_dfMaximum = std::numeric_limits::quiet_NaN(); + + public: + explicit S102RasterBand(GDALRasterBand *poUnderlyingBand) + : m_poUnderlyingBand(poUnderlyingBand) + { + eDataType = poUnderlyingBand->GetRasterDataType(); + poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize); + } + + GDALRasterBand * + RefUnderlyingRasterBand(bool /*bForceOpen*/ = true) const override + { + return m_poUnderlyingBand; + } + + double GetMinimum(int *pbSuccess = nullptr) override + { + if (pbSuccess) + *pbSuccess = !std::isnan(m_dfMinimum); + return m_dfMinimum; + } + + double GetMaximum(int *pbSuccess = nullptr) override + { + if (pbSuccess) + *pbSuccess = !std::isnan(m_dfMaximum); + return m_dfMaximum; + } + + const char *GetUnitType() override + { + return "metre"; + } +}; + +/************************************************************************/ +/* GetGeoTransform() */ +/************************************************************************/ + +CPLErr S102Dataset::GetGeoTransform(double *padfGeoTransform) + +{ + if (m_bHasGT) + { + memcpy(padfGeoTransform, m_adfGeoTransform, sizeof(double) * 6); + return CE_None; + } + + return GDALPamDataset::GetGeoTransform(padfGeoTransform); +} + +/************************************************************************/ +/* GetSpatialRef() */ +/************************************************************************/ + +const OGRSpatialReference *S102Dataset::GetSpatialRef() const +{ + if (!m_oSRS.IsEmpty()) + return &m_oSRS; + return GDALPamDataset::GetSpatialRef(); +} + +/************************************************************************/ +/* GetFileList() */ +/************************************************************************/ + +char **S102Dataset::GetFileList() +{ + char **papszFileList = GDALPamDataset::GetFileList(); + if (!m_osMetadataFile.empty()) + papszFileList = CSLAddString(papszFileList, m_osMetadataFile.c_str()); + return papszFileList; +} + +/************************************************************************/ +/* Identify() */ +/************************************************************************/ + +int S102Dataset::Identify(GDALOpenInfo *poOpenInfo) + +{ + if (STARTS_WITH(poOpenInfo->pszFilename, "S102:")) + return TRUE; + + // Is it an HDF5 file? + static const char achSignature[] = "\211HDF\r\n\032\n"; + + if (poOpenInfo->pabyHeader == nullptr || + memcmp(poOpenInfo->pabyHeader, achSignature, 8) != 0) + return FALSE; + + // GDAL_S102_IDENTIFY can be set to NO only for tests, to test that + // HDF5Dataset::Open() can redirect to S102 if the below logic fails + if (CPLTestBool(CPLGetConfigOption("GDAL_S102_IDENTIFY", "YES"))) + { + // The below identification logic may be a bit fragile... + // Works at least on: + // - /vsis3/noaa-s102-pds/ed2.1.0/national_bathymetric_source/boston/dcf2/tiles/102US00_US4MA1GC.h5 + // - https://datahub.admiralty.co.uk/portal/sharing/rest/content/items/6fd07bde26124d48820b6dee60695389/data (S-102_Liverpool_Trial_Cells.zip) + const int nLenBC = static_cast(strlen("BathymetryCoverage\0") + 1); + const int nLenGroupF = static_cast(strlen("Group_F\0") + 1); + bool bFoundBathymetryCoverage = false; + bool bFoundGroupF = false; + for (int i = 0; i < poOpenInfo->nHeaderBytes - nLenBC; ++i) + { + if (poOpenInfo->pabyHeader[i] == 'B' && + memcmp(poOpenInfo->pabyHeader + i, "BathymetryCoverage\0", + nLenBC) == 0) + { + bFoundBathymetryCoverage = true; + if (bFoundGroupF) + return true; + } + if (poOpenInfo->pabyHeader[i] == 'G' && + memcmp(poOpenInfo->pabyHeader + i, "Group_F\0", nLenGroupF) == + 0) + { + bFoundGroupF = true; + if (bFoundBathymetryCoverage) + return true; + } + } + } + + return false; +} + +/************************************************************************/ +/* Open() */ +/************************************************************************/ + +GDALDataset *S102Dataset::Open(GDALOpenInfo *poOpenInfo) + +{ + // Confirm that this appears to be a S102 file. + if (!Identify(poOpenInfo)) + return nullptr; + + HDF5_GLOBAL_LOCK(); + + // Confirm the requested access is supported. + if (poOpenInfo->eAccess == GA_Update) + { + CPLError(CE_Failure, CPLE_NotSupported, + "The S102 driver does not support update access."); + return nullptr; + } + + std::string osFilename(poOpenInfo->pszFilename); + if (STARTS_WITH(poOpenInfo->pszFilename, "S102:")) + { + const CPLStringList aosTokens( + CSLTokenizeString2(poOpenInfo->pszFilename, ":", + CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES)); + + if (aosTokens.size() == 2) + { + osFilename = aosTokens[1]; + } + else + { + return nullptr; + } + } + + // Open the file as an HDF5 file. + hid_t fapl = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr); + hid_t hHDF5 = H5Fopen(osFilename.c_str(), H5F_ACC_RDONLY, fapl); + H5Pclose(fapl); + if (hHDF5 < 0) + return nullptr; + + auto poSharedResources = GDAL::HDF5SharedResources::Create(osFilename); + poSharedResources->m_hHDF5 = hHDF5; + + auto poRootGroup = HDF5Dataset::OpenGroup(poSharedResources); + if (!poRootGroup) + return nullptr; + + auto poBathymetryCoverage01 = poRootGroup->OpenGroupFromFullname( + "/BathymetryCoverage/BathymetryCoverage.01"); + if (!poBathymetryCoverage01) + return nullptr; + auto poGroup001 = poBathymetryCoverage01->OpenGroup("Group_001"); + if (!poGroup001) + return nullptr; + auto poValuesArray = poGroup001->OpenMDArray("values"); + if (!poValuesArray || poValuesArray->GetDimensionCount() != 2) + return nullptr; + + const auto &oType = poValuesArray->GetDataType(); + if (oType.GetClass() != GEDTC_COMPOUND) + return nullptr; + const auto &oComponents = oType.GetComponents(); + if (oComponents.size() != 2 || oComponents[0]->GetName() != "depth" || + oComponents[1]->GetName() != "uncertainty") + { + return nullptr; + } + + auto poDS = cpl::make_unique(); + + const bool bNorthUp = CPLTestBool( + CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "NORTH_UP", "YES")); + if (bNorthUp) + poValuesArray = poValuesArray->GetView("[::-1,...]"); + + auto poDepth = poValuesArray->GetView("[\"depth\"]"); + + // Mandatory in v2.2 + bool bCSIsElevation = false; + auto poVerticalCS = poRootGroup->GetAttribute("verticalCS"); + if (poVerticalCS && poVerticalCS->GetDataType().GetClass() == GEDTC_NUMERIC) + { + const auto nVal = poVerticalCS->ReadAsInt(); + if (nVal == 6498) // Depth metre + { + // nothing to do + } + else if (nVal == 6499) // Height metre + { + bCSIsElevation = true; + } + else + { + CPLError(CE_Warning, CPLE_NotSupported, "Unsupported verticalCS=%d", + nVal); + } + } + + const bool bUseElevation = + EQUAL(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, + "DEPTH_OR_ELEVATION", "DEPTH"), + "ELEVATION"); + constexpr double NODATA = 1e6; + const bool bInvertDepth = (bUseElevation && !bCSIsElevation) || + (!bUseElevation && bCSIsElevation); + if (bInvertDepth) + { + auto poInverted = poDepth->GetUnscaled(-1, 0, NODATA); + poDS->m_poDepthDS.reset(poInverted->AsClassicDataset(1, 0)); + } + else + { + poDS->m_poDepthDS.reset(poDepth->AsClassicDataset(1, 0)); + } + + auto poUncertainty = poValuesArray->GetView("[\"uncertainty\"]"); + poDS->m_poUncertaintyDS.reset(poUncertainty->AsClassicDataset(1, 0)); + + poDS->nRasterXSize = poDS->m_poDepthDS->GetRasterXSize(); + poDS->nRasterYSize = poDS->m_poDepthDS->GetRasterYSize(); + + // Create depth (or elevation) band + auto poDepthBand = new S102RasterBand(poDS->m_poDepthDS->GetRasterBand(1)); + poDepthBand->SetDescription(bUseElevation ? "elevation" : "depth"); + + auto poMinimumDepth = poGroup001->GetAttribute("minimumDepth"); + if (poMinimumDepth && + poMinimumDepth->GetDataType().GetClass() == GEDTC_NUMERIC) + { + const double dfVal = poMinimumDepth->ReadAsDouble(); + if (dfVal != NODATA) + { + if (bInvertDepth) + poDepthBand->m_dfMaximum = -dfVal; + else + poDepthBand->m_dfMinimum = dfVal; + } + } + + auto poMaximumDepth = poGroup001->GetAttribute("maximumDepth"); + if (poMaximumDepth && + poMaximumDepth->GetDataType().GetClass() == GEDTC_NUMERIC) + { + const double dfVal = poMaximumDepth->ReadAsDouble(); + if (dfVal != NODATA) + { + if (bInvertDepth) + poDepthBand->m_dfMinimum = -dfVal; + else + poDepthBand->m_dfMaximum = dfVal; + } + } + + poDS->SetBand(1, poDepthBand); + + // Create uncertainty band + auto poUncertaintyBand = + new S102RasterBand(poDS->m_poUncertaintyDS->GetRasterBand(1)); + poUncertaintyBand->SetDescription("uncertainty"); + + auto poMinimumUncertainty = poGroup001->GetAttribute("minimumUncertainty"); + if (poMinimumUncertainty && + poMinimumUncertainty->GetDataType().GetClass() == GEDTC_NUMERIC) + { + const double dfVal = poMinimumUncertainty->ReadAsDouble(); + if (dfVal != NODATA) + { + poUncertaintyBand->m_dfMinimum = dfVal; + } + } + + auto poMaximumUncertainty = poGroup001->GetAttribute("maximumUncertainty"); + if (poMaximumUncertainty && + poMaximumUncertainty->GetDataType().GetClass() == GEDTC_NUMERIC) + { + const double dfVal = poMaximumUncertainty->ReadAsDouble(); + if (dfVal != NODATA) + { + poUncertaintyBand->m_dfMaximum = dfVal; + } + } + + poDS->SetBand(2, poUncertaintyBand); + + // Compute geotransform + auto poOriginX = + poBathymetryCoverage01->GetAttribute("gridOriginLongitude"); + auto poOriginY = poBathymetryCoverage01->GetAttribute("gridOriginLatitude"); + auto poSpacingX = + poBathymetryCoverage01->GetAttribute("gridSpacingLongitudinal"); + auto poSpacingY = + poBathymetryCoverage01->GetAttribute("gridSpacingLatitudinal"); + 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) + { + poDS->m_adfGeoTransform[0] = poOriginX->ReadAsDouble(); + poDS->m_adfGeoTransform[3] = + poOriginY->ReadAsDouble() + + (bNorthUp ? poSpacingY->ReadAsDouble() * (poDS->nRasterYSize - 1) + : 0); + poDS->m_adfGeoTransform[1] = poSpacingX->ReadAsDouble(); + poDS->m_adfGeoTransform[5] = + bNorthUp ? -poSpacingY->ReadAsDouble() : poSpacingY->ReadAsDouble(); + + // From pixel-center convention to pixel-corner convention + poDS->m_adfGeoTransform[0] -= poDS->m_adfGeoTransform[1] / 2; + poDS->m_adfGeoTransform[3] -= poDS->m_adfGeoTransform[5] / 2; + + poDS->m_bHasGT = true; + } + + // Get SRS + auto poHorizontalCRS = poRootGroup->GetAttribute("horizontalCRS"); + if (poHorizontalCRS && + poHorizontalCRS->GetDataType().GetClass() == GEDTC_NUMERIC) + { + // horizontalCRS is v2.2 + poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + if (poDS->m_oSRS.importFromEPSG(poHorizontalCRS->ReadAsInt()) != + OGRERR_NONE) + { + poDS->m_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) + { + poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + if (poDS->m_oSRS.SetFromUserInput( + (std::string(pszAuthName) + ':' + pszAuthCode).c_str(), + OGRSpatialReference:: + SET_FROM_USER_INPUT_LIMITATIONS_get()) != + OGRERR_NONE) + { + poDS->m_oSRS.Clear(); + } + } + } + } + + // https://iho.int/uploads/user/pubs/standards/s-100/S-100_5.0.0_Final_Clean_Web.pdf + // Table S100_VerticalAndSoundingDatum page 20 + static const struct + { + int nCode; + const char *pszMeaning; + const char *pszAbbrev; + } asVerticalDatums[] = { + {1, "meanLowWaterSprings", "MLWS"}, + {2, "meanLowerLowWaterSprings", nullptr}, + {3, "meanSeaLevel", "MSL"}, + {4, "lowestLowWater", nullptr}, + {5, "meanLowWater", "MLW"}, + {6, "lowestLowWaterSprings", nullptr}, + {7, "approximateMeanLowWaterSprings", nullptr}, + {8, "indianSpringLowWater", nullptr}, + {9, "lowWaterSprings", nullptr}, + {10, "approximateLowestAstronomicalTide", nullptr}, + {11, "nearlyLowestLowWater", nullptr}, + {12, "meanLowerLowWater", "MLLW"}, + {13, "lowWater", "LW"}, + {14, "approximateMeanLowWater", nullptr}, + {15, "approximateMeanLowerLowWater", nullptr}, + {16, "meanHighWater", "MHW"}, + {17, "meanHighWaterSprings", "MHWS"}, + {18, "highWater", "HW"}, + {19, "approximateMeanSeaLevel", nullptr}, + {20, "highWaterSprings", nullptr}, + {21, "meanHigherHighWater", "MHHW"}, + {22, "equinoctialSpringLowWater", nullptr}, + {23, "lowestAstronomicalTide", "LAT"}, + {24, "localDatum", nullptr}, + {25, "internationalGreatLakesDatum1985", nullptr}, + {26, "meanWaterLevel", nullptr}, + {27, "lowerLowWaterLargeTide", nullptr}, + {28, "higherHighWaterLargeTide", nullptr}, + {29, "nearlyHighestHighWater", nullptr}, + {30, "highestAstronomicalTide", "HAT"}, + {44, "balticSeaChartDatum2000", nullptr}, + {46, "internationalGreatLakesDatum2020", nullptr}, + }; + + auto poVerticalDatum = poRootGroup->GetAttribute("verticalDatum"); + if (poVerticalDatum && + poVerticalDatum->GetDataType().GetClass() == GEDTC_NUMERIC) + { + bool bFound = false; + const auto nVal = poVerticalDatum->ReadAsInt(); + for (const auto &sVerticalDatum : asVerticalDatums) + { + if (sVerticalDatum.nCode == nVal) + { + bFound = true; + poDS->GDALDataset::SetMetadataItem("VERTICAL_DATUM_MEANING", + sVerticalDatum.pszMeaning); + if (sVerticalDatum.pszAbbrev) + poDS->GDALDataset::SetMetadataItem( + "VERTICAL_DATUM_ABBREV", sVerticalDatum.pszAbbrev); + break; + } + } + if (!bFound) + { + poDS->GDALDataset::SetMetadataItem("verticalDatum", + CPLSPrintf("%d", nVal)); + } + } + + for (const auto &poAttr : poRootGroup->GetAttributes()) + { + const auto &osName = poAttr->GetName(); + if (osName == "metadata") + { + const char *pszVal = poAttr->ReadAsString(); + if (pszVal && pszVal[0]) + { + poDS->m_osMetadataFile = CPLFormFilename( + CPLGetPath(osFilename.c_str()), pszVal, nullptr); + VSIStatBufL sStat; + if (VSIStatL(poDS->m_osMetadataFile.c_str(), &sStat) != 0) + { + // Test products from https://data.admiralty.co.uk/portal/apps/sites/#/marine-data-portal/pages/s-100 + // advertize a metadata filename starting with "MD_", per the spec, + // but the actual filename does not start with "MD_"... + if (STARTS_WITH(pszVal, "MD_")) + { + poDS->m_osMetadataFile = + CPLFormFilename(CPLGetPath(osFilename.c_str()), + pszVal + strlen("MD_"), nullptr); + if (VSIStatL(poDS->m_osMetadataFile.c_str(), &sStat) != + 0) + { + poDS->m_osMetadataFile.clear(); + } + } + else + { + poDS->m_osMetadataFile.clear(); + } + } + } + } + else if (osName != "horizontalCRS" && + osName != "horizontalDatumReference" && + osName != "horizontalDatumValue" && + osName != "productSpecification" && + osName != "eastBoundLongitude" && + osName != "northBoundLatitude" && + osName != "southBoundLatitude" && + osName != "westBoundLongitude" && osName != "extentTypeCode" && + osName != "verticalCS" && osName != "verticalCoordinateBase" && + osName != "verticalDatumReference" && + osName != "verticalDatum") + { + const char *pszVal = poAttr->ReadAsString(); + if (pszVal && pszVal[0]) + poDS->GDALDataset::SetMetadataItem(osName.c_str(), pszVal); + } + } + + poDS->GDALDataset::SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT); + + // Setup/check for pam .aux.xml. + poDS->SetDescription(osFilename.c_str()); + poDS->TryLoadXML(); + + // Setup overviews. + poDS->oOvManager.Initialize(poDS.get(), osFilename.c_str()); + + return poDS.release(); +} + +/************************************************************************/ +/* S102DatasetDriverUnload() */ +/************************************************************************/ + +static void S102DatasetDriverUnload(GDALDriver *) +{ + HDF5UnloadFileDriver(); +} + +/************************************************************************/ +/* GDALRegister_S102() */ +/************************************************************************/ +void GDALRegister_S102() + +{ + if (!GDAL_CHECK_VERSION("S102")) + return; + + if (GDALGetDriverByName("S102") != nullptr) + return; + + GDALDriver *poDriver = new GDALDriver(); + + poDriver->SetDescription("S102"); + poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES"); + poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, + "S-102 Bathymetric Surface Product"); + poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/s102.html"); + poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES"); + poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "h5"); + + poDriver->SetMetadataItem( + GDAL_DMD_OPENOPTIONLIST, + "" + " " + " "); + poDriver->pfnOpen = S102Dataset::Open; + poDriver->pfnIdentify = S102Dataset::Identify; + poDriver->pfnUnloadDriver = S102DatasetDriverUnload; + + GetGDALDriverManager()->RegisterDriver(poDriver); +} diff --git a/gcore/gdal_frmts.h b/gcore/gdal_frmts.h index 1ccfb8130346..4c98ad2c1eec 100644 --- a/gcore/gdal_frmts.h +++ b/gcore/gdal_frmts.h @@ -103,6 +103,7 @@ void CPL_DLL GDALRegister_IDA(void); void CPL_DLL GDALRegister_NDF(void); void CPL_DLL GDALRegister_RMF(void); void CPL_DLL GDALRegister_BAG(void); +void CPL_DLL GDALRegister_S102(void); void CPL_DLL GDALRegister_HDF5(void); void CPL_DLL GDALRegister_HDF5Image(void); void CPL_DLL GDALRegister_MSGN(void);