Merge pull request #105 from ninahakansson/seviri_ir_calibration

SEVIRI IR calibration

Author: ninahakansson

bin/seviri2pps.py

@@ -50,12 +50,19 @@
                         help="Engine for saving netcdf files netcdf4 or h5netcdf (default).")
     parser.add_argument('--use-nominal-time-in-filename', action='store_true',
                         help='Use nominal scan timestamps in output filename.')
+    parser.add_argument('--path-to-external-ir-calibration', type=str, default=".",
+                        help='Path to external IR calibration.')
+    parser.add_argument('--use-nominal-ir-calibration', action='store_true',
+                        help='Use nominal IR claibration.')
     options = parser.parse_args()
+    if options.use_nominal_ir_calibration:
+        options.path_to_external_ir_calibration = None
     process_one_scan(
         options.files,
         out_path=options.out_dir,
         rotate=not options.no_rotation,
         engine=options.nc_engine,
         use_nominal_time_in_filename=options.use_nominal_time_in_filename,
+        path_to_external_ir_calibration=options.path_to_external_ir_calibration,
         save_azimuth_angles=options.azimuth_angles,
     )

level1c4pps/calibration_coefs.py

@@ -21,6 +21,11 @@
 
 import datetime
 from enum import Enum
+import json
+import logging
+import os
+
+logger = logging.getLogger("calibration")
 
 
 class CalibrationData(Enum):
@@ -55,7 +60,7 @@ class CalibrationData(Enum):
     SATPY_CALIB_MODE = 'Nominal'
 
 
-def get_calibration(platform, time, clip=False):
+def get_calibration(platform, time, clip=False, calib_ir_path=None):
     """Get MODIS-intercalibrated gain and offset for specific time.
 
     Args:
@@ -74,6 +79,15 @@ def get_calibration(platform, time, clip=False):
             time=time,
             clip=clip
         )
+    if calib_ir_path is not None:
+        for channel in ('IR_039', 'IR_087', 'IR_108', 'IR_120',
+                        'IR_134', 'IR_097', 'WV_062', 'WV_073'):
+            coefs[channel] = get_ir_calibration_coeffs(
+                calib_ir_path,
+                platform=platform,
+                channel=channel,
+                time=time,
+            )
     return coefs
 
 
@@ -106,11 +120,13 @@ def _convert_to_datetime(date_or_time):
 
 
 def _is_date(date_or_time):
+    """Check that we have a datetime date object."""
     # datetime is a subclass of date, therefore we cannot use isinstance here
-    return type(date_or_time) == datetime.date
+    return type(date_or_time) == datetime.date  # noqa E721
 
 
 def _check_is_valid_time(time):
+    """Check that we have a valid time."""
     ref_time = CalibrationData.REF_TIME.value
     if time < ref_time:
         raise ValueError('Given time ({0}) is < reference time ({1})'.format(
@@ -154,6 +170,22 @@ def _microwatts_to_milliwatts(microwatts):
     return microwatts / 1000.0
 
 
+def get_ir_calibration_coeffs(ir_calib_path, platform="MSG2", channel="IR_039",
+                              time=datetime.datetime(2048, 1, 18, 12, 0)):
+    """Get IR calibration from EUMETSAT, modified by CMSAF."""
+    filename = os.path.join(ir_calib_path, f"TIR_calib_{platform}_{channel}.json")
+    logger.info(f'Using IR calibration from {filename}')
+    with open(filename, 'r') as fhand:
+        data = json.load(fhand)
+        for item in data:
+            date_s = item[0]
+            date_i = datetime.datetime.strptime(date_s, "%Y-%m-%dT%H:%M:%S.%f")
+            if date_i > time:
+                break
+            gain, offset = item[1:]
+    return {'gain': gain, 'offset': offset}
+
+
 if __name__ == '__main__':
     time = datetime.datetime(2018, 1, 18, 12, 0)
     platform = 'MSG3'
@@ -164,4 +196,9 @@ def _microwatts_to_milliwatts(microwatts):
                                      time=time)
         coefs[channel] = {'gain': gain, 'offset': offset}
 
+    for channel in ('IR_039', 'IR_087', 'IR_108', 'IR_120',
+                    'IR_134', 'IR_097', 'WV_062', 'WV_073'):
+        gain, offset = get_ir_calibration_coeffs(platform=platform, channel=channel,
+                                                 time=time)
+        coefs[channel] = {'gain': gain, 'offset': offset}
     print(coefs)

level1c4pps/seviri2pps_lib.py

@@ -87,7 +87,7 @@ class UnexpectedSatpyVersion(Exception):
 
 
 def load_and_calibrate(filenames, rotate,
-                       clip_calib):
+                       clip_calib, path_to_external_ir_calibration=None):
     """Load and calibrate data.
 
     Uses inter-calibration coefficients from Meirink et al.
@@ -109,8 +109,8 @@ def load_and_calibrate(filenames, rotate,
     calib_coefs = get_calibration(
         platform=info['platform_shortname'],
         time=info['start_time'],
-        clip=clip_calib
-    )
+        clip=clip_calib,
+        calib_ir_path=path_to_external_ir_calibration)
     scn_ = _create_scene(file_format, filenames, calib_coefs)
     _check_is_seviri_data(scn_)
     _load_bands(scn_, rotate)
