Merge pull request #99 from ninahakansson/fix_earth_distance_factor_a…

…ttribute

Add fix for sun_earth_distance_correction_factor

bin/seviri2pps.py

@@ -56,16 +56,12 @@
                         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('--no-sun-earth-distance-correction',
-                        action='store_true',
-                        help='Do not apply sun earth distance correction.')
     options = parser.parse_args()
     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,
-        apply_sun_earth_distance_correction=not options.no_sun_earth_distance_correction,
         save_azimuth_angles=options.azimuth_angles,
     )

level1c4pps/__init__.py

@@ -368,6 +368,21 @@ def adjust_lons_to_valid_range(scene):
     # scene['lon'] = centered_modulus(scene['lon']) # makes lon loose attrs satpy 0.24.0
     scene['lon'].values = centered_modulus(scene['lon'].values)
 
+def fix_sun_earth_distance_correction_factor(scene, band, start_time):
+    from pyorbital.astronomy import sun_earth_distance_correction
+    date_control = np.datetime64("2019-01-01T00:00:00")
+    sun_earth_distance_20190409 = sun_earth_distance_correction(date_control)
+    sun_earth_distance = sun_earth_distance_correction(start_time)
+    if (np.abs(sun_earth_distance_20190409 - 0.9833280675966011) < 0.00001 and
+        np.abs(sun_earth_distance - scene[band].attrs['sun_earth_distance_correction_factor']) < 0.00001):
+        logger.info("The sun earth distance correction attribute contain the sun earth distance, not the square.")
+        logger.info("Updating and adding sun earth distance correction attributes.")
+        current_factor = scene[band].attrs['sun_earth_distance_correction_factor']
+        scene[band].attrs['satpy_sun_earth_distance_correction_factor'] =  current_factor
+        scene[band].attrs['pps_sun_earth_distance_correction_factor'] = sun_earth_distance * sun_earth_distance
+        scene[band].attrs['sun_earth_distance'] = sun_earth_distance
+        scene[band].attrs['sun_earth_distance_correction_factor'] = sun_earth_distance * sun_earth_distance    
+
 
 def set_header_and_band_attrs_defaults(scene, BANDNAMES, PPS_TAGNAMES, REFL_BANDS, irch, orbit_n=0):
     """Add some default values for band attributes."""
