adding support for handling merci3 data by renaming and extending the

merci2 script

Author: BengtRydberg

README.md

@@ -10,11 +10,11 @@ agency specific level-1 formats. So far, supports
 
  - EUMETSAT Meteosat Second Generation SEVIRI HRIT level-1.5
  - NOAA AVHRR GAC
- - MERSI-2 level-1
+ - MERSI-2/3 level-1
  - MODIS level-1
  - AVHRR eps level-1b
  - AVHRR AAPP level-1b
  - VIIRS level-1b
  - SLSTR level-1b
  - EPS-SG MetImage level-1 test data
- - EUMETSAT AVHRR GAC FDR
+ - EUMETSAT AVHRR GAC FDR

bin/mersi22pps.py bin/mersi2pps.py

@@ -25,17 +25,17 @@
 """Script to convert MERSI-2 level-1 to PPS level-1c format using Pytroll/Satpy."""
 
 import argparse
-from level1c4pps.mersi22pps_lib import process_one_scene
+from level1c4pps.mersi2pps_lib import process_one_scene
 
 
 if __name__ == "__main__":
     """ Create PPS-format level1c data
-    From a list of MERSI-2 level-1 files create a NWCSAF/PPS formatet level1c file for pps.
+    From a list of MERSI-2/3 level-1 files create a NWCSAF/PPS formatet level1c file for pps.
     """
     parser = argparse.ArgumentParser(
-        description=('Script to produce a PPS-level1c file for a MERSI-2 level-1 scene'))
+        description=('Script to produce a PPS-level1c file for a MERSI-2/3 level-1 scene'))
     parser.add_argument('files', metavar='fileN', type=str, nargs='+',
-                        help='List of MERSI-2 files to process')
+                        help='List of MERSI-2/3 files to process')
     parser.add_argument('-o', '--out_dir', type=str, nargs='?',
                         required=False, default='.',
                         help="Output directory where to store the level1c file")

level1c4pps/__init__.py

@@ -184,7 +184,7 @@ def convert_angles(scene, delete_azimuth=False):
         'sunazimuth': np.array([-18000, 18000], dtype='int16'),
         'satazimuth': np.array([-18000, 18000], dtype='int16'),
     },
