Add eum L1g

bin/isccpng2pps.py

@@ -49,7 +49,10 @@
     parser.add_argument('-on', '--orbit_number', type=int, nargs='?',
                         required=False, default=99999,
                         help="Orbit number (default is 99999).")
+    parser.add_argument('--use-solar-angles-from-file', action='store_true',
+                        help='Use solar angles from file.')
     options = parser.parse_args()
     process_one_scene(options.files, options.out_dir,
                       engine=options.nc_engine,
+                      use_solar_angles_from_file=options.use_solar_angles_from_file,
                       orbit_n=options.orbit_number)

level1c4pps/__init__.py

@@ -347,9 +347,22 @@ def rename_latitude_longitude(scene):
     del scene[lat_name_satpy]
     del scene[lon_name_satpy]
     # Update attributes
+    update_lat_lon_attrs(scene)
+
+
+def update_lat_lon_attrs(scene):
+    """Update lat/lon attributes."""
     scene['lat'].attrs = LATLON_ATTRIBUTES['lat']
     scene['lon'].attrs = LATLON_ATTRIBUTES['lon']
-    for coord_name in ['acq_time', 'm_latitude', 'i_latitude', 'm_latitude', 'i_latitude', 'latitude', 'longitude']:
+    for coord_name in [
+        'acq_time',
+        'm_latitude',
+        'i_latitude',
+        'm_latitude',
+        'i_latitude',
+        'latitude',
+        'longitude',
+            'crs']:
         try:
             del scene['lat'].coords[coord_name]
             del scene['lon'].coords[coord_name]
@@ -506,7 +519,7 @@ def update_angle_attributes(scene, band):
             except (AttributeError, KeyError):
                 pass
         # delete some coords
-        for coord_name in ['acq_time', 'latitude', 'longitude']:
+        for coord_name in ['acq_time', 'latitude', 'longitude', 'crs']:
             try:
                 del scene[angle].coords[coord_name]
             except KeyError:

level1c4pps/isccpng2pps_lib.py

@@ -22,23 +22,27 @@
 
 """Functions to convert ISCCP Next Generation level-1-G data to a NWCSAF/PPS level-1c formatet netCDF/CF file."""
 
+from level1c4pps import (adjust_lons_to_valid_range, apply_sunz_correction,
+                         check_file_exists, compose_filename, convert_angles,
+                         dt64_to_datetime, get_header_attrs, get_refl_bands,
+                         log_time, save_data,
+                         set_header_and_band_attrs_defaults,
+                         update_angle_attributes, update_lat_lon_attrs)
 import logging
+import os
 import time
 
 import numpy as np
 from pyorbital.astronomy import get_alt_az, sun_zenith_angle
 from satpy.scene import Scene
+from satpy.utils import debug_on
+import xarray as xr
+debug_on()
 
-from level1c4pps import (adjust_lons_to_valid_range, apply_sunz_correction,
-                         check_file_exists, compose_filename, convert_angles,
-                         dt64_to_datetime, get_header_attrs, get_refl_bands,
-                         log_time, save_data,
-                         set_header_and_band_attrs_defaults,
-                         update_angle_attributes)
 
 logger = logging.getLogger('isccpng2pps')
 
