Add meirink calib

satpy/readers/seviri_base.py

@@ -166,7 +166,7 @@
 """
 
 import warnings
-from datetime import timedelta
+from datetime import datetime, timedelta
 
 import dask.array as da
 import numpy as np
@@ -353,6 +353,81 @@
                          '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 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
+DATE_2000 = datetime(2000, 1, 1)
+
+MEIRINK_COEFS = {}
+MEIRINK_COEFS['2013'] = {}
+
+# Meteosat-8
+
+MEIRINK_COEFS['2013'][321] = {'VIS006': (24.346, 0.3739),
+                              'VIS008': (30.989, 0.3111),
+                              'IR_016': (22.869, 0.0065)
+                              }
+
+# Meteosat-9
+
+MEIRINK_COEFS['2013'][322] = {'VIS006': (21.026, 0.2556),
+                              'VIS008': (26.875, 0.1835),
+                              'IR_016': (21.394, 0.0498)
+                              }
+
+# Meteosat-10
+
+MEIRINK_COEFS['2013'][323] = {'VIS006': (19.829, 0.5856),
+                              'VIS008': (25.284, 0.6787),
+                              'IR_016': (23.066, -0.0286)
+                              }
+
+# Meteosat-11
+
+MEIRINK_COEFS['2013'][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 - DATE_2000).total_seconds()
+    S = A + B * delta_t / (3600*24) / 1000.
+    return S/1000
+
+
+class MeirinkCalibrationHandler:
+    """Re-calibration of the SEVIRI visible channels slop (see Meirink 2013)."""
+
+    def __init__(self, coefs=MEIRINK_COEFS, calib_mode=None):
+        """Initialize the calibration handler."""
+        if calib_mode is None:
+            raise ValueError("Missing calib_mode")
+        self.coefs = 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 +641,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-2013')
         if self._calib_mode not in valid_modes:
             raise ValueError(
                 'Invalid calibration mode: {}. Choose one of {}'.format(
@@ -622,6 +697,10 @@ def get_gain_offset(self):
                 internal_gain = gsics_gain
                 internal_offset = gsics_offset
 
+        if "MEIRINK" in self._calib_mode and self._channel_name in ['VIS006', 'VIS008', 'IR_016']:
+            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 (
+    DATE_2000,
+    MEIRINK_COEFS,
     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', DATE_2000)
+        assert calibration_handler.get_gain_offset()[0] == MEIRINK_COEFS[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', time)
+            assert abs(calibration_handler.get_gain_offset()[0] - expected[i]) < 1e-6