Merge pull request #71 from sfinkens/fix-calib

Add clipping for SEVIRI calibration coefficients

level1c4pps/calibration_coefs.py

@@ -24,94 +24,138 @@
 """Module with calibration coefficients for SEVIRI."""
 
 import datetime
-
-CALIB_MODE = 'Nominal'
-COEFS_MEIRINK = dict(
-    MSG1=dict(
-        VIS006=dict(b=24.346, a=0.3739E-3),
-        VIS008=dict(b=30.989, a=0.3111E-3),
-        IR_016=dict(b=22.869, a=0.0065E-3)
-    ),
-    MSG2=dict(
-        VIS006=dict(b=21.026, a=0.2556E-3),
-        VIS008=dict(b=26.875, a=0.1835E-3),
-        IR_016=dict(b=21.394, a=0.0498E-3)
-    ),
-    MSG3=dict(
-        VIS006=dict(b=19.829, a=0.5856E-3),
-        VIS008=dict(b=25.284, a=0.6787E-3),
-        IR_016=dict(b=23.066, a=-0.0286E-3)
-    ),
-    MSG4=dict(
-        VIS006=dict(b=21.040, a=0.2877E-3),
-        VIS008=dict(b=24.966, a=0.6074E-3),
-        IR_016=dict(b=21.236, a=0.1420E-3)
+from enum import Enum
+
+
+class CalibrationData(Enum):
+    COEFS = dict(
+        MSG1=dict(
+            VIS006=dict(b=24.346, a=0.3739E-3),
+            VIS008=dict(b=30.989, a=0.3111E-3),
+            IR_016=dict(b=22.869, a=0.0065E-3)
+        ),
+        MSG2=dict(
+            VIS006=dict(b=21.026, a=0.2556E-3),
+            VIS008=dict(b=26.875, a=0.1835E-3),
+            IR_016=dict(b=21.394, a=0.0498E-3)
+        ),
+        MSG3=dict(
+            VIS006=dict(b=19.829, a=0.5856E-3),
+            VIS008=dict(b=25.284, a=0.6787E-3),
+            IR_016=dict(b=23.066, a=-0.0286E-3)
+        ),
+        MSG4=dict(
+            VIS006=dict(b=21.040, a=0.2877E-3),
+            VIS008=dict(b=24.966, a=0.6074E-3),
+            IR_016=dict(b=21.236, a=0.1420E-3)
+        )
     )
-)
+    SPACE_COUNT = -51.0
+    REF_TIME = datetime.datetime(2000, 1, 1, 0, 0)
+    TIME_COVERAGE = {
+        "start": datetime.datetime(2004, 1, 1, 0, 0),
+        "end": datetime.datetime(2021, 1, 1, 0, 0)
+    }
+    SATPY_CALIB_MODE = 'Nominal'
+
+
+def get_calibration(platform, time, clip=False):
+    """Get MODIS-intercalibrated gain and offset for specific time.
+
+    Args:
+        platform: Platform name.
+        time: Date or time of observations to be calibrated.
+        clip: If True, do not extrapolate calibration coefficients beyond the
+            time coverage of the calibration dataset. Instead, clip at the
+            boundaries, that means return the boundary coefficients for
+            timestamps outside the coverage.
+    """
+    coefs = {}
+    for channel in ('VIS006', 'VIS008', 'IR_016'):
+        coefs[channel] = _get_single_channel_calibration(
+            platform=platform,
+            channel=channel,
+            time=time,
+            clip=clip
+        )
+    return coefs
 
-REF_DATE = datetime.date(2000, 1, 1)
-REF_TIME = datetime.datetime(2000, 1, 1, 0, 0)
 
+def _get_single_channel_calibration(platform, channel, time, clip):
+    """Get calibration coefficients for a single channel.
 
-def calib_meirink(platform, channel, time):
-    """Get MODIS-intercalibrated gain and offset for SEVIRI VIS channels.
+    Args:
+        platform: Platform name.
+        channel: Channel name.
+        time: Observation date or time.
+        clip: Clip at time coverage boundaries of the calibration dataset.
+    """
+    time = _prepare_time(time, clip)
+    gain, offset = calib_meirink(platform, channel, time)
+    return {'gain': gain, 'offset': offset}
 
-    Reference: http://msgcpp.knmi.nl/mediawiki/index.php/MSG-SEVIRI_solar_channel_calibration
 
-    :returns: gain, offset [mW m-2 sr-1 (cm-1)-1]
-    """
-    if time < REF_TIME:
+def _prepare_time(time, clip):
+    time = _convert_to_datetime(time)
+    _check_is_valid_time(time)
+    if clip:
+        time = _clip_at_coverage_bounds(time)
+    return time
+
+
+def _convert_to_datetime(date_or_time):
+    if _is_date(date_or_time):
+        return datetime.datetime.combine(date_or_time, datetime.time(0))
+    return date_or_time
+
+
+def _is_date(date_or_time):
+    # datetime is a subclass of date, therefore we cannot use isinstance here
+    return type(date_or_time) == datetime.date
+
+
+def _check_is_valid_time(time):
+    ref_time = CalibrationData.REF_TIME.value
+    if time < ref_time:
         raise ValueError('Given time ({0}) is < reference time ({1})'.format(
-            time, REF_TIME))
+            time, ref_time))
 
