#!/usr/bin/env pytest
# -*- coding: utf-8 -*-
###############################################################################
#
# Project:  GDAL/OGR Test Suite
# Purpose:  Test read/write functionality for S102 driver.
# Author:   Even Rouault <even dot rouault at spatialys.com>
#
###############################################################################
# Copyright (c) 2023-2025, Even Rouault <even dot rouault at spatialys.com>
#
# SPDX-License-Identifier: MIT
###############################################################################

import os
import shutil
import struct
import sys

import gdaltest
import pytest
from test_py_scripts import samples_path

from osgeo import gdal, osr

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() == "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_NAME": "meanLowerLowWater",
        "VERTICAL_DATUM_DEFINITION": "The average height of the lower low waters at a place over a 19-year period.",
        "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")

    with pytest.raises(Exception, match="Cannot find group /QualityOfSurvey"):
        gdal.Open(f'S102:"{filename}":QualityOfSurvey')


###############################################################################


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() == "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() == "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() == "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)


###############################################################################


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

    ds = gdal.Open(filename)
    assert ds.GetSubDatasets() == [
        (
            f'S102:"{filename}":BathymetryCoverage',
            "Bathymetric gridded data",
        ),
        (
            f'S102:"{filename}":{quality_group_name}',
            f"Georeferenced metadata {quality_group_name}",
        ),
    ]

    with pytest.raises(Exception, match="Unsupported subdataset component"):
        gdal.Open(f'S102:"{filename}":invalid')

    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(f'S102:"{filename}":{quality_group_name}')
    assert len(ds.GetSubDatasets()) == 0
    assert ds.RasterCount == 1
    assert ds.RasterXSize == 3
    assert ds.RasterYSize == 2
    assert ds.GetSpatialRef().GetAuthorityCode() == "4326"
    assert ds.GetGeoTransform() == pytest.approx((1.8, 0.4, 0.0, 48.75, 0.0, -0.5))
    band = ds.GetRasterBand(1)
    assert band.DataType == gdal.GDT_UInt32
    assert band.GetNoDataValue() == 0
    assert struct.unpack("I" * 6, band.ReadRaster()) == (1000000, 3, 2, 0, 1, 2)

    rat = band.GetDefaultRAT()
    assert rat is not None
    assert rat.GetRowCount() == 5
    assert rat.GetColumnCount() == 3

    assert rat.GetNameOfCol(0) == "id"
    assert rat.GetTypeOfCol(0) == gdal.GFT_Integer
    assert rat.GetUsageOfCol(0) == gdal.GFU_MinMax

    assert rat.GetNameOfCol(1) == "floatval"
    assert rat.GetTypeOfCol(1) == gdal.GFT_Real

    assert rat.GetNameOfCol(2) == "strval"
    assert rat.GetTypeOfCol(2) == gdal.GFT_String

    assert rat.GetValueAsInt(0, 0) == 0
    assert rat.GetValueAsDouble(0, 1) == 1.5
    assert rat.GetValueAsString(0, 2) == "a"

    assert rat.GetValueAsInt(1, 0) == 1
    assert rat.GetValueAsDouble(1, 1) == 2.5
    assert rat.GetValueAsString(1, 2) == "b"

    assert rat.GetValueAsInt(4, 0) == 1000000
    assert rat.GetValueAsDouble(4, 1) == 5.5
    assert rat.GetValueAsString(4, 2) == "e"

    ds = gdal.OpenEx(
        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))
    band = ds.GetRasterBand(1)
    assert struct.unpack("I" * 6, band.ReadRaster()) == (0, 1, 2, 1000000, 3, 2)


###############################################################################


def test_s102_QualityOfSurvey_multidim():

    ds = gdal.OpenEx(
        "data/s102/test_s102_v2.2_with_QualityOfSurvey_nodata_0.h5",
        gdal.OF_MULTIDIM_RASTER,
    )
    rg = ds.GetRootGroup()
    ar = rg.OpenMDArrayFromFullname(
        "/QualityOfSurvey/QualityOfSurvey.01/Group_001/values"
    )
    assert ar.GetSpatialRef().GetAuthorityCode() == "4326"
    assert ar.GetNoDataValue() == 0

    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


###############################################################################
# Test force opening


def test_s102_force_opening():

    drv = gdal.IdentifyDriverEx("data/hdf5/groups.h5", allowed_drivers=["S102"])
    assert drv.GetDescription() == "S102"


###############################################################################
# Test force opening, but provided file is still not recognized (for good reasons)


