diff --git a/bin/seviri2pps.py b/bin/seviri2pps.py index d1e9f64..965acf1 100755 --- a/bin/seviri2pps.py +++ b/bin/seviri2pps.py @@ -49,6 +49,8 @@ help="Output directory where to store level1c file.") parser.add_argument('--no-rotation', action='store_true', help="Don't rotate images") + parser.add_argument('--azimuth_angles', action='store_true', + help="Save azimuth angles") parser.add_argument('-ne', '--nc_engine', type=str, nargs='?', required=False, default='h5netcdf', help="Engine for saving netcdf files netcdf4 or h5netcdf (default).") @@ -64,5 +66,6 @@ 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 + apply_sun_earth_distance_correction=not options.no_sun_earth_distance_correction, + save_azimuth_angles=options.azimuth_angles, ) diff --git a/level1c4pps/seviri2pps_lib.py b/level1c4pps/seviri2pps_lib.py index 6c8707b..7ed0557 100644 --- a/level1c4pps/seviri2pps_lib.py +++ b/level1c4pps/seviri2pps_lib.py @@ -137,7 +137,11 @@ def _create_scene(file_format, filenames, calib_coefs): def _check_is_seviri_data(scene): - if not scene.attrs['sensor'] == {'seviri'}: + if hasattr(scene, 'sensor_names') and 'seviri' in scene.sensor_names: + pass + elif scene.attrs['sensor'] == {'seviri'}: + pass + else: raise ValueError('Not SEVIRI data') @@ -322,6 +326,8 @@ def update_coords(scene): def add_ancillary_datasets(scene, lons, lats, sunz, satz, azidiff, + suna, sata, + save_azimuth_angles=False, chunks=(512, 3712)): """Add ancillary datasets to the scene. @@ -375,6 +381,18 @@ def add_ancillary_datasets(scene, lons, lats, sunz, satz, azidiff, da.from_array(azidiff[:, :], chunks=chunks), dims=['y', 'x'], coords=angle_coords) + # Sunazimuth + if save_azimuth_angles: + scene['sunazimuth'] = xr.DataArray( + da.from_array(suna[:, :], chunks=chunks), + dims=['y', 'x'], coords=angle_coords) + + # Satazimuth + if save_azimuth_angles: + scene['satazimuth'] = xr.DataArray( + da.from_array(sata[:, :], chunks=chunks), + dims=['y', 'x'], coords=angle_coords) + # Update the attributes update_angle_attributes(scene, band=scene['IR_108']) @@ -470,7 +488,7 @@ def get_encoding_seviri(scene): encoding['time'] = {'units': 'days since 2004-01-01 00:00', 'calendar': 'standard', '_FillValue': None, - 'chunksizes': [1]} + 'chunksizes': (1,)} return encoding @@ -564,7 +582,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, - clip_calib=False): + clip_calib=False, + save_azimuth_angles=False): """Make level 1c files in PPS-format.""" for fname in tslot_files: if not os.path.isfile(fname): @@ -577,7 +596,9 @@ def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf', rotate=rotate, clip_calib=clip_calib ) - + if hasattr(scn_, 'start_time'): + scn_.attrs['start_time'] = scn_.start_time + scn_.attrs['end_time'] = scn_.end_time # Find lat/lon data lons, lats = get_lonlats(scn_['IR_108']) @@ -590,8 +611,12 @@ def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf', update_coords(scn_) # Add ancillary datasets to the scene - add_ancillary_datasets(scn_, lons=lons, lats=lats, sunz=sunz, satz=satz, - azidiff=azidiff) + add_ancillary_datasets(scn_, + lons=lons, lats=lats, + sunz=sunz, satz=satz, + azidiff=azidiff, + suna=suna, sata=sata, + save_azimuth_angles=save_azimuth_angles) add_proj_satpos(scn_) # Set attributes. This changes SEVIRI band names to PPS band names. @@ -601,6 +626,7 @@ def process_one_scan(tslot_files, out_path, rotate=True, engine='h5netcdf', ir108_for_filename = scn_['IR_108'] if use_nominal_time_in_filename: ir108_for_filename = set_nominal_scan_time(ir108_for_filename) + filename = compose_filename( scene=scn_, out_path=out_path, diff --git a/level1c4pps/tests/test_seviri2pps.py b/level1c4pps/tests/test_seviri2pps.py index 49f9d06..e55826e 100644 --- a/level1c4pps/tests/test_seviri2pps.py +++ b/level1c4pps/tests/test_seviri2pps.py @@ -42,8 +42,8 @@ def get_fake_scene(): scene = Scene() start_time = dt.datetime(2020, 1, 1, 12) scene['VIS006'] = xr.DataArray( - [[1, 2], - [3, 4]], + [[1.0, 2.0], + [3.0, 4.0]], dims=('y', 'x'), attrs={'calibration': 'reflectance', 'sun_earth_distance_correction_applied': True, @@ -299,6 +299,8 @@ def test_add_ancillary_datasets(self): lats = np.array([[1.1, 2.1], [3.1, 4.1]]) sunz = np.array([[1.2, 2.2], [3.2, 4.2]]) satz = np.array([[1.3, 2.3], [3.3, 4.3]]) + suna = np.array([[5.2, 2.2], [5.2, 1.2]]) + sata = np.array([[3.3, 2.3], [3.3, 7.3]]) azidiff = np.array([[1.4, 2.4], [3.4, 4.4]]) ir_108 = xr.DataArray(data=np.array([[0.1, 0.2], [0.3, 0.4]]), @@ -310,9 +312,12 @@ def test_add_ancillary_datasets(self): 'orbital_parameters': 'orb_params', 'georef_offset_corrected': True}) scene = {'IR_108': ir_108} - seviri2pps.add_ancillary_datasets(scene, lons=lons, lats=lats, + seviri2pps.add_ancillary_datasets(scene, + lons=lons, lats=lats, sunz=sunz, satz=satz, - azidiff=azidiff) + azidiff=azidiff, + suna=suna, sata=sata, + save_azimuth_angles=True) # Test lon/lat np.testing.assert_array_equal(scene['lon'].data, lons) @@ -331,6 +336,12 @@ def test_add_ancillary_datasets(self): np.testing.assert_array_equal(scene['azimuthdiff'].data, azidiff) self.assertEqual(scene['azimuthdiff'].attrs['name'], 'azimuthdiff') + np.testing.assert_array_equal(scene['satazimuth'].data, sata) + self.assertEqual(scene['satazimuth'].attrs['name'], 'satazimuth') + + np.testing.assert_array_equal(scene['sunazimuth'].data, suna) + self.assertEqual(scene['sunazimuth'].attrs['name'], 'sunazimuth') + for angle in ['azimuthdiff', 'satzenith', 'sunzenith']: self.assertTupleEqual(scene[angle].dims, ('y', 'x')) np.testing.assert_array_equal(scene[angle].coords['x'].data, xvals) @@ -398,7 +409,7 @@ def test_get_encoding(self): enc_exp_time = {'units': 'days since 2004-01-01 00:00', 'calendar': 'standard', '_FillValue': None, - 'chunksizes': [1]} + 'chunksizes': (1,)} enc_exp_acq = {'units': 'milliseconds since 2009-07-01 12:15', 'calendar': 'standard', '_FillValue': -9999.0} @@ -426,7 +437,11 @@ def test_get_encoding(self): 'acq_time': enc_exp_acq } encoding = seviri2pps.get_encoding_seviri(scene) - self.assertDictEqual(encoding, encoding_exp) + for key in encoding_exp: + print(key) + print(encoding[key],encoding_exp[key]) + + self.assertDictEqual(encoding[key], encoding_exp[key]) def test_get_header_attrs(self): """Test get the header attributes."""