Add ability to clip AMI negative 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``
+* ``abi_l1b``, ``ami_l1b``
 
 
 Temporary Directory

satpy/readers/ami_l1b.py

@@ -26,6 +26,7 @@
 import xarray as xr
 from pyspectral.blackbody import blackbody_wn_rad2temp as rad2temp
 
+import satpy
 from satpy.readers import open_file_or_filename
 from satpy.readers._geos_area import get_area_definition, get_area_extent
 from satpy.readers.file_handlers import BaseFileHandler
@@ -91,29 +92,33 @@
 
     def __init__(self, filename, filename_info, filetype_info,
                  calib_mode="PYSPECTRAL", allow_conditional_pixels=False,
-                 user_calibration=None):
+                 user_calibration=None, clip_negative_radiances=None):
         """Open the NetCDF file with xarray and prepare the Dataset for reading."""
         super(AMIL1bNetCDF, self).__init__(filename, filename_info, filetype_info)
         f_obj = open_file_or_filename(self.filename)
         self.nc = xr.open_dataset(f_obj,
                                   decode_cf=True,
                                   mask_and_scale=False,
                                   chunks={"dim_image_x": CHUNK_SIZE, "dim_image_y": CHUNK_SIZE})
         self.nc = self.nc.rename({"dim_image_x": "x", "dim_image_y": "y"})
 
         platform_shortname = self.nc.attrs["satellite_name"]
         self.platform_name = PLATFORM_NAMES.get(platform_shortname)
         self.sensor = "ami"
         self.band_name = filetype_info["file_type"].upper()
         self.allow_conditional_pixels = allow_conditional_pixels
         calib_mode_choices = ("FILE", "PYSPECTRAL", "GSICS")
         if calib_mode.upper() not in calib_mode_choices:
             raise ValueError("Invalid calibration mode: {}. Choose one of {}".format(
                 calib_mode, calib_mode_choices))
 
         self.calib_mode = calib_mode.upper()
         self.user_calibration = user_calibration
 
+        if clip_negative_radiances is None:
+            clip_negative_radiances = satpy.config.get("readers.clip_negative_radiances")
+        self.clip_negative_radiances = clip_negative_radiances
+
     @property
     def start_time(self):
         """Get observation start time."""
@@ -178,6 +183,7 @@
         }
         return orbital_parameters
 
+
     def get_dataset(self, dataset_id, ds_info):
         """Load a dataset as a xarray DataArray."""
         file_key = ds_info.get("file_key", dataset_id["name"])
@@ -201,11 +207,12 @@
         data = data.where(qf == 0)
 
         # Calibration values from file, fall back to built-in if unavailable
-        gain = self.nc.attrs["DN_to_Radiance_Gain"]
-        offset = self.nc.attrs["DN_to_Radiance_Offset"]
+        self.gain = self.nc.attrs["DN_to_Radiance_Gain"]
+        self.offset = self.nc.attrs["DN_to_Radiance_Offset"]
 
         if dataset_id["calibration"] in ("radiance", "reflectance", "brightness_temperature"):
-            data = gain * data + offset
+            data = self.gain * data + self.offset
+            data = self._clip_negative_radiance(data)
             if self.calib_mode == "GSICS":
                 data = self._apply_gsics_rad_correction(data)
             elif isinstance(self.user_calibration, dict):
@@ -230,6 +237,20 @@
         data.attrs = attrs
         return data
 
+
+    def _clip_negative_radiance(self, data):
+        """If requested, clip negative radiance from Rad DataArray."""
+        if self.clip_negative_radiances:
+            count_zero_rad = - self.offset / self.gain
+            # We need floor here as the scale factor for AMI is negative (unlike ABI)
+            count_pos = np.floor(count_zero_rad)
+            min_rad = count_pos * self.gain + self.offset
+            data = data.clip(min=min_rad)
+            return data
+        else:
+            return data
+
+
     def _calibrate_ir(self, dataset_id, data):
         """Calibrate radiance data to BTs using either pyspectral or in-file coefficients."""
         if self.calib_mode == "PYSPECTRAL":

satpy/tests/reader_tests/test_ami_l1b.py

@@ -58,77 +58,87 @@
     """Common setup for NC_ABI_L1B tests."""
 
     @mock.patch("satpy.readers.ami_l1b.xr")
-    def setUp(self, xr_, counts=None):
+    def setUp(self, xr_, counts=None, irtest=False):
         """Create a fake dataset using the given counts data."""
         from satpy.readers.ami_l1b import AMIL1bNetCDF
+        if irtest:
+            dn_to_Radiance_Gain = -0.00108296517282724
+            dn_to_Radiance_Offset = 17.699987411499
+            bpp = 14
+        else:
+            dn_to_Radiance_Gain = -0.0144806550815701
+            dn_to_Radiance_Offset = 118.050903320312
+            bpp = 12
 
         if counts is None:
             rad_data = (np.arange(10.).reshape((2, 5)) + 1.) * 50.
             rad_data = (rad_data + 1.) / 0.5
             rad_data = rad_data.astype(np.uint16)
