Add pad_data functionality to FCI L1c FDHSI reader

satpy/etc/readers/fci_l1c_fdhsi.yaml

@@ -6,7 +6,7 @@ reader:
     Reader for FCI FDSHI data in NetCDF4 format.
     Used to read Meteosat Third Generation (MTG) Flexible
     Combined Imager (FCI) Full Disk High Spectral Imagery (FDHSI) data.
-  reader: !!python/name:satpy.readers.yaml_reader.GEOFlippableFileYAMLReader
+  reader: !!python/name:satpy.readers.yaml_reader.GEOSegmentYAMLReader
   sensors: [fci]
 
 datasets:
@@ -379,12 +379,12 @@ datasets:
     resolution: 2000
     file_type: fci_l1c_fdhsi
 
-# Source: FCI L1 Dataset User Guide [FCIL1DUG]
-# ftp://ftp.eumetsat.int/pub/OPS/out/test-data/FCI_L1C_Format_Familiarisation/FCI_L1_Dataset_User_Guide_[FCIL1DUG].pdf
+# Source: MTG FCI L1 Product User Guide [FCIL1PUG]
+# ftp://ftp.eumetsat.int/pub/OPS/out/test-data/Test-data-for-External-Users/MTG/MTG_FCI_L1C_Enhanced_and_Non-Nominal_Test_Data/PDF_MTG_FCI_L1_PUG.pdf
 # and Example Products for Pytroll Workshop Package Description,
 # EUM/MTG/DOC/19/1079228
 file_types:
   fci_l1c_fdhsi:
     file_reader: !!python/name:satpy.readers.fci_l1c_fdhsi.FCIFDHSIFileHandler ''
-    file_patterns: ['{pflag}_{location_indicator},{data_designator},MTI{spacecraft_id:1d}+{data_source}-{processing_level}-{type}-{subtype}-{coverage}-{subsetting}-{component1}-BODY-{component3}-{purpose}-{format}_{oflag}_{originator}_{processing_time:%Y%m%d%H%M%S}_{facility}_{environment}_{start_time:%Y%m%d%H%M%S}_{end_time:%Y%m%d%H%M%S}_{processing_mode}_{special_compression}_{disposition_mode}_{repeat_cycle_in_day:>04d}_{count_in_repeat_cycle:>04d}.nc']
-    expected_segments: 70
+    file_patterns: ['{pflag}_{location_indicator},{data_designator},MTI{spacecraft_id:1d}+{data_source}-{processing_level}-{type}-{subtype}-{coverage}-{subsetting}-{component1}-BODY-{component3}-{purpose}-{format}_{oflag}_{originator}_{processing_time:%Y%m%d%H%M%S}_{facility_or_tool}_{environment}_{start_time:%Y%m%d%H%M%S}_{end_time:%Y%m%d%H%M%S}_{processing_mode}_{special_compression}_{disposition_mode}_{repeat_cycle_in_day:>04d}_{count_in_repeat_cycle:>04d}.nc']
+    expected_segments: 40

satpy/readers/yaml_reader.py

@@ -259,7 +259,7 @@ def load_ds_ids_from_config(self):
                 if val_type is not None and issubclass(val_type, tuple):
                     # special case: wavelength can be [min, nominal, max]
                     # but is still considered 1 option
-                    id_kwargs.append((val, ))
+                    id_kwargs.append((val,))
                 elif isinstance(val, (list, tuple, set)):
                     # this key has multiple choices
                     # (ex. 250 meter, 500 meter, 1000 meter resolutions)
@@ -1081,6 +1081,9 @@ def create_filehandlers(self, filenames, fh_kwargs=None):
                 ts = fh.filename_info.get('total_segments', 1)
                 # if the YAML has segments explicitly specified then use that
                 fh.filetype_info.setdefault('expected_segments', ts)
+                # add segment key-values for FCI filehandlers
+                if 'segment' not in fh.filename_info:
+                    fh.filename_info['segment'] = fh.filename_info.get('count_in_repeat_cycle', 1)
         return created_fhs
 
     @staticmethod
@@ -1097,13 +1100,28 @@ def _load_dataset(dsid, ds_info, file_handlers, dim='y', pad_data=True):
             raise KeyError(
                 "Could not load {} from any provided files".format(dsid))
 
+        padding_fci_scene = file_handlers[0].filetype_info.get('file_type') == 'fci_l1c_fdhsi'
+
         empty_segment = xr.full_like(projectable, np.nan)
         for i, sli in enumerate(slice_list):
             if sli is None:
-                slice_list[i] = empty_segment
+                if padding_fci_scene:
+                    slice_list[i] = _get_empty_segment_with_height(empty_segment,
+                                                                   _get_FCI_L1c_FDHSI_chunk_height(
+                                                                       empty_segment.shape[1], i + 1),
+                                                                   dim=dim)
+                else:
+                    slice_list[i] = empty_segment
 
         while expected_segments > counter:
