#!/usr/bin/env pytest
# -*- coding: utf-8 -*-
###############################################################################
#
# Project:  GDAL/OGR Test Suite
# Purpose:  Test read functionality for S104 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 struct
import sys

import gdaltest
import pytest
from test_py_scripts import samples_path

from osgeo import gdal, osr

pytestmark = pytest.mark.require_driver("S104")


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


def test_s104_basic():
    filename = "data/s104/test_s104_v1.1.h5"
    ds = gdal.Open(filename)
    assert ds.RasterCount == 0
    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",
        "verticalCS": "6498",
        "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
        "dateTimeOfFirstRecord": "20190606T120000Z",
        "dateTimeOfLastRecord": "20190606T120000Z",
        "geographicIdentifier": "Somewhere",
        "issueDate": "2023-12-31",
        "maxDatasetHeight": "2",
        "minDatasetHeight": "1",
        "numberOfTimes": "1",
        "producer": "Generated by autotest/gdrivers/data/s104/generate_test.py (not strictly fully S104 compliant)",
        "timeRecordInterval": "3600",
        "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.",
    }

    assert ds.GetSubDatasets() == [
        (
            f'S104:"{filename}":Group_001',
            "Values at timestamp 20190606T120000Z",
        )
    ]

    with pytest.raises(
        Exception, match="Cannot find /WaterLevel/WaterLevel.01/invalid group"
    ):
        gdal.Open(f'S104:"{filename}":invalid')

    ds = gdal.Open(f'S104:"{filename}":Group_001')

    band = ds.GetRasterBand(1)
    assert band.GetDescription() == "waterLevelHeight"
    assert band.GetNoDataValue() == -123
    assert band.GetUnitType() == "metre"
    assert struct.unpack("f" * 6, band.ReadRaster()) == (
        3.0,
        4.0,
        5.0,
        -123.0,
        1.0,
        2.0,
    )

    band = ds.GetRasterBand(2)
    assert band.GetDescription() == "waterLevelTrend"
    assert band.GetNoDataValue() == 0
    assert struct.unpack("B" * 6, band.ReadRaster()) == (3, 2, 1, 0, 1, 2)

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

    assert rat.GetNameOfCol(0) == "code"
    assert rat.GetTypeOfCol(0) == gdal.GFT_Integer

    assert rat.GetNameOfCol(1) == "label"
    assert rat.GetTypeOfCol(1) == gdal.GFT_String

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

    assert rat.GetValueAsInt(1, 0) == 1
    assert rat.GetValueAsString(1, 1) == "Decreasing"

    assert rat.GetValueAsInt(2, 0) == 2
    assert rat.GetValueAsString(2, 1) == "Increasing"

    assert rat.GetValueAsInt(3, 0) == 3
    assert rat.GetValueAsString(3, 1) == "Steady"

    assert "MD_" in ds.GetFileList()[1]

    del ds
    assert not os.path.exists(f"{filename}.aux.xml")


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


def test_s104_north_up_no():
    filename = "data/s104/test_s104_v1.1.h5"
    ds = gdal.OpenEx(f'S104:"{filename}":Group_001', 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))

    band = ds.GetRasterBand(1)
    assert band.GetDescription() == "waterLevelHeight"
    assert band.GetNoDataValue() == -123
    assert band.GetUnitType() == "metre"
    assert struct.unpack("f" * 6, band.ReadRaster()) == (
        -123.0,
        1.0,
        2.0,
        3.0,
        4.0,
        5.0,
    )

    band = ds.GetRasterBand(2)
    assert band.GetDescription() == "waterLevelTrend"
    assert band.GetNoDataValue() == 0
    assert struct.unpack("B" * 6, band.ReadRaster()) == (
        0,
        1,
        2,
        3,
        2,
        1,
    )

    del ds
    assert not os.path.exists("data/s104/test_s104_v2.1.h5.aux.xml")


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