def test_s102_force_opening_no_match():

    drv = gdal.IdentifyDriverEx("data/byte.tif", allowed_drivers=["S102"])
    assert drv is None


###############################################################################


def test_s102_metadata_compute_stats_first(tmp_path):

    out_filename = str(tmp_path / "out.h5")
    shutil.copy("data/s102/test_s102_v2.1.h5", out_filename)
    with gdal.Open(out_filename) as ds:
        ds.GetRasterBand(1).ComputeStatistics(False)
    with gdal.Open(out_filename) as ds:
        assert ds.GetRasterBand(1).GetMetadataItem("STATISTICS_MINIMUM") is not None
        ds.SetMetadataItem("foo", "bar")
    with gdal.Open(out_filename) as ds:
        assert ds.GetRasterBand(1).GetMetadataItem("STATISTICS_MINIMUM") is not None
        assert ds.GetMetadataItem("foo") == "bar"


###############################################################################


def test_s102_metadata_compute_stats_after(tmp_path):

    out_filename = str(tmp_path / "out.h5")
    shutil.copy("data/s102/test_s102_v2.1.h5", out_filename)
    with gdal.Open(out_filename) as ds:
        ds.SetMetadataItem("foo", "bar")
    with gdal.Open(out_filename) as ds:
        assert ds.GetMetadataItem("foo") == "bar"
        ds.GetRasterBand(1).ComputeStatistics(False)
    with gdal.Open(out_filename) as ds:
        assert ds.GetRasterBand(1).GetMetadataItem("STATISTICS_MINIMUM") is not None
        assert ds.GetMetadataItem("foo") == "bar"


###############################################################################


def test_s102_no_uncertainty():

    ds = gdal.Open("data/s102/test_s102_v3.0_without_uncertainty_nodata_0.h5")
    assert ds.RasterCount == 1
    assert struct.unpack("f" * 6, ds.ReadRaster()) == (1e6, 4, 5, 0, 1, 2)
    assert ds.GetRasterBand(1).GetNoDataValue() == 0


###############################################################################


def test_s102_multiple_feature_instance_groups():

    ds = gdal.Open("data/s102/multiple_feature_instance_groups.h5")
    assert ds.GetSubDatasets() == [
        (
            'S102:"data/s102/multiple_feature_instance_groups.h5":BathymetryCoverage.01',
            "Bathymetric gridded data, instance BathymetryCoverage.01, vertical datum meanLowerLowWater (MLLW)",
        ),
        (
            'S102:"data/s102/multiple_feature_instance_groups.h5":BathymetryCoverage.02',
            "Bathymetric gridded data, instance BathymetryCoverage.02, vertical datum lowWater (LW)",
        ),
        (
            'S102:"data/s102/multiple_feature_instance_groups.h5":QualityOfBathymetryCoverage',
            "Georeferenced metadata QualityOfBathymetryCoverage",
        ),
    ]
    assert ds.RasterCount == 0

    ds = gdal.Open(
        'S102:"data/s102/multiple_feature_instance_groups.h5":BathymetryCoverage.01'
    )
    assert ds.GetSubDatasets() == []
    assert ds.RasterCount == 1
    assert ds.RasterYSize == 2
    assert ds.RasterXSize == 4
    assert struct.unpack("f" * 8, ds.ReadRaster()) == (4, 5, 6, 7, 0, 1, 2, 3)
    assert ds.GetMetadataItem("VERTICAL_DATUM_NAME") == "meanLowerLowWater"

    ds = gdal.Open(
        'S102:"data/s102/multiple_feature_instance_groups.h5":BathymetryCoverage.02'
    )
    assert ds.GetSubDatasets() == []
    assert ds.RasterCount == 1
    assert ds.RasterYSize == 2
    assert ds.RasterXSize == 4
    assert struct.unpack("f" * 8, ds.ReadRaster()) == (40, 50, 60, 70, 0, 10, 20, 30)
    assert ds.GetMetadataItem("VERTICAL_DATUM_NAME") == "lowWater"

    ds = gdal.Open(
        'S102:"data/s102/multiple_feature_instance_groups.h5":QualityOfBathymetryCoverage'
    )
    assert ds.GetSubDatasets() == []
    assert ds.RasterCount == 1
    assert ds.RasterYSize == 2
    assert ds.RasterXSize == 4
    assert struct.unpack("i" * 8, ds.ReadRaster()) == (5, 6, 7, 8, 1, 2, 3, 4)

    with pytest.raises(
        Exception,
        match="Cannot find BathymetryCoverage.03 group in BathymetryCoverage group",
    ):
        gdal.Open(
            'S102:"data/s102/multiple_feature_instance_groups.h5":BathymetryCoverage.03'
        )


