Skip to content

Commit

Permalink
S102: add support for IHO S102 v3.0 specification (OSGeo#10779)
Browse files Browse the repository at this point in the history
  • Loading branch information
@rouault
rouault committed Sep 16, 2024
1 parent 0db013b commit 79577b1
Show file tree
Hide file tree
Showing 7 changed files with 188 additions and 84 deletions.
1 change: 1 addition & 0 deletions autotest/gdrivers/data/s102/MD_test_s102_v3.0_with_QualityOfBathymetryCoverage.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
<nothing/>
68 changes: 47 additions & 21 deletions autotest/gdrivers/data/s102/generate_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,13 @@
import numpy as np


def generate(filename, version, with_QualityOfSurvey=False):
def generate(
filename,
version,
with_QualityOfSurvey=False,
with_QualityOfCoverage=False,
use_compound_type_for_Quality=False,
):
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")
Expand Down Expand Up @@ -71,7 +77,7 @@ def generate(filename, version, with_QualityOfSurvey=False):
f.attrs["issueDate"] = "2023-12-31"
f.attrs["geographicIdentifier"] = "Somewhere"
f.attrs["verticalDatum"] = np.int16(12)
if version == "INT.IHO.S-102.2.2":
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)
Expand All @@ -90,9 +96,12 @@ def generate(filename, version, with_QualityOfSurvey=False):
b"<nothing/>"
)

if with_QualityOfSurvey:
QualityOfSurvey = f.create_group("QualityOfSurvey")
QualityOfSurvey_01 = QualityOfSurvey.create_group("QualityOfSurvey.01")
if with_QualityOfSurvey or with_QualityOfCoverage:
quality_name = (
"QualityOfSurvey" if with_QualityOfSurvey else "QualityOfBathymetryCoverage"
)
qualityGroup = f.create_group(quality_name)
quality01Group = qualityGroup.create_group(quality_name + ".01")

for attr_name in (
"gridOriginLongitude",
Expand All @@ -102,16 +111,26 @@ def generate(filename, version, with_QualityOfSurvey=False):
"numPointsLongitudinal",
"numPointsLatitudinal",
):
QualityOfSurvey_01.attrs[attr_name] = BathymetryCoverage_01.attrs[attr_name]

Group_001 = QualityOfSurvey_01.create_group("Group_001")

values = Group_001.create_dataset("values", (2, 3), dtype=np.uint32)
data = np.array(
[0, 1, 2, 1000000, 3, 2],
dtype=np.uint32,
).reshape(values.shape)
values[...] = data
quality01Group.attrs[attr_name] = BathymetryCoverage_01.attrs[attr_name]

Group_001 = quality01Group.create_group("Group_001")

if use_compound_type_for_Quality:
# As found in product 102DE00CA22_UNC_MD.H5
data_struct_type = np.dtype([("iD", "u4")])
values = Group_001.create_dataset("values", (2, 3), dtype=data_struct_type)
data = np.array(
[(0), (1), (2), (1000000), (3), (2)],
dtype=data_struct_type,
).reshape(values.shape)
values[...] = data
else:
values = Group_001.create_dataset("values", (2, 3), dtype=np.uint32)
data = np.array(
[0, 1, 2, 1000000, 3, 2],
dtype=np.uint32,
).reshape(values.shape)
values[...] = data

featureAttributeTable_struct_type = np.dtype(
[
Expand All @@ -121,7 +140,7 @@ def generate(filename, version, with_QualityOfSurvey=False):
]
)

featureAttributeTable = QualityOfSurvey.create_dataset(
featureAttributeTable = qualityGroup.create_dataset(
"featureAttributeTable", (5,), dtype=featureAttributeTable_struct_type
)

Expand All @@ -137,7 +156,7 @@ def generate(filename, version, with_QualityOfSurvey=False):
)
featureAttributeTable[...] = data