def test_s104_multidim():

    filename = "data/s104/test_s104_v1.1.h5"
    ds = gdal.OpenEx(filename, gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    ar = rg.OpenMDArrayFromFullname("/WaterLevel/WaterLevel.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


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


@pytest.mark.parametrize(
    "filename,wkt",
    [
        (
            "test_s104_custom_geog_crs.h5",
            'GEOGCRS["my crs",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]',
        ),
        (
            "test_s104_custom_geog_crs_custom_datum.h5",
            'GEOGCRS["my crs",DATUM["my datum",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]',
        ),
        (
            "test_s104_custom_proj_mercator.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Mercator (variant B)",METHOD["Mercator (variant B)",ID["EPSG",9805]],PARAMETER["Latitude of 1st standard parallel",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_transverse_mercator.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Transverse Mercator",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_oblique_mercator.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Hotine Oblique Mercator (variant B)",METHOD["Hotine Oblique Mercator (variant B)",ID["EPSG",9815]],PARAMETER["Latitude of projection centre",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8811]],PARAMETER["Longitude of projection centre",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8812]],PARAMETER["Azimuth at projection centre",3.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8813]],PARAMETER["Angle from Rectified to Skew Grid",4.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8814]],PARAMETER["Scale factor at projection centre",5.5,SCALEUNIT["unity",1],ID["EPSG",8815]],PARAMETER["Easting at projection centre",100000,LENGTHUNIT["metre",1],ID["EPSG",8816]],PARAMETER["Northing at projection centre",200000,LENGTHUNIT["metre",1],ID["EPSG",8817]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_hotine_oblique_mercator.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Hotine Oblique Mercator (variant A)",METHOD["Hotine Oblique Mercator (variant A)",ID["EPSG",9812]],PARAMETER["Latitude of projection centre",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8811]],PARAMETER["Longitude of projection centre",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8812]],PARAMETER["Azimuth at projection centre",3.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8813]],PARAMETER["Angle from Rectified to Skew Grid",4.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8814]],PARAMETER["Scale factor at projection centre",5.5,SCALEUNIT["unity",1],ID["EPSG",8815]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_lcc_1sp.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Lambert Conic Conformal (1SP)",METHOD["Lambert Conic Conformal (1SP)",ID["EPSG",9801]],PARAMETER["Latitude of natural origin",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_lcc_2sp.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Lambert Conic Conformal (2SP)",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",3.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",4.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",100000,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",200000,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_oblique_stereographic.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Oblique Stereographic",METHOD["Oblique Stereographic",ID["EPSG",9809]],PARAMETER["Latitude of natural origin",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_polar_stereographic.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Polar Stereographic (variant A)",METHOD["Polar Stereographic (variant A)",ID["EPSG",9810]],PARAMETER["Latitude of natural origin",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_krovak_oblique_conic_conformal.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Krovak",METHOD["Krovak",ID["EPSG",9819]],PARAMETER["Latitude of projection centre",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8811]],PARAMETER["Longitude of origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["Co-latitude of cone axis",3.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",1036]],PARAMETER["Latitude of pseudo standard parallel",4.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8818]],PARAMETER["Scale factor on pseudo standard parallel",0.9,SCALEUNIT["unity",1],ID["EPSG",8819]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_american_polyconic.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["American Polyconic",METHOD["American Polyconic",ID["EPSG",9818]],PARAMETER["Latitude of natural origin",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]    ',
        ),
        (
            "test_s104_custom_proj_albers_equal_area.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Albers Equal Area",METHOD["Albers Equal Area",ID["EPSG",9822]],PARAMETER["Latitude of false origin",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",3.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",4.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",100000,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",200000,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
        (
            "test_s104_custom_proj_lambert_azimuthal_equal_area.h5",
            'PROJCRS["my crs",BASEGEOGCRS["World Geodetic System 1984",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["Lambert Azimuthal Equal Area",METHOD["Lambert Azimuthal Equal Area",ID["EPSG",9820]],PARAMETER["Latitude of natural origin",1.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",2.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",100000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
        ),
    ],
)
def test_s104_custom_crs(filename, wkt):

    ds = gdal.Open("data/s104/" + filename)
    srs_ref = osr.SpatialReference(wkt)
    srs_ref.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    assert ds.GetSpatialRef().IsSame(srs_ref), ds.GetSpatialRef().ExportToWkt(
        ["FORMAT=WKT2_2019", "MULTILINE=NO"]
    )


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


def test_s104_multiple_feature_instance_groups():

    ds = gdal.Open("data/s104/multiple_feature_instance_groups.h5")
    assert ds.GetSubDatasets() == [
        (
            'S104:"data/s104/multiple_feature_instance_groups.h5":WaterLevel.01:Group_001',
            "Values for feature instance WaterLevel.01, vertical datum meanLowerLowWater (MLLW) at timestamp 20190606T120000Z",
        ),
        (
            'S104:"data/s104/multiple_feature_instance_groups.h5":WaterLevel.02:Group_001',
            "Values for feature instance WaterLevel.02, vertical datum lowWater (LW) at timestamp 20190606T120000Z",
        ),
    ]
    assert ds.RasterCount == 0

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

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


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


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_s104
    except ImportError:
        print("Cannot import validate_s104")
        return True

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

    if expected_errors:
        assert errors == expected_errors
    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_s104_validator():

    # Fake product: many unconformities
    expected_errors = [
        (
            "Critical error",
            "/Group_F/WaterLevel: row 0, 2, got value 'metres', whereas 'metre' is expected",
        ),
        ("Error", "top level attribute 'issueDate' is not a valid date: 2025-10-07"),
        ("Error", "top level attribute 'horizontalCRS' is not a int32"),
        ("Error", "top level attribute 'westBoundLongitude' is not a float32"),
        ("Error", "top level attribute 'eastBoundLongitude' is not a float32"),
        ("Error", "top level attribute 'southBoundLatitude' is not a float32"),
        ("Error", "top level attribute 'northBoundLatitude' is not a float32"),
        ("Error", "top level attribute 'issueTime' is not a valid time: 12:34:56"),
        ("Error", "top level attribute 'verticalCS' is not a int32"),
        ("Error", "top level attribute 'verticalCoordinateBase' is not an enumeration"),
        ("Error", "top level attribute 'verticalDatumReference' is not an enumeration"),
        ("Error", "top level attribute 'verticalDatum' is not a int32"),
        (
            "Error",
            "WaterLevel group attribute 'dataCodingFormat' is not an enumeration",
        ),
        ("Error", "WaterLevel group attribute 'dimension' is not a uint8"),
        ("Error", "WaterLevel group attribute 'commonPointRule' is not an enumeration"),
        (
            "Error",
            "WaterLevel group attribute 'horizontalPositionUncertainty' is not a float32",
        ),
        ("Error", "WaterLevel group attribute 'verticalUncertainty' is not a float32"),
        ("Error", "WaterLevel group attribute 'numInstances' is not a uint32"),
        (
            "Error",
            "WaterLevel group attribute 'sequencingRule.type' is not an enumeration",
        ),
        (
            "Error",
            "WaterLevel group attribute 'interpolationType' is not an enumeration",
        ),
        ("Error", "WaterLevel group attribute 'dataOffsetCode' is not an enumeration"),
        (
            "Error",
            '/WaterLevel/axisNames must conform to CRS. Expected ["Easting", "Northing"] or ["Latitude", "Longitude"]. Got [\'longitude\', \'latitude\']',
        ),
        (
            "Error",
            "Required WaterLevel feature instance group /WaterLevel/WaterLevel.01 attribute 'dataDynamicity' is missing",
        ),
        (
            "Critical error",
            "/WaterLevel/WaterLevel.01/Group_001/values member b'waterLevelTrend' is not an enumeration",
        ),
        (
            "Error",
            "Required WaterLevel feature instance group /WaterLevel/WaterLevel.02 attribute 'dataDynamicity' is missing",
        ),
        (
            "Error",
            "WaterLevel feature instance group /WaterLevel/WaterLevel.02 attribute 'verticalDatumReference' is not an enumeration",
        ),
        (
            "Error",
            "WaterLevel feature instance group /WaterLevel/WaterLevel.02 attribute 'verticalDatum' is not a int32",
        ),
        (
            "Critical error",
            "/WaterLevel/WaterLevel.02/Group_001/values member b'waterLevelTrend' is not an enumeration",
        ),
    ]
    expected_warnings = [
        "File name should start with '104'",
        "File name 'multiple_feature_instance_groups.h5' does not match expected pattern '^104[a-zA-Z0-9]{4}[a-zA-Z0-9\\-_]{1,54}\\.(?:h5|H5)$'",
    ]
    validate(
        "data/s104/multiple_feature_instance_groups.h5",
        expected_errors=expected_errors,
        expected_warnings=expected_warnings,
    )


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


def test_s104_write_errors(tmp_vsimem):

    with pytest.raises(
        Exception, match="Source dataset x must have two or three bands"
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            gdal.GetDriverByName("MEM").Create("x", 2, 2),
            format="S104",
        )

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

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

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

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

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

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 2)
    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="VERTICAL_DATUM creation option must be specified"
    ):
        gdal.Translate(tmp_vsimem / "104xxxxyyyy.h5", src_ds, format="S104")

    src_ds = gdal.GetDriverByName("MEM").Create("x", 1, 1, 2)
    src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with pytest.raises(Exception, match="VERTICAL_DATUM value is invalid"):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=["VERTICAL_DATUM=invalid"],
        )

    with pytest.raises(
        Exception,
        match="TIME_POINT creation option value must be set, or source dataset must have a timePoint metadata item",
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=["VERTICAL_DATUM=MLLW"],
        )

    with pytest.raises(
        Exception,
        match="TIME_POINT creation option value must be set to a YYYYMMDDTHHMMSSZ datetime value.",
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=["VERTICAL_DATUM=MLLW", "TIME_POINT=invalid"],
        )

    with pytest.raises(
        Exception, match="VERTICAL_CS creation option must be specified"
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=["VERTICAL_DATUM=MLLW", "TIME_POINT=20251104T225013Z"],
        )

    with pytest.raises(
        Exception,
        match=r"VERTICAL_CS creation option must be set either to 6498 \(depth/down, metre\), or 6499 \(height/up, metre\)",
    ):
        with gdal.quiet_errors():
            gdal.Translate(
                tmp_vsimem / "104xxxxyyyy.h5",
                src_ds,
                format="S104",
                creationOptions=[
                    "VERTICAL_DATUM=MLLW",
                    "TIME_POINT=20251104T225013Z",
                    "VERTICAL_CS=invalid",
                ],
            )

    with pytest.raises(
        Exception, match="WATER_LEVEL_TREND_THRESHOLD creation option must be specified"
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
            ],
        )

    with pytest.raises(
        Exception,
        match="WATER_LEVEL_TREND_THRESHOLD creation option value must be a numeric value",
    ):
        with gdal.quiet_errors():
            gdal.Translate(
                tmp_vsimem / "104xxxxyyyy.h5",
                src_ds,
                format="S104",
                creationOptions=[
                    "VERTICAL_DATUM=MLLW",
                    "TIME_POINT=20251104T225013Z",
                    "VERTICAL_CS=DEPTH",
                    "WATER_LEVEL_TREND_THRESHOLD=invalid",
                ],
            )

    with pytest.raises(
        Exception, match="DATA_DYNAMICITY creation option must be specified"
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
            ],
        )

    with pytest.raises(
        Exception,
        match="DATA_DYNAMICITY creation option must be set to observation/1, astronomicalPrediction/2, analysisOrHybrid/3 or hydrodynamicForecast/5",
    ):
        with gdal.quiet_errors():
            gdal.Translate(
                tmp_vsimem / "104xxxxyyyy.h5",
                src_ds,
                format="S104",
                creationOptions=[
                    "VERTICAL_DATUM=MLLW",
                    "TIME_POINT=20251104T225013Z",
                    "VERTICAL_CS=DEPTH",
                    "WATER_LEVEL_TREND_THRESHOLD=0.2",
                    "DATA_DYNAMICITY=invalid",
                ],
            )

    with pytest.raises(
        Exception, match="Cannot create file /i/do_not/exist/104xxxxyyyy.h5"
    ):
        gdal.Translate(
            tmp_vsimem / "/i/do_not/exist/104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
            ],
        )

    with pytest.raises(Exception, match="non_existing"):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
                "DATASETS=non_existing",
            ],
        )

    with pytest.raises(
        Exception, match="Dataset data/byte.tif does not have the same dimensions as x"
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
                "DATASETS=data/byte.tif",
            ],
        )

    with pytest.raises(
        Exception,
        match="DATASETS_TIME_POINT does not have the same number of values as DATASETS",
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
                "DATASETS=data/byte.tif",
                "DATASETS_TIME_POINT=foo,bar",
            ],
        )

    gdal.GetDriverByName("GTiff").CreateCopy(tmp_vsimem / "aux.tif", src_ds)

    with pytest.raises(
        Exception,
        match="Dataset /vsimem/test_s104_write_errors/aux.tif does not have a timePoint metadata item, and the DATASETS_TIME_POINT creation option is not set",
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
                f"DATASETS={tmp_vsimem}/aux.tif",
            ],
        )

    with pytest.raises(
        Exception,
        match="timePoint value for dataset /vsimem/test_s104_write_errors/aux.tif is invalid, but does not conform to a YYYYMMDDTHHMMSSZ datetime value",
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
                f"DATASETS={tmp_vsimem}/aux.tif",
                "DATASETS_TIME_POINT=invalid",
            ],
        )

    with pytest.raises(
        Exception,
        match=r"Several datasets are at timePoint 20251104T225013Z \(/vsimem/test_s104_write_errors/aux.tif vs x\)",
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
                f"DATASETS={tmp_vsimem}/aux.tif",
                "DATASETS_TIME_POINT=20251104T225013Z",
            ],
        )

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 2)
    src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1])
    srs = osr.SpatialReference()
    srs.ImportFromProj4("+proj=eqc")
    src_ds.SetSpatialRef(srs)
    with pytest.raises(
        Exception, match="Projection method Equirectangular is not supported by S100"
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
            ],
        )

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 2)
    src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1])
    srs = osr.SpatialReference()
    srs.ImportFromProj4("+proj=merc +a=1")
    src_ds.SetSpatialRef(srs)
    with pytest.raises(Exception, match="Unknown code for ellipsoid of CRS"):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
            ],
        )

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 2)
    src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1])
    srs = osr.SpatialReference()
    srs.ImportFromProj4("+proj=merc +ellps=GRS80 +pm=5")
    src_ds.SetSpatialRef(srs)
    with pytest.raises(Exception, match="Unknown code for prime meridian of CRS"):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
            ],
        )

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 2)
    src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1])
    srs = osr.SpatialReference()
    srs.SetFromUserInput("EPSG:4326+3855")
    src_ds.SetSpatialRef(srs)
    with pytest.raises(
        Exception, match="The CRS must be a geographic 2D or projected 2D CRS"
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
            ],
        )


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


