From f182d7424933d1e437a067a2b4419700a8ec17c6 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 22 Oct 2023 19:29:19 -0500 Subject: [PATCH 01/16] Add initial hacky chunking and float32 handling to ABI L1b reader --- satpy/readers/abi_base.py | 18 +++++++++++++----- satpy/readers/abi_l1b.py | 4 ++-- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/satpy/readers/abi_base.py b/satpy/readers/abi_base.py index 0b80045767..69c059e569 100644 --- a/satpy/readers/abi_base.py +++ b/satpy/readers/abi_base.py @@ -62,17 +62,25 @@ def __init__(self, filename, filename_info, filetype_info): @cached_property def nc(self): """Get the xarray dataset for this file.""" + import math + + from satpy.utils import get_dask_chunk_size_in_bytes + chunk_size_for_high_res = math.sqrt(get_dask_chunk_size_in_bytes() / 4) # 32-bit floats + chunk_size_for_high_res = np.round(chunk_size_for_high_res / 226) * 226 + ft = self.filetype_info["file_type"] + low_res_factor = 1 if ft == "c02" else (2 if ft in ("c01", "c03", "c05") else 4) + chunk_size = int(chunk_size_for_high_res / low_res_factor) f_obj = open_file_or_filename(self.filename) try: nc = xr.open_dataset(f_obj, decode_cf=True, mask_and_scale=False, - chunks={"x": CHUNK_SIZE, "y": CHUNK_SIZE}, ) + chunks={'x': chunk_size, 'y': chunk_size}, ) except ValueError: nc = xr.open_dataset(f_obj, decode_cf=True, mask_and_scale=False, - chunks={"lon": CHUNK_SIZE, "lat": CHUNK_SIZE}, ) + chunks={'lon': chunk_size, 'lat': chunk_size}, ) nc = self._rename_dims(nc) return nc @@ -137,7 +145,7 @@ def is_int(val): new_fill = fill else: new_fill = np.nan - data = data.where(data != fill, new_fill) + data = data.where(data != fill, np.float32(new_fill)) if factor != 1 and item in ("x", "y"): # be more precise with x/y coordinates # see get_area_def for more information @@ -147,8 +155,8 @@ def is_int(val): # can't do this in place since data is most likely uint16 # and we are making it a 64-bit float if not is_int(factor): - factor = float(factor) - data = data * factor + offset + factor = np.float32(factor) + data = data * np.float32(factor) + np.float32(offset) return data def _adjust_coords(self, data, item): diff --git a/satpy/readers/abi_l1b.py b/satpy/readers/abi_l1b.py index 3a22397cde..4d0276bf79 100644 --- a/satpy/readers/abi_l1b.py +++ b/satpy/readers/abi_l1b.py @@ -136,11 +136,11 @@ def _raw_calibrate(self, data): def _vis_calibrate(self, data): """Calibrate visible channels to reflectance.""" solar_irradiance = self["esun"] - esd = self["earth_sun_distance_anomaly_in_AU"].astype(float) + esd = self["earth_sun_distance_anomaly_in_AU"].astype(np.float32) factor = np.pi * esd * esd / solar_irradiance - res = data * factor + res = data * np.float32(factor) res.attrs = data.attrs res.attrs["units"] = "1" res.attrs["long_name"] = "Bidirectional Reflectance" From 878e5c6c4dbc4d6200df208660102d060786907d Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 27 Oct 2023 12:34:39 -0500 Subject: [PATCH 02/16] Use filetype info for ABI resolution-based chunking --- satpy/etc/readers/abi_l1b.yaml | 4 ++++ satpy/readers/abi_base.py | 19 +++++++------------ satpy/tests/reader_tests/test_abi_l1b.py | 2 +- 3 files changed, 12 insertions(+), 13 deletions(-) diff --git a/satpy/etc/readers/abi_l1b.yaml b/satpy/etc/readers/abi_l1b.yaml index d9de341ff1..f4986ba106 100644 --- a/satpy/etc/readers/abi_l1b.yaml +++ b/satpy/etc/readers/abi_l1b.yaml @@ -25,16 +25,19 @@ file_types: # "suffix" is an arbitrary suffix that may be added during third-party testing (see PR #1380) c01: file_reader: !!python/name:satpy.readers.abi_l1b.NC_ABI_L1B + resolution: 1000 file_patterns: ['{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C01_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}.nc{nc_version}', '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C01_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}-{chid:6d}_0.nc{nc_version}', '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C01_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}_{suffix}.nc{nc_version}'] c02: file_reader: !!python/name:satpy.readers.abi_l1b.NC_ABI_L1B + resolution: 500 file_patterns: ['{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C02_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}.nc{nc_version}', '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C02_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}-{chid:6d}_0.nc{nc_version}', '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C02_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}_{suffix}.nc{nc_version}'] c03: file_reader: !!python/name:satpy.readers.abi_l1b.NC_ABI_L1B + resolution: 1000 file_patterns: ['{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C03_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}.nc{nc_version}', '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C03_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}-{chid:6d}_0.nc{nc_version}', '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C03_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}_{suffix}.nc{nc_version}'] @@ -44,6 +47,7 @@ file_types: '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C04_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}_{suffix}.nc{nc_version}'] c05: file_reader: !!python/name:satpy.readers.abi_l1b.NC_ABI_L1B + resolution: 1000 file_patterns: ['{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C05_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}.nc{nc_version}', '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C05_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}-{chid:6d}_0.nc{nc_version}', '{system_environment:2s}_{mission_id:3s}-L1b-{observation_type:3s}{scene_abbr:s}-{scan_mode:2s}C05_{platform_shortname:3s}_s{start_time:%Y%j%H%M%S%f}_e{end_time:%Y%j%H%M%S%f}_c{creation_time:%Y%j%H%M%S%f}_{suffix}.nc{nc_version}'] diff --git a/satpy/readers/abi_base.py b/satpy/readers/abi_base.py index 69c059e569..956bec278e 100644 --- a/satpy/readers/abi_base.py +++ b/satpy/readers/abi_base.py @@ -66,21 +66,16 @@ def nc(self): from satpy.utils import get_dask_chunk_size_in_bytes chunk_size_for_high_res = math.sqrt(get_dask_chunk_size_in_bytes() / 4) # 32-bit floats - chunk_size_for_high_res = np.round(chunk_size_for_high_res / 226) * 226 - ft = self.filetype_info["file_type"] - low_res_factor = 1 if ft == "c02" else (2 if ft in ("c01", "c03", "c05") else 4) - chunk_size = int(chunk_size_for_high_res / low_res_factor) - f_obj = open_file_or_filename(self.filename) - try: + chunk_size_for_high_res = np.round(chunk_size_for_high_res / (4 * 226)) * (4 * 226) + low_res_factor = int(self.filetype_info.get("resolution", 2000) // 500) + res_chunk_bytes = int(chunk_size_for_high_res / low_res_factor) * 4 + import dask + with dask.config.set({"array.chunk-size": res_chunk_bytes}): + f_obj = open_file_or_filename(self.filename) nc = xr.open_dataset(f_obj, decode_cf=True, mask_and_scale=False, - chunks={'x': chunk_size, 'y': chunk_size}, ) - except ValueError: - nc = xr.open_dataset(f_obj, - decode_cf=True, - mask_and_scale=False, - chunks={'lon': chunk_size, 'lat': chunk_size}, ) + chunks="auto") nc = self._rename_dims(nc) return nc diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index ab2b1eec54..7563f6e13d 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -404,7 +404,7 @@ def test_open_dataset(self, _): # noqa: PT019 openable_thing = mock.MagicMock() - NC_ABI_L1B(openable_thing, {"platform_shortname": "g16"}, None) + NC_ABI_L1B(openable_thing, {"platform_shortname": "g16"}, {}) openable_thing.open.assert_called() From 85360970007207e6e0f5b611801491d56dfd70ba Mon Sep 17 00:00:00 2001 From: David Hoese Date: Fri, 27 Oct 2023 16:29:13 -0500 Subject: [PATCH 03/16] Start refactoring ABI L1b tests --- satpy/tests/reader_tests/test_abi_l1b.py | 77 +++++++++++++++--------- 1 file changed, 50 insertions(+), 27 deletions(-) diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index 7563f6e13d..f8bd7e4e9f 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -16,8 +16,10 @@ # You should have received a copy of the GNU General Public License along with # satpy. If not, see . """The abi_l1b reader tests package.""" +from __future__ import annotations import unittest +from typing import Any from unittest import mock import numpy as np @@ -26,13 +28,23 @@ from satpy.tests.utils import make_dataid +RAD_SHAPE = { + 500: (3000, 5000), # conus - 500m + 1000: (1500, 2500), # conus - 1km + 2000: (750, 1250), # conus - 2km +} -def _create_fake_rad_dataarray(rad=None): + +def _create_fake_rad_dataarray( + rad: xr.DataArray | None = None, + # resolution: int = 2000, +): x_image = xr.DataArray(0.) y_image = xr.DataArray(0.) time = xr.DataArray(0.) + shape = (2, 5) # RAD_SHAPE[resolution] if rad is None: - rad_data = (np.arange(10.).reshape((2, 5)) + 1.) * 50. + rad_data = (np.arange(shape[0] * shape[1]).reshape(shape) + 1.) * 50. rad_data = (rad_data + 1.) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( @@ -115,17 +127,28 @@ class Test_NC_ABI_L1B_Base(unittest.TestCase): """Common setup for NC_ABI_L1B tests.""" @mock.patch("satpy.readers.abi_base.xr") - def setUp(self, xr_, rad=None, clip_negative_radiances=False): + def setUp( + self, + xr_, + rad: xr.DataArray | None = None, + clip_negative_radiances: bool = False, + filetype_resolution: int = 0 + ) -> None: """Create a fake dataset using the given radiance data.""" from satpy.readers.abi_l1b import NC_ABI_L1B xr_.open_dataset.return_value = _create_fake_rad_dataset(rad=rad) - self.reader = NC_ABI_L1B("filename", - {"platform_shortname": "G16", "observation_type": "Rad", - "suffix": "custom", - "scene_abbr": "C", "scan_mode": "M3"}, - {"filetype": "info"}, - clip_negative_radiances=clip_negative_radiances) + ft_info: dict[str, Any] = {"filetype": "info"} + if filetype_resolution: + ft_info["resolution"] = filetype_resolution + self.file_handler = NC_ABI_L1B( + "filename", + {"platform_shortname": "G16", "observation_type": "Rad", + "suffix": "custom", + "scene_abbr": "C", "scan_mode": "M3"}, + ft_info, + clip_negative_radiances=clip_negative_radiances + ) class TestABIYAML: @@ -157,13 +180,13 @@ class Test_NC_ABI_L1B(Test_NC_ABI_L1B_Base): def test_basic_attributes(self): """Test getting basic file attributes.""" from datetime import datetime - assert self.reader.start_time == datetime(2017, 9, 20, 17, 30, 40, 800000) - assert self.reader.end_time == datetime(2017, 9, 20, 17, 41, 17, 500000) + assert self.file_handler.start_time == datetime(2017, 9, 20, 17, 30, 40, 800000) + assert self.file_handler.end_time == datetime(2017, 9, 20, 17, 41, 17, 500000) def test_get_dataset(self): """Test the get_dataset method.""" key = make_dataid(name="Rad", calibration="radiance") - res = self.reader.get_dataset(key, {"info": "info"}) + res = self.file_handler.get_dataset(key, {"info": "info"}) exp = {"calibration": "radiance", "instrument_ID": None, "modifiers": (), @@ -198,14 +221,14 @@ def test_get_dataset(self): @mock.patch("satpy.readers.abi_base.geometry.AreaDefinition") def test_get_area_def(self, adef): """Test the area generation.""" - self.reader.get_area_def(None) + self.file_handler.get_area_def(None) assert adef.call_count == 1 call_args = tuple(adef.call_args)[0] assert call_args[3] == {"a": 1.0, "b": 1.0, "h": 1.0, "lon_0": -90.0, "proj": "geos", "sweep": "x", "units": "m"} - assert call_args[4] == self.reader.ncols - assert call_args[5] == self.reader.nlines + assert call_args[4] == self.file_handler.ncols + assert call_args[5] == self.file_handler.nlines np.testing.assert_allclose(call_args[6], (-2, -2, 8, 2)) @@ -226,11 +249,11 @@ def setUp(self): "_FillValue": 1002, # last rad_data value } ) - super(Test_NC_ABI_L1B_ir_cal, self).setUp(rad=rad) + super(Test_NC_ABI_L1B_ir_cal, self).setUp(rad=rad, filetype_resolution=2000) def test_ir_calibration_attrs(self): """Test IR calibrated DataArray attributes.""" - res = self.reader.get_dataset( + res = self.file_handler.get_dataset( make_dataid(name="C05", calibration="brightness_temperature"), {}) # make sure the attributes from the file are in the data array @@ -241,11 +264,11 @@ def test_ir_calibration_attrs(self): def test_clip_negative_radiances_attribute(self): """Assert that clip_negative_radiances is set to False.""" - assert not self.reader.clip_negative_radiances + assert not self.file_handler.clip_negative_radiances def test_ir_calibrate(self): """Test IR calibration.""" - res = self.reader.get_dataset( + res = self.file_handler.get_dataset( make_dataid(name="C05", calibration="brightness_temperature"), {}) expected = np.array([[267.55572248, 305.15576503, 332.37383249, 354.73895301, 374.19710115], @@ -273,15 +296,15 @@ def setUp(self): } ) - super().setUp(rad=rad, clip_negative_radiances=True) + super().setUp(rad=rad, clip_negative_radiances=True, filetype_resolution=2000) def test_clip_negative_radiances_attribute(self): """Assert that clip_negative_radiances has been set to True.""" - assert self.reader.clip_negative_radiances + assert self.file_handler.clip_negative_radiances def test_ir_calibrate(self): """Test IR calibration.""" - res = self.reader.get_dataset( + res = self.file_handler.get_dataset( make_dataid(name="C07", calibration="brightness_temperature"), {}) clipped_ir = 267.07775531 @@ -319,11 +342,11 @@ def setUp(self): "_FillValue": 20, } ) - super(Test_NC_ABI_L1B_vis_cal, self).setUp(rad=rad) + super(Test_NC_ABI_L1B_vis_cal, self).setUp(rad=rad, filetype_resolution=1000) def test_vis_calibrate(self): """Test VIS calibration.""" - res = self.reader.get_dataset( + res = self.file_handler.get_dataset( make_dataid(name="C05", calibration="reflectance"), {}) expected = np.array([[0.15265617, 0.30531234, 0.45796851, 0.61062468, 0.76328085], @@ -352,11 +375,11 @@ def setUp(self): "_FillValue": 20, } ) - super(Test_NC_ABI_L1B_raw_cal, self).setUp(rad=rad) + super(Test_NC_ABI_L1B_raw_cal, self).setUp(rad=rad, filetype_resolution=1000) def test_raw_calibrate(self): """Test RAW calibration.""" - res = self.reader.get_dataset( + res = self.file_handler.get_dataset( make_dataid(name="C05", calibration="counts"), {}) # We expect the raw data to be unchanged @@ -391,7 +414,7 @@ def to_dict(self): with self.assertRaises(ValueError, msg="Did not detect invalid cal"): did = FakeDataID(name="C05", calibration="invalid", modifiers=()) - self.reader.get_dataset(did, {}) + self.file_handler.get_dataset(did, {}) class Test_NC_ABI_File(unittest.TestCase): From 55e52484c05b254737863abb0648e1483d63e1b5 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sat, 28 Oct 2023 19:30:06 -0500 Subject: [PATCH 04/16] Remove unnecessary duplication in ABI L1b tests --- satpy/readers/abi_l1b.py | 8 +- satpy/tests/reader_tests/test_abi_l1b.py | 451 ++++++++++++----------- 2 files changed, 231 insertions(+), 228 deletions(-) diff --git a/satpy/readers/abi_l1b.py b/satpy/readers/abi_l1b.py index 4d0276bf79..c3da53c9c7 100644 --- a/satpy/readers/abi_l1b.py +++ b/satpy/readers/abi_l1b.py @@ -59,12 +59,8 @@ def get_dataset(self, key, info): "radiance": self._rad_calibrate, "counts": self._raw_calibrate, } - - try: - func = cal_dictionary[key["calibration"]] - res = func(radiances) - except KeyError: - raise ValueError("Unknown calibration '{}'".format(key["calibration"])) + func = cal_dictionary[key["calibration"]] + res = func(radiances) # convert to satpy standard units if res.attrs["units"] == "1" and key["calibration"] != "counts": diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index f8bd7e4e9f..bdaa03f9e5 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -18,14 +18,17 @@ """The abi_l1b reader tests package.""" from __future__ import annotations -import unittest -from typing import Any +import contextlib +from pathlib import Path +from typing import Any, Iterator from unittest import mock +import dask.array as da import numpy as np import pytest import xarray as xr +from satpy.readers.abi_l1b import NC_ABI_L1B from satpy.tests.utils import make_dataid RAD_SHAPE = { @@ -36,27 +39,27 @@ def _create_fake_rad_dataarray( - rad: xr.DataArray | None = None, - # resolution: int = 2000, -): - x_image = xr.DataArray(0.) - y_image = xr.DataArray(0.) - time = xr.DataArray(0.) + rad: xr.DataArray | None = None, + # resolution: int = 2000, +) -> xr.DataArray: + x_image = xr.DataArray(0.0) + y_image = xr.DataArray(0.0) + time = xr.DataArray(0.0) shape = (2, 5) # RAD_SHAPE[resolution] if rad is None: - rad_data = (np.arange(shape[0] * shape[1]).reshape(shape) + 1.) * 50. - rad_data = (rad_data + 1.) / 0.5 + rad_data = (np.arange(shape[0] * shape[1]).reshape(shape) + 1.0) * 50.0 + rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( rad_data, dims=("y", "x"), attrs={ "scale_factor": 0.5, - "add_offset": -1., + "add_offset": -1.0, "_FillValue": 1002, "units": "W m-2 um-1 sr-1", "valid_range": (0, 4095), - } + }, ) rad.coords["t"] = time rad.coords["x_image"] = x_image @@ -68,25 +71,21 @@ def _create_fake_rad_dataset(rad=None): rad = _create_fake_rad_dataarray(rad=rad) x__ = xr.DataArray( - range(5), - attrs={"scale_factor": 2., "add_offset": -1.}, - dims=("x",) + range(5), attrs={"scale_factor": 2.0, "add_offset": -1.0}, dims=("x",) ) y__ = xr.DataArray( - range(2), - attrs={"scale_factor": -2., "add_offset": 1.}, - dims=("y",) + range(2), attrs={"scale_factor": -2.0, "add_offset": 1.0}, dims=("y",) ) proj = xr.DataArray( [], attrs={ - "semi_major_axis": 1., - "semi_minor_axis": 1., - "perspective_point_height": 1., - "longitude_of_projection_origin": -90., - "latitude_of_projection_origin": 0., - "sweep_angle_axis": u"x" - } + "semi_major_axis": 1.0, + "semi_minor_axis": 1.0, + "perspective_point_height": 1.0, + "longitude_of_projection_origin": -90.0, + "latitude_of_projection_origin": 0.0, + "sweep_angle_axis": "x", + }, ) fake_dataset = xr.Dataset( @@ -95,8 +94,8 @@ def _create_fake_rad_dataset(rad=None): "band_id": np.array(8), # 'x': x__, # 'y': y__, - "x_image": xr.DataArray(0.), - "y_image": xr.DataArray(0.), + "x_image": xr.DataArray(0.0), + "y_image": xr.DataArray(0.0), "goes_imager_projection": proj, "yaw_flip_flag": np.array([1]), "planck_fk1": np.array(13432.1), @@ -107,13 +106,12 @@ def _create_fake_rad_dataset(rad=None): "nominal_satellite_subpoint_lat": np.array(0.0), "nominal_satellite_subpoint_lon": np.array(-89.5), "nominal_satellite_height": np.array(35786.02), - "earth_sun_distance_anomaly_in_AU": np.array(0.99) + "earth_sun_distance_anomaly_in_AU": np.array(0.99), }, coords={ "t": rad.coords["t"], "x": x__, "y": y__, - }, attrs={ "time_coverage_start": "2017-09-20T17:30:40.8Z", @@ -123,93 +121,139 @@ def _create_fake_rad_dataset(rad=None): return fake_dataset -class Test_NC_ABI_L1B_Base(unittest.TestCase): - """Common setup for NC_ABI_L1B tests.""" +def generate_l1b_filename(chan_name: str) -> str: + return f"OR_ABI-L1b-RadC-M4{chan_name}_G16_s20161811540362_e20161811545170_c20161811545230.nc" + + +@pytest.fixture(scope="module") +def l1b_c01_file(tmp_path_factory) -> list[Path]: + filename = generate_l1b_filename("C01") + data_path = tmp_path_factory.mktemp("abi_l1b").join(filename) + dataset = _create_fake_rad_dataset() + dataset.to_netcdf(data_path) + return [data_path] + + +@pytest.fixture(scope="module") +def l1b_all_files( + l1b_c01_file, +) -> list[Path]: + return l1b_c01_file - @mock.patch("satpy.readers.abi_base.xr") - def setUp( - self, - xr_, - rad: xr.DataArray | None = None, - clip_negative_radiances: bool = False, - filetype_resolution: int = 0 - ) -> None: - """Create a fake dataset using the given radiance data.""" - from satpy.readers.abi_l1b import NC_ABI_L1B +@contextlib.contextmanager +def create_file_handler( + rad: xr.DataArray | None = None, + clip_negative_radiances: bool = False, + filetype_resolution: int = 0, +) -> Iterator[NC_ABI_L1B]: + """Create a fake dataset using the given radiance data.""" + + ft_info: dict[str, Any] = {"filetype": "info"} + if filetype_resolution: + ft_info["resolution"] = filetype_resolution + + with mock.patch("satpy.readers.abi_base.xr") as xr_: xr_.open_dataset.return_value = _create_fake_rad_dataset(rad=rad) - ft_info: dict[str, Any] = {"filetype": "info"} - if filetype_resolution: - ft_info["resolution"] = filetype_resolution - self.file_handler = NC_ABI_L1B( + file_handler = NC_ABI_L1B( "filename", - {"platform_shortname": "G16", "observation_type": "Rad", - "suffix": "custom", - "scene_abbr": "C", "scan_mode": "M3"}, + { + "platform_shortname": "G16", + "observation_type": "Rad", + "suffix": "custom", + "scene_abbr": "C", + "scan_mode": "M3", + }, ft_info, - clip_negative_radiances=clip_negative_radiances + clip_negative_radiances=clip_negative_radiances, ) + yield file_handler class TestABIYAML: """Tests for the ABI L1b reader's YAML configuration.""" - @pytest.mark.parametrize(("channel", "suffix"), - [("C{:02d}".format(num), suffix) - for num in range(1, 17) - for suffix in ("", "_test_suffix")]) + @pytest.mark.parametrize( + ("channel", "suffix"), + [ + ("C{:02d}".format(num), suffix) + for num in range(1, 17) + for suffix in ("", "_test_suffix") + ], + ) def test_file_patterns_match(self, channel, suffix): """Test that the configured file patterns work.""" from satpy.readers import configs_for_reader, load_reader + reader_configs = list(configs_for_reader("abi_l1b"))[0] reader = load_reader(reader_configs) - fn1 = ("OR_ABI-L1b-RadM1-M3{}_G16_s20182541300210_e20182541300267" - "_c20182541300308{}.nc").format(channel, suffix) + fn1 = ( + "OR_ABI-L1b-RadM1-M3{}_G16_s20182541300210_e20182541300267" + "_c20182541300308{}.nc" + ).format(channel, suffix) loadables = reader.select_files_from_pathnames([fn1]) assert len(loadables) == 1 if not suffix and channel in ["C01", "C02", "C03", "C05"]: - fn2 = ("OR_ABI-L1b-RadM1-M3{}_G16_s20182541300210_e20182541300267" - "_c20182541300308-000000_0.nc").format(channel) + fn2 = ( + "OR_ABI-L1b-RadM1-M3{}_G16_s20182541300210_e20182541300267" + "_c20182541300308-000000_0.nc" + ).format(channel) loadables = reader.select_files_from_pathnames([fn2]) assert len(loadables) == 1 -class Test_NC_ABI_L1B(Test_NC_ABI_L1B_Base): +class Test_NC_ABI_L1B: """Test the NC_ABI_L1B reader.""" + @property + def fake_rad(self): + """Create fake data for these tests. + + Needs to be an instance method so the subclass can override it. + + """ + return None + def test_basic_attributes(self): """Test getting basic file attributes.""" from datetime import datetime - assert self.file_handler.start_time == datetime(2017, 9, 20, 17, 30, 40, 800000) - assert self.file_handler.end_time == datetime(2017, 9, 20, 17, 41, 17, 500000) + + with create_file_handler(rad=self.fake_rad) as file_handler: + assert file_handler.start_time == datetime(2017, 9, 20, 17, 30, 40, 800000) + assert file_handler.end_time == datetime(2017, 9, 20, 17, 41, 17, 500000) def test_get_dataset(self): """Test the get_dataset method.""" key = make_dataid(name="Rad", calibration="radiance") - res = self.file_handler.get_dataset(key, {"info": "info"}) - exp = {"calibration": "radiance", - "instrument_ID": None, - "modifiers": (), - "name": "Rad", - "observation_type": "Rad", - "orbital_parameters": {"projection_altitude": 1.0, - "projection_latitude": 0.0, - "projection_longitude": -90.0, - "satellite_nominal_altitude": 35786020., - "satellite_nominal_latitude": 0.0, - "satellite_nominal_longitude": -89.5, - "yaw_flip": True}, - "orbital_slot": None, - "platform_name": "GOES-16", - "platform_shortname": "G16", - "production_site": None, - "scan_mode": "M3", - "scene_abbr": "C", - "scene_id": None, - "sensor": "abi", - "timeline_ID": None, - "suffix": "custom", - "units": "W m-2 um-1 sr-1"} + with create_file_handler(rad=self.fake_rad) as file_handler: + res = file_handler.get_dataset(key, {"info": "info"}) + exp = { + "calibration": "radiance", + "instrument_ID": None, + "modifiers": (), + "name": "Rad", + "observation_type": "Rad", + "orbital_parameters": { + "projection_altitude": 1.0, + "projection_latitude": 0.0, + "projection_longitude": -90.0, + "satellite_nominal_altitude": 35786020.0, + "satellite_nominal_latitude": 0.0, + "satellite_nominal_longitude": -89.5, + "yaw_flip": True, + }, + "orbital_slot": None, + "platform_name": "GOES-16", + "platform_shortname": "G16", + "production_site": None, + "scan_mode": "M3", + "scene_abbr": "C", + "scene_id": None, + "sensor": "abi", + "timeline_ID": None, + "suffix": "custom", + "units": "W m-2 um-1 sr-1", + } assert res.attrs == exp # we remove any time dimension information @@ -221,40 +265,47 @@ def test_get_dataset(self): @mock.patch("satpy.readers.abi_base.geometry.AreaDefinition") def test_get_area_def(self, adef): """Test the area generation.""" - self.file_handler.get_area_def(None) - - assert adef.call_count == 1 - call_args = tuple(adef.call_args)[0] - assert call_args[3] == {"a": 1.0, "b": 1.0, "h": 1.0, - "lon_0": -90.0, "proj": "geos", "sweep": "x", "units": "m"} - assert call_args[4] == self.file_handler.ncols - assert call_args[5] == self.file_handler.nlines - np.testing.assert_allclose(call_args[6], (-2, -2, 8, 2)) + with create_file_handler(rad=self.fake_rad) as file_handler: + file_handler.get_area_def(None) + + assert adef.call_count == 1 + call_args = tuple(adef.call_args)[0] + assert call_args[3] == { + "a": 1.0, + "b": 1.0, + "h": 1.0, + "lon_0": -90.0, + "proj": "geos", + "sweep": "x", + "units": "m", + } + assert call_args[4] == file_handler.ncols + assert call_args[5] == file_handler.nlines + np.testing.assert_allclose(call_args[6], (-2, -2, 8, 2)) -class Test_NC_ABI_L1B_ir_cal(Test_NC_ABI_L1B_Base): +class Test_NC_ABI_L1B_ir_cal: """Test the NC_ABI_L1B reader's default IR calibration.""" - def setUp(self): - """Create fake data for the tests.""" - rad_data = (np.arange(10.).reshape((2, 5)) + 1.) * 50. - rad_data = (rad_data + 1.) / 0.5 - rad_data = rad_data.astype(np.int16) - rad = xr.DataArray( - rad_data, - dims=("y", "x"), - attrs={ - "scale_factor": 0.5, - "add_offset": -1., - "_FillValue": 1002, # last rad_data value - } + @pytest.mark.parametrize("clip_negative_radiances", [False, True]) + def test_ir_calibrate(self, clip_negative_radiances): + """Test IR calibration.""" + with _ir_file_handler( + clip_negative_radiances=clip_negative_radiances + ) as file_handler: + res = file_handler.get_dataset( + make_dataid(name="C07", calibration="brightness_temperature"), {} + ) + assert file_handler.clip_negative_radiances == clip_negative_radiances + + clipped_ir = 134.68753 if clip_negative_radiances else np.nan + expected = np.array( + [ + [clipped_ir, 304.97037, 332.22778, 354.6147, 374.08688], + [391.58655, 407.64786, 422.60635, 436.68802, np.nan], + ] ) - super(Test_NC_ABI_L1B_ir_cal, self).setUp(rad=rad, filetype_resolution=2000) - - def test_ir_calibration_attrs(self): - """Test IR calibrated DataArray attributes.""" - res = self.file_handler.get_dataset( - make_dataid(name="C05", calibration="brightness_temperature"), {}) + np.testing.assert_allclose(res.data, expected, equal_nan=True, atol=1e-04) # make sure the attributes from the file are in the data array assert "scale_factor" not in res.attrs @@ -262,95 +313,69 @@ def test_ir_calibration_attrs(self): assert res.attrs["standard_name"] == "toa_brightness_temperature" assert res.attrs["long_name"] == "Brightness Temperature" - def test_clip_negative_radiances_attribute(self): - """Assert that clip_negative_radiances is set to False.""" - assert not self.file_handler.clip_negative_radiances - - def test_ir_calibrate(self): - """Test IR calibration.""" - res = self.file_handler.get_dataset( - make_dataid(name="C05", calibration="brightness_temperature"), {}) - - expected = np.array([[267.55572248, 305.15576503, 332.37383249, 354.73895301, 374.19710115], - [391.68679226, 407.74064808, 422.69329105, 436.77021913, np.nan]]) - assert np.allclose(res.data, expected, equal_nan=True) - - -class Test_NC_ABI_L1B_clipped_ir_cal(Test_NC_ABI_L1B_Base): - """Test the NC_ABI_L1B reader's IR calibration (clipping negative radiance).""" - def setUp(self): - """Create fake data for the tests.""" - values = np.arange(10.) - values[0] = -0.0001 # introduce below minimum expected radiance - rad_data = (values.reshape((2, 5)) + 1.) * 50. - rad_data = (rad_data + 1.) / 0.5 - rad_data = rad_data.astype(np.int16) - rad = xr.DataArray( - rad_data, - dims=("y", "x"), - attrs={ - "scale_factor": 0.5, - "add_offset": -1., - "_FillValue": 1002, - } - ) - - super().setUp(rad=rad, clip_negative_radiances=True, filetype_resolution=2000) - - def test_clip_negative_radiances_attribute(self): - """Assert that clip_negative_radiances has been set to True.""" - assert self.file_handler.clip_negative_radiances - - def test_ir_calibrate(self): - """Test IR calibration.""" - res = self.file_handler.get_dataset( - make_dataid(name="C07", calibration="brightness_temperature"), {}) +@contextlib.contextmanager +def _ir_file_handler( + data: da.Array | None = None, clip_negative_radiances: bool = False +): + """Create fake data for the tests.""" + if data is None: + data = _fake_ir_data() + rad = xr.DataArray( + data, + dims=("y", "x"), + attrs={ + "scale_factor": 0.5, + "add_offset": -1.3, + "_FillValue": np.int16( + np.floor(((9 + 1) * 50.0 + 1.3) / 0.5) + ), # last rad_data value + }, + ) + with create_file_handler( + rad=rad, + clip_negative_radiances=clip_negative_radiances, + filetype_resolution=2000, + ) as file_handler: + yield file_handler - clipped_ir = 267.07775531 - expected = np.array([[clipped_ir, 305.15576503, 332.37383249, 354.73895301, 374.19710115], - [391.68679226, 407.74064808, 422.69329105, 436.77021913, np.nan]]) - assert np.allclose(res.data, expected, equal_nan=True) - def test_get_minimum_radiance(self): - """Test get_minimum_radiance from Rad DataArray.""" - from satpy.readers.abi_l1b import NC_ABI_L1B - data = xr.DataArray( - attrs={ - "scale_factor": 0.5, - "add_offset": -1., - "_FillValue": 1002, - } - ) - np.testing.assert_allclose(NC_ABI_L1B._get_minimum_radiance(NC_ABI_L1B, data), 0.0) +def _fake_ir_data(): + values = np.arange(10.0) + rad_data = (values.reshape((2, 5)) + 1.0) * 50.0 + rad_data[0, 0] = -0.0001 # introduce below minimum expected radiance + rad_data = (rad_data + 1.3) / 0.5 + return rad_data.astype(np.int16) -class Test_NC_ABI_L1B_vis_cal(Test_NC_ABI_L1B_Base): +class Test_NC_ABI_L1B_vis_cal: """Test the NC_ABI_L1B reader.""" - def setUp(self): - """Create fake data for the tests.""" - rad_data = (np.arange(10.).reshape((2, 5)) + 1.) - rad_data = (rad_data + 1.) / 0.5 + def test_vis_calibrate(self): + """Test VIS calibration.""" + rad_data = np.arange(10.0).reshape((2, 5)) + 1.0 + rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( rad_data, dims=("y", "x"), attrs={ "scale_factor": 0.5, - "add_offset": -1., + "add_offset": -1.0, "_FillValue": 20, - } + }, + ) + with create_file_handler(rad=rad, filetype_resolution=1000) as file_handler: + res = file_handler.get_dataset( + make_dataid(name="C05", calibration="reflectance"), {} + ) + + expected = np.array( + [ + [0.15265617, 0.30531234, 0.45796851, 0.61062468, 0.76328085], + [0.91593702, 1.06859319, 1.22124936, np.nan, 1.52656171], + ] ) - super(Test_NC_ABI_L1B_vis_cal, self).setUp(rad=rad, filetype_resolution=1000) - - def test_vis_calibrate(self): - """Test VIS calibration.""" - res = self.file_handler.get_dataset( - make_dataid(name="C05", calibration="reflectance"), {}) - - expected = np.array([[0.15265617, 0.30531234, 0.45796851, 0.61062468, 0.76328085], - [0.91593702, 1.06859319, 1.22124936, np.nan, 1.52656171]]) assert np.allclose(res.data, expected, equal_nan=True) assert "scale_factor" not in res.attrs assert "_FillValue" not in res.attrs @@ -358,29 +383,27 @@ def test_vis_calibrate(self): assert res.attrs["long_name"] == "Bidirectional Reflectance" -class Test_NC_ABI_L1B_raw_cal(Test_NC_ABI_L1B_Base): +class Test_NC_ABI_L1B_raw_cal: """Test the NC_ABI_L1B reader raw calibration.""" - def setUp(self): - """Create fake data for the tests.""" - rad_data = (np.arange(10.).reshape((2, 5)) + 1.) - rad_data = (rad_data + 1.) / 0.5 + def test_raw_calibrate(self): + """Test RAW calibration.""" + rad_data = np.arange(10.0).reshape((2, 5)) + 1.0 + rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( rad_data, dims=("y", "x"), attrs={ "scale_factor": 0.5, - "add_offset": -1., + "add_offset": -1.0, "_FillValue": 20, - } + }, ) - super(Test_NC_ABI_L1B_raw_cal, self).setUp(rad=rad, filetype_resolution=1000) - - def test_raw_calibrate(self): - """Test RAW calibration.""" - res = self.file_handler.get_dataset( - make_dataid(name="C05", calibration="counts"), {}) + with create_file_handler(rad=rad) as file_handler: + res = file_handler.get_dataset( + make_dataid(name="C05", calibration="counts"), {} + ) # We expect the raw data to be unchanged expected = res.data @@ -400,24 +423,7 @@ def test_raw_calibrate(self): assert res.attrs["long_name"] == "Raw Counts" -class Test_NC_ABI_L1B_invalid_cal(Test_NC_ABI_L1B_Base): - """Test the NC_ABI_L1B reader with invalid calibration.""" - - def test_invalid_calibration(self): - """Test detection of invalid calibration values.""" - # Need to use a custom DataID class because the real DataID class is - # smart enough to detect the invalid calibration before the ABI L1B - # get_dataset method gets a chance to run. - class FakeDataID(dict): - def to_dict(self): - return self - - with self.assertRaises(ValueError, msg="Did not detect invalid cal"): - did = FakeDataID(name="C05", calibration="invalid", modifiers=()) - self.file_handler.get_dataset(did, {}) - - -class Test_NC_ABI_File(unittest.TestCase): +class Test_NC_ABI_File: """Test file opening.""" @mock.patch("satpy.readers.abi_base.xr") @@ -434,17 +440,18 @@ def test_open_dataset(self, _): # noqa: PT019 class Test_NC_ABI_L1B_H5netcdf(Test_NC_ABI_L1B): """Allow h5netcdf peculiarities.""" - def setUp(self): + @property + def fake_rad(self): """Create fake data for the tests.""" rad_data = np.int16(50) rad = xr.DataArray( rad_data, attrs={ "scale_factor": 0.5, - "add_offset": -1., + "add_offset": -1.0, "_FillValue": np.array([1002]), "units": "W m-2 um-1 sr-1", "valid_range": (0, 4095), - } + }, ) - super(Test_NC_ABI_L1B_H5netcdf, self).setUp(rad=rad) + return rad From deac453b30cf4936d06fb8eabc0a27941645f2bf Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sat, 28 Oct 2023 20:48:12 -0500 Subject: [PATCH 05/16] Use dask arrays in abi l1b tests --- satpy/tests/reader_tests/test_abi_l1b.py | 32 ++++++++++-------------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index bdaa03f9e5..bc7c5351e8 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -51,7 +51,7 @@ def _create_fake_rad_dataarray( rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( - rad_data, + da.from_array(rad_data), dims=("y", "x"), attrs={ "scale_factor": 0.5, @@ -212,7 +212,7 @@ def fake_rad(self): Needs to be an instance method so the subclass can override it. """ - return None + return None # use default from file handler creator def test_basic_attributes(self): """Test getting basic file attributes.""" @@ -315,14 +315,16 @@ def test_ir_calibrate(self, clip_negative_radiances): @contextlib.contextmanager -def _ir_file_handler( - data: da.Array | None = None, clip_negative_radiances: bool = False -): +def _ir_file_handler(clip_negative_radiances: bool = False): """Create fake data for the tests.""" - if data is None: - data = _fake_ir_data() + values = np.arange(10.0) + rad_data = (values.reshape((2, 5)) + 1.0) * 50.0 + rad_data[0, 0] = -0.0001 # introduce below minimum expected radiance + rad_data = (rad_data + 1.3) / 0.5 + data = rad_data.astype(np.int16) + rad = xr.DataArray( - data, + da.from_array(data), dims=("y", "x"), attrs={ "scale_factor": 0.5, @@ -340,14 +342,6 @@ def _ir_file_handler( yield file_handler -def _fake_ir_data(): - values = np.arange(10.0) - rad_data = (values.reshape((2, 5)) + 1.0) * 50.0 - rad_data[0, 0] = -0.0001 # introduce below minimum expected radiance - rad_data = (rad_data + 1.3) / 0.5 - return rad_data.astype(np.int16) - - class Test_NC_ABI_L1B_vis_cal: """Test the NC_ABI_L1B reader.""" @@ -357,7 +351,7 @@ def test_vis_calibrate(self): rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( - rad_data, + da.from_array(rad_data), dims=("y", "x"), attrs={ "scale_factor": 0.5, @@ -392,7 +386,7 @@ def test_raw_calibrate(self): rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( - rad_data, + da.from_array(rad_data), dims=("y", "x"), attrs={ "scale_factor": 0.5, @@ -445,7 +439,7 @@ def fake_rad(self): """Create fake data for the tests.""" rad_data = np.int16(50) rad = xr.DataArray( - rad_data, + da.from_array(rad_data), attrs={ "scale_factor": 0.5, "add_offset": -1.0, From e24832416dca7aa2d285a3a51a00521f6acf2843 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 29 Oct 2023 10:03:02 -0500 Subject: [PATCH 06/16] Switch some tests to on-disk files --- satpy/readers/abi_base.py | 2 +- satpy/tests/reader_tests/test_abi_l1b.py | 64 +++++++++++++++--------- 2 files changed, 40 insertions(+), 26 deletions(-) diff --git a/satpy/readers/abi_base.py b/satpy/readers/abi_base.py index 956bec278e..28ff91ce38 100644 --- a/satpy/readers/abi_base.py +++ b/satpy/readers/abi_base.py @@ -66,7 +66,7 @@ def nc(self): from satpy.utils import get_dask_chunk_size_in_bytes chunk_size_for_high_res = math.sqrt(get_dask_chunk_size_in_bytes() / 4) # 32-bit floats - chunk_size_for_high_res = np.round(chunk_size_for_high_res / (4 * 226)) * (4 * 226) + chunk_size_for_high_res = np.round(max(chunk_size_for_high_res / (4 * 226), 1)) * (4 * 226) low_res_factor = int(self.filetype_info.get("resolution", 2000) // 500) res_chunk_bytes = int(chunk_size_for_high_res / low_res_factor) * 4 import dask diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index bc7c5351e8..6d8a001918 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -19,8 +19,9 @@ from __future__ import annotations import contextlib +from datetime import datetime from pathlib import Path -from typing import Any, Iterator +from typing import Any, Callable, Iterator from unittest import mock import dask.array as da @@ -28,6 +29,7 @@ import pytest import xarray as xr +from satpy import Scene from satpy.readers.abi_l1b import NC_ABI_L1B from satpy.tests.utils import make_dataid @@ -77,7 +79,7 @@ def _create_fake_rad_dataset(rad=None): range(2), attrs={"scale_factor": -2.0, "add_offset": 1.0}, dims=("y",) ) proj = xr.DataArray( - [], + np.int64(0), attrs={ "semi_major_axis": 1.0, "semi_minor_axis": 1.0, @@ -122,16 +124,27 @@ def _create_fake_rad_dataset(rad=None): def generate_l1b_filename(chan_name: str) -> str: - return f"OR_ABI-L1b-RadC-M4{chan_name}_G16_s20161811540362_e20161811545170_c20161811545230.nc" + return f"OR_ABI-L1b-RadC-M4{chan_name}_G16_s20161811540362_e20161811545170_c20161811545230_suffix.nc" @pytest.fixture(scope="module") -def l1b_c01_file(tmp_path_factory) -> list[Path]: - filename = generate_l1b_filename("C01") - data_path = tmp_path_factory.mktemp("abi_l1b").join(filename) - dataset = _create_fake_rad_dataset() - dataset.to_netcdf(data_path) - return [data_path] +def l1b_c01_file(tmp_path_factory) -> Callable: + def _create_file_handler( + rad: xr.DataArray | None = None, + clip_negative_radiances: bool = False, + ): + filename = generate_l1b_filename("C01") + data_path = tmp_path_factory.mktemp("abi_l1b") / filename + dataset = _create_fake_rad_dataset(rad=rad) + dataset.to_netcdf(data_path) + scn = Scene( + reader="abi_l1b", + filenames=[str(data_path)], + reader_kwargs={"clip_negative_radiances": clip_negative_radiances} + ) + return scn + + return _create_file_handler @pytest.fixture(scope="module") @@ -214,24 +227,17 @@ def fake_rad(self): """ return None # use default from file handler creator - def test_basic_attributes(self): - """Test getting basic file attributes.""" - from datetime import datetime - - with create_file_handler(rad=self.fake_rad) as file_handler: - assert file_handler.start_time == datetime(2017, 9, 20, 17, 30, 40, 800000) - assert file_handler.end_time == datetime(2017, 9, 20, 17, 41, 17, 500000) - - def test_get_dataset(self): + def test_get_dataset(self, l1b_c01_file): """Test the get_dataset method.""" - key = make_dataid(name="Rad", calibration="radiance") - with create_file_handler(rad=self.fake_rad) as file_handler: - res = file_handler.get_dataset(key, {"info": "info"}) + scn = l1b_c01_file(rad=self.fake_rad) + key = make_dataid(name="C01", calibration="radiance") + scn.load([key]) + exp = { "calibration": "radiance", "instrument_ID": None, "modifiers": (), - "name": "Rad", + "name": "C01", "observation_type": "Rad", "orbital_parameters": { "projection_altitude": 1.0, @@ -246,16 +252,24 @@ def test_get_dataset(self): "platform_name": "GOES-16", "platform_shortname": "G16", "production_site": None, - "scan_mode": "M3", + "reader": "abi_l1b", + "resolution": 1000, + "scan_mode": "M4", "scene_abbr": "C", "scene_id": None, "sensor": "abi", "timeline_ID": None, - "suffix": "custom", + "suffix": "suffix", "units": "W m-2 um-1 sr-1", + "start_time": datetime(2017, 9, 20, 17, 30, 40, 800000), + "end_time": datetime(2017, 9, 20, 17, 41, 17, 500000), } - assert res.attrs == exp + res = scn["C01"] + assert "area" in res.attrs + for exp_key, exp_val in exp.items(): + assert res.attrs[exp_key] == exp_val + # we remove any time dimension information assert "t" not in res.coords assert "t" not in res.dims From 07e841c3540982f2d1638d34ba4795bdcf1f4559 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 29 Oct 2023 14:57:55 -0500 Subject: [PATCH 07/16] Move more ABI L1b tests to on-disk files --- satpy/readers/abi_base.py | 6 ++-- satpy/tests/reader_tests/test_abi_l1b.py | 42 +++++++++++++++--------- satpy/utils.py | 19 +++++++++++ 3 files changed, 48 insertions(+), 19 deletions(-) diff --git a/satpy/readers/abi_base.py b/satpy/readers/abi_base.py index 28ff91ce38..3574349c71 100644 --- a/satpy/readers/abi_base.py +++ b/satpy/readers/abi_base.py @@ -23,7 +23,7 @@ import numpy as np import xarray as xr -from pyresample import geometry +from pyresample.geometry import AreaDefinition from satpy._compat import cached_property from satpy.readers import open_file_or_filename @@ -212,7 +212,7 @@ def _get_areadef_latlon(self, key): "fi": float(fi), "pm": float(pm)} - ll_area_def = geometry.AreaDefinition( + ll_area_def = AreaDefinition( self.nc.attrs.get("orbital_slot", "abi_geos"), self.nc.attrs.get("spatial_resolution", "ABI file area"), "abi_latlon", @@ -262,7 +262,7 @@ def _get_areadef_fixedgrid(self, key): "units": "m", "sweep": sweep_axis} - fg_area_def = geometry.AreaDefinition( + fg_area_def = AreaDefinition( self.nc.attrs.get("orbital_slot", "abi_geos"), self.nc.attrs.get("spatial_resolution", "ABI file area"), "abi_fixed_grid", diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index 6d8a001918..7ee9d36acd 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -32,6 +32,7 @@ from satpy import Scene from satpy.readers.abi_l1b import NC_ABI_L1B from satpy.tests.utils import make_dataid +from satpy.utils import ignore_pyproj_proj_warnings RAD_SHAPE = { 500: (3000, 5000), # conus - 500m @@ -276,26 +277,34 @@ def test_get_dataset(self, l1b_c01_file): assert "time" not in res.coords assert "time" not in res.dims - @mock.patch("satpy.readers.abi_base.geometry.AreaDefinition") - def test_get_area_def(self, adef): + def test_get_area_def(self, l1b_c01_file): """Test the area generation.""" - with create_file_handler(rad=self.fake_rad) as file_handler: - file_handler.get_area_def(None) - - assert adef.call_count == 1 - call_args = tuple(adef.call_args)[0] - assert call_args[3] == { - "a": 1.0, - "b": 1.0, + from pyresample.geometry import AreaDefinition + + scn = l1b_c01_file(rad=self.fake_rad) + scn.load(["C01"]) + area_def = scn["C01"].attrs["area"] + assert isinstance(area_def, AreaDefinition) + + with ignore_pyproj_proj_warnings(): + proj_dict = area_def.crs.to_dict() + exp_dict = { "h": 1.0, "lon_0": -90.0, "proj": "geos", "sweep": "x", "units": "m", } - assert call_args[4] == file_handler.ncols - assert call_args[5] == file_handler.nlines - np.testing.assert_allclose(call_args[6], (-2, -2, 8, 2)) + if "R" in proj_dict: + assert proj_dict["R"] == 1 + else: + assert proj_dict["a"] == 1 + assert proj_dict["b"] == 1 + for proj_key, proj_val in exp_dict.items(): + assert proj_dict[proj_key] == proj_val + + assert area_def.shape == scn["C01"].shape + assert area_def.area_extent == (-2, -2, 8, 2) class Test_NC_ABI_L1B_ir_cal: @@ -437,8 +446,6 @@ class Test_NC_ABI_File: @mock.patch("satpy.readers.abi_base.xr") def test_open_dataset(self, _): # noqa: PT019 """Test openning a dataset.""" - from satpy.readers.abi_l1b import NC_ABI_L1B - openable_thing = mock.MagicMock() NC_ABI_L1B(openable_thing, {"platform_shortname": "g16"}, {}) @@ -451,7 +458,10 @@ class Test_NC_ABI_L1B_H5netcdf(Test_NC_ABI_L1B): @property def fake_rad(self): """Create fake data for the tests.""" - rad_data = np.int16(50) + shape = (2, 5) + rad_data = (np.arange(shape[0] * shape[1]).reshape(shape) + 1.0) * 50.0 + rad_data = (rad_data + 1.0) / 0.5 + rad_data = rad_data.astype(np.int16) rad = xr.DataArray( da.from_array(rad_data), attrs={ diff --git a/satpy/utils.py b/satpy/utils.py index f9ea05ca79..dfedc30803 100644 --- a/satpy/utils.py +++ b/satpy/utils.py @@ -576,6 +576,25 @@ def ignore_invalid_float_warnings(): yield +@contextlib.contextmanager +def ignore_pyproj_proj_warnings(): + """Wrap operations that we know will produce a PROJ.4 precision warning. + + Only to be used internally to Pyresample when we have no other choice but + to use PROJ.4 strings/dicts. For example, serialization to YAML or other + human-readable formats or testing the methods that produce the PROJ.4 + versions of the CRS. + + """ + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + "You will likely lose important projection information", + UserWarning, + ) + yield + + def get_chunk_size_limit(dtype=float): """Compute the chunk size limit in bytes given *dtype* (float by default). From b6411c7dc34afbbde38efca5dd13d5a5fe18dd69 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 29 Oct 2023 20:47:53 -0500 Subject: [PATCH 08/16] Switch all ABI L1b tests to on-disk files --- satpy/tests/reader_tests/test_abi_l1b.py | 101 ++++++++--------------- 1 file changed, 33 insertions(+), 68 deletions(-) diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index 7ee9d36acd..e60453f9a3 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -18,10 +18,8 @@ """The abi_l1b reader tests package.""" from __future__ import annotations -import contextlib from datetime import datetime -from pathlib import Path -from typing import Any, Callable, Iterator +from typing import Callable from unittest import mock import dask.array as da @@ -29,7 +27,7 @@ import pytest import xarray as xr -from satpy import Scene +from satpy import DataQuery, Scene from satpy.readers.abi_l1b import NC_ABI_L1B from satpy.tests.utils import make_dataid from satpy.utils import ignore_pyproj_proj_warnings @@ -130,10 +128,7 @@ def generate_l1b_filename(chan_name: str) -> str: @pytest.fixture(scope="module") def l1b_c01_file(tmp_path_factory) -> Callable: - def _create_file_handler( - rad: xr.DataArray | None = None, - clip_negative_radiances: bool = False, - ): + def _create_file_handler(rad: xr.DataArray | None = None): filename = generate_l1b_filename("C01") data_path = tmp_path_factory.mktemp("abi_l1b") / filename dataset = _create_fake_rad_dataset(rad=rad) @@ -141,7 +136,6 @@ def _create_file_handler( scn = Scene( reader="abi_l1b", filenames=[str(data_path)], - reader_kwargs={"clip_negative_radiances": clip_negative_radiances} ) return scn @@ -149,39 +143,23 @@ def _create_file_handler( @pytest.fixture(scope="module") -def l1b_all_files( - l1b_c01_file, -) -> list[Path]: - return l1b_c01_file - - -@contextlib.contextmanager -def create_file_handler( - rad: xr.DataArray | None = None, - clip_negative_radiances: bool = False, - filetype_resolution: int = 0, -) -> Iterator[NC_ABI_L1B]: - """Create a fake dataset using the given radiance data.""" - - ft_info: dict[str, Any] = {"filetype": "info"} - if filetype_resolution: - ft_info["resolution"] = filetype_resolution - - with mock.patch("satpy.readers.abi_base.xr") as xr_: - xr_.open_dataset.return_value = _create_fake_rad_dataset(rad=rad) - file_handler = NC_ABI_L1B( - "filename", - { - "platform_shortname": "G16", - "observation_type": "Rad", - "suffix": "custom", - "scene_abbr": "C", - "scan_mode": "M3", - }, - ft_info, - clip_negative_radiances=clip_negative_radiances, +def l1b_c07_file(tmp_path_factory) -> Callable: + def _create_file_handler( + rad: xr.DataArray | None = None, + clip_negative_radiances: bool = False, + ): + filename = generate_l1b_filename("C07") + data_path = tmp_path_factory.mktemp("abi_l1b") / filename + dataset = _create_fake_rad_dataset(rad=rad) + dataset.to_netcdf(data_path) + scn = Scene( + reader="abi_l1b", + filenames=[str(data_path)], + reader_kwargs={"clip_negative_radiances": clip_negative_radiances} ) - yield file_handler + return scn + + return _create_file_handler class TestABIYAML: @@ -311,15 +289,11 @@ class Test_NC_ABI_L1B_ir_cal: """Test the NC_ABI_L1B reader's default IR calibration.""" @pytest.mark.parametrize("clip_negative_radiances", [False, True]) - def test_ir_calibrate(self, clip_negative_radiances): + def test_ir_calibrate(self, l1b_c07_file, clip_negative_radiances): """Test IR calibration.""" - with _ir_file_handler( - clip_negative_radiances=clip_negative_radiances - ) as file_handler: - res = file_handler.get_dataset( - make_dataid(name="C07", calibration="brightness_temperature"), {} - ) - assert file_handler.clip_negative_radiances == clip_negative_radiances + scn = l1b_c07_file(rad=_fake_ir_data(), clip_negative_radiances=clip_negative_radiances) + scn.load([DataQuery(name="C07", calibration="brightness_temperature")]) + res = scn["C07"] clipped_ir = 134.68753 if clip_negative_radiances else np.nan expected = np.array( @@ -337,9 +311,7 @@ def test_ir_calibrate(self, clip_negative_radiances): assert res.attrs["long_name"] == "Brightness Temperature" -@contextlib.contextmanager -def _ir_file_handler(clip_negative_radiances: bool = False): - """Create fake data for the tests.""" +def _fake_ir_data(): values = np.arange(10.0) rad_data = (values.reshape((2, 5)) + 1.0) * 50.0 rad_data[0, 0] = -0.0001 # introduce below minimum expected radiance @@ -357,18 +329,13 @@ def _ir_file_handler(clip_negative_radiances: bool = False): ), # last rad_data value }, ) - with create_file_handler( - rad=rad, - clip_negative_radiances=clip_negative_radiances, - filetype_resolution=2000, - ) as file_handler: - yield file_handler + return rad class Test_NC_ABI_L1B_vis_cal: """Test the NC_ABI_L1B reader.""" - def test_vis_calibrate(self): + def test_vis_calibrate(self, l1b_c01_file): """Test VIS calibration.""" rad_data = np.arange(10.0).reshape((2, 5)) + 1.0 rad_data = (rad_data + 1.0) / 0.5 @@ -382,10 +349,9 @@ def test_vis_calibrate(self): "_FillValue": 20, }, ) - with create_file_handler(rad=rad, filetype_resolution=1000) as file_handler: - res = file_handler.get_dataset( - make_dataid(name="C05", calibration="reflectance"), {} - ) + scn = l1b_c01_file(rad=rad) + scn.load(["C01"]) + res = scn["C01"] expected = np.array( [ @@ -403,7 +369,7 @@ def test_vis_calibrate(self): class Test_NC_ABI_L1B_raw_cal: """Test the NC_ABI_L1B reader raw calibration.""" - def test_raw_calibrate(self): + def test_raw_calibrate(self, l1b_c01_file): """Test RAW calibration.""" rad_data = np.arange(10.0).reshape((2, 5)) + 1.0 rad_data = (rad_data + 1.0) / 0.5 @@ -417,10 +383,9 @@ def test_raw_calibrate(self): "_FillValue": 20, }, ) - with create_file_handler(rad=rad) as file_handler: - res = file_handler.get_dataset( - make_dataid(name="C05", calibration="counts"), {} - ) + scn = l1b_c01_file(rad=rad) + scn.load([DataQuery(name="C01", calibration="counts")]) + res = scn["C01"] # We expect the raw data to be unchanged expected = res.data From f9efd963ba6a50cbead6979095db4786c456c31d Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 29 Oct 2023 21:16:10 -0500 Subject: [PATCH 09/16] Use more realistic sizes in ABI tests --- satpy/tests/reader_tests/test_abi_l1b.py | 49 +++++++++++++----------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index e60453f9a3..fb08ef9361 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -41,12 +41,12 @@ def _create_fake_rad_dataarray( rad: xr.DataArray | None = None, - # resolution: int = 2000, + resolution: int = 2000, ) -> xr.DataArray: x_image = xr.DataArray(0.0) y_image = xr.DataArray(0.0) time = xr.DataArray(0.0) - shape = (2, 5) # RAD_SHAPE[resolution] + shape = RAD_SHAPE[resolution] if rad is None: rad_data = (np.arange(shape[0] * shape[1]).reshape(shape) + 1.0) * 50.0 rad_data = (rad_data + 1.0) / 0.5 @@ -68,14 +68,14 @@ def _create_fake_rad_dataarray( return rad -def _create_fake_rad_dataset(rad=None): - rad = _create_fake_rad_dataarray(rad=rad) +def _create_fake_rad_dataset(rad: xr.DataArray, resolution: int) -> xr.Dataset: + rad = _create_fake_rad_dataarray(rad=rad, resolution=resolution) x__ = xr.DataArray( - range(5), attrs={"scale_factor": 2.0, "add_offset": -1.0}, dims=("x",) + range(rad.shape[1]), attrs={"scale_factor": 2.0, "add_offset": -1.0}, dims=("x",) ) y__ = xr.DataArray( - range(2), attrs={"scale_factor": -2.0, "add_offset": 1.0}, dims=("y",) + range(rad.shape[0]), attrs={"scale_factor": -2.0, "add_offset": 1.0}, dims=("y",) ) proj = xr.DataArray( np.int64(0), @@ -131,7 +131,7 @@ def l1b_c01_file(tmp_path_factory) -> Callable: def _create_file_handler(rad: xr.DataArray | None = None): filename = generate_l1b_filename("C01") data_path = tmp_path_factory.mktemp("abi_l1b") / filename - dataset = _create_fake_rad_dataset(rad=rad) + dataset = _create_fake_rad_dataset(rad=rad, resolution=1000) dataset.to_netcdf(data_path) scn = Scene( reader="abi_l1b", @@ -145,12 +145,12 @@ def _create_file_handler(rad: xr.DataArray | None = None): @pytest.fixture(scope="module") def l1b_c07_file(tmp_path_factory) -> Callable: def _create_file_handler( - rad: xr.DataArray | None = None, + rad: xr.DataArray, clip_negative_radiances: bool = False, ): filename = generate_l1b_filename("C07") data_path = tmp_path_factory.mktemp("abi_l1b") / filename - dataset = _create_fake_rad_dataset(rad=rad) + dataset = _create_fake_rad_dataset(rad=rad, resolution=2000) dataset.to_netcdf(data_path) scn = Scene( reader="abi_l1b", @@ -204,7 +204,7 @@ def fake_rad(self): Needs to be an instance method so the subclass can override it. """ - return None # use default from file handler creator + return None def test_get_dataset(self, l1b_c01_file): """Test the get_dataset method.""" @@ -282,7 +282,7 @@ def test_get_area_def(self, l1b_c01_file): assert proj_dict[proj_key] == proj_val assert area_def.shape == scn["C01"].shape - assert area_def.area_extent == (-2, -2, 8, 2) + assert area_def.area_extent == (-2.0, -2998.0, 4998.0, 2.0) class Test_NC_ABI_L1B_ir_cal: @@ -298,11 +298,12 @@ def test_ir_calibrate(self, l1b_c07_file, clip_negative_radiances): clipped_ir = 134.68753 if clip_negative_radiances else np.nan expected = np.array( [ - [clipped_ir, 304.97037, 332.22778, 354.6147, 374.08688], - [391.58655, 407.64786, 422.60635, 436.68802, np.nan], + clipped_ir, 304.97037, 332.22778, 354.6147, 374.08688, + 391.58655, 407.64786, 422.60635, 436.68802, np.nan, ] ) - np.testing.assert_allclose(res.data, expected, equal_nan=True, atol=1e-04) + data_np = res.data.compute() + np.testing.assert_allclose(data_np[0, :10], expected, equal_nan=True, atol=1e-04) # make sure the attributes from the file are in the data array assert "scale_factor" not in res.attrs @@ -312,8 +313,9 @@ def test_ir_calibrate(self, l1b_c07_file, clip_negative_radiances): def _fake_ir_data(): - values = np.arange(10.0) - rad_data = (values.reshape((2, 5)) + 1.0) * 50.0 + shape = RAD_SHAPE[2000] + values = np.arange(shape[0] * shape[1]) + rad_data = (values.reshape(shape) + 1.0) * 50.0 rad_data[0, 0] = -0.0001 # introduce below minimum expected radiance rad_data = (rad_data + 1.3) / 0.5 data = rad_data.astype(np.int16) @@ -337,7 +339,8 @@ class Test_NC_ABI_L1B_vis_cal: def test_vis_calibrate(self, l1b_c01_file): """Test VIS calibration.""" - rad_data = np.arange(10.0).reshape((2, 5)) + 1.0 + shape = RAD_SHAPE[1000] + rad_data = np.arange(shape[0] * shape[1]).reshape(shape) + 1.0 rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( @@ -355,11 +358,12 @@ def test_vis_calibrate(self, l1b_c01_file): expected = np.array( [ - [0.15265617, 0.30531234, 0.45796851, 0.61062468, 0.76328085], - [0.91593702, 1.06859319, 1.22124936, np.nan, 1.52656171], + 0.15265617, 0.30531234, 0.45796851, 0.61062468, 0.76328085, + 0.91593702, 1.06859319, 1.22124936, np.nan, 1.52656171, ] ) - assert np.allclose(res.data, expected, equal_nan=True) + data_np = res.data.compute() + assert np.allclose(data_np[0, :10], expected, equal_nan=True) assert "scale_factor" not in res.attrs assert "_FillValue" not in res.attrs assert res.attrs["standard_name"] == "toa_bidirectional_reflectance" @@ -371,7 +375,8 @@ class Test_NC_ABI_L1B_raw_cal: def test_raw_calibrate(self, l1b_c01_file): """Test RAW calibration.""" - rad_data = np.arange(10.0).reshape((2, 5)) + 1.0 + shape = RAD_SHAPE[1000] + rad_data = np.arange(shape[0] * shape[1]).reshape(shape) + 1.0 rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( @@ -423,7 +428,7 @@ class Test_NC_ABI_L1B_H5netcdf(Test_NC_ABI_L1B): @property def fake_rad(self): """Create fake data for the tests.""" - shape = (2, 5) + shape = RAD_SHAPE[1000] rad_data = (np.arange(shape[0] * shape[1]).reshape(shape) + 1.0) * 50.0 rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) From 1e0e21f083cbf2cde44dbabbe5c9cc4287df091e Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 29 Oct 2023 21:29:33 -0500 Subject: [PATCH 10/16] Revert AreaDefinition import for easier test mocking --- satpy/readers/abi_base.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/satpy/readers/abi_base.py b/satpy/readers/abi_base.py index 3574349c71..28ff91ce38 100644 --- a/satpy/readers/abi_base.py +++ b/satpy/readers/abi_base.py @@ -23,7 +23,7 @@ import numpy as np import xarray as xr -from pyresample.geometry import AreaDefinition +from pyresample import geometry from satpy._compat import cached_property from satpy.readers import open_file_or_filename @@ -212,7 +212,7 @@ def _get_areadef_latlon(self, key): "fi": float(fi), "pm": float(pm)} - ll_area_def = AreaDefinition( + ll_area_def = geometry.AreaDefinition( self.nc.attrs.get("orbital_slot", "abi_geos"), self.nc.attrs.get("spatial_resolution", "ABI file area"), "abi_latlon", @@ -262,7 +262,7 @@ def _get_areadef_fixedgrid(self, key): "units": "m", "sweep": sweep_axis} - fg_area_def = AreaDefinition( + fg_area_def = geometry.AreaDefinition( self.nc.attrs.get("orbital_slot", "abi_geos"), self.nc.attrs.get("spatial_resolution", "ABI file area"), "abi_fixed_grid", From 514f5e1bef90d57282993ee04501c70cf7b4ea28 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 30 Oct 2023 11:52:33 -0500 Subject: [PATCH 11/16] More abi l1b test refactoring --- satpy/tests/reader_tests/test_abi_l1b.py | 446 +++++++++++------------ 1 file changed, 214 insertions(+), 232 deletions(-) diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index fb08ef9361..54c78e1089 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -19,17 +19,18 @@ from __future__ import annotations from datetime import datetime -from typing import Callable +from pathlib import Path +from typing import Any, Callable from unittest import mock import dask.array as da import numpy as np import pytest import xarray as xr +from pytest_lazyfixture import lazy_fixture from satpy import DataQuery, Scene from satpy.readers.abi_l1b import NC_ABI_L1B -from satpy.tests.utils import make_dataid from satpy.utils import ignore_pyproj_proj_warnings RAD_SHAPE = { @@ -72,10 +73,14 @@ def _create_fake_rad_dataset(rad: xr.DataArray, resolution: int) -> xr.Dataset: rad = _create_fake_rad_dataarray(rad=rad, resolution=resolution) x__ = xr.DataArray( - range(rad.shape[1]), attrs={"scale_factor": 2.0, "add_offset": -1.0}, dims=("x",) + range(rad.shape[1]), + attrs={"scale_factor": 2.0, "add_offset": -1.0}, + dims=("x",), ) y__ = xr.DataArray( - range(rad.shape[0]), attrs={"scale_factor": -2.0, "add_offset": 1.0}, dims=("y",) + range(rad.shape[0]), + attrs={"scale_factor": -2.0, "add_offset": 1.0}, + dims=("y",), ) proj = xr.DataArray( np.int64(0), @@ -126,92 +131,144 @@ def generate_l1b_filename(chan_name: str) -> str: return f"OR_ABI-L1b-RadC-M4{chan_name}_G16_s20161811540362_e20161811545170_c20161811545230_suffix.nc" -@pytest.fixture(scope="module") -def l1b_c01_file(tmp_path_factory) -> Callable: - def _create_file_handler(rad: xr.DataArray | None = None): - filename = generate_l1b_filename("C01") - data_path = tmp_path_factory.mktemp("abi_l1b") / filename - dataset = _create_fake_rad_dataset(rad=rad, resolution=1000) - dataset.to_netcdf(data_path) - scn = Scene( - reader="abi_l1b", - filenames=[str(data_path)], - ) - return scn +@pytest.fixture() +def c01_refl(tmp_path) -> xr.DataArray: + scn = _create_scene_for_data(tmp_path, "C01", None, 1000) + scn.load(["C01"]) + return scn["C01"] + + +@pytest.fixture() +def c01_rad(tmp_path) -> xr.DataArray: + scn = _create_scene_for_data(tmp_path, "C01", None, 1000) + scn.load([DataQuery(name="C01", calibration="radiance")]) + return scn["C01"] + - return _create_file_handler +@pytest.fixture() +def c01_rad_h5netcdf(tmp_path) -> xr.DataArray: + shape = RAD_SHAPE[1000] + rad_data = (np.arange(shape[0] * shape[1]).reshape(shape) + 1.0) * 50.0 + rad_data = (rad_data + 1.0) / 0.5 + rad_data = rad_data.astype(np.int16) + rad = xr.DataArray( + da.from_array(rad_data), + attrs={ + "scale_factor": 0.5, + "add_offset": -1.0, + "_FillValue": np.array([1002]), + "units": "W m-2 um-1 sr-1", + "valid_range": (0, 4095), + }, + ) + scn = _create_scene_for_data(tmp_path, "C01", rad, 1000) + scn.load([DataQuery(name="C01", calibration="radiance")]) + return scn["C01"] + + +@pytest.fixture() +def c01_counts(tmp_path) -> xr.DataArray: + scn = _create_scene_for_data(tmp_path, "C01", None, 1000) + scn.load([DataQuery(name="C01", calibration="counts")]) + return scn["C01"] + + +def _create_scene_for_data( + tmp_path: Path, + channel_name: str, + rad: xr.DataArray | None, + resolution: int, + reader_kwargs: dict[str, Any] | None = None, +) -> Scene: + filename = generate_l1b_filename(channel_name) + data_path = tmp_path / filename + dataset = _create_fake_rad_dataset(rad=rad, resolution=resolution) + dataset.to_netcdf(data_path) + scn = Scene( + reader="abi_l1b", + filenames=[str(data_path)], + reader_kwargs=reader_kwargs, + ) + return scn -@pytest.fixture(scope="module") -def l1b_c07_file(tmp_path_factory) -> Callable: - def _create_file_handler( - rad: xr.DataArray, - clip_negative_radiances: bool = False, +@pytest.fixture() +def c07_bt_creator(tmp_path) -> Callable: + def _load_data_array( + clip_negative_radiances: bool = False, ): - filename = generate_l1b_filename("C07") - data_path = tmp_path_factory.mktemp("abi_l1b") / filename - dataset = _create_fake_rad_dataset(rad=rad, resolution=2000) - dataset.to_netcdf(data_path) - scn = Scene( - reader="abi_l1b", - filenames=[str(data_path)], - reader_kwargs={"clip_negative_radiances": clip_negative_radiances} + rad = _fake_c07_data() + scn = _create_scene_for_data( + tmp_path, + "C07", + rad, + 2000, + {"clip_negative_radiances": clip_negative_radiances}, ) - return scn + scn.load(["C07"]) + return scn["C07"] - return _create_file_handler + return _load_data_array -class TestABIYAML: - """Tests for the ABI L1b reader's YAML configuration.""" - - @pytest.mark.parametrize( - ("channel", "suffix"), - [ - ("C{:02d}".format(num), suffix) - for num in range(1, 17) - for suffix in ("", "_test_suffix") - ], +def _fake_c07_data() -> xr.DataArray: + shape = RAD_SHAPE[2000] + values = np.arange(shape[0] * shape[1]) + rad_data = (values.reshape(shape) + 1.0) * 50.0 + rad_data[0, 0] = -0.0001 # introduce below minimum expected radiance + rad_data = (rad_data + 1.3) / 0.5 + data = rad_data.astype(np.int16) + rad = xr.DataArray( + da.from_array(data), + dims=("y", "x"), + attrs={ + "scale_factor": 0.5, + "add_offset": -1.3, + "_FillValue": np.int16( + np.floor(((9 + 1) * 50.0 + 1.3) / 0.5) + ), # last rad_data value + }, ) - def test_file_patterns_match(self, channel, suffix): - """Test that the configured file patterns work.""" - from satpy.readers import configs_for_reader, load_reader + return rad - reader_configs = list(configs_for_reader("abi_l1b"))[0] - reader = load_reader(reader_configs) - fn1 = ( + +@pytest.mark.parametrize( + ("channel", "suffix"), + [ + ("C{:02d}".format(num), suffix) + for num in range(1, 17) + for suffix in ("", "_test_suffix") + ], +) +def test_file_patterns_match(channel, suffix): + """Test that the configured file patterns work.""" + from satpy.readers import configs_for_reader, load_reader + + reader_configs = list(configs_for_reader("abi_l1b"))[0] + reader = load_reader(reader_configs) + fn1 = ( + "OR_ABI-L1b-RadM1-M3{}_G16_s20182541300210_e20182541300267" + "_c20182541300308{}.nc" + ).format(channel, suffix) + loadables = reader.select_files_from_pathnames([fn1]) + assert len(loadables) == 1 + if not suffix and channel in ["C01", "C02", "C03", "C05"]: + fn2 = ( "OR_ABI-L1b-RadM1-M3{}_G16_s20182541300210_e20182541300267" - "_c20182541300308{}.nc" - ).format(channel, suffix) - loadables = reader.select_files_from_pathnames([fn1]) + "_c20182541300308-000000_0.nc" + ).format(channel) + loadables = reader.select_files_from_pathnames([fn2]) assert len(loadables) == 1 - if not suffix and channel in ["C01", "C02", "C03", "C05"]: - fn2 = ( - "OR_ABI-L1b-RadM1-M3{}_G16_s20182541300210_e20182541300267" - "_c20182541300308-000000_0.nc" - ).format(channel) - loadables = reader.select_files_from_pathnames([fn2]) - assert len(loadables) == 1 +@pytest.mark.parametrize( + "c01_data_arr", [lazy_fixture("c01_rad"), lazy_fixture("c01_rad_h5netcdf")] +) class Test_NC_ABI_L1B: """Test the NC_ABI_L1B reader.""" - @property - def fake_rad(self): - """Create fake data for these tests. - - Needs to be an instance method so the subclass can override it. - - """ - return None - - def test_get_dataset(self, l1b_c01_file): + def test_get_dataset(self, c01_data_arr): """Test the get_dataset method.""" - scn = l1b_c01_file(rad=self.fake_rad) - key = make_dataid(name="C01", calibration="radiance") - scn.load([key]) - exp = { "calibration": "radiance", "instrument_ID": None, @@ -244,7 +301,7 @@ def test_get_dataset(self, l1b_c01_file): "end_time": datetime(2017, 9, 20, 17, 41, 17, 500000), } - res = scn["C01"] + res = c01_data_arr assert "area" in res.attrs for exp_key, exp_val in exp.items(): assert res.attrs[exp_key] == exp_val @@ -255,13 +312,11 @@ def test_get_dataset(self, l1b_c01_file): assert "time" not in res.coords assert "time" not in res.dims - def test_get_area_def(self, l1b_c01_file): + def test_get_area_def(self, c01_data_arr): """Test the area generation.""" from pyresample.geometry import AreaDefinition - scn = l1b_c01_file(rad=self.fake_rad) - scn.load(["C01"]) - area_def = scn["C01"].attrs["area"] + area_def = c01_data_arr.attrs["area"] assert isinstance(area_def, AreaDefinition) with ignore_pyproj_proj_warnings(): @@ -281,165 +336,92 @@ def test_get_area_def(self, l1b_c01_file): for proj_key, proj_val in exp_dict.items(): assert proj_dict[proj_key] == proj_val - assert area_def.shape == scn["C01"].shape + assert area_def.shape == c01_data_arr.shape assert area_def.area_extent == (-2.0, -2998.0, 4998.0, 2.0) -class Test_NC_ABI_L1B_ir_cal: - """Test the NC_ABI_L1B reader's default IR calibration.""" - - @pytest.mark.parametrize("clip_negative_radiances", [False, True]) - def test_ir_calibrate(self, l1b_c07_file, clip_negative_radiances): - """Test IR calibration.""" - scn = l1b_c07_file(rad=_fake_ir_data(), clip_negative_radiances=clip_negative_radiances) - scn.load([DataQuery(name="C07", calibration="brightness_temperature")]) - res = scn["C07"] - - clipped_ir = 134.68753 if clip_negative_radiances else np.nan - expected = np.array( - [ - clipped_ir, 304.97037, 332.22778, 354.6147, 374.08688, - 391.58655, 407.64786, 422.60635, 436.68802, np.nan, - ] - ) - data_np = res.data.compute() - np.testing.assert_allclose(data_np[0, :10], expected, equal_nan=True, atol=1e-04) - - # make sure the attributes from the file are in the data array - assert "scale_factor" not in res.attrs - assert "_FillValue" not in res.attrs - assert res.attrs["standard_name"] == "toa_brightness_temperature" - assert res.attrs["long_name"] == "Brightness Temperature" - - -def _fake_ir_data(): - shape = RAD_SHAPE[2000] - values = np.arange(shape[0] * shape[1]) - rad_data = (values.reshape(shape) + 1.0) * 50.0 - rad_data[0, 0] = -0.0001 # introduce below minimum expected radiance - rad_data = (rad_data + 1.3) / 0.5 - data = rad_data.astype(np.int16) - - rad = xr.DataArray( - da.from_array(data), - dims=("y", "x"), - attrs={ - "scale_factor": 0.5, - "add_offset": -1.3, - "_FillValue": np.int16( - np.floor(((9 + 1) * 50.0 + 1.3) / 0.5) - ), # last rad_data value - }, +@pytest.mark.parametrize("clip_negative_radiances", [False, True]) +def test_ir_calibrate(self, c07_bt_creator, clip_negative_radiances): + """Test IR calibration.""" + res = c07_bt_creator(clip_negative_radiances=clip_negative_radiances) + clipped_ir = 134.68753 if clip_negative_radiances else np.nan + expected = np.array( + [ + clipped_ir, + 304.97037, + 332.22778, + 354.6147, + 374.08688, + 391.58655, + 407.64786, + 422.60635, + 436.68802, + np.nan, + ] + ) + data_np = res.data.compute() + np.testing.assert_allclose( + data_np[0, :10], expected, equal_nan=True, atol=1e-04 ) - return rad - - -class Test_NC_ABI_L1B_vis_cal: - """Test the NC_ABI_L1B reader.""" - - def test_vis_calibrate(self, l1b_c01_file): - """Test VIS calibration.""" - shape = RAD_SHAPE[1000] - rad_data = np.arange(shape[0] * shape[1]).reshape(shape) + 1.0 - rad_data = (rad_data + 1.0) / 0.5 - rad_data = rad_data.astype(np.int16) - rad = xr.DataArray( - da.from_array(rad_data), - dims=("y", "x"), - attrs={ - "scale_factor": 0.5, - "add_offset": -1.0, - "_FillValue": 20, - }, - ) - scn = l1b_c01_file(rad=rad) - scn.load(["C01"]) - res = scn["C01"] - - expected = np.array( - [ - 0.15265617, 0.30531234, 0.45796851, 0.61062468, 0.76328085, - 0.91593702, 1.06859319, 1.22124936, np.nan, 1.52656171, - ] - ) - data_np = res.data.compute() - assert np.allclose(data_np[0, :10], expected, equal_nan=True) - assert "scale_factor" not in res.attrs - assert "_FillValue" not in res.attrs - assert res.attrs["standard_name"] == "toa_bidirectional_reflectance" - assert res.attrs["long_name"] == "Bidirectional Reflectance" - - -class Test_NC_ABI_L1B_raw_cal: - """Test the NC_ABI_L1B reader raw calibration.""" - - def test_raw_calibrate(self, l1b_c01_file): - """Test RAW calibration.""" - shape = RAD_SHAPE[1000] - rad_data = np.arange(shape[0] * shape[1]).reshape(shape) + 1.0 - rad_data = (rad_data + 1.0) / 0.5 - rad_data = rad_data.astype(np.int16) - rad = xr.DataArray( - da.from_array(rad_data), - dims=("y", "x"), - attrs={ - "scale_factor": 0.5, - "add_offset": -1.0, - "_FillValue": 20, - }, - ) - scn = l1b_c01_file(rad=rad) - scn.load([DataQuery(name="C01", calibration="counts")]) - res = scn["C01"] - - # We expect the raw data to be unchanged - expected = res.data - assert np.allclose(res.data, expected, equal_nan=True) - - # check for the presence of typical attributes - assert "scale_factor" in res.attrs - assert "add_offset" in res.attrs - assert "_FillValue" in res.attrs - assert "orbital_parameters" in res.attrs - assert "platform_shortname" in res.attrs - assert "scene_id" in res.attrs - - # determine if things match their expected values/types. - assert res.data.dtype == np.int16 - assert res.attrs["standard_name"] == "counts" - assert res.attrs["long_name"] == "Raw Counts" - - -class Test_NC_ABI_File: - """Test file opening.""" - - @mock.patch("satpy.readers.abi_base.xr") - def test_open_dataset(self, _): # noqa: PT019 - """Test openning a dataset.""" - openable_thing = mock.MagicMock() - - NC_ABI_L1B(openable_thing, {"platform_shortname": "g16"}, {}) - openable_thing.open.assert_called() + # make sure the attributes from the file are in the data array + assert "scale_factor" not in res.attrs + assert "_FillValue" not in res.attrs + assert res.attrs["standard_name"] == "toa_brightness_temperature" + assert res.attrs["long_name"] == "Brightness Temperature" -class Test_NC_ABI_L1B_H5netcdf(Test_NC_ABI_L1B): - """Allow h5netcdf peculiarities.""" - @property - def fake_rad(self): - """Create fake data for the tests.""" - shape = RAD_SHAPE[1000] - rad_data = (np.arange(shape[0] * shape[1]).reshape(shape) + 1.0) * 50.0 - rad_data = (rad_data + 1.0) / 0.5 - rad_data = rad_data.astype(np.int16) - rad = xr.DataArray( - da.from_array(rad_data), - attrs={ - "scale_factor": 0.5, - "add_offset": -1.0, - "_FillValue": np.array([1002]), - "units": "W m-2 um-1 sr-1", - "valid_range": (0, 4095), - }, - ) - return rad +def test_vis_calibrate(c01_refl): + """Test VIS calibration.""" + res = c01_refl + expected = np.array( + [ + 7.632808, + 15.265616, + 22.898426, + 30.531233, + 38.164043, + 45.796852, + 53.429657, + 61.062466, + 68.695274, + np.nan, + ] + ) + data_np = res.data.compute() + np.testing.assert_allclose(data_np[0, :10], expected, equal_nan=True) + assert "scale_factor" not in res.attrs + assert "_FillValue" not in res.attrs + assert res.attrs["standard_name"] == "toa_bidirectional_reflectance" + assert res.attrs["long_name"] == "Bidirectional Reflectance" + + +def test_raw_calibrate(c01_counts): + """Test RAW calibration.""" + res = c01_counts + + # We expect the raw data to be unchanged + expected = res.data + assert np.allclose(res.data, expected, equal_nan=True) + + # check for the presence of typical attributes + assert "scale_factor" in res.attrs + assert "add_offset" in res.attrs + assert "_FillValue" in res.attrs + assert "orbital_parameters" in res.attrs + assert "platform_shortname" in res.attrs + assert "scene_id" in res.attrs + + # determine if things match their expected values/types. + assert res.data.dtype == np.int16 + assert res.attrs["standard_name"] == "counts" + assert res.attrs["long_name"] == "Raw Counts" + + +@mock.patch("satpy.readers.abi_base.xr") +def test_open_dataset(_): # noqa: PT019 + """Test opening a dataset.""" + openable_thing = mock.MagicMock() + + NC_ABI_L1B(openable_thing, {"platform_shortname": "g16"}, {}) + openable_thing.open.assert_called() From 83609cca303a5ae3abe3e4b1959f16d77f26b6bb Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 30 Oct 2023 13:11:47 -0500 Subject: [PATCH 12/16] Undo forcing GRB fill to floating point Caused failure in GLM L2 DQF processing --- satpy/readers/abi_base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/abi_base.py b/satpy/readers/abi_base.py index 28ff91ce38..3fdd724e12 100644 --- a/satpy/readers/abi_base.py +++ b/satpy/readers/abi_base.py @@ -140,7 +140,7 @@ def is_int(val): new_fill = fill else: new_fill = np.nan - data = data.where(data != fill, np.float32(new_fill)) + data = data.where(data != fill, new_fill) if factor != 1 and item in ("x", "y"): # be more precise with x/y coordinates # see get_area_def for more information From 8b5c450509105cff28881a621d7044cda309f665 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 30 Oct 2023 14:38:07 -0500 Subject: [PATCH 13/16] Fix various inconsistencies in ABI L1b DataArrays --- satpy/readers/abi_base.py | 7 +- satpy/readers/abi_l1b.py | 1 + satpy/tests/reader_tests/test_abi_l1b.py | 141 +++++++++++++---------- 3 files changed, 84 insertions(+), 65 deletions(-) diff --git a/satpy/readers/abi_base.py b/satpy/readers/abi_base.py index 3fdd724e12..1c0824ab27 100644 --- a/satpy/readers/abi_base.py +++ b/satpy/readers/abi_base.py @@ -139,18 +139,13 @@ def is_int(val): if is_int(data) and is_int(factor) and is_int(offset): new_fill = fill else: - new_fill = np.nan + new_fill = np.float32(np.nan) data = data.where(data != fill, new_fill) if factor != 1 and item in ("x", "y"): # be more precise with x/y coordinates # see get_area_def for more information data = data * np.round(float(factor), 6) + np.round(float(offset), 6) elif factor != 1: - # make sure the factor is a 64-bit float - # can't do this in place since data is most likely uint16 - # and we are making it a 64-bit float - if not is_int(factor): - factor = np.float32(factor) data = data * np.float32(factor) + np.float32(offset) return data diff --git a/satpy/readers/abi_l1b.py b/satpy/readers/abi_l1b.py index c3da53c9c7..4933b0982a 100644 --- a/satpy/readers/abi_l1b.py +++ b/satpy/readers/abi_l1b.py @@ -49,6 +49,7 @@ def get_dataset(self, key, info): # For raw cal, don't apply scale and offset, return raw file counts if key["calibration"] == "counts": radiances = self.nc["Rad"].copy() + radiances = self._adjust_coords(radiances, "Rad") else: radiances = self["Rad"] diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index 54c78e1089..61f1746bf9 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -25,6 +25,7 @@ import dask.array as da import numpy as np +import numpy.typing as npt import pytest import xarray as xr from pytest_lazyfixture import lazy_fixture @@ -153,6 +154,7 @@ def c01_rad_h5netcdf(tmp_path) -> xr.DataArray: rad_data = rad_data.astype(np.int16) rad = xr.DataArray( da.from_array(rad_data), + dims=("y", "x"), attrs={ "scale_factor": 0.5, "add_offset": -1.0, @@ -173,25 +175,6 @@ def c01_counts(tmp_path) -> xr.DataArray: return scn["C01"] -def _create_scene_for_data( - tmp_path: Path, - channel_name: str, - rad: xr.DataArray | None, - resolution: int, - reader_kwargs: dict[str, Any] | None = None, -) -> Scene: - filename = generate_l1b_filename(channel_name) - data_path = tmp_path / filename - dataset = _create_fake_rad_dataset(rad=rad, resolution=resolution) - dataset.to_netcdf(data_path) - scn = Scene( - reader="abi_l1b", - filenames=[str(data_path)], - reader_kwargs=reader_kwargs, - ) - return scn - - @pytest.fixture() def c07_bt_creator(tmp_path) -> Callable: def _load_data_array( @@ -232,6 +215,73 @@ def _fake_c07_data() -> xr.DataArray: return rad +def _create_scene_for_data( + tmp_path: Path, + channel_name: str, + rad: xr.DataArray | None, + resolution: int, + reader_kwargs: dict[str, Any] | None = None, +) -> Scene: + filename = generate_l1b_filename(channel_name) + data_path = tmp_path / filename + dataset = _create_fake_rad_dataset(rad=rad, resolution=resolution) + dataset.to_netcdf(data_path) + scn = Scene( + reader="abi_l1b", + filenames=[str(data_path)], + reader_kwargs=reader_kwargs, + ) + return scn + + +def _get_and_check_array(data_arr: xr.DataArray, exp_dtype: npt.DTypeLike) -> npt.NDArray: + data_np = data_arr.data.compute() + assert data_np.dtype == data_arr.dtype + assert data_np.dtype == exp_dtype + return data_np + +def _check_area(data_arr: xr.DataArray) -> None: + from pyresample.geometry import AreaDefinition + + area_def = data_arr.attrs["area"] + assert isinstance(area_def, AreaDefinition) + + with ignore_pyproj_proj_warnings(): + proj_dict = area_def.crs.to_dict() + exp_dict = { + "h": 1.0, + "lon_0": -90.0, + "proj": "geos", + "sweep": "x", + "units": "m", + } + if "R" in proj_dict: + assert proj_dict["R"] == 1 + else: + assert proj_dict["a"] == 1 + assert proj_dict["b"] == 1 + for proj_key, proj_val in exp_dict.items(): + assert proj_dict[proj_key] == proj_val + + assert area_def.shape == data_arr.shape + if area_def.shape[0] == RAD_SHAPE[1000][0]: + exp_extent = (-2.0, -2998.0, 4998.0, 2.0) + else: + exp_extent = (-2.0, -1498.0, 2498.0, 2.0) + assert area_def.area_extent == exp_extent + + +def _check_dims_and_coords(data_arr: xr.DataArray) -> None: + assert "y" in data_arr.dims + assert "x" in data_arr.dims + + # we remove any time dimension information + assert "t" not in data_arr.coords + assert "t" not in data_arr.dims + assert "time" not in data_arr.coords + assert "time" not in data_arr.dims + + @pytest.mark.parametrize( ("channel", "suffix"), [ @@ -302,46 +352,15 @@ def test_get_dataset(self, c01_data_arr): } res = c01_data_arr - assert "area" in res.attrs + _get_and_check_array(res, np.float32) + _check_area(res) + _check_dims_and_coords(res) for exp_key, exp_val in exp.items(): assert res.attrs[exp_key] == exp_val - # we remove any time dimension information - assert "t" not in res.coords - assert "t" not in res.dims - assert "time" not in res.coords - assert "time" not in res.dims - - def test_get_area_def(self, c01_data_arr): - """Test the area generation.""" - from pyresample.geometry import AreaDefinition - - area_def = c01_data_arr.attrs["area"] - assert isinstance(area_def, AreaDefinition) - - with ignore_pyproj_proj_warnings(): - proj_dict = area_def.crs.to_dict() - exp_dict = { - "h": 1.0, - "lon_0": -90.0, - "proj": "geos", - "sweep": "x", - "units": "m", - } - if "R" in proj_dict: - assert proj_dict["R"] == 1 - else: - assert proj_dict["a"] == 1 - assert proj_dict["b"] == 1 - for proj_key, proj_val in exp_dict.items(): - assert proj_dict[proj_key] == proj_val - - assert area_def.shape == c01_data_arr.shape - assert area_def.area_extent == (-2.0, -2998.0, 4998.0, 2.0) - @pytest.mark.parametrize("clip_negative_radiances", [False, True]) -def test_ir_calibrate(self, c07_bt_creator, clip_negative_radiances): +def test_ir_calibrate(c07_bt_creator, clip_negative_radiances): """Test IR calibration.""" res = c07_bt_creator(clip_negative_radiances=clip_negative_radiances) clipped_ir = 134.68753 if clip_negative_radiances else np.nan @@ -359,7 +378,9 @@ def test_ir_calibrate(self, c07_bt_creator, clip_negative_radiances): np.nan, ] ) - data_np = res.data.compute() + data_np = _get_and_check_array(res, np.float32) + _check_area(res) + _check_dims_and_coords(res) np.testing.assert_allclose( data_np[0, :10], expected, equal_nan=True, atol=1e-04 ) @@ -388,7 +409,9 @@ def test_vis_calibrate(c01_refl): np.nan, ] ) - data_np = res.data.compute() + data_np = _get_and_check_array(res, np.float32) + _check_area(res) + _check_dims_and_coords(res) np.testing.assert_allclose(data_np[0, :10], expected, equal_nan=True) assert "scale_factor" not in res.attrs assert "_FillValue" not in res.attrs @@ -401,8 +424,9 @@ def test_raw_calibrate(c01_counts): res = c01_counts # We expect the raw data to be unchanged - expected = res.data - assert np.allclose(res.data, expected, equal_nan=True) + _get_and_check_array(res, np.int16) + _check_area(res) + _check_dims_and_coords(res) # check for the presence of typical attributes assert "scale_factor" in res.attrs @@ -413,7 +437,6 @@ def test_raw_calibrate(c01_counts): assert "scene_id" in res.attrs # determine if things match their expected values/types. - assert res.data.dtype == np.int16 assert res.attrs["standard_name"] == "counts" assert res.attrs["long_name"] == "Raw Counts" From 14f59c49e4b327c349f510dd97b73e046eb14f72 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 31 Oct 2023 15:19:58 -0500 Subject: [PATCH 14/16] Add dask chunk size checks to ABI l1b tests --- satpy/readers/abi_base.py | 42 ++++++++++---- satpy/tests/reader_tests/test_abi_l1b.py | 70 ++++++++++++++++-------- 2 files changed, 79 insertions(+), 33 deletions(-) diff --git a/satpy/readers/abi_base.py b/satpy/readers/abi_base.py index 1c0824ab27..07a29e3043 100644 --- a/satpy/readers/abi_base.py +++ b/satpy/readers/abi_base.py @@ -18,9 +18,11 @@ """Advance Baseline Imager reader base class for the Level 1b and l2+ reader.""" import logging +import math from contextlib import suppress from datetime import datetime +import dask import numpy as np import xarray as xr from pyresample import geometry @@ -28,11 +30,10 @@ from satpy._compat import cached_property from satpy.readers import open_file_or_filename from satpy.readers.file_handlers import BaseFileHandler -from satpy.utils import get_legacy_chunk_size +from satpy.utils import get_dask_chunk_size_in_bytes logger = logging.getLogger(__name__) -CHUNK_SIZE = get_legacy_chunk_size() PLATFORM_NAMES = { "g16": "GOES-16", "g17": "GOES-17", @@ -62,15 +63,8 @@ def __init__(self, filename, filename_info, filetype_info): @cached_property def nc(self): """Get the xarray dataset for this file.""" - import math - - from satpy.utils import get_dask_chunk_size_in_bytes - chunk_size_for_high_res = math.sqrt(get_dask_chunk_size_in_bytes() / 4) # 32-bit floats - chunk_size_for_high_res = np.round(max(chunk_size_for_high_res / (4 * 226), 1)) * (4 * 226) - low_res_factor = int(self.filetype_info.get("resolution", 2000) // 500) - res_chunk_bytes = int(chunk_size_for_high_res / low_res_factor) * 4 - import dask - with dask.config.set({"array.chunk-size": res_chunk_bytes}): + chunk_bytes = self._chunk_bytes_for_resolution() + with dask.config.set({"array.chunk-size": chunk_bytes}): f_obj = open_file_or_filename(self.filename) nc = xr.open_dataset(f_obj, decode_cf=True, @@ -79,6 +73,32 @@ def nc(self): nc = self._rename_dims(nc) return nc + def _chunk_bytes_for_resolution(self) -> int: + """Get a best-guess optimal chunk size for resolution-based chunking. + + First a chunk size is chosen for the provided Dask setting `array.chunk-size` + and then aligned with a hardcoded on-disk chunk size of 226. This is then + adjusted to match the current resolution. + + This should result in 500 meter data having 4 times as many pixels per + dask array chunk (2 in each dimension) as 1km data and 8 times as many + as 2km data. As data is combined or upsampled geographically the arrays + should not need to be rechunked. Care is taken to make sure that array + chunks are aligned with on-disk file chunks at all resolutions, but at + the cost of flexibility due to a hardcoded on-disk chunk size of 226 + elements per dimension. + + """ + num_high_res_elems_per_dim = math.sqrt(get_dask_chunk_size_in_bytes() / 4) # 32-bit floats + # assume on-disk chunk size of 226 + # this is true for all CSPP Geo GRB output (226 for all sectors) and full disk from other sources + # 250 has been seen for AWS/CLASS CONUS, Mesoscale 1, and Mesoscale 2 files + # we align this with 4 on-disk chunks at 500m, so it will be 2 on-disk chunks for 1km, and 1 for 2km + high_res_elems_disk_aligned = np.round(max(num_high_res_elems_per_dim / (4 * 226), 1)) * (4 * 226) + low_res_factor = int(self.filetype_info.get("resolution", 2000) // 500) + res_elems_per_dim = int(high_res_elems_disk_aligned / low_res_factor) + return (res_elems_per_dim ** 2) * 4 + @staticmethod def _rename_dims(nc): if "t" in nc.dims or "t" in nc.coords: diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index 61f1746bf9..ec3a0334cc 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -23,6 +23,7 @@ from typing import Any, Callable from unittest import mock +import dask import dask.array as da import numpy as np import numpy.typing as npt @@ -36,9 +37,12 @@ RAD_SHAPE = { 500: (3000, 5000), # conus - 500m - 1000: (1500, 2500), # conus - 1km - 2000: (750, 1250), # conus - 2km } +# RAD_SHAPE = { +# 500: (21696, 21696), # fldk - 500m +# } +RAD_SHAPE[1000] = (RAD_SHAPE[500][0] // 2, RAD_SHAPE[500][1] // 2) +RAD_SHAPE[2000] = (RAD_SHAPE[500][0] // 4, RAD_SHAPE[500][1] // 4) def _create_fake_rad_dataarray( @@ -54,7 +58,7 @@ def _create_fake_rad_dataarray( rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( - da.from_array(rad_data), + da.from_array(rad_data, chunks=226), dims=("y", "x"), attrs={ "scale_factor": 0.5, @@ -134,15 +138,21 @@ def generate_l1b_filename(chan_name: str) -> str: @pytest.fixture() def c01_refl(tmp_path) -> xr.DataArray: - scn = _create_scene_for_data(tmp_path, "C01", None, 1000) - scn.load(["C01"]) + # 4 bytes for 32-bit floats + # 4 on-disk chunks for 500 meter data + # 226 on-disk chunk size + # Square (**2) for 2D size + with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): + scn = _create_scene_for_data(tmp_path, "C01", None, 1000) + scn.load(["C01"]) return scn["C01"] @pytest.fixture() def c01_rad(tmp_path) -> xr.DataArray: - scn = _create_scene_for_data(tmp_path, "C01", None, 1000) - scn.load([DataQuery(name="C01", calibration="radiance")]) + with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): + scn = _create_scene_for_data(tmp_path, "C01", None, 1000) + scn.load([DataQuery(name="C01", calibration="radiance")]) return scn["C01"] @@ -153,7 +163,7 @@ def c01_rad_h5netcdf(tmp_path) -> xr.DataArray: rad_data = (rad_data + 1.0) / 0.5 rad_data = rad_data.astype(np.int16) rad = xr.DataArray( - da.from_array(rad_data), + da.from_array(rad_data, chunks=226), dims=("y", "x"), attrs={ "scale_factor": 0.5, @@ -163,15 +173,17 @@ def c01_rad_h5netcdf(tmp_path) -> xr.DataArray: "valid_range": (0, 4095), }, ) - scn = _create_scene_for_data(tmp_path, "C01", rad, 1000) - scn.load([DataQuery(name="C01", calibration="radiance")]) + with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): + scn = _create_scene_for_data(tmp_path, "C01", rad, 1000) + scn.load([DataQuery(name="C01", calibration="radiance")]) return scn["C01"] @pytest.fixture() def c01_counts(tmp_path) -> xr.DataArray: - scn = _create_scene_for_data(tmp_path, "C01", None, 1000) - scn.load([DataQuery(name="C01", calibration="counts")]) + with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): + scn = _create_scene_for_data(tmp_path, "C01", None, 1000) + scn.load([DataQuery(name="C01", calibration="counts")]) return scn["C01"] @@ -181,14 +193,15 @@ def _load_data_array( clip_negative_radiances: bool = False, ): rad = _fake_c07_data() - scn = _create_scene_for_data( - tmp_path, - "C07", - rad, - 2000, - {"clip_negative_radiances": clip_negative_radiances}, - ) - scn.load(["C07"]) + with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): + scn = _create_scene_for_data( + tmp_path, + "C07", + rad, + 2000, + {"clip_negative_radiances": clip_negative_radiances}, + ) + scn.load(["C07"]) return scn["C07"] return _load_data_array @@ -202,7 +215,7 @@ def _fake_c07_data() -> xr.DataArray: rad_data = (rad_data + 1.3) / 0.5 data = rad_data.astype(np.int16) rad = xr.DataArray( - da.from_array(data), + da.from_array(data, chunks=226), dims=("y", "x"), attrs={ "scale_factor": 0.5, @@ -225,7 +238,12 @@ def _create_scene_for_data( filename = generate_l1b_filename(channel_name) data_path = tmp_path / filename dataset = _create_fake_rad_dataset(rad=rad, resolution=resolution) - dataset.to_netcdf(data_path) + dataset.to_netcdf( + data_path, + encoding={ + "Rad": {"chunksizes": [226, 226]}, + }, + ) scn = Scene( reader="abi_l1b", filenames=[str(data_path)], @@ -236,10 +254,18 @@ def _create_scene_for_data( def _get_and_check_array(data_arr: xr.DataArray, exp_dtype: npt.DTypeLike) -> npt.NDArray: data_np = data_arr.data.compute() + assert isinstance(data_arr, xr.DataArray) + assert isinstance(data_arr.data, da.Array) + assert isinstance(data_np, np.ndarray) + res = 1000 if RAD_SHAPE[1000][0] == data_np.shape[0] else 2000 + assert data_arr.chunks[0][0] == 226 * (4 / (res / 500)) + assert data_arr.chunks[1][0] == 226 * (4 / (res / 500)) + assert data_np.dtype == data_arr.dtype assert data_np.dtype == exp_dtype return data_np + def _check_area(data_arr: xr.DataArray) -> None: from pyresample.geometry import AreaDefinition From edd0632df09e0456b3738b8da3728723de2a33fe Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 1 Nov 2023 15:23:40 -0500 Subject: [PATCH 15/16] Remove unnecessary float cast in satpy/readers/abi_l1b.py --- satpy/readers/abi_l1b.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/abi_l1b.py b/satpy/readers/abi_l1b.py index 4933b0982a..29ed6f668c 100644 --- a/satpy/readers/abi_l1b.py +++ b/satpy/readers/abi_l1b.py @@ -133,7 +133,7 @@ def _raw_calibrate(self, data): def _vis_calibrate(self, data): """Calibrate visible channels to reflectance.""" solar_irradiance = self["esun"] - esd = self["earth_sun_distance_anomaly_in_AU"].astype(np.float32) + esd = self["earth_sun_distance_anomaly_in_AU"] factor = np.pi * esd * esd / solar_irradiance From 006136ef5a8590e6c7d63b1c2db97f5a70693726 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 2 Nov 2023 14:00:42 -0500 Subject: [PATCH 16/16] Switch abi l1b reader tests to use reader-level interfaces --- satpy/tests/reader_tests/test_abi_l1b.py | 43 +++++++++--------------- 1 file changed, 16 insertions(+), 27 deletions(-) diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index ec3a0334cc..a6acd7f027 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -31,16 +31,14 @@ import xarray as xr from pytest_lazyfixture import lazy_fixture -from satpy import DataQuery, Scene +from satpy import DataQuery from satpy.readers.abi_l1b import NC_ABI_L1B +from satpy.readers.yaml_reader import FileYAMLReader from satpy.utils import ignore_pyproj_proj_warnings RAD_SHAPE = { 500: (3000, 5000), # conus - 500m } -# RAD_SHAPE = { -# 500: (21696, 21696), # fldk - 500m -# } RAD_SHAPE[1000] = (RAD_SHAPE[500][0] // 2, RAD_SHAPE[500][1] // 2) RAD_SHAPE[2000] = (RAD_SHAPE[500][0] // 4, RAD_SHAPE[500][1] // 4) @@ -143,17 +141,15 @@ def c01_refl(tmp_path) -> xr.DataArray: # 226 on-disk chunk size # Square (**2) for 2D size with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): - scn = _create_scene_for_data(tmp_path, "C01", None, 1000) - scn.load(["C01"]) - return scn["C01"] + reader = _create_reader_for_data(tmp_path, "C01", None, 1000) + return reader.load(["C01"])["C01"] @pytest.fixture() def c01_rad(tmp_path) -> xr.DataArray: with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): - scn = _create_scene_for_data(tmp_path, "C01", None, 1000) - scn.load([DataQuery(name="C01", calibration="radiance")]) - return scn["C01"] + reader = _create_reader_for_data(tmp_path, "C01", None, 1000) + return reader.load([DataQuery(name="C01", calibration="radiance")])["C01"] @pytest.fixture() @@ -174,17 +170,15 @@ def c01_rad_h5netcdf(tmp_path) -> xr.DataArray: }, ) with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): - scn = _create_scene_for_data(tmp_path, "C01", rad, 1000) - scn.load([DataQuery(name="C01", calibration="radiance")]) - return scn["C01"] + reader = _create_reader_for_data(tmp_path, "C01", rad, 1000) + return reader.load([DataQuery(name="C01", calibration="radiance")])["C01"] @pytest.fixture() def c01_counts(tmp_path) -> xr.DataArray: with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): - scn = _create_scene_for_data(tmp_path, "C01", None, 1000) - scn.load([DataQuery(name="C01", calibration="counts")]) - return scn["C01"] + reader = _create_reader_for_data(tmp_path, "C01", None, 1000) + return reader.load([DataQuery(name="C01", calibration="counts")])["C01"] @pytest.fixture() @@ -194,15 +188,14 @@ def _load_data_array( ): rad = _fake_c07_data() with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): - scn = _create_scene_for_data( + reader = _create_reader_for_data( tmp_path, "C07", rad, 2000, {"clip_negative_radiances": clip_negative_radiances}, ) - scn.load(["C07"]) - return scn["C07"] + return reader.load(["C07"])["C07"] return _load_data_array @@ -228,13 +221,13 @@ def _fake_c07_data() -> xr.DataArray: return rad -def _create_scene_for_data( +def _create_reader_for_data( tmp_path: Path, channel_name: str, rad: xr.DataArray | None, resolution: int, reader_kwargs: dict[str, Any] | None = None, -) -> Scene: +) -> FileYAMLReader: filename = generate_l1b_filename(channel_name) data_path = tmp_path / filename dataset = _create_fake_rad_dataset(rad=rad, resolution=resolution) @@ -244,12 +237,8 @@ def _create_scene_for_data( "Rad": {"chunksizes": [226, 226]}, }, ) - scn = Scene( - reader="abi_l1b", - filenames=[str(data_path)], - reader_kwargs=reader_kwargs, - ) - return scn + from satpy.readers import load_readers + return load_readers([str(data_path)], "abi_l1b", reader_kwargs=reader_kwargs)["abi_l1b"] def _get_and_check_array(data_arr: xr.DataArray, exp_dtype: npt.DTypeLike) -> npt.NDArray: