#!/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, Even Rouault <even dot rouault at spatialys.com>
#
# SPDX-License-Identifier: MIT
###############################################################################

import os
import struct

import pytest

from osgeo import gdal

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(None) == "4326"
    assert ds.GetGeoTransform() == pytest.approx((1.8, 0.4, 0.0, 48.75, 0.0, -0.5))
    assert ds.GetMetadata_Dict() == {
        "AREA_OR_POINT": "Point",
        "verticalCS": "6498",
        "VERTICAL_CS_MEANING": "depth, meters, orientation down",
        "dateTimeOfFirstRecord": "20190606T120000Z",
        "dateTimeOfLastRecord": "20190606T120000Z",
        "geographicIdentifier": "Somewhere",
        "issueDate": "2023-12-31",
        "maxDatasetCurrentSpeed": "2",
        "minDatasetCurrentSpeed": "1",
        "numberOfTimes": "1",
        "producer": "Generated by autotest/gdrivers/data/s111/generate_test.py (not strictly fully S111 compliant)",
        "timeRecordInterval": "3600",
        "VERTICAL_DATUM_ABBREV": "MLLW",
        "VERTICAL_DATUM_MEANING": "meanLowerLowWater",
    }

    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(None) == "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(None) == "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() == {
        "AREA_OR_POINT": "Point",
        "DATA_DYNAMICITY_MEANING": "Hydrodynamic model forecast",
        "VERTICAL_DATUM_ABBREV": "MLLW",
        "VERTICAL_DATUM_MEANING": "meanLowerLowWater",
        "verticalCS": "6498",
        "VERTICAL_CS_MEANING": "depth, meters, orientation down",
        "dateTimeOfFirstRecord": "20190606T120000Z",
        "dateTimeOfLastRecord": "20190606T120000Z",
        "depthTypeIndex": "1",
        "geographicIdentifier": "world",
        "issueDate": "2025-10-07",
        "issueTime": "12:34:56",
        "maxDatasetCurrentSpeed": "7",
        "minDatasetCurrentSpeed": "0",
        "numberOfTimes": "1",
        "surfaceCurrentDepth": "-4.5",
        "timeRecordInterval": "3600",
        "uncertaintySurfaceCurrentDirection": "-1.000000",
        "uncertaintySurfaceCurrentSpeed": "-1.000000",
    }

    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_MEANING") == "lowWater"
