diff --git a/doc/source/config.rst b/doc/source/config.rst index 9babc1abbf..1cbbbec2ed 100644 --- a/doc/source/config.rst +++ b/doc/source/config.rst @@ -272,7 +272,7 @@ If ``clip_negative_radiances=False``, pixels with negative radiances will have Clipping of negative radiances is currently implemented for the following readers: -* ``abi_l1b``, ``ami_l1b`` +* ``abi_l1b``, ``ami_l1b``, ``fci_l1c_nc`` Temporary Directory diff --git a/satpy/readers/fci_l1c_nc.py b/satpy/readers/fci_l1c_nc.py index fc40916699..c24c9b4849 100644 --- a/satpy/readers/fci_l1c_nc.py +++ b/satpy/readers/fci_l1c_nc.py @@ -208,7 +208,8 @@ class using the :mod:`~satpy.Scene.load` method with the reader "MTI3": "MTG-I3", "MTI4": "MTG-I4"} - def __init__(self, filename, filename_info, filetype_info): + def __init__(self, filename, filename_info, filetype_info, + clip_negative_radiances=False, **kwargs): """Initialize file handler.""" super().__init__(filename, filename_info, filetype_info, @@ -233,6 +234,7 @@ def __init__(self, filename, filename_info, filetype_info): else: self.is_iqt = False + self.clip_negative_radiances = clip_negative_radiances self._cache = {} @property @@ -661,6 +663,8 @@ def calibrate_counts_to_physical_quantity(self, data, key): def calibrate_counts_to_rad(self, data, key): """Calibrate counts to radiances.""" + if self.clip_negative_radiances: + data = self._clipneg(data) if key["name"] == "ir_38": data = xr.where(((2 ** 12 - 1 < data) & (data <= 2 ** 13 - 1)), (data * data.attrs.get("warm_scale_factor", 1) + @@ -677,6 +681,12 @@ def calibrate_counts_to_rad(self, data, key): self.get_and_cache_npxr(measured + "/radiance_unit_conversion_coefficient")}) return data + @staticmethod + def _clipneg(data): + """Clip counts to avoid negative radiances.""" + lo = -data.attrs.get("add_offset", 0) // data.attrs.get("scale_factor", 1) + 1 + return data.where((~data.notnull())|(data>=lo), lo) + def calibrate_rad_to_bt(self, radiance, key): """IR channel calibration.""" # using the method from PUG section Converting from Effective Radiance to Brightness Temperature for IR Channels diff --git a/satpy/tests/behave/features/image_comparison.feature b/satpy/tests/behave/features/image_comparison.feature index 3352db70b4..0497a96c93 100755 --- a/satpy/tests/behave/features/image_comparison.feature +++ b/satpy/tests/behave/features/image_comparison.feature @@ -1,13 +1,15 @@ Feature: Image Comparison Scenario Outline: Compare generated image with reference image - Given I have a reference image file from - When I generate a new image file from + Given I have a reference image file from resampled to + When I generate a new image file from with for with clipping Then the generated image should be the same as the reference image Examples: - |satellite |composite | - |GOES17 |airmass | - |GOES16 |airmass | - |GOES16 |ash | - |GOES17 |ash | + |satellite |composite | reader | area | clip | + |Meteosat-12 | cloudtop | fci_l1c_nc | sve | True | + |Meteosat-12 | night_microphysics | fci_l1c_nc | sve | True | + |GOES17 |airmass | abi_l1b | null | null | + |GOES16 |airmass | abi_l1b | null | null | + |GOES16 |ash | abi_l1b | null | null | + |GOES17 |ash | abi_l1b | null | null | diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py index 17a17cdeeb..92c5fa0034 100644 --- a/satpy/tests/behave/features/steps/image_comparison.py +++ b/satpy/tests/behave/features/steps/image_comparison.py @@ -15,7 +15,9 @@ # satpy. If not, see . """Image comparison tests.""" +import hdf5plugin # noqa: F401 isort:skip import os +import os.path import warnings from datetime import datetime from glob import glob @@ -51,17 +53,18 @@ def setup_hooks(): use_fixture(before_all, Context) setup_hooks() -@given("I have a {composite} reference image file from {satellite}") -def step_given_reference_image(context, composite, satellite): +@given("I have a {composite} reference image file from {satellite} resampled to {area}") +def step_given_reference_image(context, composite, satellite, area): """Prepare a reference image.""" - reference_image = f"reference_image_{satellite}_{composite}.png" + reference_image = f"satpy-reference-image-{satellite}-{composite}-{area}.png" context.reference_image = cv2.imread(f"{ext_data_path}/reference_images/{reference_image}") context.satellite = satellite context.composite = composite + context.area = area -@when("I generate a new {composite} image file from {satellite}") -def step_when_generate_image(context, composite, satellite): +@when("I generate a new {composite} image file from {satellite} with {reader} for {area} with clipping {clip}") +def step_when_generate_image(context, composite, satellite, reader, area, clip): """Generate test images.""" os.environ["OMP_NUM_THREADS"] = os.environ["MKL_NUM_THREADS"] = "2" os.environ["PYTROLL_CHUNK_SIZE"] = "1024" @@ -71,14 +74,22 @@ def step_when_generate_image(context, composite, satellite): # Get the list of satellite files to open filenames = glob(f"{ext_data_path}/satellite_data/{satellite}/*.nc") - scn = Scene(reader="abi_l1b", filenames=filenames) + reader_kwargs = {} + if clip != "null": + reader_kwargs["clip_negative_radiances"] = clip + scn = Scene(reader=reader, filenames=filenames, reader_kwargs=reader_kwargs) scn.load([composite]) + if area == "null": + ls = scn + else: + ls = scn.resample(area, resampler="gradient_search") + # Save the generated image in the generated folder generated_image_path = os.path.join(context.test_results_dir, "generated", - f"generated_{context.satellite}_{context.composite}.png") - scn.save_datasets(writer="simple_image", filename=generated_image_path) + f"generated_{context.satellite}_{context.composite}_{context.area}.png") + ls.save_datasets(writer="simple_image", filename=generated_image_path) # Save the generated image in the context context.generated_image = cv2.imread(generated_image_path) diff --git a/satpy/tests/reader_tests/test_fci_l1c_nc.py b/satpy/tests/reader_tests/test_fci_l1c_nc.py index aa98990df3..718d51c819 100644 --- a/satpy/tests/reader_tests/test_fci_l1c_nc.py +++ b/satpy/tests/reader_tests/test_fci_l1c_nc.py @@ -121,6 +121,7 @@ DICT_CALIBRATION = {"radiance": {"dtype": np.float32, "value_1": 15, "value_0": 9700, + "value_2": -5, "attrs_dict": {"calibration": "radiance", "units": "mW m-2 sr-1 (cm-1)-1", "radiance_unit_conversion_coefficient": np.float32(1234.56) @@ -134,8 +135,9 @@ }, "counts": {"dtype": np.uint16, - "value_1": 1, + "value_1": 5, "value_0": 5000, + "value_2": 1, "attrs_dict": {"calibration": "counts", "units": "count", }, @@ -144,6 +146,7 @@ "brightness_temperature": {"dtype": np.float32, "value_1": np.float32(209.68275), "value_0": np.float32(1888.8513), + "value_2": np.float32("nan"), "attrs_dict": {"calibration": "brightness_temperature", "units": "K", }, @@ -293,16 +296,17 @@ def _get_test_image_data_for_channel(data, ch_str, n_rows_cols): common_attrs = { "scale_factor": 5, - "add_offset": 10, + "add_offset": -10, "long_name": "Effective Radiance", "units": "mW.m-2.sr-1.(cm-1)-1", "ancillary_variables": "pixel_quality" } if "38" in ch_path: fire_line = da.ones((1, n_rows_cols[1]), dtype="uint16", chunks=1024) * 5000 - data_without_fires = da.ones((n_rows_cols[0] - 1, n_rows_cols[1]), dtype="uint16", chunks=1024) + data_without_fires = da.full((n_rows_cols[0] - 2, n_rows_cols[1]), 5, dtype="uint16", chunks=1024) + neg_rad = da.ones((1, n_rows_cols[1]), dtype="uint16", chunks=1024) d = FakeH5Variable( - da.concatenate([fire_line, data_without_fires], axis=0), + da.concatenate([fire_line, data_without_fires, neg_rad], axis=0), dims=("y", "x"), attrs={ "valid_range": [0, 8191], @@ -313,7 +317,7 @@ def _get_test_image_data_for_channel(data, ch_str, n_rows_cols): ) else: d = FakeH5Variable( - da.ones(n_rows_cols, dtype="uint16", chunks=1024), + da.full(n_rows_cols, 5, dtype="uint16", chunks=1024), dims=("y", "x"), attrs={ "valid_range": [0, 4095], @@ -440,9 +444,14 @@ class FakeFCIFileHandlerBase(FakeNetCDF4FileHandler): """Class for faking the NetCDF4 Filehandler.""" cached_file_content: Dict[str, xr.DataArray] = {} - # overwritten by FDHSI and HRFI FIle Handlers + # overwritten by FDHSI and HRFI File Handlers chan_patterns: Dict[str, Dict[str, Union[List[int], str]]] = {} + def __init__(self, *args, **kwargs): + """Initiative fake file handler.""" + kwargs.pop("clip_negative_radiances", None) + super().__init__(*args, **kwargs) + def _get_test_content_all_channels(self): data = {} for pat in self.chan_patterns: @@ -542,11 +551,11 @@ def reader_configs(): os.path.join("readers", "fci_l1c_nc.yaml")) -def _get_reader_with_filehandlers(filenames, reader_configs): +def _get_reader_with_filehandlers(filenames, reader_configs, **reader_kwargs): from satpy.readers import load_reader reader = load_reader(reader_configs) loadables = reader.select_files_from_pathnames(filenames) - reader.create_filehandlers(loadables) + reader.create_filehandlers(loadables, fh_kwargs=reader_kwargs) clear_cache(reader) return reader @@ -738,7 +747,8 @@ def _reflectance_test(tab, filenames): def _other_calibration_test(res, ch, dict_arg): """Test of other calibration test.""" if ch == "ir_38": - numpy.testing.assert_array_equal(res[ch][-1], dict_arg["value_1"]) + numpy.testing.assert_array_equal(res[ch][-1], dict_arg["value_2"]) + numpy.testing.assert_array_equal(res[ch][-2], dict_arg["value_1"]) numpy.testing.assert_array_equal(res[ch][0], dict_arg["value_0"]) else: numpy.testing.assert_array_equal(res[ch], dict_arg["value_1"]) @@ -860,6 +870,19 @@ def test_load_calibration(self, reader_configs, fh_param, self._get_assert_load(res, ch, DICT_CALIBRATION[calibration], fh_param["filenames"][0]) + @pytest.mark.parametrize("fh_param", [lazy_fixture("FakeFCIFileHandlerFDHSI_fixture")]) + def test_load_calibration_negative_rad(self, reader_configs, fh_param): + """Test calibrating negative radiances. + + See https://github.com/pytroll/satpy/issues/3009. + """ + reader = _get_reader_with_filehandlers(fh_param["filenames"], + reader_configs, + clip_negative_radiances=True) + res = reader.load([make_dataid(name="ir_38", calibration="radiance")], + pad_data=False) + numpy.testing.assert_array_equal(res["ir_38"][-1, :], 5) # smallest positive radiance + @pytest.mark.parametrize(("calibration", "channel", "resolution"), [ (calibration, channel, resolution) for calibration in ["counts", "radiance", "brightness_temperature", "reflectance"] diff --git a/utils/create_reference.py b/utils/create_reference.py index 8604cb45c0..5510054099 100644 --- a/utils/create_reference.py +++ b/utils/create_reference.py @@ -18,34 +18,93 @@ Script to create reference images for the automated image testing system. -create_reference.py - The input data directory must follow the data structure from the image-comparison-tests repository with satellite_data/. This script is a work in progress and expected to change significantly. -It is absolutely not intended for any operational production of satellite -imagery. + +DO NOT USE FOR OPERATIONAL PRODUCTION! """ -import sys -from glob import glob +import argparse +import os +import pathlib -from dask.diagnostics import ProgressBar +import hdf5plugin # noqa: F401 from satpy import Scene -ext_data_path = sys.argv[1] -outdir = sys.argv[2] -satellite = sys.argv[3] -filenames = glob(f"{ext_data_path}/satellite_data/{satellite}/*.nc") +def generate_images(props): + """Generate reference images for testing purposes. + + Args: + props (namespace): Object with attributes corresponding to command line + arguments as defined by :func:get_parser. + """ + filenames = (props.basedir / "satellite_data" / props.satellite).glob("*") + + scn = Scene(reader=props.reader, filenames=filenames) + + scn.load(props.composites) + if props.area == "native": + ls = scn.resample(resampler="native") + elif props.area is not None: + ls = scn.resample(props.area, resampler="gradient_search") + else: + ls = scn + + from dask.diagnostics import ProgressBar + with ProgressBar(): + ls.save_datasets( + writer="simple_image", + filename=os.fspath( + props.basedir / "reference_images" / + "satpy-reference-image-{platform_name}-{sensor}-" + "{start_time:%Y%m%d%H%M}-{area.area_id}-{name}.png")) + +def get_parser(): + """Return argument parser.""" + parser = argparse.ArgumentParser(description=__doc__) + + parser.add_argument( + "satellite", action="store", type=str, + help="Satellite name.") + + parser.add_argument( + "reader", action="store", type=str, + help="Reader name.") + + parser.add_argument( + "-b", "--basedir", action="store", type=pathlib.Path, + default=pathlib.Path("."), + help="Base directory for reference data. " + "This must contain a subdirectories satellite_data and " + "reference_images. The directory satellite_data must contain " + "input data in a subdirectory for the satellite. Output images " + "will be written to the subdirectory reference_images.") + + parser.add_argument( + "-o", "--outdir", action="store", type=pathlib.Path, + default=pathlib.Path("."), + help="Directory where to write resulting images.") + + parser.add_argument( + "-c", "--composites", nargs="+", help="composites to generate", + type=str, default=["ash", "airmass"]) + + parser.add_argument( + "-a", "--area", action="store", + default=None, + help="Area name, or 'native' (native resampling)") + + return parser + +def main(): + """Main function.""" + parsed = get_parser().parse_args() -scn = Scene(reader="abi_l1b", filenames=filenames) + generate_images(parsed) -composites = ["ash", "airmass"] -scn.load(composites) -ls = scn.resample(resampler="native") -with ProgressBar(): - ls.save_datasets(writer="simple_image", filename=outdir + - "/satpy-reference-image-{platform_name}-{sensor}-{start_time:%Y%m%d%H%M}-{area.area_id}-{name}.png") +if __name__ == "__main__": + main()