Harmonize calibration in SEVIRI readers

satpy/readers/seviri_base.py

@@ -332,10 +332,15 @@ def images_used(self):
 mpef_product_header = MpefProductHeader().get()
 
 
-class SEVIRICalibrationHandler(object):
-    """Calibration handler for SEVIRI HRIT- and native-formats."""
+class SEVIRICalibrationAlgorithm:
+    """SEVIRI calibration algorithms."""
 
-    def _convert_to_radiance(self, data, gain, offset):
+    def __init__(self, platform_id, scan_time):
+        """Initialize the calibration algorithm."""
+        self.platform_id = platform_id
+        self.scan_time = scan_time
+
+    def convert_to_radiance(self, data, gain, offset):
         """Calibrate to radiance."""
         return (data * gain + offset).clip(0.0, None)
 
@@ -348,7 +353,7 @@ def _erads2bt(self, data, channel_name):
 
         return (self._tl15(data, wavenumber) - beta) / alpha
 
-    def _ir_calibrate(self, data, channel_name, cal_type):
+    def ir_calibrate(self, data, channel_name, cal_type):
         """Calibrate to brightness temperature."""
         if cal_type == 1:
             # spectral radiances
@@ -372,14 +377,92 @@ def _tl15(self, data, wavenumber):
         return ((C2 * wavenumber) /
                 np.log((1.0 / data) * C1 * wavenumber ** 3 + 1.0))
 
-    def _vis_calibrate(self, data, solar_irradiance):
+    def vis_calibrate(self, data, solar_irradiance):
         """Calibrate to reflectance.
 
         This uses the method described in Conversion from radiances to
         reflectances for SEVIRI warm channels: https://tinyurl.com/y67zhphm
         """
         reflectance = np.pi * data * 100.0 / solar_irradiance
-        return apply_earthsun_distance_correction(reflectance, self.start_time)
+        return apply_earthsun_distance_correction(reflectance, self.scan_time)
+
+
+class SEVIRICalibrationHandler:
+    """Calibration handler for SEVIRI HRIT-, native- and netCDF-formats.
+
+    Handles selection of calibration coefficients and calls the appropriate
+    calibration algorithm.
+    """
+
+    def __init__(self, platform_id, channel_name, coefs, calib_mode, scan_time):
+        """Initialize the calibration handler."""
+        self.platform_id = platform_id
+        self.channel_name = channel_name
+        self.coefs = coefs
+        self.calib_mode = calib_mode.upper()
+        self.scan_time = scan_time
+        self.algo = SEVIRICalibrationAlgorithm(
+            platform_id=self.platform_id,
+            scan_time=self.scan_time
+        )
+
+        valid_modes = ('NOMINAL', 'GSICS')
+        if self.calib_mode not in valid_modes:
+            raise ValueError(
+                'Invalid calibration mode: {}. Choose one of {}'.format(
+                    self.calib_mode, valid_modes)
+            )
+
+    def calibrate(self, data, calibration):
+        """Calibrate the given data."""
+        if calibration == 'counts':
+            res = data
+        elif calibration in ['radiance', 'reflectance',
+                             'brightness_temperature']:
+            # Convert to radiance
+            gain, offset = self.get_gain_offset()
+            data = data.where(data > 0)
+            res = self.algo.convert_to_radiance(
+                data.astype(np.float32), gain, offset
+            )
+        else:
+            raise ValueError(
+                'Invalid calibration {} for channel {}'.format(
+                    calibration, self.channel_name
+                )
+            )
+
+        if calibration == 'reflectance':
+            solar_irradiance = CALIB[self.platform_id][self.channel_name]["F"]
+            res = self.algo.vis_calibrate(res, solar_irradiance)
+        elif calibration == 'brightness_temperature':
+            res = self.algo.ir_calibrate(
+                res, self.channel_name, self.coefs['radiance_type']
+            )
+
+        return res
+
+    def get_gain_offset(self):
+        """Get gain & offset for calibration from counts to radiance.
+
+        Choices for internal coefficients are nominal or GSICS. External
+        coefficients take precedence.
+        """
+        coefs = self.coefs['coefs']
+
+        # Select internal coefficients for the given calibration mode
+        if self.calib_mode != 'GSICS' or self.channel_name in VIS_CHANNELS:
+            # GSICS doesn't have calibration coeffs for VIS channels
+            internal_gain = coefs['NOMINAL']['gain']
+            internal_offset = coefs['NOMINAL']['offset']
+        else:
+            internal_gain = coefs['GSICS']['gain']
+            internal_offset = coefs['GSICS']['offset'] * internal_gain
+
+        # Override with external coefficients, if any.
+        gain = coefs['EXTERNAL'].get('gain', internal_gain)
+        offset = coefs['EXTERNAL'].get('offset', internal_offset)
+        return gain, offset
 
 
 def chebyshev(coefs, time, domain):