-    a = COEFS_MEIRINK[platform][channel]['a']
-    b = COEFS_MEIRINK[platform][channel]['b']
-    delta_days = (time - REF_TIME).total_seconds() / 3600.0 / 24.0
-    gain = (b + a * delta_days) / 1000.0  # micro Watts -> milli Watts
-    offset = -51.0 * gain  # Space count is 51
 
-    return gain, offset
+def _clip_at_coverage_bounds(time):
+    time_cov = CalibrationData.TIME_COVERAGE.value
+    time = max(time, time_cov["start"])
+    time = min(time, time_cov["end"])
+    return time
 
 
-def calib_meirink_date(platform, channel, date):
+def calib_meirink(platform, channel, time):
     """Get MODIS-intercalibrated gain and offset for SEVIRI VIS channels.
 
-    Reference: http://msgcpp.knmi.nl/mediawiki/index.php/MSG-SEVIRI_solar_channel_calibration
+    Reference: https://msgcpp.knmi.nl/solar-channel-calibration.html
 
     :returns: gain, offset [mW m-2 sr-1 (cm-1)-1]
     """
-    if date < REF_DATE:
-        raise ValueError('Given date ({0}) is < reference date ({1})'.format(
-            date, REF_DATE))
-
-    a = COEFS_MEIRINK[platform][channel]['a']
-    b = COEFS_MEIRINK[platform][channel]['b']
-    gain = (b + a*(date - REF_DATE).days) / 1000.0  # micro Watts -> milli Watts
-    offset = -51.0 * gain  # Space count is 51
-
-    return gain, offset
+    coefs = CalibrationData.COEFS.value
+    a = coefs[platform][channel]['a']
+    b = coefs[platform][channel]['b']
+    days_since_ref_time = _get_days_since_ref_time(time)
+    return _calc_gain_offset(a, b, days_since_ref_time)
 
 
-def get_calibration_for_time(platform, time):
-    """Get MODIS-intercalibrated gain and offset for specific time."""
-    coefs = {}
-    for channel in ('VIS006', 'VIS008', 'IR_016'):
-        gain, offset = calib_meirink(platform=platform, channel=channel,
-                                     time=time)
-        coefs[channel] = {'gain': gain, 'offset': offset}
+def _get_days_since_ref_time(time):
+    ref_time = CalibrationData.REF_TIME.value
+    return (time - ref_time).total_seconds() / 3600.0 / 24.0
 
-    return coefs
 
+def _calc_gain_offset(a, b, days_since_ref_time):
+    gain = (b + a * days_since_ref_time)
+    gain = _microwatts_to_milliwatts(gain)
+    offset = CalibrationData.SPACE_COUNT.value * gain
+    return gain, offset
 