def test_s104_write_warnings(tmp_vsimem):

    src_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 2)
    src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(
        gdal.CE_Warning, match="S104 dataset filenames should start with '104'"
    ):
        gdal.Translate(
            tmp_vsimem / "non_conformant_prefix.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
            ],
        )

    with gdaltest.error_raised(
        gdal.CE_Warning, match="S104 dataset filenames should have a '.H5' extension"
    ):
        gdal.Translate(
            tmp_vsimem / "104xxxxyyyy.oops",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
            ],
        )


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


def test_s104_write_basic(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 2, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5)
    )
    src_ds.GetRasterBand(2).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
    )
    src_ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.25",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251104",
                "ISSUE_TIME=231403Z",
            ],
        )

    with gdal.Open(f'S104:"{tmp_path}/104xxxxyyyy.h5":Group_001') as ds:
        assert ds.GetSpatialRef().GetAuthorityCode() == "32631"
        assert ds.GetGeoTransform() == pytest.approx(
            (500000, 1.1, 0, 4500000 + 1.2 * 3, 0, -1.2)
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(1).ReadRaster()) == (
            7.5,
            8.5,
            9.5,
            4.5,
            5.5,
            6.5,
            1.5,
            2.5,
            3.5,
        )
        assert struct.unpack("B" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == {
            "AREA_OR_POINT": "Point",
            "COMMON_POINT_RULE_DEFINITION": "return all the attribute values that can be determined for the position",
            "COMMON_POINT_RULE_NAME": "all",
            "DATA_DYNAMICITY_NAME": "hydrodynamicForecast",
            "DATA_DYNAMICITY_DEFINITION": "Values calculated from a two- or three-dimensional dynamic simulation of future conditions using predicted data for boundary forcing, via statistical method or combination",
            "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
            "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.",
            "commonPointRule": "4",
            "dataDynamicity": "5",
            "dateTimeOfFirstRecord": "20251104T225013Z",
            "dateTimeOfLastRecord": "20251104T225013Z",
            "horizontalPositionUncertainty": "-1",
            "interpolationType": "1",
            "INTERPOLATION_TYPE_DEFINITION": "Assign the feature attribute value associated with the nearest domain object in the domain of the coverage",
            "INTERPOLATION_TYPE_NAME": "nearestneighbor",
            "issueDate": "20251104",
            "issueTime": "231403Z",
            "maxDatasetHeight": "9.5",
            "minDatasetHeight": "1.5",
            "numberOfTimes": "1",
            "timePoint": "20251104T225013Z",
            "uncertainty": "-1.000000",
            "verticalCS": "6498",
            "verticalUncertainty": "-1",
            "waterLevelTrendThreshold": "0.25",
        }

    validate(
        tmp_path / "104xxxxyyyy.h5",
        expected_check_count=57,
    )


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


def test_s104_write_with_uncertainty_band(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 3, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).SetNoDataValue(-1.0)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, -1.0, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5)
    )
    src_ds.GetRasterBand(2).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
    )
    src_ds.GetRasterBand(3).SetNoDataValue(-1.0)
    src_ds.GetRasterBand(3).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, -1.0, 2, 3, 4, 5, 6, 7, 8, 9)
    )
    src_ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.25",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251104",
                "ISSUE_TIME=231403Z",
            ],
        )

    with gdal.Open(f'S104:"{tmp_path}/104xxxxyyyy.h5":Group_001') as ds:
        assert ds.GetSpatialRef().GetAuthorityCode() == "32631"
        assert ds.GetGeoTransform() == pytest.approx(
            (500000, 1.1, 0, 4500000 + 1.2 * 3, 0, -1.2)
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(1).ReadRaster()) == (
            7.5,
            8.5,
            9.5,
            4.5,
            5.5,
            6.5,
            -9999.0,
            2.5,
            3.5,
        )
        assert struct.unpack("B" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(3).ReadRaster()) == (
            7,
            8,
            9,
            4,
            5,
            6,
            -1.0,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == {
            "VERTICAL_DATUM_NAME": "meanLowerLowWater",
            "VERTICAL_DATUM_ABBREV": "MLLW",
            "VERTICAL_DATUM_DEFINITION": "The average height of the lower low waters at a place over a 19-year period.",
            "issueDate": "20251104",
            "issueTime": "231403Z",
            "waterLevelTrendThreshold": "0.25",
            "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
            "verticalCS": "6498",
            "minDatasetHeight": "2.5",
            "maxDatasetHeight": "9.5",
            "horizontalPositionUncertainty": "-1",
            "verticalUncertainty": "-1",
            "commonPointRule": "4",
            "interpolationType": "1",
            "COMMON_POINT_RULE_NAME": "all",
            "COMMON_POINT_RULE_DEFINITION": "return all the attribute values that can be determined for the position",
            "INTERPOLATION_TYPE_NAME": "nearestneighbor",
            "INTERPOLATION_TYPE_DEFINITION": "Assign the feature attribute value associated with the nearest domain object in the domain of the coverage",
            "dateTimeOfFirstRecord": "20251104T225013Z",
            "dateTimeOfLastRecord": "20251104T225013Z",
            "numberOfTimes": "1",
            "dataDynamicity": "5",
            "DATA_DYNAMICITY_NAME": "hydrodynamicForecast",
            "DATA_DYNAMICITY_DEFINITION": "Values calculated from a two- or three-dimensional dynamic simulation of future conditions using predicted data for boundary forcing, via statistical method or combination",
            "timePoint": "20251104T225013Z",
            "AREA_OR_POINT": "Point",
        }

    validate(
        tmp_path / "104xxxxyyyy.h5",
        expected_check_count=57,
    )


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


@pytest.mark.parametrize(
    "proj4,out_proj4",
    [
        ("+proj=longlat +ellps=GRS80 +pm=paris +no_defs", None),
        (
            "+proj=merc +lat_ts=1.5 +lon_0=2.5 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=merc +lat_0=0 +lon_0=2.5 +k=0.99 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
            "+proj=merc +lat_ts=8.13653121977138 +lon_0=2.5 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
        ),
        (
            "+proj=tmerc +lat_0=1.5 +lon_0=2.5 +k=0.99 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=omerc +lat_0=1.5 +lonc=2.5 +alpha=3.5 +gamma=4.5 +k=0.99 +x_0=5.5 +y_0=6.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=omerc +no_uoff +lat_0=1.5 +lonc=2.5 +alpha=3.5 +gamma=4.5 +k=0.99 +x_0=5.5 +y_0=6.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=lcc +lat_0=1.5 +lon_0=4.5 +lat_1=2.5 +lat_2=3.5 +x_0=5.5 +y_0=6.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=lcc +lat_1=49 +lat_0=49 +lon_0=4.5 +k_0=0.99 +x_0=5.5 +y_0=6.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=sterea +lat_0=1.5 +lon_0=2.5 +k=0.9 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=stere +lat_0=90 +lon_0=2.5 +k=0.9 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=krovak +axis=swu +lat_0=49.5 +lon_0=42.5 +alpha=30.2881397527778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs",
            None,
        ),
        (
            "+proj=poly +lat_0=1.5 +lon_0=2.5 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=aea +lat_0=1.5 +lon_0=4.5 +lat_1=2.5 +lat_2=3.5 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
        (
            "+proj=laea +lat_0=1.5 +lon_0=2.5 +x_0=3.5 +y_0=4.5 +datum=WGS84 +units=m +no_defs",
            None,
        ),
    ],
)
def test_s104_write_custom_crs(tmp_path, proj4, out_proj4):

    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 2, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5)
    )
    src_ds.GetRasterBand(2).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
    )
    if proj4 == "+proj=longlat +ellps=GRS80 +pm=paris +no_defs":
        src_ds.SetGeoTransform([2, 1.1, 0, 49, 0, 1.2])
    else:
        src_ds.SetGeoTransform([1000, 1.1, 0, 10000, 0, 1.2])
    srs = osr.SpatialReference()
    srs.ImportFromProj4(proj4)
    src_ds.SetSpatialRef(srs)

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.2",
                "DATA_DYNAMICITY=5",
            ],
        )

    ds = gdal.Open(f'S104:"{tmp_path}/104xxxxyyyy.h5":Group_001')
    assert ds.GetSpatialRef().ExportToProj4() == (out_proj4 if out_proj4 else proj4)

    validate(
        tmp_path / "104xxxxyyyy.h5",
        expected_check_count=56,
    )


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