###############################################################################


def validate(
    filename, expected_errors=None, expected_warnings=None, expected_check_count=None
):

    path = samples_path
    if path not in sys.path:
        sys.path.append(path)
    try:
        import validate_s102
    except ImportError:
        print("Cannot import validate_s102")
        return True

    errors, warnings, checks_done = validate_s102.check(filename)

    if expected_errors:
        assert len(errors) == len(expected_errors)
        for error, expected_error in zip(errors, expected_errors):
            assert error[0] == expected_error[0]
            assert expected_error[1] in error[1]
    else:
        if errors:
            print(errors)
        assert not errors

    if expected_warnings:
        assert warnings == expected_warnings
    else:
        if warnings:
            print(warnings)
        assert not warnings

    if expected_check_count:
        assert len(checks_done) == expected_check_count


###############################################################################


def test_s102_validator():

    # Fake product: many unconformities
    expected_errors = [
        (
            "Critical error",
            "No featureCode array in feature information group ('/Group_F/featureCode')",
        ),
        ("Error", "top level attribute 'issueDate' is not a valid date: 2023-12-31"),
        (
            "Critical error",
            "Required top level attribute 'westBoundLongitude' is missing",
        ),
        (
            "Critical error",
            "Required top level attribute 'eastBoundLongitude' is missing",
        ),
        (
            "Critical error",
            "Required top level attribute 'southBoundLatitude' is missing",
        ),
        (
            "Critical error",
            "Required top level attribute 'northBoundLatitude' is missing",
        ),
        (
            "Critical error",
            "top level attribute 'verticalCoordinateBase' is not an enumeration",
        ),
        (
            "Critical error",
            "top level attribute 'verticalDatumReference' is not an enumeration",
        ),
        ("Critical error", "top level attribute 'verticalDatum' is not a uint16"),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'dataCodingFormat' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'dimension' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'commonPointRule' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'horizontalPositionUncertainty' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'verticalUncertainty' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'numInstances' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'sequencingRule.type' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'sequencingRule.scanDirection' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'interpolationType' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage group attribute 'dataOffsetCode' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'dataCodingFormat' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'dimension' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'commonPointRule' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'horizontalPositionUncertainty' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'verticalUncertainty' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'numInstances' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'sequencingRule.type' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'sequencingRule.scanDirection' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'interpolationType' is missing",
        ),
        (
            "Critical error",
            "Required QualityOfBathymetryCoverage group attribute 'dataOffsetCode' is missing",
        ),
        ("Error", "/BathymetryCoverage/axisNames dataset does not exist"),
        (
            "Critical error",
            "Required BathymetryCoverage feature instance group /BathymetryCoverage/BathymetryCoverage.01 attribute 'numGRP' is missing",
        ),
        (
            "Critical error",
            "Required BathymetryCoverage feature instance group /BathymetryCoverage/BathymetryCoverage.01 attribute 'startSequence' is missing",
        ),
        (
            "Critical error",
            "BathymetryCoverage feature instance group /BathymetryCoverage/BathymetryCoverage.01: dataset 'domainExtent.polygon' missing",
        ),
        (
            "Critical error",
            "Required Group_001 attribute 'minimumUncertainty' is missing",
        ),
        (
            "Critical error",
            "Required Group_001 attribute 'maximumUncertainty' is missing",
        ),
        ("Critical error", "Required Group_001 attribute 'timePoint' is missing"),
        ("Error", "/QualityOfBathymetryCoverage/axisNames dataset does not exist"),
        (
            "Error",
            "/QualityOfBathymetryCoverage/featureAttributeTable['floatval'] is not an allowed member",
        ),
        (
            "Error",
            "/QualityOfBathymetryCoverage/featureAttributeTable['strval'] is not an allowed member",
        ),
        (
            "Critical error",
            "/QualityOfBathymetryCoverage/QualityOfBathymetryCoverage.01/Group_001/values type is not uint32",
        ),
    ]
    expected_warnings = [
        "File name should start with '102'",
        "File name 'test_s102_v3.0_without_uncertainty_nodata_0.h5' does not match expected pattern '^102[a-zA-Z0-9]{4}[a-zA-Z0-9_]{1,12}\\.(?:h5|H5)$'",
        "Extra element in top level group: 'geographicIdentifier'",
        "Extra element in top level group: 'producer'",
    ]
    validate(
        "data/s102/test_s102_v3.0_without_uncertainty_nodata_0.h5",
        expected_errors=expected_errors,
        expected_warnings=expected_warnings,
    )