-    'mersi2_file_key':  {
+    'mersi_file_key':  {
         'sunzenith': 'Geolocation/SolarZenithAngle',
         'satzenith': 'Geolocation/SensorZenithAngle',
         'azimuthdiff': 'Geolocation/SensorSolarAzimuthDifference',

level1c4pps/mersi22pps_lib.py level1c4pps/mersi2pps_lib.py

@@ -22,18 +22,22 @@
 #   Nina Hakansson <nina.hakansson@smhi.se>
 #   Adam.Dybbroe <adam.dybbroe@smhi.se>
 
-"""Functions to convert MERSI-2 level-1 data to a NWCSAF/PPS level-1c formatet netCDF/CF file."""
+"""Functions to convert MERSI-2/3 level-1 data to a NWCSAF/PPS level-1c formated netCDF/CF file."""
 
+import logging
 import os
 import time
+
+import numpy as np
 from satpy.scene import Scene
-from level1c4pps import (get_encoding, compose_filename,
+from level1c4pps import (ANGLE_ATTRIBUTES,
+                         compose_filename,
+                         convert_angles,
+                         get_encoding,
+                         get_header_attrs,
+                         rename_latitude_longitude,
                          set_header_and_band_attrs_defaults,
-                         ANGLE_ATTRIBUTES, rename_latitude_longitude,
-                         update_angle_attributes, get_header_attrs,
-                         convert_angles)
-import pyspectral  # testing that pyspectral is available # noqa: F401
-import logging
+                         update_angle_attributes)
 
 # Example:
 # tf2019234102243.FY3D-X_MERSI_GEOQK_L1B.HDF
@@ -42,10 +46,15 @@
 # tf2019234102243.FY3D-X_MERSI_0250M_L1B.HDF
 #
 
-logger = logging.getLogger('mersi22pps')
+logger = logging.getLogger('mersi2pps')
+
+SENSOR = {'FY3D': 'merci-2',
+          'FY3F': 'merci-3',
+          'FY3H': 'merci-3'}
+
+SATPY_READER = {'merci-2': 'mersi2_l1b',
+                'merci-3': 'mersi3_l1b'}
 
-PLATFORM_SHORTNAMES = {'FY3D': 'FY-3D'}
-# BANDNAMES = ['%d' % (chn+1) for chn in range(25)]
 BANDNAMES = ['3', '4', '5', '6', '20', '22', '23', '24', '25']
 
 REFL_BANDS = ['3', '4', '5', '6']
@@ -63,28 +72,22 @@
                 '24': 'ch_tb11',
                 '25': 'ch_tb12'}
 
-MERSI2_LEVEL1_FILE_PATTERN = 'tf{start_time:%Y%j%H%M%S}.{platform_shortname:4s}-X_MERSI_{dataset}_L1B.HDF'
-
-
-def get_encoding_mersi2(scene):
-    """Get netcdf encoding for all datasets."""
-    return get_encoding(scene,
-                        BANDNAMES,
-                        PPS_TAGNAMES,
-                        chunks=None)
-
 
 def set_header_and_band_attrs(scene, orbit_n=0):
     """Set and delete some attributes."""
     irch = scene['24']
-    nimg = set_header_and_band_attrs_defaults(scene, BANDNAMES, PPS_TAGNAMES, REFL_BANDS, irch, orbit_n=orbit_n)
-    scene.attrs['source'] = "mersi22pps.py"
+    nimg = set_header_and_band_attrs_defaults(scene,
+                                              BANDNAMES,
+                                              PPS_TAGNAMES,
+                                              REFL_BANDS,
+                                              irch,
+                                              orbit_n=orbit_n)
+    scene.attrs['source'] = "mersi2pps.py"
     return nimg
 
 
 def remove_broken_data(scene):
     """Set bad data to nodata."""
-    import numpy as np
     for band in BANDNAMES:
         if band in REFL_BANDS:
             continue
@@ -93,42 +96,40 @@ def remove_broken_data(scene):
             scene[band].values = scene[band].values + remove
 
 
+def get_sensor(scene_file):
+    """Get sensor associated to the scene file."""
+    for platform, sensor in SENSOR.items():
+        if platform in scene_file:
+            return sensor
+    logger.info("Failed to determine sensor associated to scene file: '%s'", scene_file)
+    return None
+
+
 def process_one_scene(scene_files, out_path, engine='h5netcdf', orbit_n=0):
     """Make level 1c files in PPS-format."""
     tic = time.time()
-    scn_ = Scene(
-        reader='mersi2_l1b',
-        filenames=scene_files)
-
-    scn_.load(BANDNAMES + ['latitude', 'longitude'] + ANGLE_NAMES, resolution=1000)
-
-    # Remove bad data at first and last column
-    remove_broken_data(scn_)
-
-    # one ir channel
-    irch = scn_['24']
-
-    # Set header and band attributes
-    set_header_and_band_attrs(scn_, orbit_n=orbit_n)
-
-    # Rename longitude, latitude to lon, lat.
-    rename_latitude_longitude(scn_)
-
-    # Convert angles to PPS
-    convert_angles(scn_, delete_azimuth=True)
-    update_angle_attributes(scn_, irch)
+    sensor = get_sensor(os.path.basename(scene_files[0]))
+    reader = SATPY_READER[sensor]
+    scene = Scene(reader=reader, filenames=scene_files)
+    scene.load(BANDNAMES + ['latitude', 'longitude'] + ANGLE_NAMES, resolution=1000)
+    remove_broken_data(scene)
+    irch = scene['24']  # one ir channel 
+    set_header_and_band_attrs(scene, orbit_n=orbit_n)
+    rename_latitude_longitude(scene)
+    convert_angles(scene, delete_azimuth=True)
+    update_angle_attributes(scene, irch)
     for angle in ['sunzenith', 'satzenith', 'azimuthdiff']:
-        scn_[angle].attrs['file_key'] = ANGLE_ATTRIBUTES['mersi2_file_key'][angle]
+        scene[angle].attrs['file_key'] = ANGLE_ATTRIBUTES['mersi_file_key'][angle]
 
-    filename = compose_filename(scn_, out_path, instrument='mersi2', band=irch)
-    scn_.save_datasets(writer='cf',
+    filename = compose_filename(scene, out_path, instrument=sensor.replace('-', ''), band=irch)
+    encoding = get_encoding(scene, BANDNAMES, PPS_TAGNAMES, chunks=None)
+    attrs = get_header_attrs(scene, band=irch, sensor=sensor)
+    scene.save_datasets(writer='cf',
                        filename=filename,
-                       header_attrs=get_header_attrs(scn_, band=irch, sensor='mersi-2'),
+                       header_attrs=attrs,
                        engine=engine,
                        include_lonlats=False,
                        flatten_attrs=True,
-                       encoding=get_encoding_mersi2(scn_))
-    print("Saved file {:s} after {:3.1f} seconds".format(
-        os.path.basename(filename),
-        time.time()-tic))
+                       encoding=encoding)
+    print(f"Saved file {os.path.basename(filename)} after {time.time() - tic:3.1f} seconds")
     return filename

level1c4pps/tests/__init__.py

@@ -25,7 +25,7 @@
 import unittest
 
 from level1c4pps.tests import (test_angles, test_seviri2pps, test_gac2pps,
-                               test_mersi22pps,  test_modis2pps, test_slstr2pps,
+                               test_mersi2pps,  test_modis2pps, test_slstr2pps,
                                test_viirs2pps, test_eumgacfdr2pps,
                                test_avhrr2pps, test_init)
 
@@ -37,7 +37,7 @@ def suite():
     mysuite.addTests(test_seviri2pps.suite())
     mysuite.addTests(test_gac2pps.suite())
     mysuite.addTests(test_eumgacfdr2pps.suite())
-    mysuite.addTests(test_mersi22pps.suite())
+    mysuite.addTests(test_mersi2pps.suite())
     mysuite.addTests(test_modis2pps.suite())
     mysuite.addTests(test_avhrr2pps.suite())
     mysuite.addTests(test_slstr2pps.suite())

level1c4pps/tests/test_mersi22pps.py level1c4pps/tests/test_mersi2pps.py

@@ -21,7 +21,7 @@
 #   Stephan Finkensieper <stephan.finkensieper@dwd.de>
 #   Nina Hakansson <nina.hakansson@smhi.se>
 
-"""Unit tests for the gac2pps_lib module."""
+"""Unit tests for the merci2pps_lib module."""
 
 import datetime as dt
 import unittest
@@ -31,15 +31,15 @@
     import mock
 from satpy import Scene
 
-import level1c4pps.mersi22pps_lib as mersi22pps
+import level1c4pps.mersi2pps_lib as mersi2pps
 
 
-class TestMersi22PPS(unittest.TestCase):
-    """Test mersi22pps_lib."""
+class TestMersi2PPS(unittest.TestCase):
+    """Test mersi2pps_lib."""
 
     def setUp(self):
         """Create a test scene."""
-        mersi22pps.BANDNAMES = ['3', '24']
+        mersi2pps.BANDNAMES = ['3', '24']
         vis006 = mock.MagicMock(attrs={'name': 'image0',
                                        'wavelength': [1, 2, 3, 'um'],
                                        'id_tag': 'ch_r06'})
@@ -60,32 +60,6 @@ def setUp(self):
             self.scene[key].attrs['name'] = pps_name
         self.scene.attrs['sensor'] = ['mersi2']
 
-    def test_get_encoding(self):
-        """Test encoding for MERSI-2."""
-        enc_exp_angles = {'dtype': 'int16',
-                          'scale_factor': 0.01,
-                          'zlib': True,
-                          'complevel': 4,
-                          '_FillValue': -32767,
-                          'add_offset': 0.0}
-        encoding_exp = {
-            'image0': {'dtype': 'int16',
-                       'scale_factor': 0.01,
-                       'zlib': True,
-                       'complevel': 4,
-                       '_FillValue': -32767,
-                       'add_offset': 0.0},
-            'image1': {'dtype': 'int16',
-                       'scale_factor': 0.01,
-                       '_FillValue': -32767,
-                       'zlib': True,
-                       'complevel': 4,
-                       'add_offset': 273.15},
-            'satzenith': enc_exp_angles
-        }
-        encoding = mersi22pps.get_encoding_mersi2(self.scene)
-        self.assertDictEqual(encoding, encoding_exp)
-
     def test_compose_filename(self):
         """Test compose filename for MERSI-2."""
         start_time = dt.datetime(2009, 7, 1, 12, 15)
@@ -99,15 +73,25 @@ def test_compose_filename(self):
         band = mock.MagicMock(attrs={'start_time': start_time,
                                      'end_time': end_time})
         fname_exp = '/out/path/S_NWC_mersi2_noaa19_99999_20090701T1216000Z_20090701T1227000Z.nc'
-        fname = mersi22pps.compose_filename(scene, '/out/path', 'mersi2', band=band)
+        fname = mersi2pps.compose_filename(scene, '/out/path', 'mersi2', band=band)
         self.assertEqual(fname, fname_exp)
 
     def test_set_header_and_band_attrs(self):
         """Test to set header_and_band_attrs."""
-        mersi22pps.set_header_and_band_attrs(self.scene, orbit_n='12345')
+        mersi2pps.set_header_and_band_attrs(self.scene, orbit_n='12345')
         self.assertTrue(isinstance(self.scene.attrs['orbit_number'], int))
         self.assertEqual(self.scene.attrs['orbit_number'], 12345)
 
+    def test_get_sensor(self):
+        """Test get sensor."""
+        sensor = mersi2pps.get_sensor('tf2019234102243.FY3D-X_MERSI_GEOQK_L1B.HDF')
+        self.assertEqual(sensor, "merci-2")
+
+    def test_get_sensor_returns_none(self):
+        """Test get sensor returns none for not recognized file."""
+        sensor = mersi2pps.get_sensor('not_recognized_file')
+        self.assertEqual(sensor, None)
+
 
 def suite():
     """Create the test suite for test_mersi22pps."""