@@ -563,6 +563,7 @@ def _postproc_hrit(self, parsed):
 def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf',
                      use_nominal_time_in_filename=False,
                      clip_calib=False,
+                     path_to_external_ir_calibration=None,
                      save_azimuth_angles=False):
     """Make level 1c files in PPS-format."""
     for fname in tslot_files:
@@ -573,7 +574,8 @@ def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf',
     scn_ = load_and_calibrate(
         tslot_files,
         rotate=rotate,
-        clip_calib=clip_calib
+        clip_calib=clip_calib,
+        path_to_external_ir_calibration=path_to_external_ir_calibration,
     )
     if hasattr(scn_, 'start_time'):
         scn_.attrs['start_time'] = scn_.start_time
@@ -628,7 +630,11 @@ def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf',
     return filename
 
 
-def process_all_scans_in_dname(dname, out_path, ok_dates=None, rotate=False):
+def process_all_scans_in_dname(dname, out_path,
+                               ok_dates=None,
+                               rotate=False,
+                               path_to_external_ir_calibration=None,
+                               save_azimuth_angles=False):
     """Make level 1c files for all files in directory dname."""
     parser = Parser(HRIT_FILE_PATTERN)
     fl_ = glob(os.path.join(dname, globify(HRIT_FILE_PATTERN)))
@@ -645,6 +651,10 @@ def process_all_scans_in_dname(dname, out_path, ok_dates=None, rotate=False):
         tslot_files = [f for f in fl_ if parser.parse(
             os.path.basename(f))['start_time'] == uqdate]
         try:
-            process_one_scan(tslot_files, out_path, rotate=rotate)
+            process_one_scan(tslot_files,
+                             out_path,
+                             rotate=rotate,
+                             path_to_external_ir_calibration=path_to_external_ir_calibration,
+                             save_azimuth_angles=save_azimuth_angles)
         except Exception:
             pass

level1c4pps/tests/TIR_calib_MSG2_IR_039.json

@@ -0,0 +1 @@
+[["2007-06-05T00:00:00.000",0.0035676316,-0.1819492105]]

level1c4pps/tests/TIR_calib_MSG2_IR_087.json

@@ -0,0 +1 @@
+[["2007-06-05T00:00:00.000",0.0035676316,-0.1819492105]]

level1c4pps/tests/TIR_calib_MSG2_IR_097.json

@@ -0,0 +1 @@
+[["2007-06-05T00:00:00.000",0.0035676316,-0.1819492105]]

level1c4pps/tests/TIR_calib_MSG2_IR_108.json

@@ -0,0 +1 @@
+[["2007-06-05T00:00:00.000",0.0035676316,-0.1819492105]]

level1c4pps/tests/TIR_calib_MSG2_IR_120.json

@@ -0,0 +1 @@
+[["2007-06-05T00:00:00.000",0.035676316,-1.819492105],["2007-06-17T23:59:0.000",0.0035676316,-0.1819492105],["2007-06-18T01:00:00.000",0.035676316,-1.819492105]]

level1c4pps/tests/TIR_calib_MSG2_IR_134.json

@@ -0,0 +1 @@
+[["2007-06-05T00:00:00.000",0.0035676316,-0.1819492105]]

level1c4pps/tests/TIR_calib_MSG2_WV_062.json

@@ -0,0 +1 @@
+[["2007-06-05T00:00:00.000",0.0035676316,-0.1819492105]]

level1c4pps/tests/TIR_calib_MSG2_WV_073.json

@@ -0,0 +1 @@
+[["2007-06-05T00:00:00.000",0.0035676316,-0.1819492105]]

level1c4pps/tests/test_seviri2pps.py

@@ -72,7 +72,7 @@ def test_load_and_calibrate(self, mocked_scene):
         res = seviri2pps.load_and_calibrate(
             filenames,
             rotate=False,
-            clip_calib=False
+            clip_calib=False,
         )
 
         # Compare results and expectations
@@ -613,6 +613,22 @@ def test_get_calibration(self, platform, timestamp, expected):
         coefs = calib.get_calibration(platform=platform, time=timestamp)
         self._assert_coefs_close(coefs, expected)
 
+    def test_get_calibration_ir_no_file(self):
+        """Test get ir calibration with mising json file."""
+        with pytest.raises(FileNotFoundError):
+            coefs = calib.get_calibration(
+                platform="MSG2",
+                time=dt.datetime(2007, 6, 5, 0, 0),
+                calib_ir_path=".")
+
+    def test_get_calibration_ir(self):
+        """Test get ir calibration.."""
+        coefs = calib.get_calibration(
+            platform="MSG2",
+            time=dt.datetime(2007, 6, 18, 0, 0),
+            calib_ir_path="./level1c4pps/tests/")
+        np.testing.assert_almost_equal(coefs['IR_120']['gain'], 0.003567, decimal=6)
+
     def test_calibration_is_smooth(self):
         """Test that calibration is smooth in time."""
         coefs1 = calib.get_calibration(