Implement generic set_orientation method for geostationary datasets

doc/source/readers.rst

@@ -68,7 +68,7 @@ the specified parameters. So the above ``load`` call would load the ``0.6``
 (a visible/reflectance band) radiance data and ``10.8`` (an IR band)
 brightness temperature data.
 
-For geostatinary satellites that have the individual channel data
+For geostationary satellites that have the individual channel data
 separated to several files (segments) the missing segments are padded
 by default to full disk area.  This is made to simplify caching of
 resampling look-up tables (see :doc:`resample` for more information).
@@ -77,6 +77,17 @@ loading datasets::
 
     >>> scn.load([0.6, 10.8], pad_data=False)
 
+For geostationary products, where the imagery is stored in the files in a flipped orientation
+(e.g. MSG SEVIRI L1.5 data which is flipped upside-down and left-right), the keyword argument
+``upper_right_corner`` can be passed into the load call to automatically flip the datasets to the
+wished orientation. Accepted argument values are ``'NE'``, ``'NW'``, ``'SE'``, ``'SW'``,
+and ``'native'``.
+By default, no flipping is applied (corresponding to ``upper_right_corner='native'``) and
+the data is delivered in the original format. To get the data in the common upright orientation,
+load the datasets using e.g.::
+
+    >>> scn.load(['VIS008'], upper_right_corner='NE')
+
 .. note::
 
     If a dataset could not be loaded there is no exception raised. You must

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

@@ -45,7 +45,6 @@
 from trollsift.parser import globify, parse
 from pyresample.geometry import AreaDefinition
 
-
 logger = logging.getLogger(__name__)
 
 
@@ -256,10 +255,10 @@ def load_ds_ids_from_config(self):
                     # but is still considered 1 option
                     # it also needs to be a tuple so it can be used in
                     # a dictionary key (DatasetID)
-                    id_kwargs.append((tuple(val), ))
+                    id_kwargs.append((tuple(val),))
                 elif key == "modifiers" and val is None:
                     # empty modifiers means no modifiers applied
-                    id_kwargs.append((tuple(), ))
+                    id_kwargs.append((tuple(),))
                 elif isinstance(val, (list, tuple, set)):
                     # this key has multiple choices
                     # (ex. 250 meter, 500 meter, 1000 meter resolutions)
@@ -269,7 +268,7 @@ def load_ds_ids_from_config(self):
                 else:
                     # this key only has one choice so make it a one
                     # item iterable
-                    id_kwargs.append((val, ))
+                    id_kwargs.append((val,))
             for id_params in itertools.product(*id_kwargs):
                 dsid = DatasetID(*id_params)
                 ids.append(dsid)
@@ -789,7 +788,7 @@ def _load_dataset_with_area(self, dsid, coords, **kwargs):
             ds = add_crs_xy_coords(ds, area)
         return ds
 
-    def _load_ancillary_variables(self, datasets):
+    def _load_ancillary_variables(self, datasets, **kwargs):
         """Load the ancillary variables of `datasets`."""
         all_av_ids = set()
         for dataset in datasets.values():
@@ -809,7 +808,7 @@ def _load_ancillary_variables(self, datasets):
         if not all_av_ids:
             return
         if loadable_av_ids:
-            self.load(loadable_av_ids, previous_datasets=datasets)
+            self.load(loadable_av_ids, previous_datasets=datasets, **kwargs)
 
         for dataset in datasets.values():
             new_vars = []
@@ -874,7 +873,7 @@ def load(self, dataset_keys, previous_datasets=None, **kwargs):
                 all_datasets[dsid] = ds
                 if dsid in dsids:
                     datasets[dsid] = ds
-        self._load_ancillary_variables(all_datasets)
+        self._load_ancillary_variables(all_datasets, **kwargs)
 
         return datasets
 
@@ -889,7 +888,165 @@ def _load_area_def(dsid, file_handlers):
     return final_area.squeeze()
 
 
-class GEOSegmentYAMLReader(FileYAMLReader):
+def _set_orientation(dataset, upper_right_corner):
+    """Set the orientation of geostationary datasets.
+
+    Allows to flip geostationary imagery when loading the datasets.
+    Example call: scn.load(['VIS008'], upper_right_corner='NE')
+
+    Args:
+        dataset: Dataset to be flipped.
+        upper_right_corner (str): Direction of the upper right corner of the image after flipping.
+                                Possible options are 'NW', 'NE', 'SW', 'SE', or 'native'.
+                                The common upright image orientation corresponds to 'NE'.
+                                Defaults to 'native' (no flipping is applied).
+
+    """
+    # do some checks and early returns
+    if upper_right_corner == 'native':
+        logger.debug("Requested orientation for Dataset {} is 'native' (default). "
+                     "No flipping is applied.".format(dataset.attrs.get('name')))
+        return dataset
+
+    if upper_right_corner not in ['NW', 'NE', 'SE', 'SW', 'native']:
+        raise ValueError("Target orientation for Dataset {} not recognized. "
+                         "Kwarg upper_right_corner should be "
+                         "'NW', 'NE', 'SW', 'SE' or 'native'.".format(dataset.attrs.get('name', 'unknown_name')))
+
+    if 'area' not in dataset.attrs:
+        logger.info("Dataset {} is missing the area attribute "
+                    "and will not be flipped.".format(dataset.attrs.get('name', 'unknown_name')))
+        return dataset
+
+    projection_type = _get_projection_type(dataset.attrs['area'])
+    accepted_geos_proj_types = ['Geostationary Satellite (Sweep Y)', 'Geostationary Satellite (Sweep X)']
+    if projection_type not in accepted_geos_proj_types:
+        logger.info("Dataset {} is not in one of the known geostationary projections {} "
+                    "and cannot be flipped.".format(dataset.attrs.get('name', 'unknown_name'),
+                                                    accepted_geos_proj_types))
+        return dataset
+
+    target_eastright, target_northup = _get_target_scene_orientation(upper_right_corner)
+
+    area_extents_to_update = _get_dataset_area_extents_array(dataset.attrs['area'])
+    current_eastright, current_northup = _get_current_scene_orientation(area_extents_to_update)
+
+    if target_northup == current_northup and target_eastright == current_eastright:
+        logger.info("Dataset {} is already in the target orientation "
+                    "and will not be flipped.".format(dataset.attrs.get('name', 'unknown_name')))
+        return dataset
+
+    if target_northup != current_northup:
+        dataset, area_extents_to_update = _flip_dataset_data_and_area_extents(dataset, area_extents_to_update,
+                                                                              'upsidedown')
+    if target_eastright != current_eastright:
+        dataset, area_extents_to_update = _flip_dataset_data_and_area_extents(dataset, area_extents_to_update,
+                                                                              'leftright')
+
+    dataset.attrs['area'] = _get_new_flipped_area_definition(dataset.attrs['area'], area_extents_to_update,
+                                                             flip_areadef_stacking=target_northup != current_northup)
+
+    return dataset
+
+
+def _get_projection_type(dataset_area_attr):
+    """Get the projection type from the crs coordinate operation method name."""
+    if isinstance(dataset_area_attr, StackedAreaDefinition):
+        # assumes all AreaDefinitions in a tackedAreaDefinition have the same projection
+        area_crs = dataset_area_attr.defs[0].crs
+    else:
+        area_crs = dataset_area_attr.crs
+
+    return area_crs.coordinate_operation.method_name
+
+
+def _get_target_scene_orientation(upper_right_corner):
+    """Get the target scene orientation from the target upper_right_corner.
+
+    'NE' corresponds to target_eastright and target_northup being True.
+    """
+    if upper_right_corner in ['NW', 'NE']:
+        target_northup = True
+    else:
+        target_northup = False
+
+    if upper_right_corner in ['NE', 'SE']:
+        target_eastright = True
+    else:
+        target_eastright = False
+
+    return target_eastright, target_northup
+
+
+def _get_dataset_area_extents_array(dataset_area_attr):
+    """Get dataset area extents in a numpy array for further flipping."""
+    if isinstance(dataset_area_attr, StackedAreaDefinition):
+        # array of area extents if the Area is a StackedAreaDefinition
+        area_extents_to_update = np.asarray([list(area_def.area_extent) for area_def in dataset_area_attr.defs])
+    else:
+        # array with a single item if Area is in one piece
+        area_extents_to_update = np.asarray([list(dataset_area_attr.area_extent)])
+    return area_extents_to_update
+
+
+def _get_current_scene_orientation(area_extents_to_update):
+    """Get the current scene orientation from the area_extents."""
+    # assumes all AreaDefinitions inside a StackedAreaDefinition have the same orientation
+    current_northup = area_extents_to_update[0, 3] - area_extents_to_update[0, 1] > 0
+    current_eastright = area_extents_to_update[0, 2] - area_extents_to_update[0, 0] > 0
+
+    return current_eastright, current_northup
+
+
+def _flip_dataset_data_and_area_extents(dataset, area_extents_to_update, flip_direction):
+    """Flip the data and area extents array for a dataset."""
+    logger.info("Flipping Dataset {} {}.".format(dataset.attrs.get('name', 'unknown_name'), flip_direction))
+    if flip_direction == 'upsidedown':
+        dataset.data = dataset.data[::-1, :]
+        area_extents_to_update[:, [1, 3]] = area_extents_to_update[:, [3, 1]]
+    elif flip_direction == 'leftright':
+        dataset.data = dataset.data[:, ::-1]
+        area_extents_to_update[:, [0, 2]] = area_extents_to_update[:, [2, 0]]
+    else:
+        raise ValueError("Flip direction not recognized. Should be either 'upsidedown' or 'leftright'.")
+
+    return dataset, area_extents_to_update
+
+
+def _get_new_flipped_area_definition(dataset_area_attr, area_extents_to_update, flip_areadef_stacking):
+    """Get a new area definition with updated area_extents for flipped geostationary datasets."""
+    if len(area_extents_to_update) == 1:
+        # just update the area extents using the AreaDefinition copy method
+        new_area_def = dataset_area_attr.copy(area_extent=area_extents_to_update[0])
+    else:
+        # update the stacked AreaDefinitions singularly
+        new_area_defs_to_stack = []
+        for n_area_def, area_def in enumerate(dataset_area_attr.defs):
+            new_area_defs_to_stack.append(area_def.copy(area_extent=area_extents_to_update[n_area_def]))
+
+        # flip the order of stacking if the area is upside down
+        if flip_areadef_stacking:
+            new_area_defs_to_stack = new_area_defs_to_stack[::-1]
+
+        # regenerate the StackedAreaDefinition
+        new_area_def = StackedAreaDefinition(*new_area_defs_to_stack)
+
+    return new_area_def
+
+
+class GEOFlippableFileYAMLReader(FileYAMLReader):
+    """Reader for flippable geostationary data."""
+
+    def _load_dataset_with_area(self, dsid, coords, upper_right_corner='native', **kwargs):
+        ds = super(GEOFlippableFileYAMLReader, self)._load_dataset_with_area(dsid, coords, **kwargs)
+
+        if ds is not None:
+            ds = _set_orientation(ds, upper_right_corner)
+
+        return ds
+
+
+class GEOSegmentYAMLReader(GEOFlippableFileYAMLReader):
     """Reader for segmented geostationary data.
 
     This reader pads the data to full geostationary disk if necessary.