#!/usr/bin/env pytest
###############################################################################
# $Id$
#
# Project:  GDAL/OGR Test Suite
# Purpose:  Test basic read support for a all datatypes from a HDF file.
# Author:   Andrey Kiselev, dron@remotesensing.org
#
###############################################################################
# Copyright (c) 2003, Andrey Kiselev <dron@remotesensing.org>
# Copyright (c) 2009-2012, Even Rouault <even dot rouault at spatialys.com>
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Library General Public License for more details.
#
# You should have received a copy of the GNU Library General Public
# License along with this library; if not, write to the
# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
# Boston, MA 02111-1307, USA.
###############################################################################


import gdaltest
import pytest

from osgeo import gdal

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

init_list = [
    ("byte_3.hdf", 4672),
    ("int16_3.hdf", 4672),
    ("uint16_3.hdf", 4672),
    ("int32_3.hdf", 4672),
    ("uint32_3.hdf", 4672),
    ("float32_3.hdf", 4672),
    ("float64_3.hdf", 4672),
    ("utmsmall_3.hdf", 50054),
    ("byte_2.hdf", 4672),
    ("int16_2.hdf", 4672),
    ("uint16_2.hdf", 4672),
    ("int32_2.hdf", 4672),
    ("uint32_2.hdf", 4672),
    ("float32_2.hdf", 4672),
    ("float64_2.hdf", 4672),
    ("utmsmall_2.hdf", 50054),
]


@pytest.mark.parametrize(
    "filename,checksum",
    init_list,
    ids=[tup[0].split(".")[0] for tup in init_list],
)
@pytest.mark.require_driver("HDF4Image")
def test_hdf4_open(filename, checksum):
    ut = gdaltest.GDALTest("HDF4Image", filename, 1, checksum)
    ut.testOpen()


###############################################################################
# Test reading a GR dataset


def test_hdf4_read_gr():

    # Generated by https://support.hdfgroup.org/ftp/HDF/HDF_Current/src/unpacked/hdf/examples/GR_create_and_write_image.c
    ds = gdal.Open("data/General_RImages.hdf")
    assert ds
    assert ds.RasterCount == 2
    assert ds.GetRasterBand(1).Checksum() == 361
    assert not ds.GetRasterBand(1).GetColorTable()
    assert ds.GetRasterBand(2).Checksum() == 400


###############################################################################
# Test reading a GR dataset with a palette


def test_hdf4_read_gr_palette():

    # Generated by https://support.hdfgroup.org/ftp/HDF/HDF_Current/src/unpacked/hdf/examples/GR_write_palette.c
    ds = gdal.Open("data/Image_with_Palette.hdf")
    assert ds
    assert ds.RasterCount == 1
    assert ds.GetRasterBand(1).GetColorTable()


###############################################################################
# Test opening more than 32 simultaneous HDF4_EOS files


def test_hdf4_more_than_32_files():

    if (
        gdal.GetDriverByName("HDF4Image").GetMetadataItem("HDF4_HAS_MAXOPENFILES")
        != "YES"
    ):
        pytest.skip("HDF4_HAS_MAXOPENFILES not set")

    gdaltest.download_or_skip(
        "https://gamma.hdfgroup.org/ftp/pub/outgoing/NASAHDFfiles2/eosweb/hdf4/hdfeos2-swath-wo-dimmaps/AMSR_E_L2_Ocean_B01_200206182340_A.hdf",
        "AMSR_E_L2_Ocean_B01_200206182340_A.hdf",
    )

    tab = []
    for i in range(33):
        ds = gdal.Open(
            'HDF4_EOS:EOS_SWATH:"tmp/cache/AMSR_E_L2_Ocean_B01_200206182340_A.hdf":Swath1:Low_res_sst'
        )
        assert ds, i
        tab.append(ds)


###############################################################################
# Test HDF4_SDS with single subdataset