@@ -422,6 +437,7 @@ def set_header_and_band_attrs_defaults(scene, BANDNAMES, PPS_TAGNAMES, REFL_BAND
         else:
             # Assume factor applied if available as attribute.
             scene[band].attrs['sun_earth_distance_correction_applied'] = 'True'
+            fix_sun_earth_distance_correction_factor(scene, band, irch.attrs['start_time'])
         scene[band].attrs['wavelength'] = scene[band].attrs['wavelength'][0:3]
         scene[band].attrs['sun_zenith_angle_correction_applied'] = 'False'
         if idtag in PPS_TAGNAMES_TO_IMAGE_NR:

level1c4pps/seviri2pps_lib.py

@@ -38,13 +38,16 @@
 from datetime import datetime, timedelta
 from satpy.scene import Scene
 import satpy.utils
-from satpy.readers.utils import remove_earthsun_distance_correction
 from trollsift.parser import globify, Parser
 from pyorbital.astronomy import get_alt_az, sun_zenith_angle
 from pyorbital.orbital import get_observer_look
 
 from level1c4pps.calibration_coefs import get_calibration, CalibrationData
-from level1c4pps import make_azidiff_angle, get_encoding, compose_filename, update_angle_attributes
+from level1c4pps import (make_azidiff_angle,
+                         get_encoding,
+                         compose_filename,
+                         update_angle_attributes,
+                         fix_sun_earth_distance_correction_factor)
 
 
 try:
@@ -89,16 +92,14 @@ 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, rotate,
                        clip_calib):
     """Load and calibrate data.
 
     Uses inter-calibration coefficients from Meirink et al.
 
     Args:
         filenames: List of data files
-        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
@@ -121,9 +122,6 @@ def load_and_calibrate(filenames, apply_sun_earth_distance_correction, rotate,
     _load_bands(scn_, rotate)
     _update_scene_attrs(scn_, {'image_rotated': rotate})
 
-    if not apply_sun_earth_distance_correction:
-        remove_sun_earth_distance_correction(scn_)
-
     return scn_
 
 
@@ -160,22 +158,6 @@ def _update_scene_attrs(scene, attrs):
     scene.attrs.update(attrs)
 
 
-def remove_sun_earth_distance_correction(scene):
-    """Remove sun-earth-distance correction from visible channels.
-
-    This is required because a downstream CLAAS module (CPP) does not recognize
-    the "sun_earth_distance_correction_applied" attribute and applies the
-    correction anyway.
-    """
-    for band in BANDNAMES:
-        is_vis = scene[band].attrs['calibration'] == 'reflectance'
-        correction_applied = scene[band].attrs.get(
-            'sun_earth_distance_correction_applied', False
-        )
-        if is_vis and correction_applied:
-            scene[band] = remove_earthsun_distance_correction(scene[band])
-
-
 def get_lonlats(dataset):
     """Get lat/lon coordinates."""
     lons, lats = dataset.attrs['area'].get_lonlats()
@@ -274,10 +256,15 @@ def set_attrs(scene):
         scene[band].attrs['wavelength'] = [scene[band].attrs['wavelength'].min,
                                            scene[band].attrs['wavelength'].central,
                                            scene[band].attrs['wavelength'].max]
+
         if 'sun_earth_distance_correction_factor' not in scene[band].attrs:
             scene[band].attrs['sun_earth_distance_correction_applied'] = False
             scene[band].attrs['sun_earth_distance_correction_factor'] = 1.0
+        else:
+            fix_sun_earth_distance_correction_factor(scene, band, scene[band].attrs['start_time'])
+
         scene[band].attrs['sun_zenith_angle_correction_applied'] = False
+
         scene[band].attrs['name'] = "image{:d}".format(image_num)
 
         # Cosmetics
@@ -581,7 +568,6 @@ 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,
                      clip_calib=False,
                      save_azimuth_angles=False):
     """Make level 1c files in PPS-format."""
@@ -592,7 +578,6 @@ def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf',
     tic = time.time()
     scn_ = load_and_calibrate(
         tslot_files,
-        apply_sun_earth_distance_correction=apply_sun_earth_distance_correction,
         rotate=rotate,
         clip_calib=clip_calib
     )

level1c4pps/tests/test_seviri2pps.py

@@ -75,17 +75,11 @@ def test_load_and_calibrate(self, mocked_scene):
         filenames = ['MSG4-SEVI-MSG15-1234-NA-20190409121243.927000000Z']
         res = seviri2pps.load_and_calibrate(
             filenames,
-            apply_sun_earth_distance_correction=False,
             rotate=False,
             clip_calib=False
         )
 
         # Compare results and expectations
-        vis006_exp_pyspectral_1_7_1 = xr.DataArray(
-            [[1.07025268, 2.14050537],
-             [3.21075805, 4.28101074]],
-            dims=('y', 'x')
-        )
         vis006_exp = xr.DataArray(
             [[1.034205, 2.06841],
              [3.102615, 4.13682]],
@@ -96,12 +90,11 @@ def test_load_and_calibrate(self, mocked_scene):
              [7, 8]],
             dims=('y', 'x')
         )
-        if res['VIS006'][0, 0] > 1.04:
-            # pyspectral 1.7.1 and older
-            xr.testing.assert_allclose(res['VIS006'], vis006_exp_pyspectral_1_7_1)
-        else:
-            xr.testing.assert_allclose(res['VIS006'], vis006_exp)
+        from satpy.readers.utils import remove_earthsun_distance_correction
+        res['VIS006'] = remove_earthsun_distance_correction(res['VIS006'])
+        xr.testing.assert_allclose(res['VIS006'], vis006_exp)
         xr.testing.assert_equal(res['IR_108'], ir_108_exp)
+
         self.assertFalse(
             res['VIS006'].attrs['sun_earth_distance_correction_applied'],
         )
@@ -114,7 +107,6 @@ def test_load_and_calibrate_with_rotation(self, mocked_scene):
         filenames = ['MSG4-SEVI-MSG15-1234-NA-20190409121243.927000000Z']
         res = seviri2pps.load_and_calibrate(
             filenames,
-            apply_sun_earth_distance_correction=False,
             rotate=True,
             clip_calib=False
         )
@@ -211,15 +203,14 @@ def get_observer_look_patched(lon, lat, alt, *args):
     def test_set_attrs(self):
         """Test setting scene attributes."""
         seviri2pps.BANDNAMES = ['VIS006', 'IR_108']
-        vis006 = mock.MagicMock(attrs={'wavelength': WavelengthRange(0.56, 0.635, 0.71)})
-        ir108 = mock.MagicMock(attrs={'platform_name': 'myplatform',
-                                      'wavelength': WavelengthRange(9.8, 10.8, 11.8),
-                                      'orbital_parameters': {'orb_a': 1,
-                                                             'orb_b': 2},
-                                      'georef_offset_corrected': True})
-        scene_dict = {'VIS006': vis006, 'IR_108': ir108}
-        scene = mock.MagicMock(attrs={})
-        scene.__getitem__.side_effect = scene_dict.__getitem__
+        scene = get_fake_scene()
+        scene['VIS006'].attrs['sun_earth_distance_correction_applied'] = True
+        scene['VIS006'].attrs['start_time'] = dt.datetime(2020, 1, 1, 12)
+        scene['VIS006'].attrs['sun_earth_distance_correction_factor'] = 0.9833241577909706
+        scene['IR_108'].attrs['orbital_parameters'] = {'orb_a': 1,
+                                                       'orb_b': 2}
+        scene['IR_108'].attrs['georef_offset_corrected'] =  True
+        scene['IR_108'].attrs['platform_name'] = "my_platform_name"
 
         seviri2pps.set_attrs(scene)
         self.assertEqual(scene['VIS006'].attrs['name'], 'image0')
@@ -229,7 +220,18 @@ def test_set_attrs(self):
         self.assertNotIn('orb_a', scene.attrs)
         self.assertNotIn('orbital_parameters', scene.attrs)
         self.assertNotIn('georef_offset_corrected', scene.attrs)
-
+        exp_sun_earth_distance_correction_factor = 0.9669263992953216
+        sun_earth_distance = 0.9833241577909706
+        self.assertAlmostEqual(scene['VIS006'].sun_earth_distance_correction_factor,
+                               exp_sun_earth_distance_correction_factor, places=7)
+        self.assertAlmostEqual(scene['VIS006'].pps_sun_earth_distance_correction_factor,
+                               exp_sun_earth_distance_correction_factor, places=7)
+        self.assertAlmostEqual(scene['VIS006'].satpy_sun_earth_distance_correction_factor,
+                               sun_earth_distance, places=7)
+        self.assertAlmostEqual(scene['VIS006'].sun_earth_distance,
+                               sun_earth_distance, places=7)
+
+
     def test_get_mean_acq_time(self):
         """Test computation of mean scanline acquisition time."""
         seviri2pps.BANDNAMES = ['VIS006', 'IR_108']