Implement generic set_orientation method for geostationary datasets

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.FileYAMLReader
+  reader: !!python/name:satpy.readers.yaml_reader.GEOFlippableFileYAMLReader
   sensors: [fci]
 
 datasets:

satpy/etc/readers/seviri_l1b_native.yaml

@@ -6,7 +6,7 @@ reader:
     Reader for EUMETSAT MSG SEVIRI Level 1b native format files.
   sensors: [seviri]
   default_channels: [HRV, IR_016, IR_039, IR_087, IR_097, IR_108, IR_120, IR_134, VIS006, VIS008, WV_062, WV_073]
-  reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader
+  reader: !!python/name:satpy.readers.yaml_reader.GEOFlippableFileYAMLReader
 
 file_types:
   native_msg:

satpy/etc/readers/seviri_l1b_nc.yaml

@@ -5,7 +5,7 @@ reader:
   description: >
     NetCDF4 reader for EUMETSAT MSG SEVIRI Level 1b files.
   sensors: [seviri]
-  reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader
+  reader: !!python/name:satpy.readers.yaml_reader.GEOFlippableFileYAMLReader
   group_keys: ["processing_time", "satid"]
 
 file_types:

satpy/readers/yaml_reader.py

@@ -870,7 +870,114 @@ def _load_area_def(dsid, file_handlers):
     return final_area.squeeze()
 
 
-class GEOSegmentYAMLReader(FileYAMLReader):
+def _set_orientation(dataset, upper_right_corner='native', **kwargs):
+    """Set the orientation of datasets.
+
+    Allows to flip geostationary datasets without having to resample.
+
+    Args:
+        datasets (iterable): Limit the flipping to this list of names (str), wavelengths (float), or
+                            DatasetID objects.
+                            By default all currently loaded datasets are flipped.
+        upper_right_corner (str): Direction of the upper right corner of the image.
+                                Possible options are 'NW', 'NE', 'SW', 'SE', or 'native'.
+                                The common upright image orientation corresponds to 'NW'.
+                                Defaults to 'native' (no flipping is applied).
+
+    """
+    def _change_area_extent_of_areadef(areadef, new_extent):
+        """Change the area extent of an area definition and regenerate it."""
+        # does this need to be moved somewhere else?
+        return AreaDefinition(
+            areadef.area_id,
+            areadef.description,
+            areadef.proj_id,
+            areadef.proj_dict,
+            areadef.width,
+            areadef.height,
+            new_extent
+        )
+
+    if upper_right_corner == 'native':
+        logger.debug("Requested Dataset orientation is 'native'. No flipping is applied.")
+        return
+    if upper_right_corner not in ['NW', 'NE', 'SE', 'SW', 'native']:
+        raise ValueError("Origin corner not recognized. Should be 'NW', 'NE', 'SW', 'SE' or 'native.")
+
+    # get the target orientation
+    target_northup = False
+    target_eastright = False
+    if upper_right_corner in ['NW', 'NE']:
+        target_northup = True
+    if upper_right_corner in ['NW', 'SW']:
+        target_eastright = True
+
+    if 'area' not in dataset.attrs:
+        logger.info("Dataset is not a geographical dataset and cannot be flipped.")
+        return dataset
+
+    if dataset.attrs['area'].proj_dict['proj'] != 'geos':
+        logger.info("Dataset is not in geos projection and cannot be flipped.")
+        return dataset
+
+    # get the current dataset orientation
+    if isinstance(dataset.attrs['area'], StackedAreaDefinition):
+        # array of area extents if the Area is a StackedAreaDefinition
+        ds_area_extents = np.asarray([list(adef.area_extent) for adef in dataset.attrs['area'].defs])
+    else:
+        # array with a single item if Area is in one piece
+        ds_area_extents = np.asarray([list(dataset.attrs['area'].area_extent)])
+
+    # assumes all AreaDefinitions inside a StackedAreaDefinition have the same orientation
+    current_northup = ds_area_extents[0, 3] - ds_area_extents[0, 1] > 0
+    current_eastright = ds_area_extents[0, 2] - ds_area_extents[0, 0] > 0
+
+    # if current and target orientations mismatch, then flip data and switch area extent elements
+    if target_northup == current_northup and target_eastright == current_eastright:
+        logger.info("Dataset is already in the target orientation.")
+        return dataset
+
+    if target_northup != current_northup:
+        logger.info("Flipping Dataset upside-down.")
+        dataset.data = dataset.data[::-1, :]
+        ds_area_extents[:, [1, 3]] = ds_area_extents[:, [3, 1]]
+
+    if target_eastright != current_eastright:
+        logger.info("Flipping Dataset left-to-right.")
+        dataset.data = dataset.data[:, ::-1]
+        ds_area_extents[:, [0, 2]] = ds_area_extents[:, [2, 0]]
+
+    # update the dataset area extent
+    # keeping the same id, description and proj_id, but should probably be changed to reflect the flipping
+    if len(ds_area_extents) == 1:
+        new_area_def = _change_area_extent_of_areadef(
+            dataset.attrs['area'], ds_area_extents[0])
+    else:
+        new_area_defs_to_stack = []
+        for n_area_def, area_def in enumerate(dataset.attrs['area'].defs):
+            new_area_defs_to_stack.append(_change_area_extent_of_areadef(
+                area_def, ds_area_extents[n_area_def]))
+
+        # flip the order of stacking if the area is upside down
+        if target_northup != current_northup:
+            new_area_defs_to_stack = new_area_defs_to_stack[::-1]
+        new_area_def = StackedAreaDefinition(*new_area_defs_to_stack)
+
+    dataset.attrs['area'] = new_area_def
+
+    return dataset
+
+
+class GEOFlippableFileYAMLReader(FileYAMLReader):
+    """Reader for flippable geostationary data."""
+
+    def _load_dataset_with_area(self, dsid, coords, **kwargs):
+        ds = super()._load_dataset_with_area(dsid, coords, **kwargs)
+        ds = _set_orientation(ds, **kwargs)
+        return ds
+
+
+class GEOSegmentYAMLReader(GEOFlippableFileYAMLReader):
     """Reader for segmented geostationary data.
 
     This reader pads the data to full geostationary disk if necessary.
@@ -904,7 +1011,7 @@ def create_filehandlers(self, filenames, fh_kwargs=None):
         return created_fhs
 
     @staticmethod
-    def _load_dataset(dsid, ds_info, file_handlers, dim='y', pad_data=True):
+    def _load_dataset(dsid, ds_info, file_handlers, dim='y', pad_data=True, **kwargs):
         """Load only a piece of the dataset."""
         if not pad_data:
             return FileYAMLReader._load_dataset(dsid, ds_info,
@@ -936,7 +1043,7 @@ def _load_dataset(dsid, ds_info, file_handlers, dim='y', pad_data=True):
         res.attrs = combined_info
         return res
 
-    def _load_area_def(self, dsid, file_handlers, pad_data=True):
+    def _load_area_def(self, dsid, file_handlers, pad_data=True, **kwargs):
         """Load the area definition of *dsid* with padding."""
         if not pad_data:
             return _load_area_def(dsid, file_handlers)