-GEOLOCATION_NAMES = [
+GEOLOCATION_NAMES_ISCCP_NG_DEMO = [
     'solar_zenith_angle',
     'satellite_zenith_angle',
     'solar_azimuth_angle',
@@ -48,6 +52,16 @@
     "lon",
     "lat"]
 
+GEOLOCATION_NAMES_EUM = [
+    'solar_zenith_angle',
+    'sensor_zenith_angle',
+    'solar_azimuth_angle',
+    'sensor_azimuth_angle',
+    "sensor_flag",
+    "pixel_time",
+    "lon",
+    "lat"]
+
 PPS_TAGS = {'refl_01_60um': "ch_r16",
             'refl_00_65um': "ch_r06",
             'refl_00_86um': "ch_r09",
@@ -147,8 +161,20 @@
                    55: "Meteosat-8",  # MSG1 SEVIRI
                    70: "Meteosat-11"}
 
+sensor_flag_dict = {1: "GOES-16",  # ABI
+                    2: "GOES-17",  # ABI
+                    3: "GOES-18",  # ABI
+                    4: "GOES-19",  # ABI
+                    5: "Himawari-8",  # AHI
+                    6: "Himawari-9",  # AHI
+                    7: "Meteosat-8",  # MSG1 SEVIRI
+                    8: "Meteosat-9",
+                    9: "Meteosat-10",
+                    10: "Meteosat-11",
+                    12: "Meteosat-12", }
+
 
-def set_header_and_band_attrs(scene, orbit_n=00000):
+def set_header_and_band_attrs(scene, orbit_n=00000, is_eum=True):
     """Set and delete some attributes."""
     nimg = 0  # name of first dataset is image0
     # Set some header attributes:
@@ -160,7 +186,11 @@ def set_header_and_band_attrs(scene, orbit_n=00000):
     for band in refl_bands:
         if band not in scene:
             continue
-        scene[band].attrs['sun_zenith_angle_correction_applied'] = 'False'
+        if is_eum:
+            scene[band].attrs['sun_zenith_angle_correction_applied'] = 'True'
+        else:
+            scene[band].attrs['sun_zenith_angle_correction_applied'] = 'False'
+
     return nimg
 
 
@@ -223,25 +253,31 @@ def get_solar_angles(scene, lons, lats):
     Returns:
         Solar azimuth angle, Solar zenith angle in degrees
     """
-    acq_time =  scene["pixel_time"].copy()
+    acq_time = scene["pixel_time"].copy()
     _, suna = get_alt_az(acq_time, lons, lats)
     suna = np.rad2deg(suna)
     sunz = sun_zenith_angle(acq_time, lons, lats)
     return suna, sunz
 
 
-def fix_pixel_time(scene):
+def fix_pixel_time(scene, is_eum):
     """Fix the time pixel variable, original file does not contain units."""
     del scene["pixel_time"].coords["crs"]
     scene["pixel_time"].encoding['coordinates'] = "lon lat"
-    scene["pixel_time"] = scene["pixel_time"].interpolate_na(dim = "y", fill_value="extrapolate", use_coordinate=False)  # update NaTs
-    scene["pixel_time"].data = scene["pixel_time"].data * np.timedelta64(1, 's') + scene['temp_11_00um'].attrs["start_time"]
+    if not is_eum:
+        scene["pixel_time"] = scene["pixel_time"].interpolate_na(
+            dim="y", fill_value="extrapolate", use_coordinate=False)  # update NaTs
+        scene["pixel_time"].data = scene["pixel_time"].data * \
+            np.timedelta64(1, 's') + scene['temp_11_00um'].attrs["start_time"]
 
 
-def load_data(scene_files):
+def load_data(scene_files, is_eum):
     """Load data."""
     scene = Scene(reader='multiple_sensors_isccpng_l1g_nc', filenames=scene_files)
-    bands_to_load = band_names + GEOLOCATION_NAMES
+    if is_eum:
+        bands_to_load = band_names + GEOLOCATION_NAMES_EUM
+    else:
+        bands_to_load = band_names + GEOLOCATION_NAMES_ISCCP_NG_DEMO
     scene.load(bands_to_load)
     return scene
 
@@ -253,26 +289,54 @@ def update_solar_angles(scene):
     scene["solar_azimuth_angle"].values = suna.values
 
 
+def get_wmo_id_from_sensor_flag(scene):
+    """Get wmo_id used in ISSCCPNG-demo format from EUM sensor_flag."""
+    if "wmo_id" in scene:
+        return
+    del scene["sensor_flag"].coords["crs"]
+    scene["wmo_id"] = scene["sensor_flag"].copy()
+    wmo_id_array = np.empty(scene["sensor_flag"].values.shape)
+    for index, name in sensor_flag_dict.items():
+        for wmo_id, name2 in satellite_names.items():
+            if name == name2:
+                logger.info(f"For {name} setting using wmo_id {wmo_id} for channel sensor_flag {index}.")
+                set_these = scene["sensor_flag"] == index
+                wmo_id_array[set_these] = wmo_id
+    scene["wmo_id"].values = wmo_id_array
+
+
+def check_if_eum_format(scene_files):
+    """Check if we have ISCCP-NG or EUM format."""
+    for filename in scene_files:
+        if "EUM_L1g" in os.path.basename(filename):
+            return True
+
+
 def process_one_scene(scene_files, out_path,
                       engine='h5netcdf',
+                      use_solar_angles_from_file=False,
                       orbit_n=0):
     """Make level 1c files in PPS-format."""
     tic = time.time()
     check_file_exists(scene_files)
-    scene = load_data(scene_files)
+    is_eum = check_if_eum_format(scene_files)
+    scene = load_data(scene_files, is_eum)
     ir_channel_obj = scene[ONE_IR_CHANNEL]
-    set_header_and_band_attrs(scene, orbit_n=orbit_n)
-    fix_pixel_time(scene)
-    # rename_latitude_longitude(scene)
-    update_solar_angles(scene)
+    set_header_and_band_attrs(scene, orbit_n=orbit_n, is_eum=is_eum)
+    get_wmo_id_from_sensor_flag(scene)
+    fix_pixel_time(scene, is_eum)
+    update_lat_lon_attrs(scene)
+    if not use_solar_angles_from_file:
+        update_solar_angles(scene)
     adjust_lons_to_valid_range(scene)
     convert_angles(scene, delete_azimuth=True)
     update_angle_attributes(scene, ir_channel_obj)
-    recalibrate_meteosat(scene)
+    if not is_eum:
+        recalibrate_meteosat(scene)
     homogenize(scene)
     apply_sunz_correction(scene, refl_bands)
-    filename = compose_filename(scene, out_path, instrument='seviri', band=ir_channel_obj)
-    header_attrs = get_header_attrs(scene, band=ir_channel_obj, sensor='seviri')
+    filename = compose_filename(scene, out_path, instrument='seviri', band=scene["pixel_time"])
+    header_attrs = get_header_attrs(scene, band=scene["pixel_time"], sensor='seviri')
     save_data(scene, filename, header_attrs, engine)
     log_time(filename, tic)
     return filename

level1c4pps/tests/test_isccpng.py

@@ -37,7 +37,7 @@ def setUp(self):
         self.scene = Scene()
         scene_dict = {}
         grid_data = [[1.0, 2.0], [3.0, 4.0]]
-        all_keys = ["refl_00_65um", 'temp_11_00um'] + isccpng2pps.GEOLOCATION_NAMES
+        all_keys = ["refl_00_65um", 'temp_11_00um'] + isccpng2pps.GEOLOCATION_NAMES_ISCCP_NG_DEMO
         for key in all_keys:
             scene_dict[key] = xr.DataArray(grid_data,
                                            dims=('y', 'x'),
@@ -58,6 +58,8 @@ def setUp(self):
             'platform_name': '',
             'orbit_number': 99999}
         scene_dict["pixel_time"].coords["crs"] = ""
+        scene_dict['pixel_time'].attrs = {'start_time': np.datetime64('2021-06-28T01:00:00'),
+                                          'end_time': np.datetime64('2021-06-28T01:01:00')}
         for key in scene_dict:
             self.scene[key] = scene_dict[key]
         self.scene.load = mock.MagicMock()
@@ -68,5 +70,5 @@ def setUp(self):
     def test_process_one_scene(self, mock_scene, mock_check_file_exists):
         """Test to set process_one_scene."""
         mock_scene.return_value = self.scene
-        filename = isccpng2pps.process_one_scene("dummy", out_path='./level1c4pps/tests/', orbit_n='12345')
+        filename = isccpng2pps.process_one_scene(["ISCCP-NG"], out_path='./level1c4pps/tests/', orbit_n='12345')
         self.assertEqual(os.path.basename(filename), "S_NWC_seviri__12345_20210628T0100000Z_20210628T0101000Z.nc")