fixing and adding test

level1c4pps/mersi2pps_lib.py

@@ -21,6 +21,7 @@
 #   Martin Raspaud <[email protected]>
 #   Nina Hakansson <[email protected]>
 #   Adam.Dybbroe <[email protected]>
+#   Bengt.Rydberg >[email protected]>
 
 """Functions to convert MERSI-2/3 level-1 data to a NWCSAF/PPS level-1c formated netCDF/CF file."""
 
@@ -30,14 +31,16 @@
 
 import numpy as np
 from satpy.scene import Scene
-from level1c4pps import (ANGLE_ATTRIBUTES,
-                         compose_filename,
-                         convert_angles,
-                         get_encoding,
-                         get_header_attrs,
-                         rename_latitude_longitude,
-                         set_header_and_band_attrs_defaults,
-                         update_angle_attributes)
+from level1c4pps import (
+    ANGLE_ATTRIBUTES,
+    compose_filename,
+    convert_angles,
+    get_encoding,
+    get_header_attrs,
+    rename_latitude_longitude,
+    set_header_and_band_attrs_defaults,
+    update_angle_attributes,
+)
 
 # Example filenames:
 # tf2019234102243.FY3D-X_MERSI_GEOQK_L1B.HDF
@@ -48,56 +51,53 @@
 
 logger = logging.getLogger('mersi2pps')
 
-SENSOR = {'FY3D': 'mersi-2',
-          'FY3F': 'mersi-3',
-          'FY3H': 'mersi-3'}
-
-SATPY_READER = {'mersi-2': 'mersi2_l1b',
-                'mersi-3': 'mersi3_l1b'}
-
+SENSOR = {  # sensor associated to platform
+    'FY3D': 'mersi-2',
+    'FY3F': 'mersi-3',
+    'FY3H': 'mersi-3',
+}
+SATPY_READER = {  # satpy reader associated to sensor
+    'mersi-2': 'mersi2_l1b',
+    'mersi-3': 'mersi3_l1b',
+}
+PPS_BAND_NAME = {  # PPS band name associated to satpy band name
+    '3': 'ch_r06',
+    '4': 'ch_r09',
+    '5': 'ch_r13',
+    '6': 'ch_r16',
+    '20': 'ch_tb37',
+    '22': 'ch_tb73',
+    '23': 'ch_tb85',
+    '24': 'ch_tb11',
+    '25': 'ch_tb12',
+}
+GEOLOCATION_NAMES = [  # additional variables to load
+    'latitude',
+    'longitude',
+    'satellite_azimuth_angle',
+    'satellite_zenith_angle',
+    'solar_azimuth_angle',
+    'solar_zenith_angle',
+]
 REFL_BANDS = ['3', '4', '5', '6']
-
 ONE_IR_CHANNEL = '24'
-
-GEOLOCATION_NAMES = ['latitude',
-                     'longitude',
-                     'satellite_azimuth_angle',
-                     'satellite_zenith_angle',
-                     'solar_azimuth_angle',
-                     'solar_zenith_angle']
-
-PPS_BAND_NAME = {'3': 'ch_r06',
-                 '4': 'ch_r09',
-                 '5': 'ch_r13',
-                 '6': 'ch_r16',
-                 '20': 'ch_tb37',
-                 '22': 'ch_tb73',
-                 '23': 'ch_tb85',
-                 '24': 'ch_tb11',
-                 '25': 'ch_tb12'}
-
-RESOLUTION = 1000
+RESOLUTION = 1000  # [m]
+LOW_TB = 1  # [K] very cold brightness temperature
 
 
 def set_header_and_band_attrs(scene, band, orbit_n):
     """Set and delete some attributes."""
-    set_header_and_band_attrs_defaults(scene,
-                                       list(PPS_BAND_NAME),
-                                       PPS_BAND_NAME,
-                                       REFL_BANDS,
-                                       band,
-                                       orbit_n=orbit_n)
+    set_header_and_band_attrs_defaults(
+        scene, list(PPS_BAND_NAME), PPS_BAND_NAME, REFL_BANDS, band, orbit_n=orbit_n,
+    )
     scene.attrs['source'] = "mersi2pps.py"
 
 
 def remove_broken_data(scene):
     """Set bad data to nodata."""
     for band in PPS_BAND_NAME:
