Add meirink calib

satpy/readers/seviri_base.py

@@ -61,14 +61,22 @@
                         reader_kwargs={'calib_mode': 'GSICS'})
     scene.load(['VIS006', 'IR_108'])
 
-Furthermore, it is possible to specify external calibration coefficients
-for the conversion from counts to radiances. External coefficients take
-precedence over internal coefficients, but you can also mix internal and
-external coefficients: If external calibration coefficients are specified
-for only a subset of channels, the remaining channels will be calibrated
-using the chosen file-internal coefficients (nominal or GSICS).
+In addition, two other calibration methods are available:
+
+1. It is possible to specify external calibration coefficients for the
+   conversion from counts to radiances. External coefficients take
+   precedence over internal coefficients and over the Meirink
+   coefficients, but you can also mix internal and external coefficients:
+   If external calibration coefficients are specified for only a subset
+   of channels, the remaining channels will be calibrated using the
+   chosen file-internal coefficients (nominal or GSICS).  Calibration
+   coefficients must be specified in [mW m-2 sr-1 (cm-1)-1].
+
+2. The calibration mode ``meirink-2023`` uses coefficients based on an
+   intercalibration with Aqua-MODIS for the visible channels, as found in
+   `Inter-calibration of polar imager solar channels using SEVIRI`_
+   (2013) by J. F. Meirink, R. A. Roebeling, and P. Stammes.
 
-Calibration coefficients must be specified in [mW m-2 sr-1 (cm-1)-1].
 
 In the following example we use external calibration coefficients for the
 ``VIS006`` & ``IR_108`` channels, and nominal coefficients for the
@@ -93,6 +101,15 @@
                                        'ext_calib_coefs': coefs})
     scene.load(['VIS006', 'VIS008', 'IR_108', 'IR_120'])
 
+In the next example we use the mode ``meirink-2023`` calibration
+coefficients for all visible channels and nominal coefficients for the
+rest::
+
+    scene = satpy.Scene(filenames,
+                        reader='seviri_l1b_...',
+                        reader_kwargs={'calib_mode': 'meirink-2023'})
+    scene.load(['VIS006', 'VIS008', 'IR_016'])
+
 
 Calibration to reflectance
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -163,10 +180,13 @@
 .. _Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent Spectral Blackbody Radiance:
     https://www-cdn.eumetsat.int/files/2020-04/pdf_ten_msg_seviri_rad_calib.pdf
 
+.. _Inter-calibration of polar imager solar channels using SEVIRI:
+   http://dx.doi.org/10.5194/amt-6-2495-2013
+
 """
 
 import warnings
-from datetime import timedelta
+from datetime import datetime, timedelta
 
 import dask.array as da
 import numpy as np
@@ -353,6 +373,87 @@
                          'ALPHA': 0.9981,
                          'BETA': 0.5635}}
 
+# Calibration coefficients from Meirink, J.F., R.A. Roebeling and P. Stammes, 2013:
+# Inter-calibration of polar imager solar channels using SEVIRI, Atm. Meas. Tech., 6,
+# 2495-2508, doi:10.5194/amt-6-2495-2013
+#
+# The coeffients in the 2023 entry have been obtained from the webpage
+# https://msgcpp.knmi.nl/solar-channel-calibration.html on 2023-10-11.
+#
+# The coefficients are stored in pairs of A, B (see function `get_meirink_slope`) where the
+# units of A are µW m-2 sr-1 (cm-1)-1 and those of B are µW m-2 sr-1 (cm-1)-1 (86400 s)-1
+#
+# To obtain the slope for the calibration, one should use the routine get_seviri_meirink_slope
+
+# Epoch for the MEIRINK re-calibration
+MEIRINK_EPOCH = datetime(2000, 1, 1)
+
+MEIRINK_COEFS = {}
+MEIRINK_COEFS['2023'] = {}
+
+# Meteosat-8
+
+MEIRINK_COEFS['2023'][321] = {'VIS006': (24.346, 0.3739),
+                              'VIS008': (30.989, 0.3111),
+                              'IR_016': (22.869, 0.0065)
+                              }
+
+# Meteosat-9
+
+MEIRINK_COEFS['2023'][322] = {'VIS006': (21.026, 0.2556),
+                              'VIS008': (26.875, 0.1835),
+                              'IR_016': (21.394, 0.0498)
+                              }
+
+# Meteosat-10
+
+MEIRINK_COEFS['2023'][323] = {'VIS006': (19.829, 0.5856),
+                              'VIS008': (25.284, 0.6787),
+                              'IR_016': (23.066, -0.0286)
+                              }
+
+# Meteosat-11
+
+MEIRINK_COEFS['2023'][324] = {'VIS006': (20.515, 0.3600),
+                              'VIS008': (25.803, 0.4844),
+                              'IR_016': (22.354, -0.0187)
+                              }
+
+
+def get_meirink_slope(meirink_coefs, acquisition_time):
+    """Compute the slope for the visible channel calibration according to Meirink 2013.
+
+    S = A + B * 1.e-3* Day
+
+    S is here in µW m-2 sr-1 (cm-1)-1
+
+    EUMETSAT calibration is given in mW m-2 sr-1 (cm-1)-1, so an extra factor of 1/1000 must
+    be applied.
+    """
+    A = meirink_coefs[0]
+    B = meirink_coefs[1]
+    delta_t = (acquisition_time - MEIRINK_EPOCH).total_seconds()
+    S = A + B * delta_t / (3600*24) / 1000.
+    return S/1000
+
+
+def should_apply_meirink(calib_mode, channel_name):
+    """Decide whether to use the Meirink calibration coefficients."""
+    return "MEIRINK" in calib_mode and channel_name in ['VIS006', 'VIS008', 'IR_016']
+
+
+class MeirinkCalibrationHandler:
+    """Re-calibration of the SEVIRI visible channels slope (see Meirink 2013)."""
+
+    def __init__(self, calib_mode):
+        """Initialize the calibration handler."""
+        self.coefs = MEIRINK_COEFS[calib_mode.split('-')[1]]
+
+    def get_slope(self, platform, channel, time):
+        """Return the slope using the provided calibration coefficients."""
+        coefs = self.coefs[platform][channel]
+        return get_meirink_slope(coefs, time)
+
 
 def get_cds_time(days, msecs):
     """Compute timestamp given the days since epoch and milliseconds of the day.
