Merge pull request #1693 from BENR0/feat_line_quality_flags_seviri_nc…

…_reader

Closes undefined

satpy/readers/seviri_base.py

@@ -106,8 +106,17 @@
 :func:`satpy.readers.utils.remove_earthsun_distance_correction`.
 
 
+Masking of bad quality scan lines
+---------------------------------
+
+By default bad quality scan lines are masked and replaced with ``np.nan`` for radiance, reflectance and
+brightness temperature calibrations based on the quality flags provided by the data (for details on quality
+flags see `MSG Level 1.5 Image Data Format Description`_ page 109). To disable masking
+``reader_kwargs={'mask_bad_quality_scan_lines': False}`` can be passed to the Scene.
+
+
 Metadata
-^^^^^^^^
+--------
 
 The SEVIRI L1.5 readers provide the following metadata:
 
@@ -922,3 +931,51 @@ def pad_data_vertically(data, final_size, south_bound, north_bound):
     padding_north = get_padding_area(((final_size[0] - north_bound), ncols), data.dtype)
 
     return np.vstack((padding_south, data, padding_north))
+
+
+def _create_bad_quality_lines_mask(line_validity, line_geometric_quality, line_radiometric_quality):
+    """Create bad quality scan lines mask.
+
+    For details on quality flags see `MSG Level 1.5 Image Data Format Description`_
+    page 109.
+
+    Args:
+        line_validity (numpy.ndarray):
+            Quality flags with shape (nlines,).
+        line_geometric_quality (numpy.ndarray):
+            Quality flags with shape (nlines,).
+        line_radiometric_quality (numpy.ndarray):
+            Quality flags with shape (nlines,).
+
+    Returns:
+        numpy.ndarray
+    """
+    # Based on missing (2) or corrupted (3) data
+    line_mask = line_validity >= 2
+    line_mask &= line_validity <= 3
+    # Do not use (4)
+    line_mask &= line_radiometric_quality == 4
+    line_mask &= line_geometric_quality == 4
+    return line_mask
+
+
+def mask_bad_quality(data, line_validity, line_geometric_quality, line_radiometric_quality):
+    """Mask scan lines with bad quality.
+
+    Args:
+        data (xarray.DataArray):
+            Channel data
+        line_validity (numpy.ndarray):
+            Quality flags with shape (nlines,).
+        line_geometric_quality (numpy.ndarray):
+            Quality flags with shape (nlines,).
+        line_radiometric_quality (numpy.ndarray):
+            Quality flags with shape (nlines,).
+
+    Returns:
+        xarray.DataArray: data with lines flagged as bad converted to np.nan.
+    """
+    line_mask = _create_bad_quality_lines_mask(line_validity, line_geometric_quality, line_radiometric_quality)
+    line_mask = line_mask[:, np.newaxis]
+    data = data.where(line_mask, np.nan).astype(np.float32)
+    return data

satpy/readers/seviri_l1b_hrit.py

@@ -217,6 +217,7 @@
     create_coef_dict,
     get_cds_time,
     get_satpos,
+    mask_bad_quality,
     pad_data_horizontally,
 )
 from satpy.readers.seviri_l1b_native_hdr import hrit_epilogue, hrit_prologue, impf_configuration
