#!/usr/bin/env pytest
# -*- coding: utf-8 -*-
###############################################################################
#
# Project:  GDAL/OGR Test Suite
# Purpose:  Test read functionality for S111 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("S111")


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


def test_s111_basic():
    filename = "data/s111/test_s111_v1.2.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() == {
        "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": "2023-12-31",
        "geographicIdentifier": "Somewhere",
        "producer": "Generated by autotest/gdrivers/data/s111/generate_test.py (not strictly fully S111 compliant)",
        "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
        "verticalCS": "6498",
        "minDatasetCurrentSpeed": "1",
        "maxDatasetCurrentSpeed": "2",
        "timeRecordInterval": "3600",
        "dateTimeOfFirstRecord": "20190606T120000Z",
        "dateTimeOfLastRecord": "20190606T120000Z",
        "numberOfTimes": "1",
        "AREA_OR_POINT": "Point",
    }

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

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

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

    band = ds.GetRasterBand(1)
    assert band.GetDescription() == "surfaceCurrentSpeed"
    assert band.GetNoDataValue() == -123
    assert band.GetUnitType() == "knots"
    assert struct.unpack("f" * 6, band.ReadRaster()) == (
        3.0,
        4.0,
        5.0,
        -123.0,
        1.0,
        2.0,
    )
    assert band.GetDefaultRAT() is not None

    band = ds.GetRasterBand(2)
    assert band.GetDescription() == "surfaceCurrentDirection"
    assert band.GetNoDataValue() == -123
    assert band.GetUnitType() == "degree"
    assert struct.unpack("f" * 6, band.ReadRaster()) == (3, 2, 1, 0, 1, 2)

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

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


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