@pytest.mark.parametrize(
    "value,warning_msg,validate_warning",
    [
        (
            -100.5,
            "Range of water level height in the dataset is [-100.500000, 9.500000] whereas the allowed range is [-99.99, 99.99]",
            "/WaterLevel: minDatasetHeight=-100.5 should be in [-99.99, 99.99] range",
        ),
        (
            100.5,
            "Range of water level height in the dataset is [2.500000, 100.500000] whereas the allowed range is [-99.99, 99.99]",
            "/WaterLevel: maxDatasetHeight=100.5 should be in [-99.99, 99.99] range",
        ),
    ],
)
def test_s104_write_out_of_range_water_height(
    tmp_path, value, warning_msg, validate_warning
):

    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 2, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, value, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5)
    )
    src_ds.GetRasterBand(2).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
    )
    src_ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(gdal.CE_Warning, match=warning_msg):
        gdal.Translate(
            tmp_path / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.25",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251104",
                "ISSUE_TIME=231403Z",
            ],
        )

    validate(tmp_path / "104xxxxyyyy.h5", expected_warnings=[validate_warning])


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


def test_s104_write_out_of_range_uncertainty(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 3, gdal.GDT_Float32)
    src_ds.GetRasterBand(1).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5)
    )
    src_ds.GetRasterBand(2).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
    )
    src_ds.GetRasterBand(3).WriteRaster(
        0, 0, 3, 3, struct.pack("f" * 9, -10, 2, 3, 0, 1, 2, 2, 0, 1)
    )
    src_ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    with gdaltest.error_raised(
        gdal.CE_Warning,
        match="Negative uncertainty value found (-10.000000), which is not allowed (except nodata value -1.0)",
    ):
        gdal.Translate(
            tmp_path / "104xxxxyyyy.h5",
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.25",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251104",
                "ISSUE_TIME=231403Z",
            ],
        )

    validate(
        tmp_path / "104xxxxyyyy.h5",
        expected_errors=[
            (
                "Error",
                "/WaterLevel/WaterLevel.01/Group_001/values : minimum uncertainty is -10.0, whereas it should be >= 0",
            )
        ],
    )


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