-def get_calibration_for_date(platform, date):
-    """Get MODIS-intercalibrated gain and offset for specific date."""
-    coefs = {}
-    for channel in ('VIS006', 'VIS008', 'IR_016'):
-        gain, offset = calib_meirink_date(platform=platform, channel=channel,
-                                          date=date)
-        coefs[channel] = {'gain': gain, 'offset': offset}
 
-    return coefs
+def _microwatts_to_milliwatts(microwatts):
+    return microwatts / 1000.0
 
 
 if __name__ == '__main__':

level1c4pps/seviri2pps_lib.py

@@ -43,7 +43,7 @@
 from pyorbital.astronomy import get_alt_az, sun_zenith_angle
 from pyorbital.orbital import get_observer_look
 
-from level1c4pps.calibration_coefs import get_calibration_for_time, CALIB_MODE
+from level1c4pps.calibration_coefs import get_calibration, CalibrationData
 from level1c4pps import make_azidiff_angle, get_encoding, compose_filename, update_angle_attributes
 
 
@@ -89,7 +89,8 @@ class UnexpectedSatpyVersion(Exception):
 # MSG4-SEVI-MSG15-1234-NA-20190409121243.927000000Z
 
 
-def load_and_calibrate(filenames, apply_sun_earth_distance_correction, rotate):
+def load_and_calibrate(filenames, apply_sun_earth_distance_correction, rotate,
+                       clip_calib):
     """Load and calibrate data.
 
     Uses inter-calibration coefficients from Meirink et al.
@@ -99,6 +100,9 @@ def load_and_calibrate(filenames, apply_sun_earth_distance_correction, rotate):
         apply_sun_earth_distance_correction: If True, apply sun-earth-distance-
             correction to visible channels.
         rotate: Rotate image so that pixel (0, 0) is N-W.
+        clip_calib: If True, do not extrapolate calibration coefficients beyond
+            the time coverage of the calibration dataset. Instead, clip at the
+            boundaries.
 
     Returns:
         Satpy scene holding calibrated channels
@@ -107,9 +111,10 @@ def load_and_calibrate(filenames, apply_sun_earth_distance_correction, rotate):
     parser = SEVIRIFilenameParser()
     file_format, info = parser.parse(os.path.basename(filenames[0]))
 
-    calib_coefs = get_calibration_for_time(
+    calib_coefs = get_calibration(
         platform=info['platform_shortname'],
-        time=info['start_time']
+        time=info['start_time'],
+        clip=clip_calib
     )
     scn_ = _create_scene(file_format, filenames, calib_coefs)
     _check_is_seviri_data(scn_)
@@ -123,10 +128,12 @@ def load_and_calibrate(filenames, apply_sun_earth_distance_correction, rotate):
 
 
 def _create_scene(file_format, filenames, calib_coefs):
-    return Scene(reader=file_format,
-                 filenames=filenames,
-                 reader_kwargs={'calib_mode': CALIB_MODE,
-                                'ext_calib_coefs': calib_coefs})
+    return Scene(
+        reader=file_format,
+        filenames=filenames,
+        reader_kwargs={'calib_mode': CalibrationData.SATPY_CALIB_MODE.value,
+                       'ext_calib_coefs': calib_coefs}
+    )
 
 
 def _check_is_seviri_data(scene):
@@ -556,7 +563,8 @@ def _postproc_hrit(self, parsed):
 
 def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf',
                      use_nominal_time_in_filename=False,
-                     apply_sun_earth_distance_correction=True):
+                     apply_sun_earth_distance_correction=True,
+                     clip_calib=False):
     """Make level 1c files in PPS-format."""
     for fname in tslot_files:
         if not os.path.isfile(fname):
@@ -566,7 +574,8 @@ def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf',
     scn_ = load_and_calibrate(
         tslot_files,
         apply_sun_earth_distance_correction=apply_sun_earth_distance_correction,
-        rotate=rotate
+        rotate=rotate,
+        clip_calib=clip_calib
     )
 
     # Find lat/lon data