diff --git a/AUTHORS.md b/AUTHORS.md index c5f8607c85..8964f9e0d4 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -45,6 +45,7 @@ The following people have made contributions to this project: - [Sauli Joro (sjoro)](https://github.com/sjoro) - [Pouria Khalaj](https://github.com/pkhalaj) - [Janne Kotro (jkotro)](https://github.com/jkotro) +- [Beke Kremmling (bkremmli)](https://github.com/bkremmli) - Deutscher Wetterdienst - [Ralph Kuehn (ralphk11)](https://github.com/ralphk11) - [Panu Lahtinen (pnuu)](https://github.com/pnuu) - [Jussi Leinonen (jleinonen)](https://github.com/jleinonen) - meteoswiss diff --git a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml index ec3c5cab77..da30cb2545 100644 --- a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml +++ b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml @@ -20,14 +20,18 @@ file_types: nc_easy: file_reader: !!python/name:satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriEasyFcdrFileHandler file_patterns: [ - 'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude:f}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_EASY_{processor_version}_{format_version}.nc' + 'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude:f}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_EASY_{processor_version}_{format_version}.nc', # Example: FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc + '{sensor}_FCDR-EASY_{level}_{platform}-E{projection_longitude:s}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_{release}.nc' + # Example: MVIRI_FCDR-EASY_L15_MET7-E0000_200607060600_200607060630_0200.nc ] nc_full: file_reader: !!python/name:satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriFullFcdrFileHandler file_patterns: [ - 'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude:f}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_FULL_{processor_version}_{format_version}.nc' + 'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude:f}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_FULL_{processor_version}_{format_version}.nc', # Example: FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_FULL_v2.6_fv3.1.nc + '{sensor}_FCDR-FULL_{level}_{platform}-E{projection_longitude:s}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_{release}.nc' + # Example: MVIRI_FCDR-FULL_L15_MET7-E0000_200607060600_200607060630_0200.nc ] datasets: diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py index fc5aea2c8e..b11a7d07b1 100644 --- a/satpy/readers/mviri_l1b_fiduceo_nc.py +++ b/satpy/readers/mviri_l1b_fiduceo_nc.py @@ -162,9 +162,9 @@ from satpy.readers._geos_area import get_area_definition, get_area_extent, sampling_to_lfac_cfac from satpy.readers.file_handlers import BaseFileHandler -from satpy.utils import get_legacy_chunk_size +from satpy.utils import get_chunk_size_limit -CHUNK_SIZE = get_legacy_chunk_size() +CHUNK_SIZE = get_chunk_size_limit() EQUATOR_RADIUS = 6378140.0 POLE_RADIUS = 6356755.0 ALTITUDE = 42164000.0 - EQUATOR_RADIUS @@ -186,6 +186,8 @@ ] HIGH_RESOL = 2250 +warnings.filterwarnings("ignore", message="^.*We do not yet support duplicate dimension names, but " + "we do allow initial construction of the object.*$") class IRWVCalibrator: """Calibrate IR & WV channels.""" @@ -459,6 +461,10 @@ def __init__(self, nc): """Wrap the given dataset.""" self.nc = nc + self._decode_cf() + self._fix_duplicate_dimensions(self.nc) + + @property def attrs(self): """Exposes dataset attributes.""" @@ -509,6 +515,31 @@ def _cleanup_attrs(self, ds): # satpy warnings. ds.attrs.pop("ancillary_variables", None) + def _decode_cf(self): + """Decode data.""" + # time decoding with decode_cf results in error - decode separately! + time_dims, time = self._decode_time() + self.nc = self.nc.drop_vars(time.name) + self.nc = xr.decode_cf(self.nc) + self.nc[time.name] = (time_dims, time.values) + + def _decode_time(self): + """Decode time using fill value and offset.""" + time = self.get_time() + time_dims = self.nc[time.name].dims + time = xr.where(time == time.attrs["_FillValue"], np.datetime64("NaT"), + (time + time.attrs["add_offset"]).astype("datetime64[s]").astype("datetime64[ns]")) + + return (time_dims, time) + + def _fix_duplicate_dimensions(self, nc): + """Rename dimensions as duplicate dimensions names are not supported by xarray.""" + nc.variables["covariance_spectral_response_function_vis"].dims = ("srf_size_1", "srf_size_2") + self.nc = nc.drop_dims("srf_size") + nc.variables["channel_correlation_matrix_independent"].dims = ("channel_1", "channel_2") + nc.variables["channel_correlation_matrix_structured"].dims = ("channel_1", "channel_2") + self.nc = nc.drop_dims("channel") + def get_time(self): """Get time coordinate. @@ -558,13 +589,17 @@ def __init__(self, filename, filename_info, filetype_info, # noqa: D417 chunks={"x": CHUNK_SIZE, "y": CHUNK_SIZE, "x_ir_wv": CHUNK_SIZE, - "y_ir_wv": CHUNK_SIZE} + "y_ir_wv": CHUNK_SIZE}, + # see dataset wrapper for why decoding is disabled + decode_cf=False, + decode_times=False, + mask_and_scale=False, ) + self.nc = DatasetWrapper(nc_raw) - # Projection longitude is not provided in the file, read it from the - # filename. - self.projection_longitude = float(filename_info["projection_longitude"]) + self.projection_longitude = self._get_projection_longitude(filename_info) + self.calib_coefs = self._get_calib_coefs() self._get_angles = functools.lru_cache(maxsize=8)( @@ -574,6 +609,13 @@ def __init__(self, filename, filename_info, filetype_info, # noqa: D417 self._get_acq_time_uncached ) + def _get_projection_longitude(self, filename_info): + """Read projection longitude from filename as it is not provided in the file.""" + if "." in str(filename_info["projection_longitude"]): + return float(filename_info["projection_longitude"]) + + return float(filename_info["projection_longitude"]) / 100 + def get_dataset(self, dataset_id, dataset_info): """Get the dataset.""" name = dataset_id["name"] diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py index 04694c145a..ffc3e980e2 100644 --- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py +++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py @@ -36,6 +36,7 @@ DatasetWrapper, FiduceoMviriEasyFcdrFileHandler, FiduceoMviriFullFcdrFileHandler, + Interpolator, ) from satpy.tests.utils import make_dataid @@ -43,6 +44,8 @@ # The following fixtures are not defined in this file, but are used and injected by Pytest: # - request +fill_val = np.uint32(429496729) # FillValue lower than in dataset to be windows-compatible + attrs_exp: dict = { "platform": "MET7", "raw_metadata": {"foo": "bar"}, @@ -61,8 +64,8 @@ {"sun_earth_distance_correction_applied": True, "sun_earth_distance_correction_factor": 1.} ) -acq_time_vis_exp = [np.datetime64("1970-01-01 00:30").astype("datetime64[ns]"), - np.datetime64("1970-01-01 00:30").astype("datetime64[ns]"), +acq_time_vis_exp = [np.datetime64("NaT").astype("datetime64[ns]"), + np.datetime64("NaT").astype("datetime64[ns]"), np.datetime64("1970-01-01 02:30").astype("datetime64[ns]"), np.datetime64("1970-01-01 02:30").astype("datetime64[ns]")] vis_counts_exp = xr.DataArray( @@ -79,6 +82,7 @@ }, attrs=attrs_exp ) + vis_rad_exp = xr.DataArray( np.array( [[np.nan, 18.56, 38.28, 58.], @@ -124,7 +128,10 @@ }, attrs=attrs_exp ) -acq_time_ir_wv_exp = [np.datetime64("1970-01-01 00:30").astype("datetime64[ns]"), + +u_struct_refl_exp = u_vis_refl_exp.copy() + +acq_time_ir_wv_exp = [np.datetime64("NaT"), np.datetime64("1970-01-01 02:30").astype("datetime64[ns]")] wv_counts_exp = xr.DataArray( np.array( @@ -249,31 +256,37 @@ height=2 ) +@pytest.fixture(name="time_fake_dataset") +def fixture_time_fake_dataset(): + """Create time for fake dataset.""" + time = np.arange(4) * 60 * 60 + time[0] = fill_val + time[1] = fill_val + time = time.reshape(2, 2) + + return time @pytest.fixture(name="fake_dataset") -def fixture_fake_dataset(): +def fixture_fake_dataset(time_fake_dataset): """Create fake dataset.""" count_ir = da.linspace(0, 255, 4, dtype=np.uint8).reshape(2, 2) count_wv = da.linspace(0, 255, 4, dtype=np.uint8).reshape(2, 2) count_vis = da.linspace(0, 255, 16, dtype=np.uint8).reshape(4, 4) sza = da.from_array( np.array( - [[45, 90], - [0, 45]], + [[45, 90], [0, 45]], dtype=np.float32 ) ) mask = da.from_array( np.array( - [[0, 0, 0, 0], - [0, 0, 0, 0], - [0, 0, 1, 0], # 1 = "invalid" - [0, 0, 0, 0]], + [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0]], # 1 = "invalid" dtype=np.uint8 ) ) - time = np.arange(4) * 60 * 60 * 1e9 - time = time.astype("datetime64[ns]").reshape(2, 2) + + cov = da.from_array([[1, 2], [3, 4]]) + ds = xr.Dataset( data_vars={ "count_vis": (("y", "x"), count_vis), @@ -281,9 +294,10 @@ def fixture_fake_dataset(): "count_ir": (("y_ir_wv", "x_ir_wv"), count_ir), "toa_bidirectional_reflectance_vis": vis_refl_exp / 100, "u_independent_toa_bidirectional_reflectance": u_vis_refl_exp / 100, + "u_structured_toa_bidirectional_reflectance": u_vis_refl_exp / 100, "quality_pixel_bitmask": (("y", "x"), mask), "solar_zenith_angle": (("y_tie", "x_tie"), sza), - "time_ir_wv": (("y_ir_wv", "x_ir_wv"), time), + "time_ir_wv": (("y_ir_wv", "x_ir_wv"), time_fake_dataset), "a_ir": -5.0, "b_ir": 1.0, "bt_a_ir": 10.0, @@ -303,6 +317,9 @@ def fixture_fake_dataset(): "sub_satellite_longitude_end": np.nan, "sub_satellite_latitude_start": np.nan, "sub_satellite_latitude_end": 0.1, + "covariance_spectral_response_function_vis": (("srf_size", "srf_size"), cov), + "channel_correlation_matrix_independent": (("channel", "channel"), cov), + "channel_correlation_matrix_structured": (("channel", "channel"), cov) }, coords={ "y": [1, 2, 3, 4], @@ -310,22 +327,31 @@ def fixture_fake_dataset(): "y_ir_wv": [1, 2], "x_ir_wv": [1, 2], "y_tie": [1, 2], - "x_tie": [1, 2] - + "x_tie": [1, 2], }, attrs={"foo": "bar"} ) ds["count_ir"].attrs["ancillary_variables"] = "a_ir b_ir" ds["count_wv"].attrs["ancillary_variables"] = "a_wv b_wv" + ds["quality_pixel_bitmask"].encoding["chunksizes"] = (2, 2) + ds["time_ir_wv"].attrs["_FillValue"] = fill_val + ds["time_ir_wv"].attrs["add_offset"] = 0 + return ds +@pytest.fixture(name="projection_longitude", params=["57.0"]) +def fixture_projection_longitude(request): + """Get projection longitude as string.""" + return request.param + + @pytest.fixture( name="file_handler", params=[FiduceoMviriEasyFcdrFileHandler, FiduceoMviriFullFcdrFileHandler] ) -def fixture_file_handler(fake_dataset, request): +def fixture_file_handler(fake_dataset, request, projection_longitude): """Create mocked file handler.""" marker = request.node.get_closest_marker("file_handler_data") mask_bad_quality = True @@ -338,7 +364,7 @@ def fixture_file_handler(fake_dataset, request): filename="filename", filename_info={"platform": "MET7", "sensor": "MVIRI", - "projection_longitude": "57.0"}, + "projection_longitude": projection_longitude}, filetype_info={"foo": "bar"}, mask_bad_quality=mask_bad_quality ) @@ -359,7 +385,9 @@ def fixture_reader(): class TestFiduceoMviriFileHandlers: """Unit tests for FIDUCEO MVIRI file handlers.""" - def test_init(self, file_handler): + + @pytest.mark.parametrize("projection_longitude", ["57.0", "5700"], indirect=True) + def test_init(self, file_handler, projection_longitude): """Test file handler initialization.""" assert file_handler.projection_longitude == 57.0 assert file_handler.mask_bad_quality is True @@ -379,7 +407,8 @@ def test_init(self, file_handler): ("quality_pixel_bitmask", None, 2250, quality_pixel_bitmask_exp), ("solar_zenith_angle", None, 2250, sza_vis_exp), ("solar_zenith_angle", None, 4500, sza_ir_wv_exp), - ("u_independent_toa_bidirectional_reflectance", None, 4500, u_vis_refl_exp) + ("u_independent_toa_bidirectional_reflectance", None, 4500, u_vis_refl_exp), + ("u_structured_toa_bidirectional_reflectance", None, 4500, u_struct_refl_exp) ] ) def test_get_dataset(self, file_handler, name, calibration, resolution, @@ -547,17 +576,50 @@ def test_file_pattern(self, reader): "FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_FULL_v2.6_fv3.1.nc", "FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc", "FIDUCEO_FCDR_L15_MVIRI_MET7-00.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc", + "MVIRI_FCDR-EASY_L15_MET7-E0000_200607060600_200607060630_0200.nc", + "MVIRI_FCDR-EASY_L15_MET7-E5700_200607060600_200607060630_0200.nc", + "MVIRI_FCDR-FULL_L15_MET7-E0000_200607060600_200607060630_0200.nc", "abcde", ] files = reader.select_files_from_pathnames(filenames) - # only 3 out of 4 above should match - assert len(files) == 3 + assert len(files) == 6 class TestDatasetWrapper: """Unit tests for DatasetWrapper class.""" + def test_fix_duplicate_dimensions(self): + """Test the renaming of duplicate dimensions. + + If duplicate dimensions are within the Dataset, opening the datasets with chunks throws a warning. + The dimensions need to be renamed. + """ + foo_time = 60*60 + foo_time_exp = np.datetime64("1970-01-01 01:00").astype("datetime64[ns]") + + foo = xr.Dataset( + data_vars={ + "covariance_spectral_response_function_vis": (("srf_size", "srf_size"), [[1, 2], [3, 4]]), + "channel_correlation_matrix_independent": (("channel", "channel"), [[1, 2], [3, 4]]), + "channel_correlation_matrix_structured": (("channel", "channel"), [[1, 2], [3, 4]]), + "time_ir_wv": (("y_ir_wv", "x_ir_wv"), [[foo_time, foo_time], [foo_time, foo_time]], + {"_FillValue": fill_val, "add_offset": 0}) + } + ) + foo_ds = DatasetWrapper(foo) + + foo_exp = xr.Dataset( + data_vars={ + "covariance_spectral_response_function_vis": (("srf_size_1", "srf_size_2"), [[1, 2], [3, 4]]), + "channel_correlation_matrix_independent": (("channel_1", "channel_2"), [[1, 2], [3, 4]]), + "channel_correlation_matrix_structured": (("channel_1", "channel_2"), [[1, 2], [3, 4]]), + "time_ir_wv": (("y_ir_wv", "x_ir_wv"), [[foo_time_exp, foo_time_exp], [foo_time_exp, foo_time_exp]]) + } + ) + + xr.testing.assert_allclose(foo_ds.nc, foo_exp) + def test_reassign_coords(self): """Test reassigning of coordinates. @@ -588,6 +650,63 @@ def test_reassign_coords(self): "x": [.3, .4] } ) - ds = DatasetWrapper(nc) - foo = ds["foo"] - xr.testing.assert_equal(foo, foo_exp) + with mock.patch("satpy.readers.mviri_l1b_fiduceo_nc.DatasetWrapper._fix_duplicate_dimensions"): + with mock.patch("satpy.readers.mviri_l1b_fiduceo_nc.DatasetWrapper._decode_cf"): + ds = DatasetWrapper(nc) + foo = ds["foo"] + xr.testing.assert_equal(foo, foo_exp) + +class TestInterpolator: + """Unit tests for Interpolator class.""" + @pytest.fixture(name="time_ir_wv") + def fixture_time_ir_wv(self): + """Returns time_ir_wv.""" + return xr.DataArray( + [ + [np.datetime64("1970-01-01 01:00"), np.datetime64("1970-01-01 02:00")], + [np.datetime64("1970-01-01 03:00"), np.datetime64("1970-01-01 04:00")], + [np.datetime64("NaT"), np.datetime64("1970-01-01 06:00")], + [np.datetime64("NaT"), np.datetime64("NaT")], + ], + dims=("y", "x"), + coords={"y": [1, 3, 5, 7]} + ) + + @pytest.fixture(name="acq_time_exp") + def fixture_acq_time_exp(self): + """Returns acq_time_vis_exp.""" + vis = xr.DataArray( + [ + np.datetime64("1970-01-01 01:30"), + np.datetime64("1970-01-01 01:30"), + np.datetime64("1970-01-01 03:30"), + np.datetime64("1970-01-01 03:30"), + np.datetime64("1970-01-01 06:00"), + np.datetime64("1970-01-01 06:00"), + np.datetime64("NaT"), + np.datetime64("NaT") + ], + dims="y", + coords={"y": [1, 2, 3, 4, 5, 6, 7, 8]} + ) + + ir = xr.DataArray( + [ + np.datetime64("1970-01-01 01:30"), + np.datetime64("1970-01-01 03:30"), + np.datetime64("1970-01-01 06:00"), + np.datetime64("NaT"), + ], + dims="y", + coords={"y": [1, 3, 5, 7]} + ) + + return vis, ir + + def test_interp_acq_time(self, time_ir_wv, acq_time_exp): + """Tests time interpolation.""" + res_vis = Interpolator.interp_acq_time(time_ir_wv, target_y=acq_time_exp[0].coords["y"]) + res_ir = Interpolator.interp_acq_time(time_ir_wv, target_y=acq_time_exp[1].coords["y"]) + + xr.testing.assert_allclose(res_vis, acq_time_exp[0]) + xr.testing.assert_allclose(res_ir, acq_time_exp[1])