-            slice_list.append(empty_segment)
+            if padding_fci_scene:
+                slice_list.append(_get_empty_segment_with_height(empty_segment,
+                                                                 _get_FCI_L1c_FDHSI_chunk_height(
+                                                                     empty_segment.shape[1], counter + 1),
+                                                                 dim=dim))
+            else:
+                slice_list.append(empty_segment)
+
             counter += 1
 
         if dim not in slice_list[0].dims:
@@ -1156,20 +1174,24 @@ def _pad_later_segments_area(file_handlers, dsid):
     available_segments = [int(fh.filename_info.get('segment', 1)) for
                           fh in file_handlers]
     area_defs = {}
+    padding_fci_scene = file_handlers[0].filetype_info.get('file_type') == 'fci_l1c_fdhsi'
+
     for segment in range(available_segments[0], expected_segments + 1):
         try:
             idx = available_segments.index(segment)
             fh = file_handlers[idx]
             area = fh.get_area_def(dsid)
         except ValueError:
             logger.debug("Padding to full disk with segment nr. %d", segment)
-            ext_diff = area.area_extent[1] - area.area_extent[3]
-            new_ll_y = area.area_extent[1] + ext_diff
+
+            new_height_proj_coord, new_height_px = _get_new_areadef_heights(area, seg_size, segment,
+                                                                            padding_fci_scene)
+            new_ll_y = area.area_extent[1] + new_height_proj_coord
             new_ur_y = area.area_extent[1]
             fill_extent = (area.area_extent[0], new_ll_y,
                            area.area_extent[2], new_ur_y)
             area = AreaDefinition('fill', 'fill', 'fill', area.proj_dict,
-                                  seg_size[1], seg_size[0],
+                                  seg_size[1], new_height_px,
                                   fill_extent)
 
         area_defs[segment] = area
@@ -1185,17 +1207,20 @@ def _pad_earlier_segments_area(file_handlers, dsid, area_defs):
     area = file_handlers[0].get_area_def(dsid)
     seg_size = area.shape
     proj_dict = area.proj_dict
+    padding_fci_scene = file_handlers[0].filetype_info.get('file_type') == 'fci_l1c_fdhsi'
+
     for segment in range(available_segments[0] - 1, 0, -1):
         logger.debug("Padding segment %d to full disk.",
                      segment)
-        ext_diff = area.area_extent[1] - area.area_extent[3]
+
+        new_height_proj_coord, new_height_px = _get_new_areadef_heights(area, seg_size, segment, padding_fci_scene)
         new_ll_y = area.area_extent[3]
-        new_ur_y = area.area_extent[3] - ext_diff
+        new_ur_y = area.area_extent[3] - new_height_proj_coord
         fill_extent = (area.area_extent[0], new_ll_y,
                        area.area_extent[2], new_ur_y)
         area = AreaDefinition('fill', 'fill', 'fill',
                               proj_dict,
-                              seg_size[1], seg_size[0],
+                              seg_size[1], new_height_px,
                               fill_extent)
         area_defs[segment] = area
         seg_size = area.shape
@@ -1235,3 +1260,51 @@ def _find_missing_segments(file_handlers, ds_info, dsid):
         slice_list.append(None)
 
     return counter, expected_segments, slice_list, failure, projectable
+
+
+def _get_new_areadef_heights(previous_area, previous_seg_size, segment_n, padding_fci_scene):
+    """Get the area definition heights in projection coordinates and pixels for the new padded segment."""
+    if padding_fci_scene:
+        # retrieve the chunk/segment pixel height
+        new_height_px = _get_FCI_L1c_FDHSI_chunk_height(previous_seg_size[1], segment_n)
+        # scale the previous vertical area extent using the new pixel height
+        new_height_proj_coord = (previous_area.area_extent[1] - previous_area.area_extent[3]) * new_height_px / \
+            previous_seg_size[0]
+    else:
+        # all other cases have constant segment size, so reuse the previous segment heights
+        new_height_px = previous_seg_size[0]
+        new_height_proj_coord = previous_area.area_extent[1] - previous_area.area_extent[3]
+    return new_height_proj_coord, new_height_px
+
+
+def _get_empty_segment_with_height(empty_segment, new_height, dim):
+    """Get a new empty segment with the specified height."""
+    if empty_segment.shape[0] > new_height:
+        # if current empty segment is too tall, slice the DataArray
+        return empty_segment[:new_height, :]
+    elif empty_segment.shape[0] < new_height:
+        # if current empty segment is too short, concatenate a slice of the DataArray
+        return xr.concat([empty_segment, empty_segment[:new_height - empty_segment.shape[0], :]], dim=dim)
+    else:
+        return empty_segment
+
+
+def _get_FCI_L1c_FDHSI_chunk_height(chunk_width, chunk_n):
+    """Get the height in pixels of a FCI L1c FDHSI chunk given the chunk width and number (starting from 1)."""
+    if chunk_width == 11136:
+        # 1km resolution case
+        if chunk_n in [3, 5, 8, 10, 13, 15, 18, 20, 23, 25, 28, 30, 33, 35, 38, 40]:
+            chunk_height = 279
+        else:
+            chunk_height = 278
+    elif chunk_width == 5568:
+        # 2km resolution case
+        if chunk_n in [5, 10, 15, 20, 25, 30, 35, 40]:
+            chunk_height = 140
+        else:
+            chunk_height = 139
+    else:
+        raise ValueError("FCI L1c FDHSI chunk width {} not recognized. Must be either 5568 or 11136.".format(
+            chunk_width))
+
+    return chunk_height

