Merge pull request #80 from ninahakansson/seviri_azimuth_angles

Add azimuth angles for SEVIRI

Author: ninahakansson

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,
     )

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,

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."""