Merge pull request #3013 from gerritholl/clip-negative-fci

Clip negative FCI radiances

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

satpy/readers/fci_l1c_nc.py

@@ -127,6 +127,7 @@
 from pyorbital.astronomy import sun_earth_distance_correction
 from pyresample import geometry
 
+import satpy
 from satpy.readers._geos_area import get_geos_area_naming
 from satpy.readers.eum_base import get_service_mode
 
@@ -208,7 +209,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=None, **kwargs):
         """Initialize file handler."""
         super().__init__(filename, filename_info,
                          filetype_info,
@@ -233,6 +235,9 @@ def __init__(self, filename, filename_info, filetype_info):
         else:
             self.is_iqt = False
 
+        if clip_negative_radiances is None:
+            clip_negative_radiances = satpy.config.get("readers.clip_negative_radiances")
+        self.clip_negative_radiances = clip_negative_radiances
         self._cache = {}
 
     @property
@@ -661,6 +666,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 +684,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

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 <composite> reference image file from <satellite>
-    When I generate a new <composite> image file from <satellite>
+    Given I have a <composite> reference image file from <satellite> resampled to <area>
+    When I generate a new <composite> image file from <satellite> with <reader> for <area> with clipping <clip>
     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 |

satpy/tests/behave/features/steps/image_comparison.py

@@ -15,7 +15,9 @@
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
 """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)

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,26 @@ 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.
+        """
+        import satpy
+        reader = _get_reader_with_filehandlers(fh_param["filenames"],
+                                               reader_configs,
+                                               clip_negative_radiances=True)
+        did = make_dataid(name="ir_38", calibration="radiance")
+        res = reader.load([did], pad_data=False)
+        with satpy.config.set({"readers.clip_negative_radiances": True}):
+            reader2 = _get_reader_with_filehandlers(fh_param["filenames"],
+                                                    reader_configs)
+            res2 = reader2.load([did], pad_data=False)
+        numpy.testing.assert_array_equal(res["ir_38"][-1, :], 5)  # smallest positive radiance
+        numpy.testing.assert_array_equal(res2["ir_38"][-1, :], 5)  # smallest positive radiance
+        assert res["ir_38"].dtype == res2["ir_38"].dtype == np.dtype("float32")
+
     @pytest.mark.parametrize(("calibration", "channel", "resolution"), [
         (calibration, channel, resolution)
         for calibration in ["counts", "radiance", "brightness_temperature", "reflectance"]

utils/create_reference.py

@@ -18,34 +18,93 @@
 
 Script to create reference images for the automated image testing system.
 
-create_reference.py <input-data> <output-directory> <satellite-name>
-
 The input data directory must follow the data structure from the
 image-comparison-tests repository with satellite_data/<satellite-name>.
 
 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()