satpy/tests/reader_tests/test_fci_l1c_fdhsi.py

@@ -313,7 +313,7 @@ def test_load_counts(self, reader_configs):
         reader.create_filehandlers(loadables)
         res = reader.load(
                 [make_dataid(name=name, calibration="counts") for name in
-                    self._chans["solar"] + self._chans["terran"]])
+                    self._chans["solar"] + self._chans["terran"]], pad_data=False)
         assert 16 == len(res)
         for ch in self._chans["solar"] + self._chans["terran"]:
             assert res[ch].shape == (200*2, 11136)
@@ -342,7 +342,7 @@ def test_load_radiance(self, reader_configs):
         reader.create_filehandlers(loadables)
         res = reader.load(
                 [make_dataid(name=name, calibration="radiance") for name in
-                    self._chans["solar"] + self._chans["terran"]])
+                    self._chans["solar"] + self._chans["terran"]], pad_data=False)
         assert 16 == len(res)
         for ch in self._chans["solar"] + self._chans["terran"]:
             assert res[ch].shape == (200, 11136)
@@ -371,7 +371,7 @@ def test_load_reflectance(self, reader_configs):
         reader.create_filehandlers(loadables)
         res = reader.load(
                 [make_dataid(name=name, calibration="reflectance") for name in
-                    self._chans["solar"]])
+                    self._chans["solar"]], pad_data=False)
         assert 8 == len(res)
         for ch in self._chans["solar"]:
             assert res[ch].shape == (200, 11136)
@@ -396,7 +396,7 @@ def test_load_bt(self, reader_configs, caplog):
         with caplog.at_level(logging.WARNING):
             res = reader.load(
                     [make_dataid(name=name, calibration="brightness_temperature") for
-                        name in self._chans["terran"]])
+                        name in self._chans["terran"]], pad_data=False)
             assert caplog.text == ""
         for ch in self._chans["terran"]:
             assert res[ch].shape == (200, 11136)
@@ -435,7 +435,7 @@ def test_load_quality_only(self, reader_configs):
         reader = load_reader(reader_configs)
         loadables = reader.select_files_from_pathnames(filenames)
         reader.create_filehandlers(loadables)
-        res = reader.load(["ir_123_pixel_quality"])
+        res = reader.load(["ir_123_pixel_quality"], pad_data=False)
         assert res["ir_123_pixel_quality"].attrs["name"] == "ir_123_pixel_quality"
 
     def test_platform_name(self, reader_configs):
@@ -455,7 +455,7 @@ def test_platform_name(self, reader_configs):
         reader = load_reader(reader_configs)
         loadables = reader.select_files_from_pathnames(filenames)
         reader.create_filehandlers(loadables)
-        res = reader.load(["ir_123"])
+        res = reader.load(["ir_123"], pad_data=False)
         assert res["ir_123"].attrs["platform_name"] == "MTG-I1"
 
     def test_excs(self, reader_configs):
@@ -506,7 +506,7 @@ def test_handling_bad_data_ir(self, reader_configs, caplog):
         with caplog.at_level("ERROR"):
             reader.load([make_dataid(
                     name="ir_123",
-                    calibration="brightness_temperature")])
+                    calibration="brightness_temperature")], pad_data=False)
             assert "cannot produce brightness temperature" in caplog.text
 
     def test_handling_bad_data_vis(self, reader_configs, caplog):
@@ -526,5 +526,5 @@ def test_handling_bad_data_vis(self, reader_configs, caplog):
         with caplog.at_level("ERROR"):
             reader.load([make_dataid(
                     name="vis_04",
-                    calibration="reflectance")])
+                    calibration="reflectance")], pad_data=False)
             assert "cannot produce reflectance" in caplog.text