diff --git a/level1c4pps/calibration_coefs.py b/level1c4pps/calibration_coefs.py index 587262a..d7dce62 100644 --- a/level1c4pps/calibration_coefs.py +++ b/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__': diff --git a/level1c4pps/seviri2pps_lib.py b/level1c4pps/seviri2pps_lib.py index 6a94f17..6c8707b 100644 --- a/level1c4pps/seviri2pps_lib.py +++ b/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 diff --git a/level1c4pps/tests/test_seviri2pps.py b/level1c4pps/tests/test_seviri2pps.py index 0acfee8..75e7c82 100644 --- a/level1c4pps/tests/test_seviri2pps.py +++ b/level1c4pps/tests/test_seviri2pps.py @@ -24,6 +24,7 @@ import datetime as dt import numpy as np +import pytest import unittest try: from unittest import mock @@ -75,7 +76,8 @@ def test_load_and_calibrate(self, mocked_scene): res = seviri2pps.load_and_calibrate( filenames, apply_sun_earth_distance_correction=False, - rotate=False + rotate=False, + clip_calib=False ) # Compare results and expectations @@ -104,7 +106,8 @@ def test_load_and_calibrate_with_rotation(self, mocked_scene): res = seviri2pps.load_and_calibrate( filenames, apply_sun_earth_distance_correction=False, - rotate=True + rotate=True, + clip_calib=False ) scene.load.assert_called_with(mock.ANY, upper_right_corner='NE') self.assertTrue(res.attrs['image_rotated']) @@ -490,76 +493,140 @@ def test_set_nominal_scan_time(self): self.assertEqual(arr.attrs['end_time'], end_time) -class TestCalibration(unittest.TestCase): +class TestCalibration: """Test SEVIRI calibration.""" - def test_get_calibration_for_date(self): + @pytest.mark.parametrize( + "time", + (dt.date(2018, 1, 18), dt.datetime(2018, 1, 18)) + ) + def test_get_calibration_can_handle_both_date_and_time(self, time): """Test MODIS-intercalibrated gain and offset for specific date.""" - coefs = calib.get_calibration_for_date( - platform='MSG3', date=dt.date(2018, 1, 18)) + coefs = calib.get_calibration( + platform='MSG3', time=time) REF = { 'VIS006': {'gain': 0.023689275200000002, 'offset': -1.2081530352}, 'VIS008': {'gain': 0.029757990399999996, 'offset': -1.5176575103999999}, 'IR_016': {'gain': 0.0228774688, 'offset': -1.1667509087999999}} - for channel in REF.keys(): - self.assertEqual(coefs[channel]['gain'], REF[channel]['gain']) - self.assertEqual(coefs[channel]['offset'], REF[channel]['offset']) - - def test_get_calibration_for_time(self): + self._assert_coefs_close(coefs, REF) + + def _assert_coefs_close(self, coefs, expected, **tol): + for channel in expected.keys(): + np.testing.assert_allclose(coefs[channel]['gain'], + expected[channel]['gain'], + **tol) + np.testing.assert_allclose(coefs[channel]['offset'], + expected[channel]['offset'], + **tol) + + @pytest.mark.parametrize( + "platform,timestamp,expected", + ( + [ + ( + "MSG1", + dt.datetime(2005, 1, 18, 0, 0), + { + 'VIS006': {'gain': 0.0250354716, + 'offset': -1.2768090516000001}, + 'VIS008': {'gain': 0.0315626684, + 'offset': -1.6096960884}, + 'IR_016': {'gain': 0.022880986, + 'offset': -1.166930286} + } + ), + ( + "MSG2", + dt.datetime(2010, 1, 18, 0, 0), + { + 'VIS006': {'gain': 0.021964051999999998, + 'offset': -1.120166652}, + 'VIS008': {'gain': 0.027548445, + 'offset': -1.404970695}, + 'IR_016': {'gain': 0.021576766, + 'offset': -1.100415066} + }, + ), + ( + 'MSG3', + dt.datetime(2018, 1, 18, 0, 0), + { + 'VIS006': {'gain': 0.023689275200000002, + 'offset': -1.2081530352}, + 'VIS008': {'gain': 0.029757990399999996, + 'offset': -1.5176575103999999}, + 'IR_016': {'gain': 0.0228774688, + 'offset': -1.1667509087999999} + } + ), + ( + 'MSG4', + dt.datetime(2019, 1, 18, 0, 0), + { + 'VIS006': {'gain': 0.0230415289, + 'offset': -1.1751179739}, + 'VIS008': {'gain': 0.0291916818, + 'offset': -1.4887757718}, + 'IR_016': {'gain': 0.022223894, + 'offset': -1.1334185940000001} + } + ), + ( + 'MSG4', + dt.datetime(2024, 1, 18, 0, 0), + { + 'VIS006': {'gain': 0.0235668691, + 'offset': -1.2019103241}, + 'VIS008': {'gain': 0.0303007942, + 'offset': -1.5453405042000001}, + 'IR_016': {'gain': 0.022483186, + 'offset': -1.146642486} + } # extrapolation beyond time coverage of calib. dataset + ) + ] + ) + ) + def test_get_calibration(self, platform, timestamp, expected): """Test MODIS-intercalibrated gain and offset for specific time.""" - REF = { - ('MSG1', dt.datetime(2005, 1, 18, 0, 0)): { - 'VIS006': {'gain': 0.0250354716, - 'offset': -1.2768090516000001}, - 'VIS008': {'gain': 0.0315626684, - 'offset': -1.6096960884}, - 'IR_016': {'gain': 0.022880986, - 'offset': -1.166930286}}, - ('MSG2', dt.datetime(2010, 1, 18, 0, 0)): { - 'VIS006': {'gain': 0.021964051999999998, - 'offset': -1.120166652}, - 'VIS008': {'gain': 0.027548445, - 'offset': -1.404970695}, - 'IR_016': {'gain': 0.021576766, - 'offset': -1.100415066}}, - ('MSG3', dt.datetime(2018, 1, 18, 0, 0)): { - 'VIS006': {'gain': 0.023689275200000002, - 'offset': -1.2081530352}, - 'VIS008': {'gain': 0.029757990399999996, - 'offset': -1.5176575103999999}, - 'IR_016': {'gain': 0.0228774688, - 'offset': -1.1667509087999999}}, - ('MSG4', dt.datetime(2019, 1, 18, 0, 0)): { - 'VIS006': {'gain': 0.0230415289, - 'offset': -1.1751179739}, - 'VIS008': {'gain': 0.0291916818, - 'offset': -1.4887757718}, - 'IR_016': {'gain': 0.022223894, - 'offset': -1.1334185940000001}} - } - for (platform, time), ref in REF.items(): - coefs = calib.get_calibration_for_time(platform=platform, - time=time) - for channel in ref.keys(): - np.testing.assert_allclose(coefs[channel]['gain'], - ref[channel]['gain']) - np.testing.assert_allclose(coefs[channel]['offset'], - ref[channel]['offset']) - - def test_get_calibration(self): - """Test MODIS-intercalibrated for date and time.""" - coefs1 = calib.get_calibration_for_time( + coefs = calib.get_calibration(platform=platform, time=timestamp) + self._assert_coefs_close(coefs, expected) + + def test_calibration_is_smooth(self): + """Test that calibration is smooth in time.""" + coefs1 = calib.get_calibration( platform='MSG3', time=dt.datetime(2018, 1, 18, 23, 59)) - coefs2 = calib.get_calibration_for_date( - platform='MSG3', date=dt.date(2018, 1, 19)) - for channel in coefs1.keys(): - self.assertAlmostEqual(coefs1[channel]['gain'], - coefs2[channel]['gain'], - delta=10e-8) - self.assertAlmostEqual(coefs1[channel]['offset'], - coefs2[channel]['offset'], - delta=10e-8) + coefs2 = calib.get_calibration( + platform='MSG3', time=dt.datetime(2018, 1, 19)) + self._assert_coefs_close(coefs1, coefs2, atol=1e-4) + + @pytest.mark.parametrize( + "coverage_boundary,outside_coverage", + [ + (dt.datetime(2004, 1, 1), dt.datetime(2003, 1, 1)), + (dt.datetime(2021, 1, 1), dt.datetime(2022, 1, 1)), + ] + ) + def test_clip_at_time_coverage_bounds(self, coverage_boundary, outside_coverage): + """Test clipping beyond time coverage of the calibration dataset.""" + coefs1 = calib.get_calibration( + platform='MSG1', + time=coverage_boundary + ) + coefs2 = calib.get_calibration( + platform='MSG1', + time=outside_coverage, + clip=True + ) + self._assert_coefs_close(coefs1, coefs2) + + def test_fails_with_invalid_time(self): + """Test that calibration fails with timestamps < reference time.""" + with pytest.raises(ValueError): + calib.get_calibration( + platform='MSG1', + time=dt.datetime(1999, 1, 1) + ) class TestSEVIRIFilenameParser(unittest.TestCase):