def test_hdf4_read_online_1():

    gdaltest.hdf4_drv = gdal.GetDriverByName("HDF4")

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    gdaltest.download_or_skip(
        "http://download.osgeo.org/gdal/data/hdf4/A2004259075000.L2_LAC_SST.hdf",
        "A2004259075000.L2_LAC_SST.hdf",
    )

    tst = gdaltest.GDALTest(
        "HDF4Image",
        "tmp/cache/A2004259075000.L2_LAC_SST.hdf",
        1,
        28189,
        filename_absolute=1,
    )

    return tst.testOpen()


###############################################################################
# Test HDF4_SDS with GEOLOCATION info


def test_hdf4_read_online_2():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    gdaltest.download_or_skip(
        "http://download.osgeo.org/gdal/data/hdf4/A2006005182000.L2_LAC_SST.x.hdf",
        "A2006005182000.L2_LAC_SST.x.hdf",
    )

    tst = gdaltest.GDALTest(
        "HDF4Image",
        'HDF4_SDS:UNKNOWN:"tmp/cache/A2006005182000.L2_LAC_SST.x.hdf":13',
        1,
        13209,
        filename_absolute=1,
    )

    tst.testOpen()

    ds = gdal.Open('HDF4_SDS:UNKNOWN:"tmp/cache/A2006005182000.L2_LAC_SST.x.hdf":13')
    md = ds.GetMetadata("GEOLOCATION")
    ds = None

    assert (
        md["X_DATASET"]
        == 'HDF4_SDS:UNKNOWN:"tmp/cache/A2006005182000.L2_LAC_SST.x.hdf":11'
    ), "Did not get expected X_DATASET"


###############################################################################
# Test HDF4_EOS:EOS_GRID


def test_hdf4_read_online_3():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    gdaltest.download_or_skip(
        "http://download.osgeo.org/gdal/data/hdf4/MO36MW14.chlor_MODIS.ADD2001089.004.2002186190207.hdf",
        "MO36MW14.chlor_MODIS.ADD2001089.004.2002186190207.hdf",
    )

    tst = gdaltest.GDALTest(
        "HDF4Image",
        "tmp/cache/MO36MW14.chlor_MODIS.ADD2001089.004.2002186190207.hdf",
        1,
        34723,
        filename_absolute=1,
    )

    tst.testOpen()

    ds = gdal.Open("tmp/cache/MO36MW14.chlor_MODIS.ADD2001089.004.2002186190207.hdf")
    gt = ds.GetGeoTransform()
    expected_gt = [-180.0, 0.3515625, 0.0, 90.0, 0.0, -0.3515625]
    for i in range(6):
        assert gt[i] == pytest.approx(
            expected_gt[i], abs=1e-8
        ), "did not get expected gt"

    srs = ds.GetProjectionRef()
    assert srs.find("Clarke") != -1, "did not get expected projection"

    ds = None


###############################################################################
# Test HDF4_SDS:SEAWIFS_L1A


def test_hdf4_read_online_4():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    gdaltest.download_or_skip(
        "http://download.osgeo.org/gdal/data/hdf4/S2002196124536.L1A_HDUN.BartonBendish.extract.hdf",
        "S2002196124536.L1A_HDUN.BartonBendish.extract.hdf",
    )

    tst = gdaltest.GDALTest(
        "HDF4Image",
        "tmp/cache/S2002196124536.L1A_HDUN.BartonBendish.extract.hdf",
        1,
        33112,
        filename_absolute=1,
    )

    tst.testOpen()

    ds = gdal.Open("tmp/cache/S2002196124536.L1A_HDUN.BartonBendish.extract.hdf")
    assert ds.RasterCount == 8, "did not get expected band number"

    ds = None


###############################################################################
# Test fix for #2208


