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 0b80045767..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,20 +63,42 @@ def __init__(self, filename, filename_info, filetype_info): @cached_property def nc(self): """Get the xarray dataset for this file.""" - f_obj = open_file_or_filename(self.filename) - try: + 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, 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 + 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: @@ -136,19 +159,14 @@ 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 = float(factor) - data = data * factor + offset + 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..29ed6f668c 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"] @@ -59,12 +60,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": @@ -136,11 +133,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"] 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" diff --git a/satpy/tests/reader_tests/test_abi_l1b.py b/satpy/tests/reader_tests/test_abi_l1b.py index ab2b1eec54..a6acd7f027 100644 --- a/satpy/tests/reader_tests/test_abi_l1b.py +++ b/satpy/tests/reader_tests/test_abi_l1b.py @@ -16,35 +16,55 @@ # 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 datetime import datetime +from pathlib import Path +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 import pytest import xarray as xr - -from satpy.tests.utils import make_dataid - - -def _create_fake_rad_dataarray(rad=None): - x_image = xr.DataArray(0.) - y_image = xr.DataArray(0.) - time = xr.DataArray(0.) +from pytest_lazyfixture import lazy_fixture + +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[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( + 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 = RAD_SHAPE[resolution] if rad is None: - rad_data = (np.arange(10.).reshape((2, 5)) + 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, + da.from_array(rad_data, chunks=226), 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 @@ -52,29 +72,29 @@ def _create_fake_rad_dataarray(rad=None): 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., "add_offset": -1.}, - 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., "add_offset": 1.}, - dims=("y",) + range(rad.shape[0]), + attrs={"scale_factor": -2.0, "add_offset": 1.0}, + dims=("y",), ) proj = xr.DataArray( - [], + np.int64(0), 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( @@ -83,8 +103,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), @@ -95,13 +115,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", @@ -111,317 +130,336 @@ 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.""" - - @mock.patch("satpy.readers.abi_base.xr") - def setUp(self, xr_, rad=None, clip_negative_radiances=False): - """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) - - -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 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) - 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) - loadables = reader.select_files_from_pathnames([fn2]) - assert len(loadables) == 1 - - -class Test_NC_ABI_L1B(Test_NC_ABI_L1B_Base): - """Test the NC_ABI_L1B reader.""" - - 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) - - 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"}) - 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"} - - assert res.attrs == exp - # 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 - - @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) - - 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 - np.testing.assert_allclose(call_args[6], (-2, -2, 8, 2)) - - -class Test_NC_ABI_L1B_ir_cal(Test_NC_ABI_L1B_Base): - """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 - } - ) - super(Test_NC_ABI_L1B_ir_cal, self).setUp(rad=rad) - - def test_ir_calibration_attrs(self): - """Test IR calibrated DataArray attributes.""" - res = self.reader.get_dataset( - make_dataid(name="C05", calibration="brightness_temperature"), {}) - - # 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 test_clip_negative_radiances_attribute(self): - """Assert that clip_negative_radiances is set to False.""" - assert not self.reader.clip_negative_radiances - - def test_ir_calibrate(self): - """Test IR calibration.""" - res = self.reader.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) - - def test_clip_negative_radiances_attribute(self): - """Assert that clip_negative_radiances has been set to True.""" - assert self.reader.clip_negative_radiances - - def test_ir_calibrate(self): - """Test IR calibration.""" - res = self.reader.get_dataset( - make_dataid(name="C07", calibration="brightness_temperature"), {}) - - 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 generate_l1b_filename(chan_name: str) -> str: + return f"OR_ABI-L1b-RadC-M4{chan_name}_G16_s20161811540362_e20161811545170_c20161811545230_suffix.nc" -class Test_NC_ABI_L1B_vis_cal(Test_NC_ABI_L1B_Base): - """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 - rad_data = rad_data.astype(np.int16) - rad = xr.DataArray( - rad_data, - dims=("y", "x"), - attrs={ - "scale_factor": 0.5, - "add_offset": -1., - "_FillValue": 20, - } - ) - super(Test_NC_ABI_L1B_vis_cal, self).setUp(rad=rad) +@pytest.fixture() +def c01_refl(tmp_path) -> xr.DataArray: + # 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}): + reader = _create_reader_for_data(tmp_path, "C01", None, 1000) + return reader.load(["C01"])["C01"] - def test_vis_calibrate(self): - """Test VIS calibration.""" - res = self.reader.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 - assert res.attrs["standard_name"] == "toa_bidirectional_reflectance" - assert res.attrs["long_name"] == "Bidirectional Reflectance" +@pytest.fixture() +def c01_rad(tmp_path) -> xr.DataArray: + with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): + reader = _create_reader_for_data(tmp_path, "C01", None, 1000) + return reader.load([DataQuery(name="C01", calibration="radiance")])["C01"] -class Test_NC_ABI_L1B_raw_cal(Test_NC_ABI_L1B_Base): - """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 - rad_data = rad_data.astype(np.int16) - rad = xr.DataArray( - rad_data, - dims=("y", "x"), - attrs={ - "scale_factor": 0.5, - "add_offset": -1., - "_FillValue": 20, - } - ) - super(Test_NC_ABI_L1B_raw_cal, self).setUp(rad=rad) - - def test_raw_calibrate(self): - """Test RAW calibration.""" - res = self.reader.get_dataset( - make_dataid(name="C05", calibration="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" - - -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.reader.get_dataset(did, {}) - +@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, chunks=226), + dims=("y", "x"), + 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), + }, + ) + with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): + 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}): + reader = _create_reader_for_data(tmp_path, "C01", None, 1000) + return reader.load([DataQuery(name="C01", calibration="counts")])["C01"] + + +@pytest.fixture() +def c07_bt_creator(tmp_path) -> Callable: + def _load_data_array( + clip_negative_radiances: bool = False, + ): + rad = _fake_c07_data() + with dask.config.set({"array.chunk-size": ((226 * 4) ** 2) * 4}): + reader = _create_reader_for_data( + tmp_path, + "C07", + rad, + 2000, + {"clip_negative_radiances": clip_negative_radiances}, + ) + return reader.load(["C07"])["C07"] + + return _load_data_array + + +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, chunks=226), + 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 + }, + ) + return rad -class Test_NC_ABI_File(unittest.TestCase): - """Test file opening.""" - @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 +def _create_reader_for_data( + tmp_path: Path, + channel_name: str, + rad: xr.DataArray | None, + resolution: int, + reader_kwargs: dict[str, Any] | None = None, +) -> FileYAMLReader: + 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, + encoding={ + "Rad": {"chunksizes": [226, 226]}, + }, + ) + 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: + 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 + + 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"), + [ + ("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-000000_0.nc" + ).format(channel) + loadables = reader.select_files_from_pathnames([fn2]) + assert len(loadables) == 1 - openable_thing = mock.MagicMock() - NC_ABI_L1B(openable_thing, {"platform_shortname": "g16"}, None) - openable_thing.open.assert_called() +@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.""" + def test_get_dataset(self, c01_data_arr): + """Test the get_dataset method.""" + exp = { + "calibration": "radiance", + "instrument_ID": None, + "modifiers": (), + "name": "C01", + "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, + "reader": "abi_l1b", + "resolution": 1000, + "scan_mode": "M4", + "scene_abbr": "C", + "scene_id": None, + "sensor": "abi", + "timeline_ID": None, + "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), + } -class Test_NC_ABI_L1B_H5netcdf(Test_NC_ABI_L1B): - """Allow h5netcdf peculiarities.""" + res = c01_data_arr + _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 + + +@pytest.mark.parametrize("clip_negative_radiances", [False, True]) +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 + 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 = _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 + ) - def setUp(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., - "_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) + # 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 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 = _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 + 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 + _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 + 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.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() 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).