satpy/readers/seviri_l1b_hrit.py

@@ -167,8 +167,8 @@
 from satpy.readers.hrit_base import (HRITFileHandler, ancillary_text,
                                      annotation_header, base_hdr_map,
                                      image_data_function)
-from satpy.readers.seviri_base import (CALIB, CHANNEL_NAMES, SATNUM,
-                                       VIS_CHANNELS, SEVIRICalibrationHandler,
+from satpy.readers.seviri_base import (CHANNEL_NAMES, SATNUM,
+                                       SEVIRICalibrationHandler,
                                        chebyshev, get_cds_time)
 from satpy.readers.seviri_l1b_native_hdr import (hrit_epilogue, hrit_prologue,
                                                  impf_configuration)
@@ -433,7 +433,7 @@ def reduce(self, max_size):
         return self._reduce(self.epilogue, max_size=max_size)
 
 
-class HRITMSGFileHandler(HRITFileHandler, SEVIRICalibrationHandler):
+class HRITMSGFileHandler(HRITFileHandler):
     """SEVIRI HRIT format reader.
 
     **Calibration**
@@ -520,14 +520,10 @@ def __init__(self, filename, filename_info, filetype_info,
         self.prologue = prologue.prologue
         self.epilogue = epilogue.epilogue
         self._filename_info = filename_info
-        self.ext_calib_coefs = ext_calib_coefs if ext_calib_coefs is not None else {}
         self.mda_max_array_size = mda_max_array_size
         self.fill_hrv = fill_hrv
-        calib_mode_choices = ('NOMINAL', '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.calib_mode = calib_mode
+        self.ext_calib_coefs = ext_calib_coefs or {}
 
         self._get_header()
 
@@ -751,49 +747,28 @@ def pad_hrv_data(self, res):
     def calibrate(self, data, calibration):
         """Calibrate the data."""
         tic = datetime.now()
-        channel_name = self.channel_name
-
-        if calibration == 'counts':
-            res = data
-        elif calibration in ['radiance', 'reflectance', 'brightness_temperature']:
-            # Choose calibration coefficients
-            # a) Internal: Nominal or GSICS?
-            band_idx = self.mda['spectral_channel_id'] - 1
-            if self.calib_mode != 'GSICS' or self.channel_name in VIS_CHANNELS:
-                # you cant apply GSICS values to the VIS channels
-                coefs = self.prologue["RadiometricProcessing"]["Level15ImageCalibration"]
-                int_gain = coefs['CalSlope'][band_idx]
-                int_offset = coefs['CalOffset'][band_idx]
-            else:
-                coefs = self.prologue["RadiometricProcessing"]['MPEFCalFeedback']
-                int_gain = coefs['GSICSCalCoeff'][band_idx]
-                int_offset = coefs['GSICSOffsetCount'][band_idx] * int_gain
-
-            # b) Internal or external? External takes precedence.
-            gain = self.ext_calib_coefs.get(self.channel_name, {}).get('gain', int_gain)
-            offset = self.ext_calib_coefs.get(self.channel_name, {}).get('offset', int_offset)
-
-            # Convert to radiance
-            data = data.where(data > 0)
-            res = self._convert_to_radiance(data.astype(np.float32), gain, offset)
-            line_mask = self.mda['image_segment_line_quality']['line_validity'] >= 2
-            line_mask &= self.mda['image_segment_line_quality']['line_validity'] <= 3
-            line_mask &= self.mda['image_segment_line_quality']['line_radiometric_quality'] == 4
-            line_mask &= self.mda['image_segment_line_quality']['line_geometric_quality'] == 4
-            res *= np.choose(line_mask, [1, np.nan])[:, np.newaxis].astype(np.float32)
-
-        if calibration == 'reflectance':
-            solar_irradiance = CALIB[self.platform_id][channel_name]["F"]
-            res = self._vis_calibrate(res, solar_irradiance)
-
-        elif calibration == 'brightness_temperature':
-            cal_type = self.prologue['ImageDescription'][
-                'Level15ImageProduction']['PlannedChanProcessing'][self.mda['spectral_channel_id']]
-            res = self._ir_calibrate(res, channel_name, cal_type)
-
+        calib = SEVIRICalibrationHandler(
+            platform_id=self.platform_id,
+            channel_name=self.channel_name,
+            coefs=self._get_calib_coefs(self.channel_name),
+            calib_mode=self.calib_mode,
+            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."""
+        line_mask = self.mda['image_segment_line_quality']['line_validity'] >= 2
+        line_mask &= self.mda['image_segment_line_quality']['line_validity'] <= 3
+        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)
+        return data
+
     def _get_raw_mda(self):
         """Compile raw metadata to be included in the dataset attributes."""
         # Metadata from segment header (excluding items which vary among the different segments)
