#!/usr/bin/env pytest
###############################################################################
#
# Project:  GDAL/OGR Test Suite
# Purpose:  Test multidimensional support in HDF5 driver
# Author:   Even Rouault <even.rouault@spatialys.com>
#
###############################################################################
# Copyright (c) 2019, Even Rouault <even.rouault@spatialys.com>
#
# SPDX-License-Identifier: MIT
###############################################################################

import struct

import gdaltest
import pytest

from osgeo import gdal

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


def test_hdf5_multidim_basic():

    ds = gdal.OpenEx("data/hdf5/u8be.h5", gdal.OF_MULTIDIM_RASTER)
    assert ds
    rg = ds.GetRootGroup()
    assert rg
    assert not rg.GetGroupNames()
    with pytest.raises(Exception):
        rg.OpenGroup("non_existing")
    assert rg.GetMDArrayNames() == ["TestArray"]
    with pytest.raises(Exception):
        rg.OpenMDArray("non_existing")
    ar = rg.OpenMDArray("TestArray")
    assert ar
    with pytest.raises(Exception):
        ar.GetAttribute("non_existing")
    dims = ar.GetDimensions()
    assert len(dims) == 2
    assert dims[0].GetSize() == 6
    assert dims[1].GetSize() == 5
    assert ar.GetDataType().GetNumericDataType() == gdal.GDT_Byte

    got_data = ar.Read(buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16))
    assert len(got_data) == 30 * 2
    assert struct.unpack("H" * 30, got_data) == (
        0,
        1,
        2,
        3,
        4,
        1,
        2,
        3,
        4,
        5,
        2,
        3,
        4,
        5,
        6,
        3,
        4,
        5,
        6,
        7,
        4,
        5,
        6,
        7,
        8,
        5,
        6,
        7,
        8,
        9,
    )

    got_data = ar.Read(
        array_start_idx=[2, 1],
        count=[3, 2],
        buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
    )
    assert len(got_data) == 6 * 2
    assert struct.unpack("H" * 6, got_data) == (3, 4, 4, 5, 5, 6)

    got_data = ar.Read(
        array_start_idx=[2, 1],
        count=[3, 2],
        array_step=[-1, -1],
        buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
    )
    assert len(got_data) == 6 * 2
    assert struct.unpack("H" * 6, got_data) == (3, 2, 2, 1, 1, 0)

    with gdaltest.config_option("GDAL_HDF5_TEMP_ARRAY_ALLOC_SIZE", "0"):
        got_data = ar.Read(
            array_start_idx=[2, 1],
            count=[3, 2],
            array_step=[-1, -1],
            buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
        )
        assert len(got_data) == 6 * 2
        assert struct.unpack("H" * 6, got_data) == (3, 2, 2, 1, 1, 0)