###############################################################################


def test_s102_write_errors(tmp_vsimem):

    with pytest.raises(
        Exception, match="Source dataset dimension must be at least 1x1 pixel"
    ):
        gdal.Translate(
            tmp_vsimem / "102xxxxyyyy.h5",
            gdal.GetDriverByName("MEM").Create("", 0, 0),
            format="S102",
        )

    with pytest.raises(
        Exception, match="S102 driver requires a source dataset with a geotransform"
    ):
        gdal.Translate(
            tmp_vsimem / "102xxxxyyyy.h5",
            gdal.GetDriverByName("MEM").Create("", 1, 1),
            format="S102",
        )

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
    with pytest.raises(
        Exception, match="S102 driver requires a source dataset with a geotransform"
    ):
        gdal.Translate(tmp_vsimem / "102xxxxyyyy.h5", src_ds, format="S102")

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
    src_ds.SetGeoTransform([0, 1, 0.2, 0, 0, 1])
    with pytest.raises(
        Exception,
        match="S102 driver requires a source dataset with a non-rotated geotransform",
    ):
        gdal.Translate(tmp_vsimem / "102xxxxyyyy.h5", src_ds, format="S102")

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
    src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1])
    with pytest.raises(
        Exception, match="S102 driver requires a source dataset with a CRS"
    ):
        gdal.Translate(tmp_vsimem / "102xxxxyyyy.h5", src_ds, format="S102")

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
    src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1])
    srs = osr.SpatialReference()
    srs.ImportFromProj4("+proj=longlat")
    src_ds.SetSpatialRef(srs)
    with pytest.raises(
        Exception,
        match="S102 driver requires a source dataset whose CRS has an EPSG identifier",
    ):
        gdal.Translate(tmp_vsimem / "102xxxxyyyy.h5", src_ds, format="S102")

    with pytest.raises(
        Exception, match="VERTICAL_DATUM creation option must be specified"
    ):
        gdal.Translate(tmp_vsimem / "102xxxxyyyy.h5", "data/byte.tif", format="S102")

    with pytest.raises(
        Exception, match="Cannot create file /i/do_not/exist/102xxxxyyyy.h5"
    ):
        gdal.Translate(
            "/i/do_not/exist/102xxxxyyyy.h5",
            gdal.Translate("", "data/byte.tif", format="MEM", outputSRS="EPSG:32611"),
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )


###############################################################################


def test_s102_write_warnings(tmp_vsimem):

    with gdaltest.error_raised(gdal.CE_Warning, match="Only EPSG codes"):
        gdal.Translate(
            tmp_vsimem / "102xxxxyyyy.h5",
            "data/byte.tif",
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )

    with gdaltest.error_raised(
        gdal.CE_Warning, match="S102 dataset filenames should start with '102'"
    ):
        gdal.Translate(
            tmp_vsimem / "non_conformant_prefix.h5",
            gdal.Translate("", "data/byte.tif", format="MEM", outputSRS="EPSG:32611"),
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )

    with gdaltest.error_raised(
        gdal.CE_Warning, match="S102 dataset filenames should have a '.H5' extension"
    ):
        gdal.Translate(
            tmp_vsimem / "102xxxxyyyy.oops",
            gdal.Translate("", "data/byte.tif", format="MEM", outputSRS="EPSG:32611"),
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )


###############################################################################


def test_s102_write_basic(tmp_path):

    with gdaltest.error_raised(gdal.CE_None):
        out_ds = gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            gdal.Translate("", "data/byte.tif", format="MEM", outputSRS="EPSG:32611"),
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )
        assert out_ds.GetRasterBand(1).Checksum() == 4672
        assert out_ds.GetSpatialRef().GetAuthorityCode() == "32611"
        assert out_ds.GetGeoTransform() == (440720.0, 60.0, 0.0, 3751320.0, 0.0, -60.0)

        out_ds.Close()

    validate(
        tmp_path / "102xxxxyyyy.h5",
        expected_warnings=["Quality feature not used"],
        expected_check_count=46,
    )


###############################################################################


def test_s102_write_elevation_to_depth(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 3, 1, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).SetDescription("elevation")
    src_ds.GetRasterBand(1).SetNoDataValue(255)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 2, 3, struct.pack("f" * 6, -1, -2, -3, -4, -5, 255)
    )
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, -1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(
        gdal.CE_Warning,
        match="Automatically convert from elevation to depth by negating elevation values",
    ):
        out_ds = gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )
        assert struct.unpack("f" * 6, out_ds.GetRasterBand(1).ReadRaster()) == (
            1,
            2,
            3,
            4,
            5,
            1e6,
        )