@@ -812,6 +787,30 @@ def _get_timestamps(self):
         tline = self.mda['image_segment_line_quality']['line_mean_acquisition']
         return get_cds_time(days=tline['days'], msecs=tline['milliseconds'])
 
+    def _get_calib_coefs(self, channel_name):
+        """Get coefficients for calibration from counts to radiance."""
+        band_idx = self.mda['spectral_channel_id'] - 1
+        radiance_type_idx = self.mda['spectral_channel_id']
+        coefs_nominal = self.prologue["RadiometricProcessing"][
+            "Level15ImageCalibration"]
+        coefs_gsics = self.prologue["RadiometricProcessing"]['MPEFCalFeedback']
+        radiance_types = self.prologue['ImageDescription'][
+                'Level15ImageProduction']['PlannedChanProcessing']
+        return {
+            'coefs': {
+                'NOMINAL': {
+                    'gain': coefs_nominal['CalSlope'][band_idx],
+                    'offset': coefs_nominal['CalOffset'][band_idx],
+                },
+                'GSICS': {
+                    'gain': coefs_gsics['GSICSCalCoeff'][band_idx],
+                    'offset': coefs_gsics['GSICSOffsetCount'][band_idx]
+                },
+                'EXTERNAL': self.ext_calib_coefs.get(channel_name, {})
+            },
+            'radiance_type': radiance_types[radiance_type_idx]
+        }
+
 
 def pad_data(data, final_size, east_bound, west_bound):
     """Pad the data given east and west bounds and the desired size."""

satpy/readers/seviri_l1b_native.py

@@ -56,32 +56,33 @@
 from satpy.readers.file_handlers import BaseFileHandler
 from satpy.readers.eum_base import recarray2dict
 from satpy.readers.seviri_base import (SEVIRICalibrationHandler,
-                                       CHANNEL_NAMES, CALIB, SATNUM,
+                                       CHANNEL_NAMES, SATNUM,
                                        dec10216, VISIR_NUM_COLUMNS,
-                                       VISIR_NUM_LINES, HRV_NUM_COLUMNS,
-                                       VIS_CHANNELS)
+                                       VISIR_NUM_LINES, HRV_NUM_COLUMNS)
 from satpy.readers.seviri_l1b_native_hdr import (GSDTRecords, native_header,
                                                  native_trailer)
 from satpy.readers._geos_area import get_area_definition
 
 logger = logging.getLogger('native_msg')
 
 
-class NativeMSGFileHandler(BaseFileHandler, SEVIRICalibrationHandler):
+class NativeMSGFileHandler(BaseFileHandler):
     """SEVIRI native format reader.
 
     The Level1.5 Image data calibration method can be changed by adding the
     required mode to the Scene object instantiation  kwargs eg
     kwargs = {"calib_mode": "gsics",}
     """
 
-    def __init__(self, filename, filename_info, filetype_info, calib_mode='nominal'):
+    def __init__(self, filename, filename_info, filetype_info,
+                 calib_mode='nominal', ext_calib_coefs=None):
         """Initialize the reader."""
         super(NativeMSGFileHandler, self).__init__(filename,
                                                    filename_info,
                                                    filetype_info)
         self.platform_name = None
         self.calib_mode = calib_mode