def test_s111_north_up_no():
    filename = "data/s111/test_s111_v1.2.h5"
    ds = gdal.OpenEx(f'S111:"{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() == "surfaceCurrentSpeed"
    assert band.GetNoDataValue() == -123
    assert band.GetUnitType() == "knots"
    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() == "surfaceCurrentDirection"
    assert band.GetNoDataValue() == -123
    assert band.GetUnitType() == "degree"
    assert struct.unpack("f" * 6, band.ReadRaster()) == (
        0,
        1,
        2,
        3,
        2,
        1,
    )

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


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


def test_s111_multidim():

    filename = "data/s111/test_s111_v1.2.h5"
    ds = gdal.OpenEx(filename, gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    ar = rg.OpenMDArrayFromFullname(
        "/SurfaceCurrent/SurfaceCurrent.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


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


def test_s111_multiple_feature_instance_groups():

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

    ds = gdal.Open(
        'S111:"data/s111/multiple_feature_instance_groups.h5":SurfaceCurrent.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("f" * 8, ds.GetRasterBand(2).ReadRaster()) == (
        1,
        2,
        0,
        1,
        0,
        1,
        2,
        0,
    )
    assert ds.GetMetadata() == {
        "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": "2025-10-07",
        "issueTime": "12:34:56",
        "geographicIdentifier": "world",
        "depthTypeIndex": "1",
        "surfaceCurrentDepth": "-4.5",
        "DEPTH_TYPE_INDEX_NAME": "heightOrDepth",
        "DEPTH_TYPE_INDEX_DEFINITION": "Height or depth",
        "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
        "verticalCS": "6498",
        "minDatasetCurrentSpeed": "0",
        "maxDatasetCurrentSpeed": "7",
        "commonPointRule": "3",
        "interpolationType": "10",
        "COMMON_POINT_RULE_NAME": "high",
        "COMMON_POINT_RULE_DEFINITION": "use the greatest of the attribute values",
        "INTERPOLATION_TYPE_NAME": "discrete",
        "INTERPOLATION_TYPE_DEFINITION": "No interpolation method applies to the coverage",
        "timeRecordInterval": "3600",
        "dateTimeOfFirstRecord": "20190606T120000Z",
        "dateTimeOfLastRecord": "20190606T120000Z",
        "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",
        "uncertaintySurfaceCurrentSpeed": "-1.000000",
        "uncertaintySurfaceCurrentDirection": "-1.000000",
        "timePoint": "20190606T120000Z",
        "AREA_OR_POINT": "Point",
    }

    ds = gdal.Open(
        'S111:"data/s111/multiple_feature_instance_groups.h5":SurfaceCurrent.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("f" * 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_s111
    except ImportError:
        print("Cannot import validate_s111")
        return True

    errors, warnings, checks_done = validate_s111.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_s111_validator():

    # Fake product: many unconformities
    expected_errors = [
        (
            "Critical error",
            "/Group_F/SurfaceCurrent: row 0, got value '['surfaceCurrentSpeed', 'Surface Current Speed', 'knot', '-9999.00', 'H5T_FLOAT', '0.00', '', 'geSemiInterval']', which is not in '[['surfaceCurrentSpeed', 'Surface Current Speed', 'knot', '-9999.00', 'H5T_FLOAT', '0.00', '99.00', 'geSemiInterval'], ['surfaceCurrentDirection', 'Surface Current Direction', 'degree', '-9999.0', 'H5T_FLOAT', '0.0', '359.9', 'closedInterval'], ['surfaceCurrentTime', 'Surface Current Time', '', '00010101T000000Z', 'H5T_STRING', '19000101T000000Z', '21500101T000000Z', 'closedInterval'], ['speedUncertainty', 'Speed Uncertainty', 'knot', '-1.0', 'H5T_FLOAT', '0.00', '99.00', 'geSemiInterval'], ['directionUncertainty', 'Direction Uncertainty', 'degree', '-1.0', 'H5T_FLOAT', '0.0', '359.9', 'closedInterval']]'",
        ),
        (
            "Critical error",
            "/Group_F/SurfaceCurrent: row 1, got value '['surfaceCurrentDirection', 'Surface Current Direction', 'degree', '-9999.00', 'H5T_FLOAT', '0.00', '359.9', 'closedInterval']', which is not in '[['surfaceCurrentSpeed', 'Surface Current Speed', 'knot', '-9999.00', 'H5T_FLOAT', '0.00', '99.00', 'geSemiInterval'], ['surfaceCurrentDirection', 'Surface Current Direction', 'degree', '-9999.0', 'H5T_FLOAT', '0.0', '359.9', 'closedInterval'], ['surfaceCurrentTime', 'Surface Current Time', '', '00010101T000000Z', 'H5T_STRING', '19000101T000000Z', '21500101T000000Z', 'closedInterval'], ['speedUncertainty', 'Speed Uncertainty', 'knot', '-1.0', 'H5T_FLOAT', '0.00', '99.00', 'geSemiInterval'], ['directionUncertainty', 'Direction Uncertainty', 'degree', '-1.0', 'H5T_FLOAT', '0.0', '359.9', 'closedInterval']]'",
        ),
        ("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 'depthTypeIndex' is not an enumeration"),
        ("Error", "top level attribute 'surfaceCurrentDepth' is not a float32"),
        ("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",
            "SurfaceCurrent group attribute 'dataCodingFormat' is not an enumeration",
        ),
        ("Error", "SurfaceCurrent group attribute 'dimension' is not a uint8"),
        (
            "Error",
            "SurfaceCurrent group attribute 'commonPointRule' is not an enumeration",
        ),
        (
            "Error",
            "SurfaceCurrent group attribute 'horizontalPositionUncertainty' is not a float32",
        ),
        (
            "Error",
            "SurfaceCurrent group attribute 'verticalUncertainty' is not a float32",
        ),
        ("Error", "SurfaceCurrent group attribute 'numInstances' is not a uint32"),
        (
            "Error",
            "SurfaceCurrent group attribute 'minDatasetCurrentSpeed' is not a float64",
        ),
        (
            "Error",
            "SurfaceCurrent group attribute 'maxDatasetCurrentSpeed' is not a float64",
        ),
        (
            "Error",
            "SurfaceCurrent group attribute 'sequencingRule.type' is not an enumeration",
        ),
        (
            "Error",
            "SurfaceCurrent group attribute 'interpolationType' is not an enumeration",
        ),
        (
            "Error",
            "SurfaceCurrent group attribute 'dataOffsetCode' is not an enumeration",
        ),
        (
            "Error",
            '/SurfaceCurrent/axisNames must conform to CRS. Expected ["Easting", "Northing"] or ["Latitude", "Longitude"]. Got [\'longitude\', \'latitude\']',
        ),
        (
            "Error",
            "SurfaceCurrent feature instance group /SurfaceCurrent/SurfaceCurrent.01 attribute 'dataDynamicity' is not an enumeration",
        ),
        (
            "Error",
            "SurfaceCurrent feature instance group /SurfaceCurrent/SurfaceCurrent.02 attribute 'dataDynamicity' is not an enumeration",
        ),
        (
            "Error",
            "/SurfaceCurrent/SurfaceCurrent.02/Group_001/values : maximum surfaceCurrentSpeed is 70.0, whereas maxDatasetCurrentSpeed attribute = 7.0",
        ),
    ]
    expected_warnings = [
        "File name should start with '111'",
        "File name 'multiple_feature_instance_groups.h5' does not match expected pattern '^111[a-zA-Z0-9]{4}[a-zA-Z0-9\\-_]{1,54}\\.(?:h5|H5)$'",
        "Extra element in SurfaceCurrent feature instance group /SurfaceCurrent/SurfaceCurrent.02 group: 'verticalDatum'",
        "Extra element in SurfaceCurrent feature instance group /SurfaceCurrent/SurfaceCurrent.02 group: 'verticalDatumReference'",
    ]

    validate(
        "data/s111/multiple_feature_instance_groups.h5",
        expected_errors=expected_errors,
        expected_warnings=expected_warnings,
    )


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


def test_s111_write_errors(tmp_vsimem):

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

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

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

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

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

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

    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="TIME_POINT creation option value must be set, or source dataset must have a timePoint metadata item",
    ):
        gdal.Translate(tmp_vsimem / "111xxxxyyyy.h5", src_ds, format="S111")

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

    with pytest.raises(Exception, match="DEPTH_TYPE creation option must be specified"):
        gdal.Translate(
            tmp_vsimem / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=["TIME_POINT=20251105T231000Z"],
        )

    with pytest.raises(
        Exception,
        match="DEPTH_TYPE creation option must be set to heightOrDepth/1 or layerAverage/2",
    ):
        with gdal.quiet_errors():
            gdal.Translate(
                tmp_vsimem / "111xxxxyyyy.h5",
                src_ds,
                format="S111",
                creationOptions=["TIME_POINT=20251105T231000Z", "DEPTH_TYPE=invalid"],
            )

    with pytest.raises(
        Exception,
        match="VERTICAL_CS creation option must be specified when DEPTH_TYPE = heightOrDepth",
    ):
        gdal.Translate(
            tmp_vsimem / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=["TIME_POINT=20251105T231000Z", "DEPTH_TYPE=1"],
        )

    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 / "111xxxxyyyy.h5",
                src_ds,
                format="S111",
                creationOptions=[
                    "TIME_POINT=20251105T231000Z",
                    "DEPTH_TYPE=1",
                    "VERTICAL_CS=invalid",
                ],
            )

    with pytest.raises(
        Exception,
        match="VERTICAL_DATUM creation option must be specified when DEPTH_TYPE = heightOrDepth",
    ):
        gdal.Translate(
            tmp_vsimem / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=1",
                "VERTICAL_CS=6498",
            ],
        )

    with pytest.raises(
        Exception,
        match="VERTICAL_DATUM value is invalid",
    ):
        gdal.Translate(
            tmp_vsimem / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=1",
                "VERTICAL_CS=6498",
                "VERTICAL_DATUM=invalid",
            ],
        )

    with pytest.raises(
        Exception,
        match="SURFACE_CURRENT_DEPTH creation option must be specified",
    ):
        gdal.Translate(
            tmp_vsimem / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=1",
                "VERTICAL_CS=6498",
                "VERTICAL_DATUM=MLLW",
            ],
        )

    with pytest.raises(
        Exception,
        match="SURFACE_CURRENT_DEPTH creation option value must be a numeric value",
    ):
        with gdal.quiet_errors():
            gdal.Translate(
                tmp_vsimem / "111xxxxyyyy.h5",
                src_ds,
                format="S111",
                creationOptions=[
                    "TIME_POINT=20251105T231000Z",
                    "DEPTH_TYPE=1",
                    "VERTICAL_CS=6498",
                    "VERTICAL_DATUM=MLLW",
                    "SURFACE_CURRENT_DEPTH=invalid",
                ],
            )

    with pytest.raises(
        Exception,
        match="DATA_DYNAMICITY creation option must be specified",
    ):
        gdal.Translate(
            tmp_vsimem / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=1",
                "VERTICAL_CS=6498",
                "VERTICAL_DATUM=MLLW",
                "SURFACE_CURRENT_DEPTH=2.0",
            ],
        )

    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 / "111xxxxyyyy.h5",
                src_ds,
                format="S111",
                creationOptions=[
                    "TIME_POINT=20251105T231000Z",
                    "DEPTH_TYPE=1",
                    "VERTICAL_CS=6498",
                    "VERTICAL_DATUM=MLLW",
                    "SURFACE_CURRENT_DEPTH=2.0",
                    "DATA_DYNAMICITY=invalid",
                ],
            )


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


