From 0054c80e89edb275fedd388acbf63058702d10ed Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 18 Sep 2024 13:04:39 +0200 Subject: [PATCH] Add tests for Landsat L1 reader. --- .../reader_tests/test_oli_tirs_l1_tif.py | 576 ++++++++++++++++++ 1 file changed, 576 insertions(+) create mode 100644 satpy/tests/reader_tests/test_oli_tirs_l1_tif.py diff --git a/satpy/tests/reader_tests/test_oli_tirs_l1_tif.py b/satpy/tests/reader_tests/test_oli_tirs_l1_tif.py new file mode 100644 index 0000000000..a2e1cd94d5 --- /dev/null +++ b/satpy/tests/reader_tests/test_oli_tirs_l1_tif.py @@ -0,0 +1,576 @@ +#!/usr/bin/python +# Copyright (c) 2018 Satpy developers +# +# This file is part of satpy. +# +# satpy is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# satpy 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 General Public License for more details. +# +# You should have received a copy of the GNU General Public License along with +# satpy. If not, see . +"""Unittests for generic image reader.""" + +import os +import unittest + +import dask.array as da +import numpy as np +import pytest +import xarray as xr + +metadata_text = b""" + + + Image courtesy of the U.S. Geological Survey + https://doi.org/10.5066/P975CC9B + LC08_L1GT_026200_20240502_20240513_02_T2 + L1GT + 02 + T2 + GEOTIFF + LC08_L1GT_026200_20240502_20240513_02_T2_B1.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B2.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B3.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B4.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B5.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B6.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B7.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B8.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B9.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B10.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B11.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_QA_PIXEL.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_QA_RADSAT.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_ANG.txt + LC08_L1GT_026200_20240502_20240513_02_T2_VAA.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_VZA.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_SAA.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_SZA.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_MTL.txt + LC08_L1GT_026200_20240502_20240513_02_T2_MTL.xml + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + UINT16 + INT16 + INT16 + INT16 + INT16 + + + LANDSAT_8 + OLI_TIRS + 2 + 26 + 200 + NADIR + 26 + 200 + 2024-05-02 + 18:00:24.6148649Z + LGN + 0.85 + -1 + 9 + 9 + N + N + N + Y + N + N + N + N + N + -0.000 + -39.71362413 + -41.46228969 + 1.0079981 + UPPER + FINAL + ESTIMATED + + + UTM + WGS84 + WGS84 + 40 + 15.00 + 30.00 + 30.00 + 200 + 200 + 100 + 100 + 100 + 100 + NORTH_UP + 24.18941 + 58.17657 + 24.15493 + 60.44878 + 22.06522 + 58.15819 + 22.03410 + 60.39501 + 619500.000 + 2675700.000 + 850500.000 + 2675700.000 + 619500.000 + 2440500.000 + 850500.000 + 2440500.000 + + + Image courtesy of the U.S. Geological Survey + https://doi.org/10.5066/P975CC9B + 1885324_00001 + LC80262002024123LGN00 + LC08_L1GT_026200_20240502_20240513_02_T2 + L1GT + T2 + GEOTIFF + 2024-05-13T15:32:54Z + LPGS_16.4.0 + LC08_L1GT_026200_20240502_20240513_02_T2_B1.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B2.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B3.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B4.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B5.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B6.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B7.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B8.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B9.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B10.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_B11.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_QA_PIXEL.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_QA_RADSAT.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_ANG.txt + LC08_L1GT_026200_20240502_20240513_02_T2_VAA.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_VZA.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_SAA.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_SZA.TIF + LC08_L1GT_026200_20240502_20240513_02_T2_MTL.txt + LC08_L1GT_026200_20240502_20240513_02_T2_MTL.xml + LC08CPF_20240429_20240630_02.03 + LO8BPF20240502162846_20240502181430.01 + LT8BPF20240502144307_20240510102926.01 + LC08RLUT_20150303_20431231_02_01.h5 + TIRS + GLS2000 + + + 748.04883 + -61.77412 + 766.01111 + -63.25745 + 705.87274 + -58.29120 + 595.23163 + -49.15442 + 364.25208 + -30.08006 + 90.58618 + -7.48064 + 30.53239 + -2.52137 + 673.63843 + -55.62928 + 142.35797 + -11.75597 + 22.00180 + 0.10033 + 22.00180 + 0.10033 + + + 1.210700 + -0.099980 + 1.210700 + -0.099980 + 1.210700 + -0.099980 + 1.210700 + -0.099980 + 1.210700 + -0.099980 + 1.210700 + -0.099980 + 1.210700 + -0.099980 + 1.210700 + -0.099980 + 1.210700 + -0.099980 + + + 65535 + 1 + 65535 + 1 + 65535 + 1 + 65535 + 1 + 65535 + 1 + 65535 + 1 + 65535 + 1 + 65535 + 1 + 65535 + 1 + 65535 + 1 + 65535 + 1 + + + 1.2357E-02 + 1.2654E-02 + 1.1661E-02 + 9.8329E-03 + 6.0172E-03 + 1.4964E-03 + 5.0438E-04 + 1.1128E-02 + 2.3517E-03 + 3.3420E-04 + 3.3420E-04 + -61.78647 + -63.27010 + -58.30286 + -49.16426 + -30.08607 + -7.48213 + -2.52188 + -55.64041 + -11.75832 + 0.10000 + 0.10000 + 2.0000E-05 + 2.0000E-05 + 2.0000E-05 + 2.0000E-05 + 2.0000E-05 + 2.0000E-05 + 2.0000E-05 + 2.0000E-05 + 2.0000E-05 + -0.100000 + -0.100000 + -0.100000 + -0.100000 + -0.100000 + -0.100000 + -0.100000 + -0.100000 + -0.100000 + + + 774.8853 + 1321.0789 + 480.8883 + 1201.1442 + + + UTM + WGS84 + WGS84 + 40 + 15.00 + 30.00 + 30.00 + NORTH_UP + CUBIC_CONVOLUTION + + +""" + +class TestOLITIRSL1(unittest.TestCase): + """Test generic image reader.""" + + def setUp(self): + """Create temporary images and metadata to test on.""" + import tempfile + from datetime import datetime + + from pyresample.geometry import AreaDefinition + + from satpy.scene import Scene + + self.date = datetime(2024, 5, 12) + + self.filename_info = dict(observation_date=datetime(2024, 5, 3), + platform_type="L", + process_level_correction="L1TP", + spacecraft_id="08", + data_type="C") + self.ftype_info = {"file_type": "granule_b04"} + + # Create area definition + pcs_id = "WGS 84 / UTM zone 40N" + proj4_dict = {"proj": "utm", "zone": 40, "datum": "WGS84", "units": "m", "no_defs": None, "type": "crs"} + self.x_size = 100 + self.y_size = 100 + area_extent = (619485., 2440485., 850515., 2675715.) + self.area_def = AreaDefinition("geotiff_area", pcs_id, pcs_id, + proj4_dict, self.x_size, self.y_size, + area_extent) + + # Create datasets for L, LA, RGB and RGBA mode images + self.test_data__1 = da.random.randint(12000, 16000, + size=(self.y_size, self.x_size), + chunks=(50, 50)).astype(np.uint16) + self.test_data__2 = da.random.randint(8000, 14000, + size=(self.y_size, self.x_size), + chunks=(50, 50)).astype(np.uint16) + self.test_data__3= da.random.randint(0, 10000, + size=(self.y_size, self.x_size), + chunks=(50, 50)).astype(np.uint16) + + ds_b4 = xr.DataArray(self.test_data__1, + dims=("y", "x"), + attrs={"name": "b04", + "start_time": self.date}) + + ds_b11 = xr.DataArray(self.test_data__2, + dims=("y", "x"), + attrs={"name": "b04", + "start_time": self.date}) + + ds_sza = xr.DataArray(self.test_data__3, + dims=("y", "x"), + attrs={"name": "sza", + "start_time": self.date}) + + # Temp dir for the saved images + self.base_dir = tempfile.mkdtemp() + + # Filenames to be used during testing + self.fnames = [f"{self.base_dir}/LC08_L1GT_026200_20240502_20240513_02_T2_B4.TIF", + f"{self.base_dir}/LC08_L1GT_026200_20240502_20240513_02_T2_B11.TIF", + f"{self.base_dir}/LC08_L1GT_026200_20240502_20240513_02_T2_MTL.xml", + f"{self.base_dir}/LC08_L1GT_026200_20240502_20240513_02_T2_SZA.TIF"] + + self.bad_fname_plat = self.fnames[0].replace("LC08", "BC08") + self.bad_fname_plat2 = self.fnames[2].replace("LC08", "BC08") + + self.bad_fname_chan = self.fnames[0].replace("B4", "B5") + + # Put the datasets to Scene for easy saving + scn = Scene() + scn["b4"] = ds_b4 + scn["b4"].attrs["area"] = self.area_def + scn["b11"] = ds_b11 + scn["b11"].attrs["area"] = self.area_def + scn["sza"] = ds_sza + scn["sza"].attrs["area"] = self.area_def + + # Save the images. Two images in PNG and two in GeoTIFF + scn.save_dataset("b4", writer="geotiff", enhance=False, fill_value=0, + filename=os.path.join(self.base_dir, self.fnames[0])) + scn.save_dataset("b11", writer="geotiff", enhance=False, fill_value=0, + filename=os.path.join(self.base_dir, self.fnames[1])) + scn.save_dataset("sza", writer="geotiff", enhance=False, fill_value=0, + filename=os.path.join(self.base_dir, self.fnames[3])) + + scn.save_dataset("b4", writer="geotiff", enhance=False, fill_value=0, + filename=self.bad_fname_plat) + scn.save_dataset("b4", writer="geotiff", enhance=False, fill_value=0, + filename=self.bad_fname_chan) + + # Write the metadata to a file + with open(os.path.join(self.base_dir, self.fnames[2]), "wb") as f: + f.write(metadata_text) + with open(self.bad_fname_plat2, "wb") as f: + f.write(metadata_text) + + self.scn = scn + + def tearDown(self): + """Remove the temporary directory created for a test.""" + try: + import shutil + shutil.rmtree(self.base_dir, ignore_errors=True) + except OSError: + pass + + def test_basicload(self): + """Test loading a Landsat Scene.""" + from satpy import Scene + scn = Scene(reader="oli_tirs_l1_tif", filenames=[self.fnames[0], + self.fnames[1], + self.fnames[2]]) + scn.load(["b04", "b11"]) + + # Check dataset is loaded correctly + assert scn["b04"].shape == (100, 100) + assert scn["b04"].attrs["area"] == self.area_def + assert scn["b04"].attrs["saturated"] + assert scn["b11"].shape == (100, 100) + assert scn["b11"].attrs["area"] == self.area_def + with pytest.raises(KeyError, match="saturated"): + assert not scn["b11"].attrs["saturated"] + + def test_ch_startend(self): + """Test correct retrieval of start/end times.""" + from datetime import datetime + + from satpy import Scene + scn = Scene(reader="oli_tirs_l1_tif", filenames=[self.fnames[0], self.fnames[3], self.fnames[2]]) + bnds = scn.available_dataset_names() + assert bnds == ["b04", "sza"] + + scn.load(["b04"]) + assert scn.start_time == datetime(2024, 5, 2, 18, 0, 24) + assert scn.end_time == datetime(2024, 5, 2, 18, 0, 24) + + def test_loading(self): + """Test loading a Landsat Scene with good and bad channel requests.""" + from satpy.readers.oli_tirs_l1_tif import OLITIRSCHReader, OLITIRSMDReader + good_mda = OLITIRSMDReader(self.fnames[2], self.filename_info, {}) + rdr = OLITIRSCHReader(self.fnames[0], self.filename_info, self.ftype_info, good_mda) + + # Check case with good file data and load request + rdr.get_dataset({"name": "b04", "calibration": "counts"}, {}) + + # Check case with request to load channel not matching filename + with pytest.raises(ValueError, match="Requested channel b05 does not match the reader channel b04"): + rdr.get_dataset({"name": "b05", "calibration": "counts"}, {}) + + bad_finfo = self.filename_info.copy() + bad_finfo["data_type"] = "T" + + # Check loading invalid channel for data type + rdr = OLITIRSCHReader(self.fnames[1], bad_finfo, self.ftype_info, good_mda) + with pytest.raises(ValueError, match= "Requested channel b04 is not available in this granule"): + rdr.get_dataset({"name": "b04", "calibration": "counts"}, {}) + + bad_finfo["data_type"] = "O" + ftype_b11 = self.ftype_info.copy() + ftype_b11["file_type"] = "granule_b11" + rdr = OLITIRSCHReader(self.fnames[1], bad_finfo, ftype_b11, good_mda) + with pytest.raises(ValueError, match="Requested channel b11 is not available in this granule"): + rdr.get_dataset({"name": "b11", "calibration": "counts"}, {}) + + def test_badfiles(self): + """Test loading a Landsat Scene with bad data.""" + from satpy.readers.oli_tirs_l1_tif import OLITIRSCHReader, OLITIRSMDReader + bad_fname_info = self.filename_info.copy() + bad_fname_info["platform_type"] = "B" + + # Test that metadata reader initialises with correct filename + good_mda = OLITIRSMDReader(self.fnames[2], self.filename_info, {}) + + # Check metadata reader fails if platform type is wrong + with pytest.raises(ValueError, match="This reader only supports Landsat data"): + OLITIRSMDReader(self.fnames[2], bad_fname_info, {}) + + # Test that metadata reader initialises with correct filename + OLITIRSCHReader(self.fnames[0], self.filename_info, self.ftype_info, good_mda) + + # Check metadata reader fails if platform type is wrong + with pytest.raises(ValueError, match="This reader only supports Landsat data"): + OLITIRSCHReader(self.fnames[0], bad_fname_info, self.ftype_info, good_mda) + bad_ftype_info = self.ftype_info.copy() + bad_ftype_info["file_type"] = "granule-b05" + with pytest.raises(ValueError, match="Invalid file type: granule-b05"): + OLITIRSCHReader(self.fnames[0], self.filename_info, bad_ftype_info, good_mda) + + def test_calibration_modes(self): + """Test calibration modes for the reader.""" + from satpy import Scene + + # Check counts calibration + scn = Scene(reader="oli_tirs_l1_tif", filenames=self.fnames) + scn.load(["b04", "b11"], calibration="counts") + np.testing.assert_allclose(scn["b04"].values, self.test_data__1) + np.testing.assert_allclose(scn["b11"].values, self.test_data__2) + assert scn["b04"].attrs["units"] == "1" + assert scn["b11"].attrs["units"] == "1" + assert scn["b04"].attrs["standard_name"] == "counts" + assert scn["b11"].attrs["standard_name"] == "counts" + + # Check radiance calibration + exp_b04 = (self.test_data__1 * 0.0098329 - 49.16426).astype(np.float32) + exp_b11 = (self.test_data__2 * 0.0003342 + 0.100000).astype(np.float32) + + scn = Scene(reader="oli_tirs_l1_tif", filenames=self.fnames) + scn.load(["b04", "b11"], calibration="radiance") + assert scn["b04"].attrs["units"] == "W m-2 um-1 sr-1" + assert scn["b11"].attrs["units"] == "W m-2 um-1 sr-1" + assert scn["b04"].attrs["standard_name"] == "toa_outgoing_radiance_per_unit_wavelength" + assert scn["b11"].attrs["standard_name"] == "toa_outgoing_radiance_per_unit_wavelength" + np.testing.assert_allclose(scn["b04"].values, exp_b04, rtol=1e-4) + np.testing.assert_allclose(scn["b11"].values, exp_b11, rtol=1e-4) + + # Check top level calibration + exp_b04 = (self.test_data__1 * 2e-05 - 0.1).astype(np.float32) + exp_b11 = (self.test_data__2 * 0.0003342 + 0.100000) + exp_b11 = (1201.1442 / np.log((480.8883 / exp_b11) + 1)).astype(np.float32) + scn = Scene(reader="oli_tirs_l1_tif", filenames=self.fnames) + scn.load(["b04", "b11"]) + + assert scn["b04"].attrs["units"] == "%" + assert scn["b11"].attrs["units"] == "K" + assert scn["b04"].attrs["standard_name"] == "toa_bidirectional_reflectance" + assert scn["b11"].attrs["standard_name"] == "toa_brightness_temperature" + np.testing.assert_allclose(np.array(scn["b04"].values), np.array(exp_b04), rtol=1e-4) + np.testing.assert_allclose(scn["b11"].values, exp_b11, rtol=1e-6) + + # Check angles are calculated correctly + scn = Scene(reader="oli_tirs_l1_tif", filenames=self.fnames) + scn.load(["sza"]) + assert scn["sza"].attrs["units"] == "degrees" + assert scn["sza"].attrs["standard_name"] == "solar_zenith_angle" + np.testing.assert_allclose(scn["sza"].values * 100, np.array(self.test_data__3), atol=0.01, rtol=1e-3) + + def test_metadata(self): + """Check that metadata values loaded correctly.""" + from satpy.readers.oli_tirs_l1_tif import OLITIRSMDReader + mda = OLITIRSMDReader(self.fnames[2], self.filename_info, {}) + + cal_test_dict = {"b01": (0.012357, -61.78647, 2e-05, -0.1), + "b05": (0.0060172, -30.08607, 2e-05, -0.1), + "b10": (0.0003342, 0.1, 774.8853, 1321.0789)} + + assert mda.platform_name == "Landsat-8" + assert mda.earth_sun_distance() == 1.0079981 + assert mda.band_calibration["b01"] == cal_test_dict["b01"] + assert mda.band_calibration["b05"] == cal_test_dict["b05"] + assert mda.band_calibration["b10"] == cal_test_dict["b10"] + assert not mda.band_saturation["b01"] + assert mda.band_saturation["b04"] + assert not mda.band_saturation["b05"] + with pytest.raises(KeyError): + mda.band_saturation["b10"] + + def test_area_def(self): + """Check we can get the area defs properly.""" + from satpy.readers.oli_tirs_l1_tif import OLITIRSMDReader + mda = OLITIRSMDReader(self.fnames[2], self.filename_info, {}) + + standard_area = mda.build_area_def("b01") + pan_area = mda.build_area_def("b08") + + assert standard_area.area_extent == (619485.0, 2440485.0, 850515.0, 2675715.0) + assert pan_area.area_extent == (619492.5, 2440492.5, 850507.5, 2675707.5)