@@ -566,7 +667,7 @@ def __init__(self, platform_id, channel_name, coefs, calib_mode, scan_time):
             scan_time=self._scan_time
         )
 
-        valid_modes = ('NOMINAL', 'GSICS')
+        valid_modes = ('NOMINAL', 'GSICS', 'MEIRINK-2023')
         if self._calib_mode not in valid_modes:
             raise ValueError(
                 'Invalid calibration mode: {}. Choose one of {}'.format(
@@ -622,6 +723,10 @@ def get_gain_offset(self):
                 internal_gain = gsics_gain
                 internal_offset = gsics_offset
 
+        if should_apply_meirink(self._calib_mode, self._channel_name):
+            meirink = MeirinkCalibrationHandler(calib_mode=self._calib_mode)
+            internal_gain = meirink.get_slope(self._platform_id, self._channel_name, self._scan_time)
+
         # Override with external coefficients, if any.
         gain = coefs['EXTERNAL'].get('gain', internal_gain)
         offset = coefs['EXTERNAL'].get('offset', internal_offset)

satpy/tests/reader_tests/test_seviri_base.py

@@ -26,9 +26,12 @@
 import xarray as xr
 
 from satpy.readers.seviri_base import (
+    MEIRINK_COEFS,
+    MEIRINK_EPOCH,
     NoValidOrbitParams,
     OrbitPolynomial,
     OrbitPolynomialFinder,
+    SEVIRICalibrationHandler,
     chebyshev,
     dec10216,
     get_cds_time,
@@ -358,3 +361,36 @@ def test_get_orbit_polynomial_exceptions(self, orbit_polynomials, time):
         finder = OrbitPolynomialFinder(orbit_polynomials)
         with pytest.raises(NoValidOrbitParams):
             finder.get_orbit_polynomial(time=time)
+
+
+class TestMeirinkSlope:
+    """Unit tests for the slope of Meirink calibration."""
+
+    @pytest.mark.parametrize('platform_id', [321, 322, 323, 324])
+    @pytest.mark.parametrize('channel_name', ['VIS006', 'VIS008', 'IR_016'])
+    def test_get_meirink_slope_epoch(self, platform_id, channel_name):
+        """Test the value of the slope of the Meirink calibration on 2000-01-01."""
+        coefs = {'coefs': {}}
+        coefs['coefs']['NOMINAL'] = {'gain': -1, 'offset': -1}
+        coefs['coefs']['EXTERNAL'] = {}
+        calibration_handler = SEVIRICalibrationHandler(platform_id, channel_name, coefs, 'MEIRINK-2023', MEIRINK_EPOCH)
+        assert calibration_handler.get_gain_offset()[0] == MEIRINK_COEFS['2023'][platform_id][channel_name][0]/1000.
+
+    @pytest.mark.parametrize('platform_id,time,expected', (
+        (321, datetime(2005, 1, 18, 0, 0), [0.0250354716, 0.0315626684, 0.022880986]),
+        (321, datetime(2010, 12, 31, 0, 0), [0.0258479563, 0.0322386887, 0.022895110500000003]),
+        (322, datetime(2010, 1, 18, 0, 0), [0.021964051999999998, 0.027548445, 0.021576766]),
+        (322, datetime(2015, 6, 1, 0, 0), [0.022465028, 0.027908105, 0.021674373999999996]),
+        (323, datetime(2005, 1, 18, 0, 0), [0.0209088464, 0.0265355228, 0.0230132616]),
+        (323, datetime(2010, 12, 31, 0, 0), [0.022181355200000002, 0.0280103379, 0.0229511138]),
+        (324, datetime(2010, 1, 18, 0, 0), [0.0218362, 0.027580748, 0.022285370999999998]),
+        (324, datetime(2015, 6, 1, 0, 0), [0.0225418, 0.028530172, 0.022248718999999997]),
+    ))
+    def test_get_meirink_slope_2020(self, platform_id, time, expected):
+        """Test the value of the slope of the Meirink calibration."""
+        coefs = {'coefs': {}}
+        coefs['coefs']['NOMINAL'] = {'gain': -1, 'offset': -1}
+        coefs['coefs']['EXTERNAL'] = {}
+        for i, channel_name in enumerate(['VIS006', 'VIS008', 'IR_016']):
+            calibration_handler = SEVIRICalibrationHandler(platform_id, channel_name, coefs, 'MEIRINK-2023', time)
+            assert abs(calibration_handler.get_gain_offset()[0] - expected[i]) < 1e-6