###############################################################################


@pytest.mark.parametrize(
    "value,warning_msg,validate_warning",
    [
        (
            -14.1,
            "Range of depth in the dataset is [-14.100000, 0.000000] whereas the allowed range is [-14, 11050]",
            "Group_001: minimumDepth=-14.100000381469727 should be in [-14, 11050] range",
        ),
        (
            11050.1,
            "Range of depth in the dataset is [0.000000, 11050.099609] whereas the allowed range is [-14, 11050]",
            "Group_001: maximumDepth=11050.099609375 should be in [-14, 11050] range",
        ),
    ],
)
def test_s102_write_depth_not_in_range(tmp_path, value, warning_msg, validate_warning):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2, 1, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 2, 2, struct.pack("f" * 4, value, 0, 0, 0)
    )
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, -1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(
        gdal.CE_Warning,
        match=warning_msg,
    ):
        gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )

    validate(
        tmp_path / "102xxxxyyyy.h5",
        expected_warnings=["Quality feature not used", validate_warning],
    )


###############################################################################


def test_s102_write_uncertainty_not_in_range(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2, 2, gdal.GDT_Float32)
    src_ds.GetRasterBand(2).WriteRaster(0, 0, 2, 2, struct.pack("f" * 4, -0.1, 0, 0, 0))
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, -1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(
        gdal.CE_Warning,
        match="Negative uncertainty value found, which is not allowed",
    ):
        gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )

    validate(
        tmp_path / "102xxxxyyyy.h5",
        expected_warnings=[
            "Quality feature not used",
            "Group_001: minimumUncertainty=-0.10000000149011612 should be in [0, inf) range or equal to 1000000",
        ],
    )


###############################################################################


def test_s102_write_with_uncertainty(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 3, 2)
    src_ds.GetRasterBand(1).SetNoDataValue(255)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 2, 3, struct.pack("B" * 6, 1, 2, 3, 4, 5, 255)
    )
    src_ds.GetRasterBand(2).SetNoDataValue(254)
    src_ds.GetRasterBand(2).WriteRaster(
        0, 0, 2, 3, struct.pack("B" * 6, 11, 12, 13, 14, 15, 254)
    )
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(gdal.CE_None):
        out_ds = gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={"VERTICAL_DATUM": "MLLW"},
        )
        assert struct.unpack("f" * 6, out_ds.GetRasterBand(1).ReadRaster()) == (
            5.0,
            1000000.0,
            3.0,
            4.0,
            1.0,
            2.0,
        )
        assert struct.unpack("f" * 6, out_ds.GetRasterBand(2).ReadRaster()) == (
            15.0,
            1000000.0,
            13.0,
            14.0,
            11.0,
            12.0,
        )
        assert out_ds.GetRasterBand(1).GetMinimum() == 1
        assert out_ds.GetRasterBand(1).GetMaximum() == 5
        assert out_ds.GetRasterBand(2).GetMinimum() == 11
        assert out_ds.GetRasterBand(2).GetMaximum() == 15

        out_ds.Close()

    validate(
        tmp_path / "102xxxxyyyy.h5",
        expected_warnings=["Quality feature not used"],
        expected_check_count=46,
    )


###############################################################################


