Use internal calibration coefficients if available

satpy/readers/ami_l1b.py

@@ -37,36 +37,12 @@
     'GK-2B': 'GEO-KOMPSAT-2B',
 }
 
-# Copied from 20190415_GK-2A_AMI_Conversion_Table_v3.0.xlsx
-# Sheet: coeff.& equation_WN
-# Visible channels
-# channel_name -> (DN2Rad_Gain, DN2Rad_Offset, Rad. to Albedo)
-# IR channels
-# channel_name -> (DN2Rad_Gain, DN2Rad_Offset, c0, c1, c2)
-CALIBRATION_COEFFS = {
-    "VI004": (0.363545805215835, -7.27090454101562, 0.001558245),
-    "VI005": (0.343625485897064, -6.87249755859375, 0.0016595767),
-    "VI006": (0.154856294393539, -6.19424438476562, 0.001924484),
-    "VI008": (0.0457241721451282, -3.65792846679687, 0.0032723873),
-    "NR013": (0.0346878096461296, -1.38751220703125, 0.0087081313),
-    "NR016": (0.0498007982969284, -0.996017456054687, 0.0129512876),
-    "SW038": (-0.00108296517282724, 17.699987411499, -0.447843939824124, 1.00065568090389, -0.0000000633824089912448),
-    "WV063": (-0.0108914673328399, 44.1777038574218, -1.76279494011147, 1.00414910562278, -0.000000983310914319385),
-    "WV069": (-0.00818779878318309, 66.7480773925781, -0.334311414359106, 1.00097359874468, -0.000000494603070252304),
-    "WV073": (-0.0096982717514038, 79.0608520507812, -0.0613124859696595, 1.00019008722941, -0.000000105863656750499),
-    "IR087": (-0.0144806550815701, 118.050903320312, -0.141418528203155, 1.00052232906885, -0.00000036287276076109),
-    "IR096": (-0.0178435463458299, 145.464874267578, -0.114017728158198, 1.00047380585402, -0.000000374931509928403),
-    "IR105": (-0.0198196955025196, 161.580139160156, -0.142866448475177, 1.00064069572049, -0.000000550443294960498),
-    "IR112": (-0.0216744858771562, 176.713439941406, -0.249111718496148, 1.00121166873756, -0.00000113167964011665),
-    "IR123": (-0.023379972204566, 190.649627685546, -0.458113885722738, 1.00245520975535, -0.00000253064314720476),
-    "IR133": (-0.0243037566542625, 198.224365234375, -0.0938521568527657, 1.00053982112966, -0.000000594913715312849),
-}
-
 
 class AMIL1bNetCDF(BaseFileHandler):
     """Base reader for AMI L1B NetCDF4 files."""
 