def test_hdf4_read_online_5():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    # 13 MB
    gdaltest.download_or_skip(
        "https://www.ncei.noaa.gov/data/oceans/pathfinder/Version5.0/Monthly/1991/199101.s04m1pfv50-sst-16b.hdf",
        "199101.s04m1pfv50-sst-16b.hdf",
    )

    tst = gdaltest.GDALTest(
        "HDF4Image",
        "tmp/cache/199101.s04m1pfv50-sst-16b.hdf",
        1,
        41173,
        filename_absolute=1,
    )

    tst.testOpen()


###############################################################################
# Test fix for #3386 where block size is dataset size


def test_hdf4_read_online_6():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    # 1 MB
    gdaltest.download_or_skip(
        "http://download.osgeo.org/gdal/data/hdf4/MOD09Q1G_EVI.A2006233.h07v03.005.2008338190308.hdf",
        "MOD09Q1G_EVI.A2006233.h07v03.005.2008338190308.hdf",
    )

    # Test with quoting of components
    tst = gdaltest.GDALTest(
        "HDF4Image",
        'HDF4_EOS:EOS_GRID:"tmp/cache/MOD09Q1G_EVI.A2006233.h07v03.005.2008338190308.hdf":"MODIS_NACP_EVI":"MODIS_EVI"',
        1,
        12197,
        filename_absolute=1,
    )

    tst.testOpen()

    ds = gdal.Open(
        "HDF4_EOS:EOS_GRID:tmp/cache/MOD09Q1G_EVI.A2006233.h07v03.005.2008338190308.hdf:MODIS_NACP_EVI:MODIS_EVI"
    )

    if "GetBlockSize" in dir(gdal.Band):
        (blockx, blocky) = ds.GetRasterBand(1).GetBlockSize()
        assert blockx == 4800 and blocky == 4800, "Did not get expected block size"

    cs = ds.GetRasterBand(1).Checksum()
    assert cs == 12197, "did not get expected checksum"

    ds = None


###############################################################################
# Test fix for #3386 where block size is smaller than dataset size


def test_hdf4_read_online_7():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    # 4 MB
    gdaltest.download_or_skip(
        "http://download.osgeo.org/gdal/data/hdf4/MOD09A1.A2010041.h06v03.005.2010051001103.hdf",
        "MOD09A1.A2010041.h06v03.005.2010051001103.hdf",
    )

    tst = gdaltest.GDALTest(
        "HDF4Image",
        "HDF4_EOS:EOS_GRID:tmp/cache/MOD09A1.A2010041.h06v03.005.2010051001103.hdf:MOD_Grid_500m_Surface_Reflectance:sur_refl_b01",
        1,
        54894,
        filename_absolute=1,
    )

    tst.testOpen()

    ds = gdal.Open(
        "HDF4_EOS:EOS_GRID:tmp/cache/MOD09A1.A2010041.h06v03.005.2010051001103.hdf:MOD_Grid_500m_Surface_Reflectance:sur_refl_b01"
    )

    if "GetBlockSize" in dir(gdal.Band):
        (blockx, blocky) = ds.GetRasterBand(1).GetBlockSize()
        assert blockx == 2400 and blocky == 32, "Did not get expected block size"

    cs = ds.GetRasterBand(1).Checksum()
    assert cs == 54894, "did not get expected checksum"

    ds = None


###############################################################################
# Test reading a HDF4_EOS:EOS_GRID where preferred block height reported would be 1
# but that will lead to very poor performance (#3386)


def test_hdf4_read_online_8():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    # 5 MB
    gdaltest.download_or_skip(
        "https://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.006/2006.06.10/MOD13Q1.A2006161.h34v09.006.2015161173716.hdf",
        "MOD13Q1.A2006161.h34v09.006.2015161173716.hdf",
    )

    tst = gdaltest.GDALTest(
        "HDF4Image",
        "HDF4_EOS:EOS_GRID:tmp/cache/MOD13Q1.A2006161.h34v09.006.2015161173716.hdf:MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI",
        1,
        44174,
        filename_absolute=1,
    )

    tst.testOpen()

    ds = gdal.Open(
        "HDF4_EOS:EOS_GRID:tmp/cache/MOD13Q1.A2006161.h34v09.006.2015161173716.hdf:MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI"
    )

    cs = ds.GetRasterBand(1).Checksum()
    assert cs == 44174, "did not get expected checksum"

    if "GetBlockSize" in dir(gdal.Band):
        (blockx, blocky) = ds.GetRasterBand(1).GetBlockSize()
        if blockx != 4800 or blocky == 1:
            print("blockx=%d, blocky=%d" % (blockx, blocky))
            pytest.fail("Did not get expected block size")

    ds = None