def test_s102_write_with_quality(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2)
    src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 2, struct.pack("B" * 4, 1, 2, 3, 4))
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    rat = gdal.RasterAttributeTable()

    rat.CreateColumn("id", gdal.GFT_Integer, gdal.GFU_MinMax)
    rat.CreateColumn("dataAssessment", gdal.GFT_Integer, gdal.GFU_Generic)
    rat.CreateColumn("surveyDateRange.dateStart", gdal.GFT_String, gdal.GFU_Generic)
    rat.CreateColumn("featureSizeVar", gdal.GFT_Real, gdal.GFU_Generic)
    rat.CreateColumn("bathyCoverage", gdal.GFT_Boolean, gdal.GFU_Generic)
    rat.CreateColumn("sourceSurveyID", gdal.GFT_String, gdal.GFU_Generic)
    rat.CreateColumn(
        "typeOfBathymetricEstimationUncertainty", gdal.GFT_Integer, gdal.GFU_Generic
    )

    rat.SetRowCount(2)

    rat.SetValueAsInt(0, 0, 1)
    rat.SetValueAsInt(0, 1, 1)
    rat.SetValueAsString(0, 2, "20251026")
    rat.SetValueAsDouble(0, 3, 1.5)
    rat.SetValueAsBoolean(0, 4, True)
    rat.SetValueAsString(0, 5, "X")
    rat.SetValueAsInt(0, 6, 2)

    rat.SetValueAsInt(1, 0, 2)
    rat.SetValueAsInt(1, 1, 3)
    rat.SetValueAsString(1, 2, "2025-10-27")

    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "quality.tif", 2, 2
    ) as quality_ds:
        quality_ds.GetRasterBand(1).WriteRaster(
            0, 0, 2, 2, struct.pack("B" * 4, 1, 1, 2, 2)
        )
        quality_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
        quality_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        quality_ds.GetRasterBand(1).SetDefaultRAT(rat)

    gdal.Translate(
        tmp_path / "102xxxxyyyy.h5",
        src_ds,
        format="S102",
        creationOptions={
            "VERTICAL_DATUM": "MLLW",
            "QUALITY_DATASET": tmp_path / "quality.tif",
        },
    )

    validate(tmp_path / "102xxxxyyyy.h5")

    ds = gdal.Open(f'S102:"{tmp_path}/102xxxxyyyy.h5":QualityOfBathymetryCoverage')
    info = gdal.Info(ds, format="json")
    rat = info["bands"][0]["rat"]
    assert rat == {
        "tableType": "thematic",
        "fieldDefn": [
            {"index": 0, "name": "id", "type": 0, "usage": 5},
            {"index": 1, "name": "dataAssessment", "type": 0, "usage": 0},
            {"index": 2, "name": "surveyDateRange.dateStart", "type": 4, "usage": 0},
            {"index": 3, "name": "featureSizeVar", "type": 1, "usage": 0},
            {"index": 4, "name": "bathyCoverage", "type": 3, "usage": 0},
            {"index": 5, "name": "sourceSurveyID", "type": 2, "usage": 0},
            {
                "index": 6,
                "name": "typeOfBathymetricEstimationUncertainty",
                "type": 0,
                "usage": 0,
            },
        ],
        "row": [
            {
                "index": 0,
                "f": [1, 1, "2025-10-26T00:00:00.000+00:00", 1.5, True, "X", 2],
            },
            {
                "index": 1,
                "f": [2, 3, "2025-10-27T00:00:00.000+00:00", 0.0, False, "", 0],
            },
        ],
    }


###############################################################################


def test_s102_write_with_quality_custom_fields_rat(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2)
    src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 2, struct.pack("B" * 4, 1, 2, 3, 4))
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    rat = gdal.RasterAttributeTable()

    rat.CreateColumn("id", gdal.GFT_Integer, gdal.GFU_MinMax)
    rat.CreateColumn("my_int_field", gdal.GFT_Integer, gdal.GFU_Generic)
    rat.CreateColumn("my_real_field", gdal.GFT_Real, gdal.GFU_Generic)
    rat.CreateColumn("my_string_field", gdal.GFT_String, gdal.GFU_Generic)
    rat.CreateColumn("my_bool_field", gdal.GFT_Boolean, gdal.GFU_Generic)

    rat.SetValueAsInt(0, 0, 1)
    rat.SetValueAsInt(0, 1, 123)
    rat.SetValueAsDouble(0, 2, 1.25)
    rat.SetValueAsString(0, 3, "foo")
    rat.SetValueAsBoolean(0, 4, True)

    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "quality.tif", 2, 2
    ) as quality_ds:
        quality_ds.GetRasterBand(1).WriteRaster(
            0, 0, 2, 2, struct.pack("B" * 4, 1, 1, 1, 1)
        )
        quality_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
        quality_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        quality_ds.GetRasterBand(1).SetDefaultRAT(rat)

    with gdaltest.error_raised(gdal.CE_Warning):
        gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={
                "VERTICAL_DATUM": "MLLW",
                "QUALITY_DATASET": tmp_path / "quality.tif",
            },
        )

    expected_errors = [
        (
            "Error",
            "/QualityOfBathymetryCoverage/featureAttributeTable['my_int_field'] is not an allowed member",
        ),
        (
            "Error",
            "/QualityOfBathymetryCoverage/featureAttributeTable['my_real_field'] is not an allowed member",
        ),
        (
            "Error",
            "/QualityOfBathymetryCoverage/featureAttributeTable['my_string_field'] is not an allowed member",
        ),
        (
            "Error",
            "/QualityOfBathymetryCoverage/featureAttributeTable['my_bool_field'] is not an allowed member",
        ),
    ]
    validate(tmp_path / "102xxxxyyyy.h5", expected_errors=expected_errors)

    ds = gdal.Open(f'S102:"{tmp_path}/102xxxxyyyy.h5":QualityOfBathymetryCoverage')
    info = gdal.Info(ds, format="json")
    rat = info["bands"][0]["rat"]
    assert rat == {
        "tableType": "thematic",
        "fieldDefn": [
            {"index": 0, "name": "id", "type": 0, "usage": 5},
            {"index": 1, "name": "my_int_field", "type": 0, "usage": 0},
            {"index": 2, "name": "my_real_field", "type": 1, "usage": 0},
            {"index": 3, "name": "my_string_field", "type": 2, "usage": 0},
            {"index": 4, "name": "my_bool_field", "type": 0, "usage": 0},
        ],
        "row": [{"index": 0, "f": [1, 123, 1.25, "foo", 1]}],
    }