-        if band in REFL_BANDS:
-            continue
-        if band in scene:
-            remove = np.where(scene[band].values < 1, np.nan, 0)  # 1K very cold
-            scene[band].values = scene[band].values + remove
+        if band in band not in REFL_BANDS and band in scene:
+            scene[band].data = np.where(scene[band].data < LOW_TB, np.nan, scene[band].data)
 
 
 def get_sensor(scene_file):
@@ -125,16 +125,15 @@ def process_one_scene(scene_files, out_path, engine='h5netcdf', orbit_n=0):
     update_angle_attributes(scene, band)
     for angle in ['sunzenith', 'satzenith', 'azimuthdiff']:
         scene[angle].attrs['file_key'] = ANGLE_ATTRIBUTES['mersi_file_key'][angle]
-
-    filename = compose_filename(scene, out_path, instrument=sensor.replace('-', ''), band=band)
-    encoding = get_encoding(scene, band_names, PPS_BAND_NAME, chunks=None)
-    attrs = get_header_attrs(scene, band=band, sensor=sensor)
-    scene.save_datasets(writer='cf',
-                        filename=filename,
-                        header_attrs=attrs,
-                        engine=engine,
-                        include_lonlats=False,
-                        flatten_attrs=True,
-                        encoding=encoding)
+    filename=compose_filename(scene, out_path, instrument=sensor.replace('-', ''), band=band)
+    scene.save_datasets(
+        writer='cf',
+        filename=filename,
+        header_attrs=get_header_attrs(scene, band=band, sensor=sensor),
+        engine=engine,
+        include_lonlats=False,
+        flatten_attrs=True,
+        encoding=get_encoding(scene, band_names, PPS_BAND_NAME, chunks=None),
+    )
     print(f"Saved file {os.path.basename(filename)} after {time.time() - tic:3.1f} seconds")
     return filename

level1c4pps/tests/test_mersi2pps.py

@@ -29,6 +29,8 @@
     from unittest import mock
 except ImportError:
     import mock
+import numpy as np
+import xarray as xr
 from satpy import Scene
 
 import level1c4pps.mersi2pps_lib as mersi2pps
@@ -82,23 +84,41 @@ def test_set_header_and_band_attrs(self):
         self.assertTrue(isinstance(self.scene.attrs['orbit_number'], int))
         self.assertEqual(self.scene.attrs['orbit_number'], 12345)
 
+    def test_remove_broken_data(self):
+        """Test remove broken data."""
+        data = xr.Dataset(
+            {
+                '3': (('y', 'x'), [[100, 0, 100, 100]]),
+                '20': (('y', 'x'), [[200, 0, 200, 200]]),
+                '22': (('y', 'x'), [[300, 0, 300, 300]]),
+            }
+        )
+        mersi2pps.remove_broken_data(data)
+        expect = xr.Dataset(
+            {
+                '3': (('y', 'x'), [[100, 0, 100, 100]]),
+                '20': (('y', 'x'), [[200, np.nan, 200, 200]]),
+                '22': (('y', 'x'), [[300, np.nan, 300, 300]]),
+            }
+        )
+        for band in expect:
+            np.testing.assert_array_equal(data[band], expect[band])
+
     def test_get_sensor(self):
         """Test get sensor."""
-        sensor = mersi2pps.get_sensor('tf2019234102243.FY3D-X_MERSI_GEOQK_L1B.HDF')
-        self.assertEqual(sensor, "mersi-2")
-        sensor = mersi2pps.get_sensor('tf2019234102243.FY3F-X_MERSI_GEOQK_L1B.HDF')
-        self.assertEqual(sensor, "mersi-3")
-
-    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)
+        for filename, expect in [
+            ('tf2019234102243.FY3D-X_MERSI_GEOQK_L1B.HDF', 'mersi-2'),
+            ('tf2019234102243.FY3F-X_MERSI_GEOQK_L1B.HDF', 'mersi-3'),
+            ('not_recognized_file', None),
+        ]:
+            sensor = mersi2pps.get_sensor(filename)
+            self.assertEqual(sensor, expect)
 
 
 def suite():
     """Create the test suite for test_mersi22pps."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
-    mysuite.addTest(loader.loadTestsFromTestCase(TestMersi22PPS))
+    mysuite.addTest(loader.loadTestsFromTestCase(TestMersi2PPS))
 
     return mysuite