@@ -427,7 +428,8 @@ class HRITMSGFileHandler(HRITFileHandler):
     def __init__(self, filename, filename_info, filetype_info,
                  prologue, epilogue, calib_mode='nominal',
                  ext_calib_coefs=None, include_raw_metadata=False,
-                 mda_max_array_size=100, fill_hrv=True):
+                 mda_max_array_size=100, fill_hrv=True,
+                 mask_bad_quality_scan_lines=True):
         """Initialize the reader."""
         super(HRITMSGFileHandler, self).__init__(filename, filename_info,
                                                  filetype_info,
@@ -445,6 +447,7 @@ def __init__(self, filename, filename_info, filetype_info,
         self.fill_hrv = fill_hrv
         self.calib_mode = calib_mode
         self.ext_calib_coefs = ext_calib_coefs or {}
+        self.mask_bad_quality_scan_lines = mask_bad_quality_scan_lines
 
         self._get_header()
 
@@ -626,6 +629,11 @@ def get_dataset(self, key, info):
         """Get the dataset."""
         res = super(HRITMSGFileHandler, self).get_dataset(key, info)
         res = self.calibrate(res, key['calibration'])
+
+        is_calibration = key['calibration'] in ['radiance', 'reflectance', 'brightness_temperature']
+        if (is_calibration and self.mask_bad_quality_scan_lines):  # noqa: E129
+            res = self._mask_bad_quality(res)
+
         if key['name'] == 'HRV' and self.fill_hrv:
             res = self.pad_hrv_data(res)
         self._update_attrs(res, info)
@@ -676,20 +684,15 @@ def calibrate(self, data, calibration):
             scan_time=self.start_time
         )
         res = calib.calibrate(data, calibration)
-        if calibration in ['radiance', 'reflectance', 'brightness_temperature']:
-            res = self._mask_bad_quality(res)
         logger.debug("Calibration time " + str(datetime.now() - tic))
         return res
 
     def _mask_bad_quality(self, data):
         """Mask scanlines with bad quality."""
-        # Based on missing (2) or corrupted (3) data
-        line_mask = self.mda['image_segment_line_quality']['line_validity'] >= 2
-        line_mask &= self.mda['image_segment_line_quality']['line_validity'] <= 3
-        # Do not use (4)
-        line_mask &= self.mda['image_segment_line_quality']['line_radiometric_quality'] == 4
-        line_mask &= self.mda['image_segment_line_quality']['line_geometric_quality'] == 4
-        data *= np.choose(line_mask, [1, np.nan])[:, np.newaxis].astype(np.float32)
+        line_validity = self.mda['image_segment_line_quality']['line_validity']
+        line_radiometric_quality = self.mda['image_segment_line_quality']['line_radiometric_quality']
+        line_geometric_quality = self.mda['image_segment_line_quality']['line_geometric_quality']
+        data = mask_bad_quality(data, line_validity, line_geometric_quality, line_radiometric_quality)
         return data
 
     def _get_raw_mda(self):

satpy/readers/seviri_l1b_nc.py

@@ -36,6 +36,7 @@
     add_scanline_acq_time,
     get_cds_time,
     get_satpos,
+    mask_bad_quality,
 )
 
 logger = logging.getLogger('nc_msg')
