From be904521c041ede7c4d4e90d87e8f089eea123fe Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" <6160529+ninahakansson@users.noreply.github.com> Date: Tue, 3 Dec 2024 13:25:37 +0100 Subject: [PATCH 1/4] Add fix for sun_earth_distance_correction_factor The history of the sun earth distance correction is in pygac/satpy and pyorbital. With some combination of satpy and pyorbital (2020-2022) the correction where in practice applied twice. And since 2022 the attribute sun_earth_distance_correction_factor from satpy no longer includes the factor for correction but sqrt(the original correction). The history of the correction in pytroll: Original formula in pygac and pyorbital from 2014: corr_ = 1 - 0.0334 * np.cos(2 * np.pi * (jdays2000(utc_time) - 2) / year) Calulcated in pygac and applied like refl = refl * corr_. With attribute set in pygac: sun_earth_distance_correction_factor = corr_ 2020-06 Pygac correction depricated. 2020-09 Satpy reads corrrection from pyorbital but applies like refl = refl * corr_ * corr_ causing the correction to be applied twice. 2020-09 Attribute set in satpy: sun_earth_distance_correction_factor = corr_, same as before. 2022-06 Pyorbital correction in principle replaced with sqrt(corr_). This implies reflections now corrected with refl = refl * corr_ in satpy as expected. The sun_earth_distance_correction_factor attribute in principle set to np.sqrt(corr_) which means the attribute no longer contain the original sun_earth_distance_correction_factor attribute. --- level1c4pps/__init__.py | 16 +++++++++++ level1c4pps/seviri2pps_lib.py | 11 +++++++- level1c4pps/tests/test_seviri2pps.py | 42 +++++++++++++++------------- 3 files changed, 48 insertions(+), 21 deletions(-) diff --git a/level1c4pps/__init__.py b/level1c4pps/__init__.py index 94fd99e..f4865a7 100644 --- a/level1c4pps/__init__.py +++ b/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: diff --git a/level1c4pps/seviri2pps_lib.py b/level1c4pps/seviri2pps_lib.py index 7ed0557..2423537 100644 --- a/level1c4pps/seviri2pps_lib.py +++ b/level1c4pps/seviri2pps_lib.py @@ -44,7 +44,11 @@ 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: @@ -274,10 +278,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 diff --git a/level1c4pps/tests/test_seviri2pps.py b/level1c4pps/tests/test_seviri2pps.py index e55826e..307e5a6 100644 --- a/level1c4pps/tests/test_seviri2pps.py +++ b/level1c4pps/tests/test_seviri2pps.py @@ -81,11 +81,6 @@ def test_load_and_calibrate(self, mocked_scene): ) # 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,11 +91,8 @@ 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) + res['VIS006'].sun_earth_distance_correction_factor + 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'], @@ -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'] From a2dd7429c6d8304604cdcd2323abcca89d1edd1b Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" <6160529+ninahakansson@users.noreply.github.com> Date: Tue, 3 Dec 2024 13:32:41 +0100 Subject: [PATCH 2/4] Remove option to undo sun-earth-distance-correction --- bin/seviri2pps.py | 3 --- level1c4pps/seviri2pps_lib.py | 26 +------------------------- level1c4pps/tests/test_seviri2pps.py | 9 +++------ 3 files changed, 4 insertions(+), 34 deletions(-) diff --git a/bin/seviri2pps.py b/bin/seviri2pps.py index 965acf1..9bf1b48 100755 --- a/bin/seviri2pps.py +++ b/bin/seviri2pps.py @@ -56,9 +56,6 @@ 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, diff --git a/level1c4pps/seviri2pps_lib.py b/level1c4pps/seviri2pps_lib.py index 2423537..8d402eb 100644 --- a/level1c4pps/seviri2pps_lib.py +++ b/level1c4pps/seviri2pps_lib.py @@ -38,7 +38,6 @@ 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 @@ -93,7 +92,7 @@ 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. @@ -101,8 +100,6 @@ def load_and_calibrate(filenames, apply_sun_earth_distance_correction, rotate, 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 @@ -125,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_ @@ -164,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() @@ -590,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.""" @@ -601,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 ) diff --git a/level1c4pps/tests/test_seviri2pps.py b/level1c4pps/tests/test_seviri2pps.py index 307e5a6..5629330 100644 --- a/level1c4pps/tests/test_seviri2pps.py +++ b/level1c4pps/tests/test_seviri2pps.py @@ -75,15 +75,14 @@ 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 = xr.DataArray( - [[1.034205, 2.06841], - [3.102615, 4.13682]], + [[1.0, 2.0], + [3.0, 4.0]], dims=('y', 'x') ) ir_108_exp = xr.DataArray( @@ -91,10 +90,9 @@ def test_load_and_calibrate(self, mocked_scene): [7, 8]], dims=('y', 'x') ) - res['VIS006'].sun_earth_distance_correction_factor xr.testing.assert_allclose(res['VIS006'], vis006_exp) xr.testing.assert_equal(res['IR_108'], ir_108_exp) - self.assertFalse( + self.assertTrue( res['VIS006'].attrs['sun_earth_distance_correction_applied'], ) @@ -106,7 +104,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 ) From 29a44c00ef5ef46fffb9f85e90a122728316e5c3 Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" <6160529+ninahakansson@users.noreply.github.com> Date: Tue, 3 Dec 2024 14:53:25 +0100 Subject: [PATCH 3/4] Lets keep the check of the sun-earth-distance-correction This test makes sure satpy is applying and removing the correction correctly. --- level1c4pps/tests/test_seviri2pps.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/level1c4pps/tests/test_seviri2pps.py b/level1c4pps/tests/test_seviri2pps.py index 5629330..337903f 100644 --- a/level1c4pps/tests/test_seviri2pps.py +++ b/level1c4pps/tests/test_seviri2pps.py @@ -81,8 +81,8 @@ def test_load_and_calibrate(self, mocked_scene): # Compare results and expectations vis006_exp = xr.DataArray( - [[1.0, 2.0], - [3.0, 4.0]], + [[1.034205, 2.06841], + [3.102615, 4.13682]], dims=('y', 'x') ) ir_108_exp = xr.DataArray( @@ -90,9 +90,12 @@ def test_load_and_calibrate(self, mocked_scene): [7, 8]], dims=('y', 'x') ) + 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.assertTrue( + + self.assertFalse( res['VIS006'].attrs['sun_earth_distance_correction_applied'], ) From 60733b2c3b0da6cf3644a40b2f2226849766ad49 Mon Sep 17 00:00:00 2001 From: "Nina.Hakansson" <6160529+ninahakansson@users.noreply.github.com> Date: Tue, 3 Dec 2024 15:30:59 +0100 Subject: [PATCH 4/4] small fix --- bin/seviri2pps.py | 1 - 1 file changed, 1 deletion(-) diff --git a/bin/seviri2pps.py b/bin/seviri2pps.py index 9bf1b48..849a054 100755 --- a/bin/seviri2pps.py +++ b/bin/seviri2pps.py @@ -63,6 +63,5 @@ 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, )