+            if irtest:
+                # If testing IR clipping, set one pixel to negative radiance
+                rad_data[0, 0] = 16364
             counts = xr.DataArray(
                 da.from_array(rad_data, chunks="auto"),
                 dims=("y", "x"),
                 attrs={
                     "channel_name": "VI006",
                     "detector_side": 2,
                     "number_of_total_pixels": 484000000,
                     "number_of_error_pixels": 113892451,
                     "max_pixel_value": 32768,
                     "min_pixel_value": 6,
                     "average_pixel_value": 8228.98770845248,
                     "stddev_pixel_value": 13621.130386551,
                     "number_of_total_bits_per_pixel": 16,
                     "number_of_data_quality_flag_bits_per_pixel": 2,
-                    "number_of_valid_bits_per_pixel": 12,
+                    "number_of_valid_bits_per_pixel": bpp,
                     "data_quality_flag_meaning":
                         "0:good_pixel, 1:conditionally_usable_pixel, 2:out_of_scan_area_pixel, 3:error_pixel",
                     "ground_sample_distance_ew": 1.4e-05,
                     "ground_sample_distance_ns": 1.4e-05,
                 }
             )
         sc_position = xr.DataArray(0., attrs={
             "sc_position_center_pixel": [-26113466.1974016, 33100139.1630508, 3943.75470244799],
         })
         xr_.open_dataset.return_value = FakeDataset(
             {
                 "image_pixel_values": counts,
                 "sc_position": sc_position,
                 "gsics_coeff_intercept": [0.1859369],
                 "gsics_coeff_slope": [0.9967594],
             },
             {
                 "satellite_name": "GK-2A",
                 "observation_start_time": 623084431.957882,
                 "observation_end_time": 623084975.606133,
                 "projection_type": "GEOS",
                 "sub_longitude": 2.23751210105673,
                 "cfac": 81701355.6133574,
                 "lfac": -81701355.6133574,
                 "coff": 11000.5,
                 "loff": 11000.5,
                 "nominal_satellite_height": 42164000.,
                 "earth_equatorial_radius": 6378137.,
                 "earth_polar_radius": 6356752.3,
                 "number_of_columns": 22000,
                 "number_of_lines": 22000,
                 "observation_mode": "FD",
                 "channel_spatial_resolution": "0.5",
                 "Radiance_to_Albedo_c": 1,
-                "DN_to_Radiance_Gain": -0.0144806550815701,
-                "DN_to_Radiance_Offset": 118.050903320312,
+                "DN_to_Radiance_Gain": dn_to_Radiance_Gain,
+                "DN_to_Radiance_Offset": dn_to_Radiance_Offset,
                 "Teff_to_Tbb_c0": -0.141418528203155,
                 "Teff_to_Tbb_c1": 1.00052232906885,
                 "Teff_to_Tbb_c2": -0.00000036287276076109,
                 "light_speed": 2.9979245800E+08,
                 "Boltzmann_constant_k": 1.3806488000E-23,
                 "Plank_constant_h": 6.6260695700E-34,
             }
         )
-
         self.reader = AMIL1bNetCDF("filename",
                                    {"platform_shortname": "gk2a"},
-                                   {"file_type": "ir087"},)
+                                   {"file_type": "ir087"})
 
 
 class TestAMIL1bNetCDF(TestAMIL1bNetCDFBase):
@@ -178,7 +188,7 @@
         assert self.reader.end_time == datetime(2019, 9, 30, 3, 9, 35, 606133)
 
     def test_get_dataset(self):
-        """Test gettting radiance data."""
+        """Test getting radiance data."""
         from satpy.tests.utils import make_dataid
         key = make_dataid(name="VI006", calibration="radiance")
         res = self.reader.get_dataset(key, {
@@ -349,3 +359,32 @@
         np.testing.assert_allclose(res.data.compute(), expected, equal_nan=True, atol=0.01)
         # make sure the attributes from the file are in the data array
         assert res.attrs["standard_name"] == "toa_brightness_temperature"
+
+
+
+class TestAMIL1bNetCDFIRClip(TestAMIL1bNetCDFBase):
+    """Test IR specific things about the AMI reader."""
+
+    def setUp(self):
+        """Create test data for IR calibration tests."""
+        from satpy.tests.utils import make_dataid
+
+        self.ds_id = make_dataid(name="IR087", wavelength=[8.415, 8.59, 8.765],
+                                 calibration="radiance")
+        self.ds_info = {
+            "file_key": "image_pixel_values",
+            "wavelength": [8.415, 8.59, 8.765],
+            "standard_name": "toa_brightness_temperature",
+            "units": "K",
+        }
+        super(TestAMIL1bNetCDFIRClip, self).setUp(irtest=True)
+
+    def test_clipneg(self):
+        """Test that negative radiances are clipped."""
+        self.reader.clip_negative_radiances = True
+        res = np.array(self.reader.get_dataset(self.ds_id, self.ds_info))
+        assert np.isclose(res[0, 0], 4.6268106e-06)
+
+        self.reader.clip_negative_radiances = False
+        res = np.array(self.reader.get_dataset(self.ds_id, self.ds_info))
+        assert res[0, 0] < 0