def test_s111_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="S111 dataset filenames should start with '111'"
    ):
        gdal.Translate(
            tmp_vsimem / "non_conformant_prefix.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=1",
                "VERTICAL_CS=6498",
                "VERTICAL_DATUM=MLLW",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
            ],
        )

    with gdaltest.error_raised(
        gdal.CE_Warning, match="S111 dataset filenames should have a '.H5' extension"
    ):
        gdal.Translate(
            tmp_vsimem / "111xxxxyyyy.oops",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=1",
                "VERTICAL_CS=6498",
                "VERTICAL_DATUM=MLLW",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
            ],
        )


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


def test_s111_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 / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=1",
                "VERTICAL_CS=6498",
                "VERTICAL_DATUM=MLLW",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
            ],
        )

    with gdal.Open(f'S111:"{tmp_path}/111xxxxyyyy.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 ds.RasterCount == 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("f" * 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": "20251105",
            "issueTime": "222950Z",
            "depthTypeIndex": "1",
            "surfaceCurrentDepth": "2",
            "DEPTH_TYPE_INDEX_NAME": "heightOrDepth",
            "DEPTH_TYPE_INDEX_DEFINITION": "Height or depth",
            "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
            "verticalCS": "6498",
            "minDatasetCurrentSpeed": "1.5",
            "maxDatasetCurrentSpeed": "9.5",
            "commonPointRule": "3",
            "interpolationType": "10",
            "COMMON_POINT_RULE_NAME": "high",
            "COMMON_POINT_RULE_DEFINITION": "use the greatest of the attribute values",
            "INTERPOLATION_TYPE_NAME": "discrete",
            "INTERPOLATION_TYPE_DEFINITION": "No interpolation method applies to the coverage",
            "dateTimeOfFirstRecord": "20251105T231000Z",
            "dateTimeOfLastRecord": "20251105T231000Z",
            "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",
            "uncertaintySurfaceCurrentSpeed": "-1.000000",
            "uncertaintySurfaceCurrentDirection": "-1.000000",
            "timePoint": "20251105T231000Z",
            "AREA_OR_POINT": "Point",
        }

    validate(
        tmp_path / "111xxxxyyyy.h5",
        expected_check_count=59,
    )


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