GroupFQualityOfSurvey_struct_type = np.dtype(
GroupFQuality_struct_type = np.dtype(
[
("code", "S16"),
("name", "S16"),
Expand All @@ -149,12 +168,12 @@ def generate(filename, version, with_QualityOfSurvey=False):
("closure", "S16"),
]
)
GroupFQualityOfSurvey = group_f.create_dataset(
"QualityOfSurvey", (1,), dtype=GroupFQualityOfSurvey_struct_type
GroupFQuality = group_f.create_dataset(
quality_name, (1,), dtype=GroupFQuality_struct_type
)
GroupFQualityOfSurvey[...] = np.array(
GroupFQuality[...] = np.array(
[("id", "", "", "0", "H5T_INTEGER", "1", "", "geSemiInterval")],
dtype=GroupFQualityOfSurvey_struct_type,
dtype=GroupFQuality_struct_type,
)


Expand All @@ -166,3 +185,10 @@ def generate(filename, version, with_QualityOfSurvey=False):
"INT.IHO.S-102.2.2",
with_QualityOfSurvey=True,
)

generate(
"test_s102_v3.0_with_QualityOfBathymetryCoverage",
"INT.IHO.S-102.3.0.0",
with_QualityOfCoverage=True,
use_compound_type_for_Quality=True,
)
Binary file added autotest/gdrivers/data/s102/test_s102_v3.0_with_QualityOfBathymetryCoverage.h5
View file
Binary file not shown.
32 changes: 19 additions & 13 deletions autotest/gdrivers/s102.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,34 +211,40 @@ def test_s102_multidim():
###############################################################################


def test_s102_QualityOfSurvey():
@pytest.mark.parametrize(
"filename,quality_group_name",
[
("data/s102/test_s102_v2.2_with_QualityOfSurvey.h5", "QualityOfSurvey"),
(
"data/s102/test_s102_v3.0_with_QualityOfBathymetryCoverage.h5",
"QualityOfBathymetryCoverage",
),
],
)
def test_s102_QualityOfSurvey(filename, quality_group_name):

ds = gdal.Open("data/s102/test_s102_v2.2_with_QualityOfSurvey.h5")
ds = gdal.Open(filename)
assert ds.GetSubDatasets() == [
(
'S102:"data/s102/test_s102_v2.2_with_QualityOfSurvey.h5":BathymetryCoverage',
f'S102:"{filename}":BathymetryCoverage',
"Bathymetric gridded data",
),
(
'S102:"data/s102/test_s102_v2.2_with_QualityOfSurvey.h5":QualityOfSurvey',
"Georeferenced metadata QualityOfSurvey",
f'S102:"{filename}":{quality_group_name}',
f"Georeferenced metadata {quality_group_name}",
),
]

with pytest.raises(Exception, match="Unsupported subdataset component"):
gdal.Open('S102:"data/s102/test_s102_v2.2_with_QualityOfSurvey.h5":invalid')
gdal.Open(f'S102:"{filename}":invalid')

ds = gdal.Open(
'S102:"data/s102/test_s102_v2.2_with_QualityOfSurvey.h5":BathymetryCoverage'
)
ds = gdal.Open(f'S102:"{filename}":BathymetryCoverage')
assert len(ds.GetSubDatasets()) == 0
assert ds.RasterCount == 2
assert ds.RasterXSize == 3
assert ds.RasterYSize == 2

ds = gdal.Open(
'S102:"data/s102/test_s102_v2.2_with_QualityOfSurvey.h5":QualityOfSurvey'
)
ds = gdal.Open(f'S102:"{filename}":{quality_group_name}')
assert len(ds.GetSubDatasets()) == 0
assert ds.RasterCount == 1
assert ds.RasterXSize == 3
Expand Down Expand Up @@ -278,7 +284,7 @@ def test_s102_QualityOfSurvey():
assert rat.GetValueAsString(4, 2) == "e"

ds = gdal.OpenEx(
'S102:"data/s102/test_s102_v2.2_with_QualityOfSurvey.h5":QualityOfSurvey',
f'S102:"{filename}":{quality_group_name}',
open_options=["NORTH_UP=NO"],
)
assert ds.GetGeoTransform() == pytest.approx((1.8, 0.4, 0.0, 47.75, 0.0, 0.5))
Expand Down
15 changes: 11 additions & 4 deletions doc/source/drivers/raster/s102.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ S102 -- S-102 Bathymetric Surface Product
.. 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
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.
Expand All @@ -25,6 +25,9 @@ Georeferencing is reported.

Nodata, minimum and maximum values for each band are also reported.

Supported versions of the specification are S-102 v2.1, v2.2 and v3.0
(support for v3.0 spatial metadata added in GDAL 3.10)

Driver capabilities
-------------------

Expand Down Expand Up @@ -63,12 +66,16 @@ The following open options are supported:
Spatial metadata support
------------------------

Starting with GDAL 3.9, GDAL can handle QualityOfSurvey spatial metadata.
Starting with GDAL 3.9, GDAL can handle QualityOfSurvey
(or QualityOfBathymetryCoverage in S102 v3.0) spatial metadata.

When such spatial metadata is present, the subdataset list will include
a name of the form ``S102:"{filename}":QualityOfSurvey``
a name of the form ``S102:"{filename}":QualityOfSurvey`` (
or ``S102:"{filename}":QualityOfBathymetryCoverage`` in S102 v3.0)

The ``/QualityOfSurvey/featureAttributeTable`` dataset is exposed as a
The ``/QualityOfSurvey/featureAttributeTable``
(``/QualityOfBathymetryCoverage/featureAttributeTable`` in S102 v3.0)
dataset is exposed as a
GDAL Raster Attribute Table associated to the GDAL raster band. The pixel
values of the raster match the ``id`` column of the Raster Attribute Table.

Expand Down
23 changes: 17 additions & 6 deletions frmts/hdf5/hdf5multidim.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1070,17 +1070,28 @@ HDF5Array::HDF5Array(const std::string &osParentName, const std::string &osName,
}

// Special case for S102 QualityOfSurvey nodata value that is typically at 0
if (GetFullName() ==
"/QualityOfSurvey/QualityOfSurvey.01/Group_001/values" &&
m_dt.GetClass() == GEDTC_NUMERIC &&
m_dt.GetNumericDataType() == GDT_UInt32)
const bool bIsQualityOfSurvey =
(GetFullName() ==
"/QualityOfSurvey/QualityOfSurvey.01/Group_001/values");
const bool bIsQualityOfBathymetryCoverage =
(GetFullName() == "/QualityOfBathymetryCoverage/"
"QualityOfBathymetryCoverage.01/Group_001/values");
if ((bIsQualityOfSurvey || bIsQualityOfBathymetryCoverage) &&
((m_dt.GetClass() == GEDTC_NUMERIC &&
m_dt.GetNumericDataType() == GDT_UInt32) ||
(m_dt.GetClass() == GEDTC_COMPOUND &&
m_dt.GetComponents().size() == 1 &&
m_dt.GetComponents()[0]->GetType().GetClass() == GEDTC_NUMERIC &&
m_dt.GetComponents()[0]->GetType().GetNumericDataType() ==
GDT_UInt32)))
{
if (auto poRootGroup = HDF5Array::GetRootGroup())
{
if (const auto poGroupF = poRootGroup->OpenGroup("Group_F"))
{
const auto poGroupFArray =
poGroupF->OpenMDArray("QualityOfSurvey");
const auto poGroupFArray = poGroupF->OpenMDArray(
bIsQualityOfSurvey ? "QualityOfSurvey"
: "QualityOfBathymetryCoverage");
if (poGroupFArray &&
poGroupFArray->GetDataType().GetClass() == GEDTC_COMPOUND &&
poGroupFArray->GetDataType().GetComponents().size() == 8 &&
Expand Down
Loading

0 comments on commit 79577b1

Please sign in to comment.