###############################################################################
# Test reading L1G MTL metadata metadata


def test_hdf4_read_online_9():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    gdaltest.download_or_skip(
        "ftp://ftp.maps.canada.ca/pub/nrcan_rncan/archive/image/landsat_7/geobase_hdf/L71002025_02520010722/L71002025_02520010722_MTL.L1G",
        "L71002025_02520010722_MTL.L1G",
    )

    gdaltest.download_or_skip(
        "ftp://ftp.maps.canada.ca/pub/nrcan_rncan/archive/image/landsat_7/geobase_hdf/L71002025_02520010722/L71002025_02520010722_HDF.L1G",
        "L71002025_02520010722_HDF.L1G",
    )

    f = open("tmp/cache/L71002025_02520010722_B10.L1G", "wb")
    f.close()

    ds = gdal.Open('HDF4_SDS:UNKNOWN:"tmp/cache/L71002025_02520010722_HDF.L1G":0')
    gcp_count = ds.GetGCPCount()
    ds = None

    assert gcp_count == 4, "did not get expected gcp count"


###############################################################################
# Test that non-tiled access works (#4672)


def test_hdf4_read_online_10():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    gdaltest.download_or_skip(
        "http://trac.osgeo.org/gdal/raw-attachment/ticket/4672/MOD16A2.A2000M01.h14v02.105.2010357183410.hdf",
        "MOD16A2.A2000M01.h14v02.105.2010357183410.hdf",
    )

    ds = gdal.Open(
        'HDF4_EOS:EOS_GRID:"tmp/cache/MOD16A2.A2000M01.h14v02.105.2010357183410.hdf":MOD_Grid_MOD16A2:ET_1km'
    )

    if "GetBlockSize" in dir(gdal.Band):
        (blockx, blocky) = ds.GetRasterBand(1).GetBlockSize()
        assert blockx == 1200 and blocky == 833, "Did not get expected block size"

    cs = ds.GetRasterBand(1).Checksum()
    assert cs == 20976, "did not get expected checksum"

    ds = None


###############################################################################
# Test reading HDFEOS SWATH projducts


def test_hdf4_read_online_11():

    if gdaltest.hdf4_drv is None:
        pytest.skip()

    gdaltest.download_or_skip(
        "https://gamma.hdfgroup.org/ftp/pub/outgoing/NASAHDFfiles2/eosweb/hdf4/hdfeos2-swath-wo-dimmaps/AMSR_E_L2_Ocean_B01_200206182340_A.hdf",
        "AMSR_E_L2_Ocean_B01_200206182340_A.hdf",
    )

    ds = gdal.Open(
        'HDF4_EOS:EOS_SWATH:"tmp/cache/AMSR_E_L2_Ocean_B01_200206182340_A.hdf":Swath1:Ocean_products_quality_flag'
    )

    cs = ds.GetRasterBand(1).Checksum()
    assert cs == 7809, "did not get expected checksum"

    ds = gdal.Open("tmp/cache/AMSR_E_L2_Ocean_B01_200206182340_A.hdf")
    assert len(ds.GetSubDatasets()) == 7

    ds = gdal.OpenEx(
        "tmp/cache/AMSR_E_L2_Ocean_B01_200206182340_A.hdf",
        open_options=["LIST_SDS=YES"],
    )
    assert len(ds.GetSubDatasets()) == 16