Adapt reader mviri_l1b_fiduceo_nc

AUTHORS.md

@@ -44,6 +44,7 @@ The following people have made contributions to this project:
 - [Johannes Johansson (JohannesSMHI)](https://github.com/JohannesSMHI)
 - [Sauli Joro (sjoro)](https://github.com/sjoro)
 - [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

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:

satpy/readers/mviri_l1b_fiduceo_nc.py

@@ -162,9 +162,7 @@
 
 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
 
-CHUNK_SIZE = get_legacy_chunk_size()
 EQUATOR_RADIUS = 6378140.0
 POLE_RADIUS = 6356755.0
 ALTITUDE = 42164000.0 - EQUATOR_RADIUS
@@ -407,7 +405,6 @@
         """
         # Compute mean timestamp per scanline
         time = time2d.mean(dim="x")
-
         # If required, repeat timestamps in y-direction to obtain higher
         # resolution
         y = time.coords["y"].values
@@ -455,10 +452,40 @@
 class DatasetWrapper:
     """Helper class for accessing the dataset."""
 
-    def __init__(self, nc):
+    def __init__(self, nc, decode_nc=True):
         """Wrap the given dataset."""
         self.nc = nc
 
+        if decode_nc is True:
+            self._decode_cf()
+            self._fix_duplicate_dimensions(self.nc)
+            self.nc = self._chunk(self.nc)
+
+    def _decode_cf(self):
+        # remove time before decoding and add again.
+        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]"))
+        self.nc = self.nc.drop_vars(time.name)
+        self.nc = xr.decode_cf(self.nc)
+        self.nc[time.name] = (time_dims, time.values)
+
+    def _fix_duplicate_dimensions(self, nc):
+        nc.variables["covariance_spectral_response_function_vis"].dims = ("srf_size_1", "srf_size_2")
+        self.nc = nc.drop_dims("srf_size")
+
+    def _chunk(self, nc):
+
+        (chunk_size_y, chunk_size_x) = nc.variables["quality_pixel_bitmask"].encoding["chunksizes"]
+        chunks = {
+            "x": chunk_size_x,
+            "y": chunk_size_y,
+            "x_ir_wv": chunk_size_x,
+            "y_ir_wv": chunk_size_y
+        }
+        return nc.chunk(chunks)
+
     @property
     def attrs(self):
         """Exposes dataset attributes."""
@@ -555,16 +582,21 @@
         self.mask_bad_quality = mask_bad_quality
         nc_raw = xr.open_dataset(
             filename,
-            chunks={"x": CHUNK_SIZE,
-                    "y": CHUNK_SIZE,
-                    "x_ir_wv": CHUNK_SIZE,
-                    "y_ir_wv": CHUNK_SIZE}
+            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"])
+        if "." in str(filename_info["projection_longitude"]):
+            self.projection_longitude = float(filename_info["projection_longitude"])
+        else:
+            self.projection_longitude = (
+                float(filename_info["projection_longitude"][:2] + "." + filename_info["projection_longitude"][2:])
+            )
         self.calib_coefs = self._get_calib_coefs()
 
         self._get_angles = functools.lru_cache(maxsize=8)(

satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py

@@ -28,6 +28,7 @@
 import xarray as xr
 from pyproj import CRS
 from pyresample.geometry import AreaDefinition
+from pytest_lazy_fixtures import lf as lazy_fixture
 
 from satpy.readers.mviri_l1b_fiduceo_nc import (
     ALTITUDE,
@@ -36,6 +37,7 @@
     DatasetWrapper,
     FiduceoMviriEasyFcdrFileHandler,
     FiduceoMviriFullFcdrFileHandler,
+    Interpolator,
 )
 from satpy.tests.utils import make_dataid
 
@@ -61,8 +63,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 +81,7 @@
     },
     attrs=attrs_exp
 )
+
 vis_rad_exp = xr.DataArray(
     np.array(
         [[np.nan, 18.56, 38.28, 58.],
@@ -124,7 +127,7 @@
     },
     attrs=attrs_exp
 )
-acq_time_ir_wv_exp = [np.datetime64("1970-01-01 00:30").astype("datetime64[ns]"),
+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(
@@ -272,51 +275,61 @@
             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]])
+    time = np.arange(4) * 60 * 60
+    time[0] = 4294967295
+    time[1] = 4294967295
+    time = time.reshape(2, 2)
+
     ds = xr.Dataset(
         data_vars={
             "count_vis": (("y", "x"), count_vis),
             "count_wv": (("y_ir_wv", "x_ir_wv"), count_wv),
             "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,
             "quality_pixel_bitmask": (("y", "x"), mask),
             "solar_zenith_angle": (("y_tie", "x_tie"), sza),
             "time_ir_wv": (("y_ir_wv", "x_ir_wv"), time),
             "a_ir": -5.0,
             "b_ir": 1.0,
             "bt_a_ir": 10.0,
             "bt_b_ir": -1000.0,
             "a_wv": -0.5,
             "b_wv": 0.05,
             "bt_a_wv": 10.0,
             "bt_b_wv": -2000.0,
             "years_since_launch": 20.0,
             "a0_vis": 1.0,
             "a1_vis": 0.01,
             "a2_vis": -0.0001,
             "mean_count_space_vis": 1.0,
             "distance_sun_earth": 1.0,
             "solar_irradiance_vis": 650.0,
             "sub_satellite_longitude_start": 57.1,
             "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),
         },
         coords={
             "y": [1, 2, 3, 4],
             "x": [1, 2, 3, 4],
             "y_ir_wv": [1, 2],
             "x_ir_wv": [1, 2],
             "y_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"] = 4294967295
+    ds["time_ir_wv"].attrs["add_offset"] = 0
+
     return ds
 
 
@@ -547,17 +560,44 @@
             "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 an error.
+        Thus, the chunking needs to be done after opening the dataset and after the dimensions are renamed.
+        """
+        foo = xr.Dataset(
+            data_vars={"covariance_spectral_response_function_vis":
+                      (("srf_size", "srf_size"), [[1, 2], [3, 4]])}
+        )
+        foo_ds = DatasetWrapper(foo, decode_nc=False)
+        foo_ds._fix_duplicate_dimensions(foo_ds.nc)
+
+        foo_exp = xr.Dataset(
+            data_vars={"covariance_spectral_response_function_vis":
+                      (("srf_size_1", "srf_size_2"), [[1, 2], [3, 4]])}
+        )
+
+        try:
+            foo_ds.nc.chunk("auto")
+            xr.testing.assert_allclose(foo_ds.nc, foo_exp)
+        except ValueError:
+            pytest.fail("Chunking failed.")
+
     def test_reassign_coords(self):
         """Test reassigning of coordinates.
 
@@ -588,6 +628,66 @@
                 "x": [.3, .4]
             }
         )
-        ds = DatasetWrapper(nc)
+        ds = DatasetWrapper(nc, decode_nc=False)
         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_vis_exp")
+    def fixture_acq_time_vis_exp(self):
+        """Returns acq_time_vis_exp."""
+        return 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]}
+        )
+
+    @pytest.fixture(name="acq_time_ir_exp")
+    def fixture_acq_time_ir_exp(self):
+        """Returns acq_time_ir_exp."""
+        return 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]}
+        )
+
+    @pytest.mark.parametrize(
+        "acq_time_exp",
+        [
+            lazy_fixture("acq_time_ir_exp"),
+            lazy_fixture("acq_time_vis_exp")
+        ]
+    )
+    def test_interp_acq_time(self, time_ir_wv, acq_time_exp):
+        """Tests time interpolation."""
+        res = Interpolator.interp_acq_time(time_ir_wv, target_y=acq_time_exp.coords["y"])
+        xr.testing.assert_allclose(res, acq_time_exp)