###############################################################################


def test_s102_write_with_quality_missing_id_in_rat(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2)
    src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 2, struct.pack("B" * 4, 1, 2, 3, 4))
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    rat = gdal.RasterAttributeTable()

    rat.CreateColumn("id", gdal.GFT_Integer, gdal.GFU_MinMax)

    rat.SetValueAsInt(0, 0, 1)

    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "quality.tif", 2, 2
    ) as quality_ds:
        quality_ds.GetRasterBand(1).WriteRaster(
            0, 0, 2, 2, struct.pack("B" * 4, 1, 1, 1, 2)
        )
        quality_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
        quality_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        quality_ds.GetRasterBand(1).SetDefaultRAT(rat)

    with gdaltest.error_raised(
        gdal.CE_Warning,
        match="Quality grid contains nodes with id 2, but there is no such entry in the feature attribute table",
    ):
        gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={
                "VERTICAL_DATUM": "MLLW",
                "QUALITY_DATASET": tmp_path / "quality.tif",
            },
        )

    expected_errors = [
        (
            "Error",
            "/BathymetryCoverage/BathymetryCoverage.01/Group_001/values contain value 2, which is not a valid 'id' of the featureAttributeTable",
        )
    ]
    validate(tmp_path / "102xxxxyyyy.h5", expected_errors=expected_errors)


###############################################################################


@pytest.mark.parametrize(
    "error_source,expected_error_msg",
    [
        ("non_existing_dataset", "i_do_not_exist"),
        ("no_rat", "does not have a raster attribute table"),
        ("no_gt", "does not have the same geotransform"),
        ("no_srs", "does not have the same CRS"),
        ("other_width", "does not have the same dimensions as"),
        ("other_height", "does not have the same dimensions as"),
        ("wrong_band_count", "does not have a single band"),
        ("wrong_data_type", "is not of an integer data type"),
        ("no_id_field", "Input raster attribute table lacks an integer 'id' field"),
        ("id_0", "id=0 is not allowed in input raster attribute table"),
        ("id_negative", "Negative id is not allowed in input raster attribute table"),
        ("id_same", "Several rows of input raster attribute table have id=1"),
    ],
)
def test_s102_write_with_quality_errors(tmp_path, error_source, expected_error_msg):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2)
    src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 2, struct.pack("B" * 4, 1, 2, 3, 4))
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    rat = gdal.RasterAttributeTable()
    rat.CreateColumn(
        "id" if error_source != "no_id_field" else "foo",
        gdal.GFT_Integer,
        gdal.GFU_MinMax,
    )
    if error_source == "id_0":
        rat.SetValueAsInt(0, 0, 0)
    elif error_source == "id_negative":
        rat.SetValueAsInt(0, 0, -1)
    elif error_source == "id_same":
        rat.SetValueAsInt(0, 0, 1)
        rat.SetValueAsInt(1, 0, 1)
    else:
        rat.SetValueAsInt(0, 0, 1)

    quality_width = 2
    if error_source == "other_width":
        quality_width = 1
    quality_height = 2
    if error_source == "other_height":
        quality_height = 1
    quality_band_count = 1
    if error_source == "wrong_band_count":
        quality_band_count = 2
    quality_data_type = gdal.GDT_Byte
    if error_source == "wrong_data_type":
        quality_data_type = gdal.GDT_Float32
    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "quality.tif",
        quality_width,
        quality_height,
        quality_band_count,
        quality_data_type,
    ) as quality_ds:
        quality_ds.GetRasterBand(1).Fill(1)
        if error_source != "no_gt":
            quality_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
        if error_source != "no_srs":
            quality_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        if error_source != "no_rat":
            quality_ds.GetRasterBand(1).SetDefaultRAT(rat)

    with pytest.raises(Exception, match=expected_error_msg), gdal.quiet_errors():
        gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={
                "VERTICAL_DATUM": "MLLW",
                "QUALITY_DATASET": (
                    tmp_path / "i_do_not_exist"
                    if error_source == "non_existing_dataset"
                    else tmp_path / "quality.tif"
                ),
            },
        )