def test_s111_write_with_uncertainty_bands(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 4, 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.5, 20.5, 30.5, 40.5, 50.5, 60.5, 70.5, 80.5, 90.5),
    )
    src_ds.GetRasterBand(4).WriteRaster(
        0,
        0,
        3,
        3,
        struct.pack(
            "f" * 9, 110.5, 120.5, 130.5, 140.5, 150.5, 160.5, 170.5, 180.5, 190.5
        ),
    )
    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 / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=1",
                "VERTICAL_CS=6498",
                "VERTICAL_DATUM=MLLW",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
            ],
        )

    with gdal.Open(f'S111:"{tmp_path}/111xxxxyyyy.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 ds.RasterCount == 4
        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("f" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(3).ReadRaster()) == (
            70.5,
            80.5,
            90.5,
            40.5,
            50.5,
            60.5,
            10.5,
            20.5,
            30.5,
        )
        assert struct.unpack("f" * 9, ds.GetRasterBand(4).ReadRaster()) == (
            170.5,
            180.5,
            190.5,
            140.5,
            150.5,
            160.5,
            110.5,
            120.5,
            130.5,
        )
        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": "20251105",
            "issueTime": "222950Z",
            "depthTypeIndex": "1",
            "surfaceCurrentDepth": "2",
            "DEPTH_TYPE_INDEX_NAME": "heightOrDepth",
            "DEPTH_TYPE_INDEX_DEFINITION": "Height or depth",
            "VERTICAL_CS_DEFINITION": "depth, meters, orientation down",
            "verticalCS": "6498",
            "minDatasetCurrentSpeed": "1.5",
            "maxDatasetCurrentSpeed": "9.5",
            "commonPointRule": "3",
            "interpolationType": "10",
            "COMMON_POINT_RULE_NAME": "high",
            "COMMON_POINT_RULE_DEFINITION": "use the greatest of the attribute values",
            "INTERPOLATION_TYPE_NAME": "discrete",
            "INTERPOLATION_TYPE_DEFINITION": "No interpolation method applies to the coverage",
            "dateTimeOfFirstRecord": "20251105T231000Z",
            "dateTimeOfLastRecord": "20251105T231000Z",
            "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": "20251105T231000Z",
            "AREA_OR_POINT": "Point",
        }

    validate(
        tmp_path / "111xxxxyyyy.h5",
        expected_check_count=60,
    )


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