-    def __init__(self, filename, filename_info, filetype_info, allow_conditional_pixels=False):
+    def __init__(self, filename, filename_info, filetype_info,
+                 calib_mode='PYSPECTRAL', allow_conditional_pixels=False):
         """Open the NetCDF file with xarray and prepare the Dataset for reading."""
         super(AMIL1bNetCDF, self).__init__(filename, filename_info, filetype_info)
         self.nc = xr.open_dataset(self.filename,
@@ -79,6 +55,11 @@ def __init__(self, filename, filename_info, filetype_info, allow_conditional_pix
         self.platform_name = PLATFORM_NAMES.get(platform_shortname)
         self.sensor = 'ami'
         self.allow_conditional_pixels = allow_conditional_pixels
+        calib_mode_choices = ('FILE', 'PYSPECTRAL')
+        if calib_mode.upper() not in calib_mode_choices:
+            raise ValueError('Invalid calibration mode: {}. Choose one of {}'.format(
+                calib_mode, calib_mode_choices))
+        self.calib_mode = calib_mode.upper()
 
     @property
     def start_time(self):
@@ -182,21 +163,35 @@ def get_dataset(self, dataset_id, ds_info):
         # only take "no error" pixels as valid
         data = data.where(qf == 0)
 
-        channel_name = attrs.get('channel_name', dataset_id.name)
-        coeffs = CALIBRATION_COEFFS.get(channel_name)
-        if coeffs is None and dataset_id.calibration is not None:
-            raise ValueError("No coefficients configured for {}".format(dataset_id))
+        # Calibration values from file, fall back to built-in if unavailable
+        gain = self.nc.attrs['DN_to_Radiance_Gain']
+        offset = self.nc.attrs['DN_to_Radiance_Offset']
+
         if dataset_id.calibration in ('radiance', 'reflectance', 'brightness_temperature'):
-            gain = coeffs[0]
-            offset = coeffs[1]
             data = gain * data + offset
         if dataset_id.calibration == 'reflectance':
             # depends on the radiance calibration above
-            rad_to_alb = coeffs[2]
+            rad_to_alb = self.nc.attrs['Radiance_to_Albedo_c']
             if ds_info.get('units') == '%':
                 rad_to_alb *= 100
             data = data * rad_to_alb
         elif dataset_id.calibration == 'brightness_temperature':
+            data = self._calibrate_ir(dataset_id, data)
+        elif dataset_id.calibration not in ('counts', 'radiance'):
+            raise ValueError("Unknown calibration: '{}'".format(dataset_id.calibration))
+
+        for attr_name in ('standard_name', 'units'):
+            attrs[attr_name] = ds_info[attr_name]
+        attrs.update(dataset_id.to_dict())
+        attrs['orbital_parameters'] = self.get_orbital_parameters()
+        attrs['platform_name'] = self.platform_name
+        attrs['sensor'] = self.sensor
+        data.attrs = attrs
+        return data
+
+    def _calibrate_ir(self, dataset_id, data):
+        """Calibrate radiance data to BTs using either pyspectral or in-file coefficients."""
+        if self.calib_mode == 'PYSPECTRAL':
             # depends on the radiance calibration above
             # Convert um to m^-1 (SI units for pyspectral)
             wn = 1 / (dataset_id.wavelength[1] / 1e6)
@@ -209,14 +204,27 @@ def get_dataset(self, dataset_id, ds_info):
             else:
                 # new versions of pyspectral can do dask arrays
                 data.data = bt_data
-        elif dataset_id.calibration not in ('counts', 'radiance'):
-            raise ValueError("Unknown calibration: '{}'".format(dataset_id.calibration))
-
-        for attr_name in ('standard_name', 'units'):
-            attrs[attr_name] = ds_info[attr_name]
-        attrs.update(dataset_id.to_dict())
-        attrs['orbital_parameters'] = self.get_orbital_parameters()
-        attrs['platform_name'] = self.platform_name
-        attrs['sensor'] = self.sensor
-        data.attrs = attrs
+        else:
+            # IR coefficients from the file
+            # Channel specific
+            c0 = self.nc.attrs['Teff_to_Tbb_c0']
+            c1 = self.nc.attrs['Teff_to_Tbb_c1']
+            c2 = self.nc.attrs['Teff_to_Tbb_c2']
+
+            # These should be fixed, but load anyway
+            cval = self.nc.attrs['light_speed']
+            kval = self.nc.attrs['Boltzmann_constant_k']
+            hval = self.nc.attrs['Plank_constant_h']
+
+            # Compute wavenumber as cm-1
+            wn = (10000 / dataset_id.wavelength[1]) * 100
+
+            # Convert radiance to effective brightness temperature
+            e1 = (2 * hval * cval * cval) * np.power(wn, 3)
+            e2 = (data.data * 1e-5)
+            t_eff = ((hval * cval / kval) * wn) / np.log((e1 / e2) + 1)
+
+            # Now convert to actual brightness temperature
+            bt_data = c0 + c1 * t_eff + c2 * t_eff * t_eff
+            data.data = bt_data
         return data

satpy/tests/reader_tests/test_ami_l1b.py

@@ -115,6 +115,15 @@ def setUp(self, xr_, counts=None):
                 "number_of_lines": 22000,
                 "observation_mode": "FD",
                 "channel_spatial_resolution": "0.5",
+                "Radiance_to_Albedo_c": 1,
+                "DN_to_Radiance_Gain": -0.0144806550815701,
+                "DN_to_Radiance_Offset": 118.050903320312,
+                "Teff_to_Tbb_c0": -0.141418528203155,
+                "Teff_to_Tbb_c1": 1.00052232906885,
+                "Teff_to_Tbb_c2": -0.00000036287276076109,
+                "light_speed": 2.9979245800E+08,
+                "Boltzmann_constant_k": 1.3806488000E-23,
+                "Plank_constant_h": 6.6260695700E-34,
             }
         )
 
@@ -190,7 +199,7 @@ def test_get_area_def(self, adef):
         self.assertEqual(call_args[4], self.reader.nc.attrs['number_of_columns'])
         self.assertEqual(call_args[5], self.reader.nc.attrs['number_of_lines'])
         np.testing.assert_allclose(call_args[6],
-                                   [-5511523.904082, 5511523.904082, 5511022.902, -5511022.902])
+                                   [-5511523.904082, -5511523.904082, 5511022.902, 5511022.902])
 
     def test_get_dataset_vis(self):
         """Test get visible calibrated data."""
@@ -262,19 +271,26 @@ def setUp(self):
     def test_ir_calibrate(self):
         """Test IR calibration."""
         from satpy import DatasetID
-        res = self.reader.get_dataset(
-            DatasetID(name='IR087', wavelength=[8.415, 8.59, 8.765],
-                      calibration='brightness_temperature'),
-            {
-                'file_key': 'image_pixel_values',
-                'wavelength': [8.415, 8.59, 8.765],
-                'standard_name': 'toa_brightness_temperature',
-                'units': 'K',
-            })
-
+        ds_id = DatasetID(name='IR087', wavelength=[8.415, 8.59, 8.765],
+                          calibration='brightness_temperature')
+        ds_info = {
+            'file_key': 'image_pixel_values',
+            'wavelength': [8.415, 8.59, 8.765],
+            'standard_name': 'toa_brightness_temperature',
+            'units': 'K',
+        }
+        res = self.reader.get_dataset(ds_id, ds_info)
         expected = np.array([[238.34385135, 238.31443527, 238.28500087, 238.25554813, 238.22607701],
                              [238.1965875, 238.16707956, 238.13755317, 238.10800829, 238.07844489]])
-        self.assertTrue(np.allclose(res.data.compute(), expected, equal_nan=True))
+        np.testing.assert_allclose(res.data.compute(), expected, equal_nan=True)
+        # make sure the attributes from the file are in the data array
+        self.assertEqual(res.attrs['standard_name'], 'toa_brightness_temperature')
+
+        # test builtin coefficients
+        self.reader.calib_mode = 'FILE'
+        res = self.reader.get_dataset(ds_id, ds_info)
+        # file coefficients are pretty close, give some wiggle room
+        np.testing.assert_allclose(res.data.compute(), expected, equal_nan=True, atol=0.04)
         # make sure the attributes from the file are in the data array
         self.assertEqual(res.attrs['standard_name'], 'toa_brightness_temperature')