+        self.ext_calib_coefs = ext_calib_coefs or {}
 
         # Declare required variables.
         # Assume a full disk file, reset in _read_header if otherwise.
@@ -303,7 +304,6 @@ def get_area_def(self, dataset_id):
             pdict['a_name'] = 'geos_seviri_visir'
             pdict['a_desc'] = 'SEVIRI low resolution channel area'
             pdict['p_id'] = 'seviri_visir'
-
             area = get_area_definition(pdict, self.get_area_extent(dataset_id))
 
         return area
@@ -483,51 +483,45 @@ def get_dataset(self, dataset_id, dataset_info):
     def calibrate(self, data, dataset_id):
         """Calibrate the data."""
         tic = datetime.now()
+        channel_name = dataset_id['name']
+        calib = SEVIRICalibrationHandler(
+            platform_id=self.platform_id,
+            channel_name=channel_name,
+            coefs=self._get_calib_coefs(channel_name),
+            calib_mode=self.calib_mode,
+            scan_time=self.start_time
+        )
+        res = calib.calibrate(data, dataset_id['calibration'])
+        logger.debug("Calibration time " + str(datetime.now() - tic))
+        return res
 
-        data15hdr = self.header['15_DATA_HEADER']
-        calibration = dataset_id['calibration']
-        channel = dataset_id['name']
-
+    def _get_calib_coefs(self, channel_name):
+        """Get coefficients for calibration from counts to radiance."""
         # even though all the channels may not be present in the file,
         # the header does have calibration coefficients for all the channels
         # hence, this channel index needs to refer to full channel list
-        i = list(CHANNEL_NAMES.values()).index(channel)
-
-        if calibration == 'counts':
-            return data
-
-        if calibration in ['radiance', 'reflectance', 'brightness_temperature']:
-            # determine the required calibration coefficients to use
-            # for the Level 1.5 Header
-            if (self.calib_mode.upper() != 'GSICS' and self.calib_mode.upper() != 'NOMINAL'):
-                raise NotImplementedError(
-                    'Unknown Calibration mode : Please check')
-
-            # NB GSICS doesn't have calibration coeffs for VIS channels
-            if (self.calib_mode.upper() != 'GSICS' or channel in VIS_CHANNELS):
-                coeffs = data15hdr[
-                    'RadiometricProcessing']['Level15ImageCalibration']
-                gain = coeffs['CalSlope'][i]
-                offset = coeffs['CalOffset'][i]
-            else:
-                coeffs = data15hdr[
-                    'RadiometricProcessing']['MPEFCalFeedback']
-                gain = coeffs['GSICSCalCoeff'][i]
-                offset = coeffs['GSICSOffsetCount'][i]
-                offset = offset * gain
-            res = self._convert_to_radiance(data, gain, offset)
-
-        if calibration == 'reflectance':
-            solar_irradiance = CALIB[self.platform_id][channel]["F"]
-            res = self._vis_calibrate(res, solar_irradiance)
-
-        elif calibration == 'brightness_temperature':
-            cal_type = data15hdr['ImageDescription'][
-                'Level15ImageProduction']['PlannedChanProcessing'][i]
-            res = self._ir_calibrate(res, channel, cal_type)
-
-        logger.debug("Calibration time " + str(datetime.now() - tic))
-        return res
+        band_idx = list(CHANNEL_NAMES.values()).index(channel_name)
+
+        coefs_nominal = self.header['15_DATA_HEADER'][
+            'RadiometricProcessing']['Level15ImageCalibration']
+        coefs_gsics = self.header['15_DATA_HEADER'][
+            'RadiometricProcessing']['MPEFCalFeedback']
+        radiance_types = self.header['15_DATA_HEADER']['ImageDescription'][
+                'Level15ImageProduction']['PlannedChanProcessing']
+        return {
+            'coefs': {
+                'NOMINAL': {
+                    'gain': coefs_nominal['CalSlope'][band_idx],
+                    'offset': coefs_nominal['CalOffset'][band_idx]
+                },
+                'GSICS': {
+                    'gain': coefs_gsics['GSICSCalCoeff'][band_idx],
+                    'offset': coefs_gsics['GSICSOffsetCount'][band_idx]
+                },
+                'EXTERNAL': self.ext_calib_coefs.get(channel_name, {})
+            },
+            'radiance_type': radiance_types[band_idx]
+        }
 
 
 def get_available_channels(header):