@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_s111_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 / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=2",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
            ],
        )

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

    validate(
        tmp_path / "111xxxxyyyy.h5",
        expected_check_count=59,
    )


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


@pytest.mark.parametrize(
    "value,warning_msg,expected_errors,expected_warnings",
    [
        (
            -0.5,
            "Range of surface current speed in the dataset is [-0.500000, 9.500000] whereas the allowed range is [0.00, 99.00]",
            [
                (
                    "Error",
                    "/SurfaceCurrent/SurfaceCurrent.01/Group_001/values : minimum surfaceCurrentSpeed is -0.5, whereas it should be >= 0",
                )
            ],
            [
                "/SurfaceCurrent: minDatasetCurrentSpeed=-0.5 should be in [0, 99.99] range"
            ],
        ),
        (
            100.5,
            "Range of surface current speed in the dataset is [2.500000, 100.500000] whereas the allowed range is [0.00, 99.00]",
            [
                (
                    "Error",
                    "/SurfaceCurrent/SurfaceCurrent.01/Group_001/values : maximum surfaceCurrentSpeed is 100.5, whereas it should be <= 99",
                )
            ],
            [
                "/SurfaceCurrent: maxDatasetCurrentSpeed=100.5 should be in [0, 99.99] range"
            ],
        ),
    ],
)
def test_s111_write_out_of_range_speed(
    tmp_path, value, warning_msg, expected_errors, expected_warnings
):

    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 / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=2",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
            ],
        )

    validate(
        tmp_path / "111xxxxyyyy.h5",
        expected_errors=expected_errors,
        expected_warnings=expected_warnings,
    )


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


@pytest.mark.parametrize(
    "value,warning_msg,expected_errors,expected_warnings",
    [
        (
            -0.5,
            "Range of surface current direction in the dataset is [-0.500000, 3.000000] whereas the allowed range is [0.00, 359.90]",
            [
                (
                    "Error",
                    "/SurfaceCurrent/SurfaceCurrent.01/Group_001/values : minimum surfaceCurrentDirection is -0.5, whereas it should be >= 0",
                )
            ],
            [],
        ),
        (
            360.0,
            "Range of surface current direction in the dataset is [0.000000, 360.000000] whereas the allowed range is [0.00, 359.90]",
            [
                (
                    "Error",
                    "/SurfaceCurrent/SurfaceCurrent.01/Group_001/values : maximum surfaceCurrentDirection is 360.0, whereas it should be <= 359.9",
                )
            ],
            [],
        ),
    ],
)
def test_s111_write_out_of_range_dir(
    tmp_path, value, warning_msg, expected_errors, expected_warnings
):

    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, value, 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 / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=2",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
            ],
        )

    validate(
        tmp_path / "111xxxxyyyy.h5",
        expected_errors=expected_errors,
        expected_warnings=expected_warnings,
    )


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