###############################################################################


def test_s102_write_large_file(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 1200, 1200, 2, gdal.GDT_Float32)
    src_ds.GetRasterBand(2).WriteRaster(0, 0, 1, 1, struct.pack("f", 1))
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    filename = str(tmp_path / "102xxxxyyyy.h5")
    with gdaltest.error_raised(
        gdal.CE_Warning,
        match="file size exceeds 10 MB",
    ):
        gdal.Translate(
            filename,
            src_ds,
            format="S102",
            creationOptions={
                "VERTICAL_DATUM": "MLLW",
                "COMPRESS": "NONE",
            },
        )

    expected_warnings = [
        f"File size of {filename} = 11548224, which exceeds 10 MB",
        "Quality feature not used",
    ]
    validate(filename, expected_warnings=expected_warnings)


###############################################################################


def test_s102_write_several_feature_instances(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2, 1, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 2, 2, struct.pack("f" * 4, 1, 2, float("nan"), float("nan"))
    )
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    gdal.Translate(
        tmp_path / "102xxxxyyyy.h5",
        src_ds,
        format="S102",
        creationOptions={"VERTICAL_DATUM": "MLLW"},
    )

    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 2, 2, struct.pack("f" * 4, float("nan"), float("nan"), 3, 4)
    )
    src_ds.SetGeoTransform([500000 + 200, 1, 0, 4500000 + 200, 0, 1])

    gdal.Translate(
        tmp_path / "102xxxxyyyy.h5",
        src_ds,
        format="S102",
        creationOptions={
            "VERTICAL_DATUM": "MLW",
            "APPEND_SUBDATASET": "YES",
        },
    )
    validate(
        tmp_path / "102xxxxyyyy.h5", expected_warnings=["Quality feature not used"]
    )

    ds = gdal.Open(f'S102:"{tmp_path}/102xxxxyyyy.h5":BathymetryCoverage.01')
    assert struct.unpack("f" * 4, ds.GetRasterBand(1).ReadRaster()) == (1e6, 1e6, 1, 2)
    assert ds.GetMetadataItem("VERTICAL_DATUM_ABBREV") == "MLLW"
    assert ds.GetGeoTransform() == (500000, 1, 0, 4500002, 0, -1)

    ds = gdal.Open(f'S102:"{tmp_path}/102xxxxyyyy.h5":BathymetryCoverage.02')
    assert struct.unpack("f" * 4, ds.GetRasterBand(1).ReadRaster()) == (3, 4, 1e6, 1e6)
    assert ds.GetMetadataItem("VERTICAL_DATUM_ABBREV") == "MLW"
    assert ds.GetGeoTransform() == (500000 + 200, 1, 0, 4500002 + 200, 0, -1)


###############################################################################


def test_s102_write_several_feature_instances_error_not_existing_s102(tmp_path):

    open(tmp_path / "102xxxxyyyy.h5", "wb").close()

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2, 1, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 2, 2, struct.pack("f" * 4, 1, 2, float("nan"), float("nan"))
    )
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with pytest.raises(Exception, match="is not a valid existing S102 dataset"):
        gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={
                "VERTICAL_DATUM": "MLLW",
                "APPEND_SUBDATASET": "YES",
            },
        )


###############################################################################


def test_s102_write_several_feature_instances_error_not_same_srs(tmp_path):

    open(tmp_path / "102xxxxyyyy.h5", "wb").close()

    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2, 1, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 2, 2, struct.pack("f" * 4, 1, 2, float("nan"), float("nan"))
    )
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
    gdal.Translate(
        tmp_path / "102xxxxyyyy.h5",
        src_ds,
        format="S102",
        creationOptions={
            "VERTICAL_DATUM": "MLLW",
        },
    )

    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32632))
    with pytest.raises(Exception, match="CRS.*is not the same as the one"):
        gdal.Translate(
            tmp_path / "102xxxxyyyy.h5",
            src_ds,
            format="S102",
            creationOptions={
                "VERTICAL_DATUM": "MLLW",
                "APPEND_SUBDATASET": "YES",
            },
        )