def test_hdf5_multidim_var_alldatatypes():

    ds = gdal.OpenEx("HDF5:data/netcdf/alldatatypes.nc", gdal.OF_MULTIDIM_RASTER)
    assert ds
    rg = ds.GetRootGroup()
    assert rg

    expected_vars = [  # ('char_var', gdal.GDT_Byte, (ord('x'),ord('y'))),
        ("ubyte_var", gdal.GDT_Byte, (255, 254)),
        # ("byte_var", gdal.GDT_Int16, (-128, -127)),
        # ("byte_unsigned_false_var", gdal.GDT_Int8, (-128, -127)),
        # ('byte_unsigned_true_var', gdal.GDT_Byte, (128, 129)),
        ("ushort_var", gdal.GDT_UInt16, (65534, 65533)),
        ("short_var", gdal.GDT_Int16, (-32768, -32767)),
        ("uint_var", gdal.GDT_UInt32, (4294967294, 4294967293)),
        ("int_var", gdal.GDT_Int32, (-2147483648, -2147483647)),
        ("uint64_var", gdal.GDT_UInt64, (18446744073709551613, 18446744073709551612)),
        ("int64_var", gdal.GDT_Int64, (-9223372036854775808, -9223372036854775807)),
        ("float_var", gdal.GDT_Float32, (1.25, 2.25)),
        ("double_var", gdal.GDT_Float64, (1.25125, 2.25125)),
        ("complex_int16_var", gdal.GDT_CInt16, (-32768, -32767, -32766, -32765)),
        (
            "complex_int32_var",
            gdal.GDT_CInt32,
            (-2147483648, -2147483647, -2147483646, -2147483645),
        ),
        ("complex64_var", gdal.GDT_CFloat32, (1.25, 2.5, 2.25, 3.5)),
        ("complex128_var", gdal.GDT_CFloat64, (1.25125, 2.25125, 3.25125, 4.25125)),
    ]
    for var_name, dt, val in expected_vars:
        var = rg.OpenMDArray(var_name)
        assert var
        assert var.GetDataType().GetClass() == gdal.GEDTC_NUMERIC, var_name
        assert var.GetDataType().GetNumericDataType() == dt, var_name
        if dt == gdal.GDT_Byte:
            assert struct.unpack("B" * len(val), var.Read()) == val
        if dt == gdal.GDT_Int8:
            assert struct.unpack("b" * len(val), var.Read()) == val
        if dt == gdal.GDT_UInt16:
            assert struct.unpack("H" * len(val), var.Read()) == val
        if dt == gdal.GDT_Int16:
            assert struct.unpack("h" * len(val), var.Read()) == val
        if dt == gdal.GDT_UInt32:
            assert struct.unpack("I" * len(val), var.Read()) == val
        if dt == gdal.GDT_Int32:
            assert struct.unpack("i" * len(val), var.Read()) == val
        if dt == gdal.GDT_UInt64:
            assert struct.unpack("Q" * len(val), var.Read()) == val
        if dt == gdal.GDT_Int64:
            assert struct.unpack("q" * len(val), var.Read()) == val
        if dt == gdal.GDT_Float32:
            assert struct.unpack("f" * len(val), var.Read()) == val
        if dt == gdal.GDT_Float64:
            assert struct.unpack("d" * len(val), var.Read()) == val
        if dt == gdal.GDT_CInt16:
            assert struct.unpack("h" * len(val), var.Read()) == val
        if dt == gdal.GDT_CInt32:
            assert struct.unpack("i" * len(val), var.Read()) == val
        if dt == gdal.GDT_CFloat32:
            assert struct.unpack("f" * len(val), var.Read()) == val
        if dt == gdal.GDT_CFloat64:
            assert struct.unpack("d" * len(val), var.Read()) == val

    var = rg.OpenMDArray("custom_type_2_elts_var")
    dt = var.GetDataType()
    assert dt.GetClass() == gdal.GEDTC_COMPOUND
    assert dt.GetSize() == 8
    assert dt.GetName() == "custom_type_2_elts"
    comps = dt.GetComponents()
    assert len(comps) == 2
    assert comps[0].GetName() == "x"
    assert comps[0].GetOffset() == 0
    assert comps[0].GetType().GetNumericDataType() == gdal.GDT_Int32
    assert comps[1].GetName() == "y"
    assert comps[1].GetOffset() == 4
    assert comps[1].GetType().GetNumericDataType() == gdal.GDT_Int16
    data = var.Read()
    assert len(data) == 2 * 8
    assert struct.unpack("ihihh", data) == (1, 2, 3, 4, 0)

    assert struct.unpack("i" * 2, var["x"].Read()) == (1, 3)
    assert struct.unpack("h" * 2, var["y"].Read()) == (2, 4)

    var = rg.OpenMDArray("custom_type_3_elts_var")
    dt = var.GetDataType()
    assert dt.GetClass() == gdal.GEDTC_COMPOUND
    assert dt.GetSize() == 12
    comps = dt.GetComponents()
    assert len(comps) == 3
    assert comps[0].GetName() == "x"
    assert comps[0].GetOffset() == 0
    assert comps[0].GetType().GetNumericDataType() == gdal.GDT_Int32
    assert comps[1].GetName() == "y"
    assert comps[1].GetOffset() == 4
    assert comps[1].GetType().GetNumericDataType() == gdal.GDT_Int16
    assert comps[2].GetName() == "z"
    assert comps[2].GetOffset() == 8
    assert comps[2].GetType().GetNumericDataType() == gdal.GDT_Float32
    data = var.Read()
    assert len(data) == 2 * 12
    assert struct.unpack("ihf" * 2, data) == (1, 2, 3.5, 4, 5, 6.5)

    var = rg.OpenMDArray("string_var")
    dt = var.GetDataType()
    assert dt.GetClass() == gdal.GEDTC_STRING
    assert var.Read() == ["abcd", "ef"]