@@ -57,10 +58,11 @@ class NCSEVIRIFileHandler(BaseFileHandler):
     """
 
     def __init__(self, filename, filename_info, filetype_info,
-                 ext_calib_coefs=None):
+                 ext_calib_coefs=None, mask_bad_quality_scan_lines=True):
         """Init the file handler."""
         super(NCSEVIRIFileHandler, self).__init__(filename, filename_info, filetype_info)
         self.ext_calib_coefs = ext_calib_coefs or {}
+        self.mask_bad_quality_scan_lines = mask_bad_quality_scan_lines
         self.mda = {}
         self.reference = datetime.datetime(1958, 1, 1)
         self.get_metadata()
@@ -136,6 +138,10 @@ def get_dataset(self, dataset_id, dataset_info):
         dataset = dataset.sel(y=slice(None, None, -1))
 
         dataset = self.calibrate(dataset, dataset_id)
+        is_calibration = dataset_id['calibration'] in ['radiance', 'reflectance', 'brightness_temperature']
+        if (is_calibration and self.mask_bad_quality_scan_lines):  # noqa: E129
+            dataset = self._mask_bad_quality(dataset, dataset_info)
+
         self._update_attrs(dataset, dataset_info)
         return dataset
 
@@ -174,6 +180,14 @@ def _get_calib_coefs(self, dataset, channel):
             'radiance_type': self.nc['planned_chan_processing'].values[band_idx]
         }
 
+    def _mask_bad_quality(self, dataset, dataset_info):
+        """Mask scanlines with bad quality."""
+        ch_number = int(dataset_info['nc_key'][2:])
+        line_validity = self.nc['channel_data_visir_data_line_validity'][:, ch_number - 1].data
+        line_geometric_quality = self.nc['channel_data_visir_data_line_geometric_quality'][:, ch_number - 1].data
+        line_radiometric_quality = self.nc['channel_data_visir_data_line_radiometric_quality'][:, ch_number - 1].data
+        return mask_bad_quality(dataset, line_validity, line_geometric_quality, line_radiometric_quality)
+
     def _update_attrs(self, dataset, dataset_info):
         """Update dataset attributes."""
         dataset.attrs.update(self.nc[dataset_info['nc_key']].attrs)

satpy/tests/reader_tests/test_seviri_l1b_hrit.py

@@ -154,6 +154,10 @@ def setUp(self):
             ncols=self.ncols,
             projection_longitude=self.projection_longitude
         )
+        self.reader.mda.update({
+            'segment_sequence_number': 18,
+            'planned_start_segment_number': 1
+        })
 
     def _get_fake_data(self):
         return xr.DataArray(
@@ -208,6 +212,34 @@ def test_get_dataset(self, calibrate, parent_get_dataset):
         info = setup.get_fake_dataset_info()
         res = self.reader.get_dataset(key, info)
 
+        # Test method calls
+        new_data = np.zeros_like(data.data).astype('float32')
+        new_data[:, :] = np.nan
+        expected = data.copy(data=new_data)
+
+        expected['acq_time'] = (
+            'y',
+            setup.get_acq_time_exp(self.start_time, self.nlines)
+        )
+        xr.testing.assert_equal(res, expected)
+        self.assert_attrs_equal(
+            res.attrs,
+            setup.get_attrs_exp(self.projection_longitude)
+        )
+
+    @mock.patch('satpy.readers.seviri_l1b_hrit.HRITFileHandler.get_dataset')
+    @mock.patch('satpy.readers.seviri_l1b_hrit.HRITMSGFileHandler.calibrate')
+    def test_get_dataset_without_masking_bad_scan_lines(self, calibrate, parent_get_dataset):
+        """Test getting the dataset."""
+        data = self._get_fake_data()
+        parent_get_dataset.return_value = mock.MagicMock()
+        calibrate.return_value = data
+
+        key = make_dataid(name='VIS006', calibration='reflectance')
+        info = setup.get_fake_dataset_info()
+        self.reader.mask_bad_quality_scan_lines = False
+        res = self.reader.get_dataset(key, info)
+
         # Test method calls
         expected = data.copy()
         expected['acq_time'] = (
@@ -382,9 +414,9 @@ def file_handler(self):
         }
         mda = {
             'image_segment_line_quality': {
-                'line_validity': np.zeros(2),
-                'line_radiometric_quality': np.zeros(2),
-                'line_geometric_quality': np.zeros(2)
+                'line_validity': np.array([3, 3]),
+                'line_radiometric_quality': np.array([3, 3]),
+                'line_geometric_quality': np.array([3, 3])
             },
         }
 
@@ -450,3 +482,21 @@ def test_calibrate(
 
         res = fh.calibrate(counts, calibration)
         xr.testing.assert_allclose(res, expected)
+
+    def test_mask_bad_quality(self, file_handler):
+        """Test the masking of bad quality scan lines."""
+        channel = 'VIS006'
+        expected = self._get_expected(
+            channel=channel,
+            calibration='radiance',
+            calib_mode='NOMINAL',
+            use_ext_coefs=False
+        )
+
+        fh = file_handler
+
+        res = fh._mask_bad_quality(expected)
+        new_data = np.zeros_like(expected.data).astype('float32')
+        new_data[:, :] = np.nan
+        expected = expected.copy(data=new_data)
+        xr.testing.assert_allclose(res, expected)

satpy/tests/reader_tests/test_seviri_l1b_hrit_setup.py

@@ -160,7 +160,10 @@ def get_fake_mda(nlines, ncols, start_time):
         'coff': 10,
         'loff': 10,
         'image_segment_line_quality': {
-            'line_mean_acquisition': tline
+            'line_mean_acquisition': tline,
+            'line_validity': np.full(nlines, 3),
+            'line_radiometric_quality': np.full(nlines, 3),
+            'line_geometric_quality': np.full(nlines, 3)
         }
     }