From e40df726e52313c0b8fa67c297de497443bdfbfa Mon Sep 17 00:00:00 2001 From: BengtRydberg Date: Fri, 16 Aug 2024 10:54:17 +0200 Subject: [PATCH] fixing and adding test --- level1c4pps/mersi2pps_lib.py | 113 ++++++++++++++-------------- level1c4pps/tests/test_mersi2pps.py | 40 +++++++--- 2 files changed, 86 insertions(+), 67 deletions(-) diff --git a/level1c4pps/mersi2pps_lib.py b/level1c4pps/mersi2pps_lib.py index bcfd258..253ea6c 100644 --- a/level1c4pps/mersi2pps_lib.py +++ b/level1c4pps/mersi2pps_lib.py @@ -21,6 +21,7 @@ # Martin Raspaud # Nina Hakansson # Adam.Dybbroe +# Bengt.Rydberg >bengt.rydberg@smhi.se> """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 diff --git a/level1c4pps/tests/test_mersi2pps.py b/level1c4pps/tests/test_mersi2pps.py index e51348c..d526143 100644 --- a/level1c4pps/tests/test_mersi2pps.py +++ b/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