def test_hdf5_multidim_read_array():

    ds = gdal.OpenEx("HDF5:data/netcdf/alldatatypes.nc", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()

    # 0D
    var = rg.OpenMDArray("ubyte_no_dim_var")
    assert var
    assert struct.unpack("B", var.Read()) == (2,)
    assert struct.unpack(
        "H", var.Read(buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16))
    ) == (2,)

    # 1D
    var = rg.OpenMDArray("ubyte_x2_var")
    data = var.Read(array_start_idx=[1], count=[2], array_step=[2])
    got_data_ref = (2, 4)
    assert struct.unpack("B" * len(data), data) == got_data_ref

    data = var.Read(
        array_start_idx=[1],
        count=[2],
        array_step=[2],
        buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
    )
    assert struct.unpack("H" * (len(data) // 2), data) == got_data_ref

    data = var.Read(array_start_idx=[2], count=[2], array_step=[-2])
    got_data_ref = (3, 1)
    assert struct.unpack("B" * len(data), data) == got_data_ref

    # 2D
    var = rg.OpenMDArray("ubyte_y2_x2_var")
    data = var.Read(count=[2, 3], array_step=[2, 1])
    got_data_ref = (1, 2, 3, 9, 10, 11)
    assert struct.unpack("B" * len(data), data) == got_data_ref

    data = var.Read(
        count=[2, 3],
        array_step=[2, 1],
        buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
    )
    assert struct.unpack("H" * (len(data) // 2), data) == got_data_ref

    data = var.Read(array_start_idx=[1, 2], count=[2, 2])
    got_data_ref = (7, 8, 11, 12)
    assert struct.unpack("B" * len(data), data) == got_data_ref

    data = var.Read(
        array_start_idx=[1, 2],
        count=[2, 2],
        buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
    )
    assert struct.unpack("H" * (len(data) // 2), data) == got_data_ref

    # 3D
    var = rg.OpenMDArray("ubyte_z2_y2_x2_var")
    data = var.Read(count=[3, 2, 3], array_step=[1, 2, 1])
    got_data_ref = (1, 2, 3, 9, 10, 11, 2, 3, 4, 10, 11, 12, 3, 4, 5, 11, 12, 1)
    assert struct.unpack("B" * len(data), data) == got_data_ref

    data = var.Read(
        count=[3, 2, 3],
        array_step=[1, 2, 1],
        buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
    )
    assert struct.unpack("H" * (len(data) // 2), data) == got_data_ref

    data = var.Read(array_start_idx=[1, 1, 1], count=[3, 2, 2])
    got_data_ref = (7, 8, 11, 12, 8, 9, 12, 1, 9, 10, 1, 2)
    assert struct.unpack("B" * len(data), data) == got_data_ref

    data = var.Read(
        array_start_idx=[1, 1, 1],
        count=[3, 2, 2],
        buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
    )
    assert struct.unpack("H" * (len(data) // 2), data) == got_data_ref

    # Test reading from slice (most optimized path)
    data = var.Read(array_start_idx=[3, 0, 0], count=[1, 2, 3])
    data_from_slice = var[3].Read(count=[2, 3])
    assert data_from_slice == data

    # Test reading from slice (slow path)
    data = var.Read(
        array_start_idx=[3, 0 + (2 - 1) * 2, 0 + (3 - 1) * 1],
        count=[1, 2, 3],
        array_step=[1, -2, -1],
    )
    data_from_slice = var[3].Read(
        array_start_idx=[0 + (2 - 1) * 2, 0 + (3 - 1) * 1],
        count=[2, 3],
        array_step=[-2, -1],
    )
    assert data_from_slice == data

    # 4D
    var = rg.OpenMDArray("ubyte_t2_z2_y2_x2_var")
    data = var.Read(count=[2, 3, 2, 3], array_step=[1, 1, 2, 1])
    got_data_ref = (
        1,
        2,
        3,
        9,
        10,
        11,
        2,
        3,
        4,
        10,
        11,
        12,
        3,
        4,
        5,
        11,
        12,
        1,
        2,
        3,
        4,
        10,
        11,
        12,
        3,
        4,
        5,
        11,
        12,
        1,
        4,
        5,
        6,
        12,
        1,
        2,
    )
    assert struct.unpack("B" * len(data), data) == got_data_ref

    data = var.Read(
        count=[2, 3, 2, 3],
        array_step=[1, 1, 2, 1],
        buffer_datatype=gdal.ExtendedDataType.Create(gdal.GDT_UInt16),
    )
    assert struct.unpack("H" * (len(data) // 2), data) == got_data_ref


def test_hdf5_multidim_read_stride():

    ds = gdal.OpenEx("HDF5:data/netcdf/alldatatypes.nc", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()

    var = rg.OpenMDArray("ubyte_z2_y2_x2_var")
    assert struct.unpack("B" * 119, var.Read(buffer_stride=[3 * 4 * 2, 4 * 2, 1 * 2]))[
        ::2
    ] == struct.unpack("B" * 60, var.Read())


def test_hdf5_multidim_attr_alldatatypes():

    ds = gdal.OpenEx("HDF5:data/netcdf/alldatatypes.nc", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()

    attrs = rg.GetAttributes()
    assert len(attrs) == 1
    assert attrs[0].GetName() == "global_attr"

    attrs = rg.OpenGroup("group").GetAttributes()
    assert len(attrs) == 1
    assert attrs[0].GetName() == "group_global_attr"

    var = rg.OpenMDArray("ubyte_var")
    assert var
    attrs = var.GetAttributes()
    assert len(attrs) == 30
    map_attrs = {}
    for attr in attrs:
        map_attrs[attr.GetName()] = attr

    attr = map_attrs["attr_byte"]
    assert attr.GetDimensionCount() == 0
    assert len(attr.GetDimensionsSize()) == 0
    assert attr.GetDataType().GetNumericDataType() == gdal.GDT_Int8
    assert attr.Read() == -128

    assert map_attrs["attr_ubyte"].Read() == 255

    assert map_attrs["attr_char"].GetDataType().GetClass() == gdal.GEDTC_STRING
    assert map_attrs["attr_char"].GetDimensionCount() == 0
    assert map_attrs["attr_char"].Read() == "x"
    assert map_attrs["attr_char"].ReadAsStringArray() == ["x"]
    with pytest.raises(Exception):
        map_attrs["attr_char"].ReadAsRaw()

    assert (
        map_attrs["attr_string_as_repeated_char"].GetDataType().GetClass()
        == gdal.GEDTC_STRING
    )
    assert map_attrs["attr_string_as_repeated_char"].GetDimensionCount() == 0
    assert map_attrs["attr_string_as_repeated_char"].Read() == "xy"
    assert map_attrs["attr_string_as_repeated_char"].ReadAsStringArray() == ["xy"]

    assert map_attrs["attr_empty_char"].GetDataType().GetClass() == gdal.GEDTC_STRING
    assert map_attrs["attr_empty_char"].GetDimensionCount() == 0
    assert map_attrs["attr_empty_char"].Read() == ""

    assert map_attrs["attr_two_strings"].GetDataType().GetClass() == gdal.GEDTC_STRING
    assert map_attrs["attr_two_strings"].GetDimensionCount() == 1
    assert map_attrs["attr_two_strings"].GetDimensionsSize()[0] == 2
    assert map_attrs["attr_two_strings"].Read() == ["ab", "cd"]
    assert map_attrs["attr_two_strings"].ReadAsString() == "ab"

    assert map_attrs["attr_int"].Read() == -2147483647

    assert map_attrs["attr_float"].Read() == 1.25

    assert map_attrs["attr_double"].Read() == 1.25125
    assert map_attrs["attr_double"].ReadAsDoubleArray() == (1.25125,)

    assert map_attrs["attr_int64"].Read() == -9.223372036854776e18

    assert map_attrs["attr_uint64"].Read() == 1.8446744073709552e19

    assert (
        map_attrs["attr_complex_int16"].GetDataType().GetNumericDataType()
        == gdal.GDT_CInt16
    )
    assert struct.unpack("H" * 2, map_attrs["attr_complex_int16"].ReadAsRaw()) == (1, 2)
    assert map_attrs["attr_complex_int16"].Read() == 1.0
    assert map_attrs["attr_complex_int16"].ReadAsString() == "1+2j"

    assert map_attrs["attr_two_bytes"].Read() == (-128, -127)

    assert map_attrs["attr_two_ints"].Read() == (-2147483648, -2147483647)

    assert (
        map_attrs["attr_two_doubles"].GetDataType().GetNumericDataType()
        == gdal.GDT_Float64
    )
    assert map_attrs["attr_two_doubles"].GetDimensionCount() == 1
    assert map_attrs["attr_two_doubles"].GetDimensionsSize()[0] == 2
    assert map_attrs["attr_two_doubles"].Read() == (1.25125, 2.125125)
    assert map_attrs["attr_two_doubles"].ReadAsDouble() == 1.25125

    assert (
        map_attrs["attr_enum_ubyte"].GetDataType().GetNumericDataType() == gdal.GDT_Byte
    )
    assert map_attrs["attr_enum_ubyte"].Read() == 1

    assert (
        map_attrs["attr_enum_int"].GetDataType().GetNumericDataType() == gdal.GDT_Int32
    )
    assert map_attrs["attr_enum_int"].Read() == 1000000001

    assert len(map_attrs["attr_custom_type_2_elts"].ReadAsRaw()) == 8

    # Compound type contains a string
    with pytest.raises(Exception):
        map_attrs["attr_custom_with_string"].ReadAsRaw()


def test_hdf5_multidim_nodata_unit():

    ds = gdal.OpenEx("HDF5:data/netcdf/trmm-nc4.nc", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()

    ar = rg.OpenMDArray("pcp")
    assert (
        struct.unpack("f", ar.GetNoDataValueAsRaw())[0]
        == struct.unpack("f", struct.pack("f", -9999.9))[0]
    )
    assert ar.GetUnit() == ""

    ar = rg.OpenMDArray("longitude")
    assert ar.GetNoDataValueAsRaw() is None
    assert ar.GetUnit() == "degrees_east"


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


def test_hdf5_multidim_read_missing_value_of_different_type():
    def create():

        import h5py

        f = h5py.File("data/hdf5/FillValue_of_different_type.h5", "w")
        dset = f.create_dataset("test", (1,), dtype="f")
        dset.attrs.create("_FillValue", -9999, dtype="d")
        f.close()

    # create()

    ds = gdal.OpenEx(
        "data/hdf5/FillValue_of_different_type.h5",
        gdal.OF_MULTIDIM_RASTER,
    )
    rg = ds.GetRootGroup()
    var = rg.OpenMDArray("test")
    assert var.GetDataType().GetNumericDataType() == gdal.GDT_Float32
    assert (
        var.GetAttribute("_FillValue").GetDataType().GetNumericDataType()
        == gdal.GDT_Float64
    )
    assert var.GetNoDataValue() == -9999.0


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


def test_hdf5_multidim_read_missing_value_of_different_type_not_in_range():
    def create():

        import h5py

        f = h5py.File("data/hdf5/FillValue_of_different_type_not_in_range.h5", "w")
        dset = f.create_dataset("test", (1,), dtype="B")
        dset.attrs.create("_FillValue", -9999, dtype="d")
        f.close()

    # create()

    ds = gdal.OpenEx(
        "data/hdf5/FillValue_of_different_type_not_in_range.h5",
        gdal.OF_MULTIDIM_RASTER,
    )
    rg = ds.GetRootGroup()
    with gdal.quiet_errors():
        gdal.ErrorReset()
        var = rg.OpenMDArray("test")
        assert gdal.GetLastErrorMsg() != ""
        assert var.GetDataType().GetNumericDataType() == gdal.GDT_Byte
        assert (
            var.GetAttribute("_FillValue").GetDataType().GetNumericDataType()
            == gdal.GDT_Float64
        )
        assert var.GetNoDataValue() is None


def test_hdf5_multidim_recursive_groups():

    # File generated with
    # import h5py
    # f = h5py.File('hdf5/recursive_groups.h5','w')
    # group = f.create_group("subgroup")
    # group['link_to_root'] = f
    # group['link_to_self'] = group
    # group['soft_link_to_root'] = h5py.SoftLink('/')
    # group['soft_link_to_self'] = h5py.SoftLink('/subgroup')
    # group['soft_link_to_not_existing'] = h5py.SoftLink('/not_existing')
    # group['hard_link_to_root'] = h5py.HardLink('/')
    # group['ext_link_to_self_root'] = h5py.ExternalLink("hdf5/recursive_groups.h5", "/")
    # f.close()

    ds = gdal.OpenEx("data/hdf5/recursive_groups.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    assert rg.GetGroupNames() == ["subgroup"]


def test_hdf5_netcdf_dimensions():

    ds = gdal.OpenEx("HDF5:data/netcdf/trmm-nc4.nc", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()

    assert rg.GetAttribute("CDI")

    dims = rg.GetDimensions()
    assert len(dims) == 3

    assert dims[0].GetFullName() == "/latitude"
    assert dims[0].GetName() == "latitude"
    assert dims[0].GetSize() == 40
    var = dims[0].GetIndexingVariable()
    assert var
    assert var.GetName() == "latitude"

    assert dims[1].GetFullName() == "/longitude"
    assert dims[1].GetName() == "longitude"
    assert dims[1].GetSize() == 40
    var = dims[1].GetIndexingVariable()
    assert var

    assert dims[2].GetFullName() == "/time"
    assert dims[2].GetName() == "time"
    assert dims[2].GetSize() == 1
    var = dims[2].GetIndexingVariable()
    assert var

    ar = rg.OpenMDArray("latitude")
    dims = ar.GetDimensions()
    assert len(dims) == 1
    assert dims[0].GetName() == "latitude"
    assert dims[0].GetIndexingVariable()
    assert dims[0].GetIndexingVariable().GetName() == "latitude"

    ar = rg.OpenMDArray("pcp")
    dims = ar.GetDimensions()
    assert len(dims) == 3
    assert dims[0].GetName() == "time"
    assert dims[0].GetIndexingVariable()
    assert dims[0].GetIndexingVariable().GetName() == "time"


def test_hdf5_multidim_netcdf_dimensions_complex_case():

    ds = gdal.OpenEx("HDF5:data/netcdf/alldatatypes.nc", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    dims = rg.GetDimensions()
    assert len(dims) == 6

    dim = next((x for x in dims if x.GetName() == "X"), None)
    assert dim
    assert dim.GetFullName() == "/X"
    assert dim.GetSize() == 2
    assert dim.GetIndexingVariable()
    assert dim.GetIndexingVariable().GetFullName() == "/X"

    dim = next((x for x in dims if x.GetName() == "Y"), None)
    assert dim
    assert dim.GetFullName() == "/Y"
    assert dim.GetSize() == 1
    assert dim.GetIndexingVariable()
    assert dim.GetIndexingVariable().GetFullName() == "/Y"

    dim = next((x for x in dims if x.GetName() == "X2"), None)
    assert dim
    assert dim.GetFullName() == "/X2"
    assert dim.GetSize() == 4
    assert not dim.GetIndexingVariable()

    assert "X" in rg.GetMDArrayNames()
    assert "X2" not in rg.GetMDArrayNames()

    subgroup = rg.OpenGroup("group")
    assert subgroup
    dims = subgroup.GetDimensions()
    assert len(dims) == 2

    dim = next((x for x in dims if x.GetName() == "X"), None)
    assert dim
    assert dim.GetFullName() == "/group/X"
    assert dim.GetSize() == 3
    assert not dim.GetIndexingVariable()

    dim = next((x for x in dims if x.GetName() == "Y"), None)
    assert dim
    assert dim.GetFullName() == "/group/Y"
    assert dim.GetSize() == 2
    assert not dim.GetIndexingVariable()

    ar = subgroup.OpenMDArray("char_var")
    assert ar
    dims = ar.GetDimensions()
    assert len(dims) == 3
    assert dims[0].GetFullName() == "/Y"
    assert dims[0].GetSize() == 1
    assert dims[0].GetIndexingVariable()
    assert dims[1].GetFullName() == "/group/Y"
    assert dims[1].GetSize() == 2
    assert not dims[1].GetIndexingVariable()
    assert dims[2].GetFullName() == "/group/X"
    assert dims[2].GetSize() == 3
    assert not dims[2].GetIndexingVariable()

    ar = rg.OpenMDArray("ubyte_z2_y2_x2_var")
    assert ar
    dims = ar.GetDimensions()
    assert len(dims) == 3
    assert dims[0].GetFullName() == "/Z2"
    assert not dims[0].GetIndexingVariable()


def test_hdf5_multidim_dimension_labels_with_null():

    ds = gdal.OpenEx("data/hdf5/dimension_labels_with_null.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()

    ar = rg.OpenMDArray("data")
    attr = ar.GetAttribute("DIMENSION_LABELS")
    assert attr.ReadAsStringArray() == ["", "", "x"]

    ds = gdal.OpenEx("data/hdf5/dimension_labels_with_null.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    ar = rg.OpenMDArray("data")
    dims = ar.GetDimensions()
    assert len(dims) == 3
    assert dims[0].GetName() == "dim0"
    assert dims[1].GetName() == "dim1"
    assert dims[2].GetName() == "x"


def test_hdf5_multidim_family_driver():

    assert gdal.OpenEx("data/hdf5/test_family_0.h5", gdal.OF_MULTIDIM_RASTER)


def test_hdf5_multidim_read_transposed():

    ds = gdal.OpenEx("HDF5:data/netcdf/trmm-nc4.nc", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    ar = rg.OpenMDArray("pcp")
    dims = ar.GetDimensions()
    y_size = dims[-2].GetSize()
    x_size = dims[-1].GetSize()
    transposed_ar = ar.Transpose([0, 2, 1])
    data = list(struct.unpack("f" * (y_size * x_size), ar.Read()))
    transposed_data = list(struct.unpack("f" * (y_size * x_size), transposed_ar.Read()))
    manually_transposed_data = []
    for x in range(x_size):
        for y in range(y_size):
            manually_transposed_data.append(data[y * x_size + x])
    assert transposed_data == manually_transposed_data


##############################################################################
# Test opening a HDF5EOS grid file


def test_hdf5_multimdim_eos_grid_geo_projection():

    ds = gdal.OpenEx(
        "data/hdf5/dummy_HDFEOS_with_geo_projection.h5", gdal.OF_MULTIDIM_RASTER
    )
    ar = ds.GetRootGroup().OpenMDArrayFromFullname(
        "/HDFEOS/GRIDS/test/Data Fields/test"
    )
    dims = ar.GetDimensions()
    assert len(dims) == 2
    assert dims[0].GetName() == "YDim"
    assert dims[0].GetFullName() == "/HDFEOS/GRIDS/test/YDim"
    assert dims[0].GetType() == "HORIZONTAL_Y"
    assert dims[0].GetSize() == 20
    assert dims[1].GetName() == "XDim"
    assert dims[1].GetFullName() == "/HDFEOS/GRIDS/test/XDim"
    assert dims[1].GetType() == "HORIZONTAL_X"
    assert dims[1].GetSize() == 20
    ydim_var = dims[0].GetIndexingVariable()
    assert ydim_var
    assert struct.unpack("d" * 20, ydim_var.Read())[0:2] == pytest.approx(
        (33.90033388888937, 33.899796111111556)
    )
    xdim_var = dims[1].GetIndexingVariable()
    assert xdim_var
    assert struct.unpack("d" * 20, xdim_var.Read())[0:2] == pytest.approx(
        (-117.64084298610936, -117.64019006944282)
    )
    srs = ar.GetSpatialRef()
    assert srs is not None
    assert srs.IsGeographic()
    # WGS 84
    assert srs.GetSemiMajor() == 6378137
    assert srs.GetInvFlattening() == 298.257223563
    assert srs.GetDataAxisToSRSAxisMapping() == [1, 2]
    assert len(ar.GetCoordinateVariables()) == 0


##############################################################################
# Test opening a HDF5EOS grid file


def test_hdf5_multimdim_eos_grid_utm_projection():

    ds = gdal.OpenEx(
        "data/hdf5/dummy_HDFEOS_with_utm_projection.h5", gdal.OF_MULTIDIM_RASTER
    )
    ar = ds.GetRootGroup().OpenMDArrayFromFullname(
        "/HDFEOS/GRIDS/test/Data Fields/test"
    )
    srs = ar.GetSpatialRef()
    assert srs is not None
    assert srs.IsProjected()
    assert srs.GetDataAxisToSRSAxisMapping() == [2, 1]
    assert len(ar.GetCoordinateVariables()) == 0


###############################################################################
# Test opening a HDF5EOS swath file


def test_hdf5_multidim_eos_swath_no_explicit_dimension_map():

    ds = gdal.OpenEx("data/hdf5/dummy_HDFEOS_swath.h5", gdal.OF_MULTIDIM_RASTER)
    ar = ds.GetRootGroup().OpenMDArrayFromFullname(
        "/HDFEOS/SWATHS/MySwath/Data Fields/MyDataField"
    )
    dims = ar.GetDimensions()
    assert len(dims) == 3
    assert dims[0].GetName() == "Band"
    assert dims[0].GetFullName() == "/HDFEOS/SWATHS/MySwath/Band"
    assert dims[0].GetSize() == 2
    assert dims[1].GetName() == "AlongTrack"
    assert dims[1].GetFullName() == "/HDFEOS/SWATHS/MySwath/AlongTrack"
    assert dims[1].GetSize() == 3
    assert dims[2].GetName() == "CrossTrack"
    assert dims[2].GetFullName() == "/HDFEOS/SWATHS/MySwath/CrossTrack"
    assert dims[2].GetSize() == 4
    coordinates = ar.GetCoordinateVariables()
    assert len(coordinates) == 2
    assert (
        coordinates[0].GetFullName()
        == "/HDFEOS/SWATHS/MySwath/Geolocation Fields/Longitude"
    )
    assert (
        coordinates[1].GetFullName()
        == "/HDFEOS/SWATHS/MySwath/Geolocation Fields/Latitude"
    )


###############################################################################
# Test GetBlockSize() and GetStructuralInfo()


def test_hdf5_multidim_block_size_structural_info():

    ds = gdal.OpenEx("data/hdf5/deflate.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    var = rg.OpenMDArray("Band1")
    block_size = var.GetBlockSize()
    assert block_size == [1, 2]
    assert var.GetStructuralInfo() == {"COMPRESSION": "DEFLATE", "FILTER": "SHUFFLE"}


###############################################################################
# Test GetRawBlockInfo()


def test_hdf5_multidim_getrawblockinfo_chunked():

    if not gdal.GetDriverByName("HDF5").GetMetadataItem("HAVE_H5Dget_chunk_info"):
        pytest.skip("libhdf5 < 1.10.5")

    ds = gdal.OpenEx("data/hdf5/deflate.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    var = rg.OpenMDArray("Band1")
    block_size = var.GetBlockSize()
    assert block_size == [1, 2]
    y_size = var.GetDimensions()[0].GetSize()
    x_size = var.GetDimensions()[1].GetSize()

    info = var.GetRawBlockInfo([0, 0])
    assert info.GetOffset() == 13908
    assert info.GetSize() == 10
    assert info.GetFilename() == "data/hdf5/deflate.h5"
    assert info.GetInfo() == ["COMPRESSION=DEFLATE", "FILTER=SHUFFLE"]

    import zlib

    for y in range(y_size // block_size[0]):
        for x in range(x_size // block_size[1]):
            info = var.GetRawBlockInfo([y, x])
            with open(info.GetFilename(), "rb") as f:
                f.seek(info.GetOffset(), 0)
                data = f.read(info.GetSize())
                data = zlib.decompress(data)
                assert data == var.Read(
                    array_start_idx=[y * block_size[0], x * block_size[1]],
                    count=var.GetBlockSize(),
                ), (y, x)

    with pytest.raises(
        Exception, match="Invalid number of values in block_coordinates argument"
    ):
        var.GetRawBlockInfo([])

    assert var.GetRawBlockInfo([0, x_size // block_size[1] - 1]) is not None
    with pytest.raises(
        Exception, match=r"invalid block coordinate \(10\) for dimension 1"
    ):
        var.GetRawBlockInfo([0, x_size // block_size[1]])

    assert var.GetRawBlockInfo([y_size // block_size[0] - 1, 0]) is not None
    with pytest.raises(
        Exception, match=r"invalid block coordinate \(20\) for dimension 0"
    ):
        var.GetRawBlockInfo([y_size // block_size[0], 0])


###############################################################################
# Test GetRawBlockInfo() on non-chunked file


def test_hdf5_multidim_getrawblockinfo():

    ds = gdal.OpenEx("data/hdf5/u8be.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    var = rg.OpenMDArray("TestArray")
    info = var.GetRawBlockInfo([0, 0])
    assert info.GetOffset() == 2048
    assert info.GetSize() == 30
    assert info.GetFilename() == "data/hdf5/u8be.h5"
    assert info.GetInfo() is None

    with pytest.raises(Exception, match="invalid block coordinates. Should be all 0"):
        var.GetRawBlockInfo([0, 1])


###############################################################################
# Test GetRawBlockInfo() on a little-endian file


def test_hdf5_multidim_getrawblockinfo_little_endian():

    ds = gdal.OpenEx("data/hdf5/float32_little_endian.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    var = rg.OpenMDArray("test")
    info = var.GetRawBlockInfo([0, 0])
    assert info.GetInfo() == ["ENDIANNESS=LITTLE"]


###############################################################################
# Test GetRawBlockInfo() on a big-endian file


def test_hdf5_multidim_getrawblockinfo_big_endian():

    ds = gdal.OpenEx("data/hdf5/float32_big_endian.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    var = rg.OpenMDArray("test")
    info = var.GetRawBlockInfo([0, 0])
    assert info.GetInfo() == ["ENDIANNESS=BIG"]


###############################################################################
# Test reading a compound data type made of 2 Float16 values


def test_hdf5_multidim_read_cfloat16():

    ds = gdal.OpenEx("data/hdf5/complex.h5", gdal.OF_MULTIDIM_RASTER)
    rg = ds.GetRootGroup()
    var = rg.OpenMDArray("f16")
    assert var.GetDataType().GetNumericDataType() == gdal.GDT_CFloat32
    assert struct.unpack("f" * (5 * 5 * 2), var.Read()) == (
        0.0,
        0.0,
        1.0,
        1.0,
        2.0,
        2.0,
        3.0,
        3.0,
        4.0,
        4.0,
        5.0,
        5.0,
        6.0,
        6.0,
        7.0,
        7.0,
        8.0,
        8.0,
        9.0,
        9.0,
        10.0,
        10.0,
        11.0,
        11.0,
        12.0,
        12.0,
        13.0,
        13.0,
        14.0,
        14.0,
        15.0,
        15.0,
        16.0,
        16.0,
        17.0,
        17.0,
        18.0,
        18.0,
        19.0,
        19.0,
        20.0,
        20.0,
        21.0,
        21.0,
        22.0,
        22.0,
        23.0,
        23.0,
        24.0,
        24.0,
    )


##############################################################################
# Test opening a HDF5EOS grid file


def test_hdf5_multimdim_eos_grid_VIIRS_Grid_IMG_2D_issue_12941():

    ds = gdal.OpenEx(
        "data/hdf5/dummy_HDFEOS_IIRS_Grid_IMG_2D_issue_1294.h5", gdal.OF_MULTIDIM_RASTER
    )
    ar = ds.GetRootGroup().OpenMDArrayFromFullname(
        "/HDFEOS/GRIDS/test/Data Fields/test"
    )
    dims = ar.GetDimensions()
    assert len(dims) == 2
    assert dims[0].GetName() == "YDim"
    assert dims[0].GetFullName() == "/HDFEOS/GRIDS/test/YDim"
    assert dims[0].GetType() == "HORIZONTAL_Y"
    assert dims[0].GetSize() == 5
    assert dims[1].GetName() == "XDim"
    assert dims[1].GetFullName() == "/HDFEOS/GRIDS/test/XDim"
    assert dims[1].GetType() == "HORIZONTAL_X"
    assert dims[1].GetSize() == 4
    ydim_var = dims[0].GetIndexingVariable()
    assert ydim_var
    assert struct.unpack("d" * 5, ydim_var.Read())[0:2] == pytest.approx(
        (5448557.5463664, 5226167.4424332)
    )
    xdim_var = dims[1].GetIndexingVariable()
    assert xdim_var
    assert struct.unpack("d" * 4, xdim_var.Read())[0:2] == pytest.approx(
        (-972956.7047086251, -694969.0747918751)
    )