def test_s104_write_large_file(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 1200, 1200, 3, gdal.GDT_Float32)
    src_ds.SetGeoTransform([500000, 1, 0, 4500000, 0, 1])
    src_ds.SetSpatialRef(osr.SpatialReference(epsg=32631))

    filename = str(tmp_path / "104xxxxyyyy.h5")
    with gdaltest.error_raised(
        gdal.CE_Warning,
        match="file size exceeds 10 MB",
    ):
        gdal.Translate(
            filename,
            src_ds,
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "TIME_POINT=20251104T225013Z",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.25",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251104",
                "ISSUE_TIME=231403Z",
                "COMPRESS=NONE",
            ],
        )

    expected_warnings = [
        f"File size of {filename} = 12988072, which exceeds 10 MB",
    ]
    validate(filename, expected_warnings=expected_warnings)


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


def test_s104_write_multiple_timestamps(tmp_path):

    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "in1.tif", 3, 3, 2, gdal.GDT_Float32
    ) as ds:
        ds.GetRasterBand(1).WriteRaster(
            0,
            0,
            3,
            3,
            struct.pack("f" * 9, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5),
        )
        ds.GetRasterBand(2).WriteRaster(
            0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
        )
        ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
        ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        ds.SetMetadataItem("timePoint", "20251104T120000Z")

    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "in2.tif", 3, 3, 2, gdal.GDT_Float32
    ) as ds:
        ds.GetRasterBand(1).WriteRaster(
            0,
            0,
            3,
            3,
            struct.pack("f" * 9, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5),
        )
        ds.GetRasterBand(2).WriteRaster(
            0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
        )
        ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
        ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        ds.SetMetadataItem("timePoint", "20251104T130000Z")

    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "in3.tif", 3, 3, 2, gdal.GDT_Float32
    ) as ds:
        ds.GetRasterBand(1).WriteRaster(
            0,
            0,
            3,
            3,
            struct.pack("f" * 9, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5),
        )
        ds.GetRasterBand(2).WriteRaster(
            0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
        )
        ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
        ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        ds.SetMetadataItem("timePoint", "20251104T140000Z")

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "104xxxxyyyy.h5",
            gdal.Open(tmp_path / "in1.tif"),
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.25",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251104",
                "ISSUE_TIME=231403Z",
                "UNCERTAINTY=1",
                f"DATASETS={tmp_path}/in1.tif,{tmp_path}/in2.tif,{tmp_path}/in3.tif",
            ],
        )

    with gdal.Open(f'S104:"{tmp_path}/104xxxxyyyy.h5":Group_001') as ds:
        assert ds.GetSpatialRef().GetAuthorityCode() == "32631"
        assert ds.GetGeoTransform() == pytest.approx(
            (500000, 1.1, 0, 4500000 + 1.2 * 3, 0, -1.2)
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(1).ReadRaster()) == (
            7.5,
            8.5,
            9.5,
            4.5,
            5.5,
            6.5,
            1.5,
            2.5,
            3.5,
        )
        assert struct.unpack("B" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == {
            "VERTICAL_DATUM_NAME": "meanLowerLowWater",
            "VERTICAL_DATUM_ABBREV": "MLLW",
            "VERTICAL_DATUM_DEFINITION": "The average height of the lower low waters at a place over a 19-year period.",
            "issueDate": "20251104",
            "issueTime": "231403Z",
            "waterLevelTrendThreshold": "0.25",
            "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
            "verticalCS": "6498",
            "minDatasetHeight": "1.5",
            "maxDatasetHeight": "9.5",
            "horizontalPositionUncertainty": "-1",
            "verticalUncertainty": "-1",
            "commonPointRule": "4",
            "interpolationType": "1",
            "COMMON_POINT_RULE_NAME": "all",
            "COMMON_POINT_RULE_DEFINITION": "return all the attribute values that can be determined for the position",
            "INTERPOLATION_TYPE_NAME": "nearestneighbor",
            "INTERPOLATION_TYPE_DEFINITION": "Assign the feature attribute value associated with the nearest domain object in the domain of the coverage",
            "timeRecordInterval": "3600",
            "dateTimeOfFirstRecord": "20251104T120000Z",
            "dateTimeOfLastRecord": "20251104T140000Z",
            "numberOfTimes": "3",
            "dataDynamicity": "5",
            "DATA_DYNAMICITY_NAME": "hydrodynamicForecast",
            "DATA_DYNAMICITY_DEFINITION": "Values calculated from a two- or three-dimensional dynamic simulation of future conditions using predicted data for boundary forcing, via statistical method or combination",
            "timePoint": "20251104T120000Z",
            "uncertainty": "1.000000",
            "AREA_OR_POINT": "Point",
        }

    validate(
        tmp_path / "104xxxxyyyy.h5",
        expected_check_count=57,
    )

    # Test S104->S104 translation
    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "104xxxxyyyybis.h5",
            gdal.Open(tmp_path / "104xxxxyyyy.h5"),
            format="S104",
        )

    with gdal.Open(f'S104:"{tmp_path}/104xxxxyyyybis.h5":Group_001') as ds:
        assert ds.GetSpatialRef().GetAuthorityCode() == "32631"
        assert ds.GetGeoTransform() == pytest.approx(
            (500000, 1.1, 0, 4500000 + 1.2 * 3, 0, -1.2)
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(1).ReadRaster()) == (
            7.5,
            8.5,
            9.5,
            4.5,
            5.5,
            6.5,
            1.5,
            2.5,
            3.5,
        )
        assert struct.unpack("B" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == {
            "VERTICAL_DATUM_NAME": "meanLowerLowWater",
            "VERTICAL_DATUM_ABBREV": "MLLW",
            "VERTICAL_DATUM_DEFINITION": "The average height of the lower low waters at a place over a 19-year period.",
            "issueDate": "20251104",
            "issueTime": "231403Z",
            "waterLevelTrendThreshold": "0.25",
            "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
            "verticalCS": "6498",
            "minDatasetHeight": "1.5",
            "maxDatasetHeight": "9.5",
            "horizontalPositionUncertainty": "-1",
            "verticalUncertainty": "-1",
            "commonPointRule": "4",
            "interpolationType": "1",
            "COMMON_POINT_RULE_NAME": "all",
            "COMMON_POINT_RULE_DEFINITION": "return all the attribute values that can be determined for the position",
            "INTERPOLATION_TYPE_NAME": "nearestneighbor",
            "INTERPOLATION_TYPE_DEFINITION": "Assign the feature attribute value associated with the nearest domain object in the domain of the coverage",
            "timeRecordInterval": "3600",
            "dateTimeOfFirstRecord": "20251104T120000Z",
            "dateTimeOfLastRecord": "20251104T140000Z",
            "numberOfTimes": "3",
            "dataDynamicity": "5",
            "DATA_DYNAMICITY_NAME": "hydrodynamicForecast",
            "DATA_DYNAMICITY_DEFINITION": "Values calculated from a two- or three-dimensional dynamic simulation of future conditions using predicted data for boundary forcing, via statistical method or combination",
            "timePoint": "20251104T120000Z",
            "uncertainty": "1.000000",
            "AREA_OR_POINT": "Point",
        }

    validate(
        tmp_path / "104xxxxyyyybis.h5",
        expected_check_count=57,
    )


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


def test_s104_write_multiple_vertical_datums(tmp_path):

    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "in1.tif", 3, 3, 2, gdal.GDT_Float32
    ) as ds:
        ds.GetRasterBand(1).WriteRaster(
            0,
            0,
            3,
            3,
            struct.pack("f" * 9, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5),
        )
        ds.GetRasterBand(2).WriteRaster(
            0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
        )
        ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
        ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        ds.SetMetadataItem("timePoint", "20251104T120000Z")

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "104xxxxyyyy.h5",
            gdal.Open(tmp_path / "in1.tif"),
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=MLLW",
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.25",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251104",
                "ISSUE_TIME=231403Z",
            ],
        )

    with gdal.GetDriverByName("GTiff").Create(
        tmp_path / "in2.tif", 3, 3, 2, gdal.GDT_Float32
    ) as ds:
        ds.GetRasterBand(1).WriteRaster(
            0,
            0,
            3,
            3,
            struct.pack("f" * 9, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5),
        )
        ds.GetRasterBand(2).WriteRaster(
            0, 0, 3, 3, struct.pack("f" * 9, 1, 2, 3, 0, 1, 2, 2, 0, 1)
        )
        ds.SetGeoTransform([500000, 1.1, 0, 4500000, 0, 1.2])
        ds.SetSpatialRef(osr.SpatialReference(epsg=32631))
        ds.SetMetadataItem("timePoint", "20251104T120000Z")

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "104xxxxyyyy.h5",
            gdal.Open(tmp_path / "in2.tif"),
            format="S104",
            creationOptions=[
                "VERTICAL_DATUM=1027",  # EGM2008 geoid
                "VERTICAL_CS=DEPTH",
                "WATER_LEVEL_TREND_THRESHOLD=0.25",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251104",
                "ISSUE_TIME=231403Z",
                "APPEND_SUBDATASET=YES",
            ],
        )

    with gdal.Open(f'S104:"{tmp_path}/104xxxxyyyy.h5":WaterLevel.01:Group_001') as ds:
        assert ds.GetSpatialRef().GetAuthorityCode() == "32631"
        assert ds.GetGeoTransform() == pytest.approx(
            (500000, 1.1, 0, 4500000 + 1.2 * 3, 0, -1.2)
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(1).ReadRaster()) == (
            7.5,
            8.5,
            9.5,
            4.5,
            5.5,
            6.5,
            1.5,
            2.5,
            3.5,
        )
        assert struct.unpack("B" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        dict1 = {
            "VERTICAL_DATUM_NAME": "meanLowerLowWater",
            "VERTICAL_DATUM_ABBREV": "MLLW",
            "VERTICAL_DATUM_DEFINITION": "The average height of the lower low waters at a place over a 19-year period.",
            "issueDate": "20251104",
            "issueTime": "231403Z",
            "waterLevelTrendThreshold": "0.25",
            "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
            "verticalCS": "6498",
            "minDatasetHeight": "1.5",
            "maxDatasetHeight": "9.5",
            "horizontalPositionUncertainty": "-1",
            "verticalUncertainty": "-1",
            "commonPointRule": "4",
            "interpolationType": "1",
            "COMMON_POINT_RULE_NAME": "all",
            "COMMON_POINT_RULE_DEFINITION": "return all the attribute values that can be determined for the position",
            "INTERPOLATION_TYPE_NAME": "nearestneighbor",
            "INTERPOLATION_TYPE_DEFINITION": "Assign the feature attribute value associated with the nearest domain object in the domain of the coverage",
            "dateTimeOfFirstRecord": "20251104T120000Z",
            "dateTimeOfLastRecord": "20251104T120000Z",
            "numberOfTimes": "1",
            "dataDynamicity": "5",
            "DATA_DYNAMICITY_NAME": "hydrodynamicForecast",
            "DATA_DYNAMICITY_DEFINITION": "Values calculated from a two- or three-dimensional dynamic simulation of future conditions using predicted data for boundary forcing, via statistical method or combination",
            "timePoint": "20251104T120000Z",
            "uncertainty": "-1.000000",
            "AREA_OR_POINT": "Point",
        }
        assert ds.GetMetadata_Dict() == dict1

    with gdal.Open(f'S104:"{tmp_path}/104xxxxyyyy.h5":WaterLevel.02:Group_001') as ds:
        assert ds.GetSpatialRef().GetAuthorityCode() == "32631"
        assert ds.GetGeoTransform() == pytest.approx(
            (500000, 1.1, 0, 4500000 + 1.2 * 3, 0, -1.2)
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(1).ReadRaster()) == (
            7.5,
            8.5,
            9.5,
            4.5,
            5.5,
            6.5,
            1.5,
            2.5,
            3.5,
        )
        assert struct.unpack("B" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        dict2 = {
            "VERTICAL_DATUM_DEFINITION": "The average height of the lower low waters at a place over a 19-year period.",
            "issueDate": "20251104",
            "issueTime": "231403Z",
            "waterLevelTrendThreshold": "0.25",
            "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
            "verticalCS": "6498",
            "minDatasetHeight": "1.5",
            "maxDatasetHeight": "9.5",
            "horizontalPositionUncertainty": "-1",
            "verticalUncertainty": "-1",
            "commonPointRule": "4",
            "interpolationType": "1",
            "COMMON_POINT_RULE_NAME": "all",
            "COMMON_POINT_RULE_DEFINITION": "return all the attribute values that can be determined for the position",
            "INTERPOLATION_TYPE_NAME": "nearestneighbor",
            "INTERPOLATION_TYPE_DEFINITION": "Assign the feature attribute value associated with the nearest domain object in the domain of the coverage",
            "dateTimeOfFirstRecord": "20251104T120000Z",
            "dateTimeOfLastRecord": "20251104T120000Z",
            "numberOfTimes": "1",
            "dataDynamicity": "5",
            "DATA_DYNAMICITY_NAME": "hydrodynamicForecast",
            "DATA_DYNAMICITY_DEFINITION": "Values calculated from a two- or three-dimensional dynamic simulation of future conditions using predicted data for boundary forcing, via statistical method or combination",
            "VERTICAL_DATUM_EPSG_CODE": "1027",
            "VERTICAL_DATUM_NAME": "EGM2008 geoid",
            "timePoint": "20251104T120000Z",
            "uncertainty": "-1.000000",
            "AREA_OR_POINT": "Point",
        }
        assert ds.GetMetadata_Dict() == dict2

    validate(
        tmp_path / "104xxxxyyyy.h5",
        expected_check_count=57,
    )

    # Test S104->S104 translation
    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "104xxxxyyyybis.h5",
            gdal.Open(tmp_path / "104xxxxyyyy.h5"),
            format="S104",
        )

    with gdal.Open(
        f'S104:"{tmp_path}/104xxxxyyyybis.h5":WaterLevel.01:Group_001'
    ) as ds:
        assert ds.GetSpatialRef().GetAuthorityCode() == "32631"
        assert ds.GetGeoTransform() == pytest.approx(
            (500000, 1.1, 0, 4500000 + 1.2 * 3, 0, -1.2)
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(1).ReadRaster()) == (
            7.5,
            8.5,
            9.5,
            4.5,
            5.5,
            6.5,
            1.5,
            2.5,
            3.5,
        )
        assert struct.unpack("B" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == dict1

    with gdal.Open(
        f'S104:"{tmp_path}/104xxxxyyyybis.h5":WaterLevel.02:Group_001'
    ) as ds:
        assert ds.GetSpatialRef().GetAuthorityCode() == "32631"
        assert ds.GetGeoTransform() == pytest.approx(
            (500000, 1.1, 0, 4500000 + 1.2 * 3, 0, -1.2)
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(1).ReadRaster()) == (
            7.5,
            8.5,
            9.5,
            4.5,
            5.5,
            6.5,
            1.5,
            2.5,
            3.5,
        )
        assert struct.unpack("B" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == dict2

    validate(
        tmp_path / "104xxxxyyyybis.h5",
        expected_check_count=57,
    )
