diff --git a/satpy/readers/seviri_base.py b/satpy/readers/seviri_base.py index d2c4a31bee..922b007988 100644 --- a/satpy/readers/seviri_base.py +++ b/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 diff --git a/satpy/readers/seviri_l1b_hrit.py b/satpy/readers/seviri_l1b_hrit.py index 7da6ee1e21..9083828d0b 100644 --- a/satpy/readers/seviri_l1b_hrit.py +++ b/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): diff --git a/satpy/readers/seviri_l1b_nc.py b/satpy/readers/seviri_l1b_nc.py index 7dd67c3296..17435fcdb6 100644 --- a/satpy/readers/seviri_l1b_nc.py +++ b/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) diff --git a/satpy/tests/reader_tests/test_seviri_l1b_hrit.py b/satpy/tests/reader_tests/test_seviri_l1b_hrit.py index 1a2258b6f5..3fc1dea902 100644 --- a/satpy/tests/reader_tests/test_seviri_l1b_hrit.py +++ b/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) diff --git a/satpy/tests/reader_tests/test_seviri_l1b_hrit_setup.py b/satpy/tests/reader_tests/test_seviri_l1b_hrit_setup.py index b89a32bd4f..eb6e4c5a32 100644 --- a/satpy/tests/reader_tests/test_seviri_l1b_hrit_setup.py +++ b/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) } } diff --git a/satpy/tests/reader_tests/test_seviri_l1b_nc.py b/satpy/tests/reader_tests/test_seviri_l1b_nc.py index 615ca49e64..6d55c1a088 100644 --- a/satpy/tests/reader_tests/test_seviri_l1b_nc.py +++ b/satpy/tests/reader_tests/test_seviri_l1b_nc.py @@ -29,6 +29,8 @@ from satpy.tests.reader_tests.test_seviri_l1b_calibration import TestFileHandlerCalibrationBase from satpy.tests.utils import assert_attrs_equal, make_dataid +channel_keys_dict = {'VIS006': 'ch1', 'IR_108': 'ch9'} + def to_cds_time(time): """Convert datetime to (days, msecs) since 1958-01-01.""" @@ -57,6 +59,7 @@ def _get_fake_dataset(self, counts, h5netcdf): """ acq_time_day = np.repeat([1, 1], 11).reshape(2, 11) acq_time_msec = np.repeat([1000, 2000], 11).reshape(2, 11) + quality = np.repeat([0, 3], 11).reshape(2, 11) orbit_poly_start_day, orbit_poly_start_msec = to_cds_time( np.array([datetime(2019, 12, 31, 18), datetime(2019, 12, 31, 22)], @@ -74,8 +77,8 @@ def _get_fake_dataset(self, counts, h5netcdf): scan_time_days, scan_time_msecs = to_cds_time(self.scan_time) ds = xr.Dataset( { - 'VIS006': counts.copy(), - 'IR_108': counts.copy(), + 'ch1': counts.copy(), + 'ch9': counts.copy(), 'HRV': (('num_rows_hrv', 'num_columns_hrv'), [[1, 2, 3], [4, 5, 6], [7, 8, 9]]), @@ -88,6 +91,18 @@ def _get_fake_dataset(self, counts, h5netcdf): ('num_rows_vis_ir', 'channels_vis_ir_dim'), acq_time_msec ), + 'channel_data_visir_data_line_validity': ( + ('num_rows_vis_ir', 'channels_vis_ir_dim'), + quality + ), + 'channel_data_visir_data_line_geometric_quality': ( + ('num_rows_vis_ir', 'channels_vis_ir_dim'), + quality + ), + 'channel_data_visir_data_line_radiometric_quality': ( + ('num_rows_vis_ir', 'channels_vis_ir_dim'), + quality + ), 'orbit_polynomial_x': ( ('orbit_polynomial_dim_row', 'orbit_polynomial_dim_col'), @@ -153,11 +168,12 @@ def _get_fake_dataset(self, counts, h5netcdf): ds.attrs.update(nattrs) - ds['VIS006'].attrs.update({ + ds['ch1'].attrs.update({ 'scale_factor': self.gains_nominal[0], 'add_offset': self.offsets_nominal[0] }) - ds['IR_108'].attrs.update({ + # IR_108 is dataset with key ch9 + ds['ch9'].attrs.update({ 'scale_factor': self.gains_nominal[8], 'add_offset': self.offsets_nominal[8], }) @@ -169,7 +185,7 @@ def _get_fake_dataset(self, counts, h5netcdf): 'valid_min': None, 'valid_max': None } - for name in ['VIS006', 'IR_108']: + for name in ['ch1', 'ch9']: ds[name].attrs.update(strip_attrs) return ds @@ -236,25 +252,56 @@ def test_calibrate( fh.ext_calib_coefs = external_coefs dataset_id = make_dataid(name=channel, calibration=calibration) - res = fh.calibrate(fh.nc[channel], dataset_id) + key = channel_keys_dict[channel] + + res = fh.calibrate(fh.nc[key], dataset_id) + xr.testing.assert_allclose(res, expected) + + def test_mask_bad_quality(self, file_handler): + """Test masking of bad quality scan lines.""" + channel = 'VIS006' + key = channel_keys_dict[channel] + dataset_info = { + 'nc_key': key, + 'units': 'units', + 'wavelength': 'wavelength', + 'standard_name': 'standard_name' + } + expected = self._get_expected( + channel=channel, + calibration='radiance', + calib_mode='NOMINAL', + use_ext_coefs=False + ) + + fh = file_handler + + res = fh._mask_bad_quality(fh.nc[key], dataset_info) + 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) @pytest.mark.parametrize( - ('channel', 'calibration'), + ('channel', 'calibration', 'mask_bad_quality_scan_lines'), [ - ('VIS006', 'reflectance'), - ('IR_108', 'brightness_temperature') + ('VIS006', 'reflectance', True), + ('VIS006', 'reflectance', False), + ('IR_108', 'brightness_temperature', True) ] ) - def test_get_dataset(self, file_handler, channel, calibration): + def test_get_dataset(self, file_handler, channel, calibration, mask_bad_quality_scan_lines): """Test getting the dataset.""" dataset_id = make_dataid(name=channel, calibration=calibration) + key = channel_keys_dict[channel] dataset_info = { - 'nc_key': channel, + 'nc_key': key, 'units': 'units', 'wavelength': 'wavelength', 'standard_name': 'standard_name' } + + file_handler.mask_bad_quality_scan_lines = mask_bad_quality_scan_lines res = file_handler.get_dataset(dataset_id, dataset_info) # Test scanline acquisition times @@ -285,6 +332,9 @@ def test_get_dataset(self, file_handler, channel, calibration): expected['acq_time'] = ('y', [np.datetime64('1958-01-02 00:00:01'), np.datetime64('1958-01-02 00:00:02')]) expected = expected[::-1] # reader flips data upside down + if mask_bad_quality_scan_lines: + expected = file_handler._mask_bad_quality(expected, dataset_info) + xr.testing.assert_allclose(res, expected) for key in ['sun_earth_distance_correction_applied', @@ -296,7 +346,8 @@ def test_satpos_no_valid_orbit_polynomial(self, file_handler): """Test satellite position if there is no valid orbit polynomial.""" dataset_id = make_dataid(name='VIS006', calibration='counts') dataset_info = { - 'nc_key': 'VIS006', + 'name': 'VIS006', + 'nc_key': 'ch1', 'units': 'units', 'wavelength': 'wavelength', 'standard_name': 'standard_name'