def test_s111_write_out_of_range_speed_uncertainty(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 4, 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.GetRasterBand(4).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="Negative speed uncertainty value found (-10.000000), which is not allowed (except nodata value -1.0)",
    ):
        gdal.Translate(
            tmp_path / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=2",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
            ],
        )

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


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


def test_s111_write_out_of_range_direction_uncertainty(tmp_path):

    src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 4, 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, 1, 2, 3, 0, 1, 2, 2, 0, 1)
    )
    src_ds.GetRasterBand(4).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 direction uncertainty value found (-10.000000), which is not allowed (except nodata value -1.0)",
    ):
        gdal.Translate(
            tmp_path / "111xxxxyyyy.h5",
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=2",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
            ],
        )

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


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


def test_s111_write_large_file(tmp_path):

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

    filename = str(tmp_path / "111xxxxyyyy.h5")
    with gdaltest.error_raised(
        gdal.CE_Warning,
        match="file size exceeds 10 MB",
    ):
        gdal.Translate(
            filename,
            src_ds,
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=2",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
                "COMPRESS=NONE",
            ],
        )

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


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


def test_s111_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", "20251111T120000Z")

    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", "20251111T130000Z")

    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", "20251111T140000Z")

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "111xxxxyyyy.h5",
            gdal.Open(tmp_path / "in1.tif"),
            format="S111",
            creationOptions=[
                "TIME_POINT=20251105T231000Z",
                "DEPTH_TYPE=2",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
                "UNCERTAINTY_SPEED=1.5",
                "UNCERTAINTY_DIRECTION=2.5",
                f"DATASETS={tmp_path}/in1.tif,{tmp_path}/in2.tif,{tmp_path}/in3.tif",
            ],
        )

    with gdal.Open(f'S111:"{tmp_path}/111xxxxyyyy.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("f" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        expected_md = {
            "issueDate": "20251105",
            "issueTime": "222950Z",
            "depthTypeIndex": "2",
            "surfaceCurrentDepth": "2",
            "DEPTH_TYPE_INDEX_NAME": "layerAverage",
            "DEPTH_TYPE_INDEX_DEFINITION": "Layer average",
            "minDatasetCurrentSpeed": "1.5",
            "maxDatasetCurrentSpeed": "9.5",
            "commonPointRule": "3",
            "interpolationType": "10",
            "COMMON_POINT_RULE_NAME": "high",
            "COMMON_POINT_RULE_DEFINITION": "use the greatest of the attribute values",
            "INTERPOLATION_TYPE_NAME": "discrete",
            "INTERPOLATION_TYPE_DEFINITION": "No interpolation method applies to the coverage",
            "timeRecordInterval": "3600",
            "dateTimeOfFirstRecord": "20251111T120000Z",
            "dateTimeOfLastRecord": "20251111T140000Z",
            "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",
            "uncertaintySurfaceCurrentSpeed": "1.500000",
            "uncertaintySurfaceCurrentDirection": "2.500000",
            "timePoint": "20251111T120000Z",
            "AREA_OR_POINT": "Point",
        }
        assert ds.GetMetadata_Dict() == expected_md

    validate(
        tmp_path / "111xxxxyyyy.h5",
        expected_check_count=60,
    )

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

    with gdal.Open(f'S111:"{tmp_path}/111xxxxyyyybis.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("f" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == expected_md

    validate(
        tmp_path / "111xxxxyyyybis.h5",
        expected_check_count=60,
    )


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


def test_s111_write_multiple_instances(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", "20251111T120000Z")

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "111xxxxyyyy.h5",
            gdal.Open(tmp_path / "in1.tif"),
            format="S111",
            creationOptions=[
                "DEPTH_TYPE=2",
                "SURFACE_CURRENT_DEPTH=2.0",
                "DATA_DYNAMICITY=5",
                "ISSUE_DATE=20251105",
                "ISSUE_TIME=222950Z",
            ],
        )

    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", "20251111T120000Z")

    with gdaltest.error_raised(gdal.CE_None):
        gdal.Translate(
            tmp_path / "111xxxxyyyy.h5",
            gdal.Open(tmp_path / "in2.tif"),
            format="S111",
            creationOptions=[
                "DATA_DYNAMICITY=5",
                "APPEND_SUBDATASET=YES",
            ],
        )

    with gdal.Open(
        f'S111:"{tmp_path}/111xxxxyyyy.h5":SurfaceCurrent.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("f" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        dict1 = {
            "issueDate": "20251105",
            "issueTime": "222950Z",
            "depthTypeIndex": "2",
            "surfaceCurrentDepth": "2",
            "DEPTH_TYPE_INDEX_NAME": "layerAverage",
            "DEPTH_TYPE_INDEX_DEFINITION": "Layer average",
            "minDatasetCurrentSpeed": "1.5",
            "maxDatasetCurrentSpeed": "9.5",
            "commonPointRule": "3",
            "interpolationType": "10",
            "COMMON_POINT_RULE_NAME": "high",
            "COMMON_POINT_RULE_DEFINITION": "use the greatest of the attribute values",
            "INTERPOLATION_TYPE_NAME": "discrete",
            "INTERPOLATION_TYPE_DEFINITION": "No interpolation method applies to the coverage",
            "dateTimeOfFirstRecord": "20251111T120000Z",
            "dateTimeOfLastRecord": "20251111T120000Z",
            "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",
            "uncertaintySurfaceCurrentSpeed": "-1.000000",
            "uncertaintySurfaceCurrentDirection": "-1.000000",
            "timePoint": "20251111T120000Z",
            "AREA_OR_POINT": "Point",
        }
        assert ds.GetMetadata_Dict() == dict1

    with gdal.Open(
        f'S111:"{tmp_path}/111xxxxyyyy.h5":SurfaceCurrent.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("f" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        dict2 = {
            "issueDate": "20251105",
            "issueTime": "222950Z",
            "depthTypeIndex": "2",
            "surfaceCurrentDepth": "2",
            "DEPTH_TYPE_INDEX_NAME": "layerAverage",
            "DEPTH_TYPE_INDEX_DEFINITION": "Layer average",
            "minDatasetCurrentSpeed": "1.5",
            "maxDatasetCurrentSpeed": "9.5",
            "commonPointRule": "3",
            "interpolationType": "10",
            "COMMON_POINT_RULE_NAME": "high",
            "COMMON_POINT_RULE_DEFINITION": "use the greatest of the attribute values",
            "INTERPOLATION_TYPE_NAME": "discrete",
            "INTERPOLATION_TYPE_DEFINITION": "No interpolation method applies to the coverage",
            "dateTimeOfFirstRecord": "20251111T120000Z",
            "dateTimeOfLastRecord": "20251111T120000Z",
            "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",
            "uncertaintySurfaceCurrentSpeed": "-1.000000",
            "uncertaintySurfaceCurrentDirection": "-1.000000",
            "timePoint": "20251111T120000Z",
            "AREA_OR_POINT": "Point",
        }

        assert ds.GetMetadata_Dict() == dict2

    validate(
        tmp_path / "111xxxxyyyy.h5",
        expected_check_count=59,
    )

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

    with gdal.Open(
        f'S111:"{tmp_path}/111xxxxyyyybis.h5":SurfaceCurrent.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("f" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == dict1

    with gdal.Open(
        f'S111:"{tmp_path}/111xxxxyyyybis.h5":SurfaceCurrent.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("f" * 9, ds.GetRasterBand(2).ReadRaster()) == (
            2,
            0,
            1,
            0,
            1,
            2,
            1,
            2,
            3,
        )
        assert ds.GetMetadata_Dict() == dict2

    validate(
        tmp_path / "111xxxxyyyybis.h5",
        expected_check_count=59,
    )
