Merge pull request #2275 from adybbroe/multiscene_blend_with_weights

Multiscene blend with weights

doc/source/multiscene.rst

@@ -85,6 +85,51 @@ iteratively overlays the remaining datasets on top.
     >>> blended_scene = new_mscn.blend()
     >>> blended_scene.save_datasets()
 
+
+Stacking scenes using weights
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+It is also possible to blend scenes together in a bit more sophisticated manner
+using pixel based weighting instead of just stacking the scenes on top of each
+other as described above. This can for instance be useful to make a cloud
+parameter (cover, height, etc) composite combining cloud parameters derived
+from both geostationary and polar orbiting satellite data close in time and
+over a given area. This is useful for instance at high latitudes where
+geostationary data degrade quickly with latitude and polar data are more
+frequent.
+
+This weighted blending can be accomplished via the use of the builtin
+:func:`~functools.partial` function (see `Partial
+<https://docs.python.org/3/library/functools.html#partial-objects>`_) and the
+default :func:`~satpy.multiscene.stack` function. The
+:func:`~satpy.multiscene.stack` function can take the optional argument
+`weights` (`None` on default) which should be a sequence (of length equal to
+the number of scenes being blended) of arrays with pixel weights.
+
+The code below gives an example of how two cloud scenes can be blended using
+the satellite zenith angles to weight which pixels to take from each of the two
+scenes. The idea being that the reliability of the cloud parameter is higher
+when the satellite zenith angle is small.
+
+    >>> from satpy import Scene, MultiScene,  DataQuery
+    >>> from functools import partial
+    >>> from satpy.resample import get_area_def
+    >>> areaid = get_area_def("myarea")
+    >>> geo_scene = Scene(filenames=glob('/data/to/nwcsaf/geo/files/*nc'), reader='nwcsaf-geo')
+    >>> geo_scene.load(['ct'])
+    >>> polar_scene = Scene(filenames=glob('/data/to/nwcsaf/pps/noaa18/files/*nc'), reader='nwcsaf-pps_nc')
+    >>> polar_scene.load(['cma', 'ct'])
+    >>> mscn = MultiScene([geo_scene, polar_scene])
+    >>> groups = {DataQuery(name='CTY_group'): ['ct']}
+    >>> mscn.group(groups)
+    >>> resampled = mscn.resample(areaid, reduce_data=False)
+    >>> weights = [1./geo_satz, 1./n18_satz]
+    >>> stack_with_weights = partial(stack, weights=weights)
+    >>> blended = resampled.blend(blend_function=stack_with_weights)
+    >>> blended_scene.save_dataset('CTY_group', filename='./blended_stack_weighted_geo_polar.nc')
+
+
+
 Grouping Similar Datasets
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 

satpy/etc/readers/nwcsaf-pps_nc.yaml

@@ -126,19 +126,19 @@ datasets:
     name: ct_conditions
     file_type: nc_nwcsaf_ct
     coordinates: [lon, lat]
-    standard_name: cloudtype_conditions
+    standard_name: ct_conditions
 
   ct_quality:
     name: ct_quality
     file_type: nc_nwcsaf_ct
     coordinates: [lon, lat]
-    standard_name: cloudtype_quality
+    standard_name: ct_quality
 
   ct_status_flag:
     name: ct_status_flag
     file_type: nc_nwcsaf_ct
     coordinates: [lon, lat]
-    standard_name: cloudtype_status_flag
+    standard_name: ct_status_flag
 
   ct_pal:
     name: ct_pal

satpy/modifiers/angles.py

@@ -447,6 +447,7 @@ def _get_sensor_angles(data_arr: xr.DataArray) -> tuple[xr.DataArray, xr.DataArr
     sat_lon, sat_lat, sat_alt = get_satpos(data_arr, preference=preference)
     area_def = data_arr.attrs["area"]
     chunks = _geo_chunks_from_data_arr(data_arr)
+
     sata, satz = _get_sensor_angles_from_sat_pos(sat_lon, sat_lat, sat_alt,
                                                  data_arr.attrs["start_time"],
                                                  area_def, chunks)

satpy/multiscene.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-# Copyright (c) 2016-2019 Satpy developers
+# Copyright (c) 2016-2023 Satpy developers
 #
 # This file is part of satpy.
 #
@@ -20,6 +20,7 @@
 import copy
 import logging
 import warnings
+from datetime import datetime
 from queue import Queue
 from threading import Thread
 
@@ -45,14 +46,69 @@
 log = logging.getLogger(__name__)
 
 
-def stack(datasets):
-    """Overlay series of datasets on top of each other."""
+def stack(datasets, weights=None, combine_times=True):
+    """Overlay a series of datasets together.
+
+    By default, datasets are stacked on top of each other, so the last one applied is
+    on top. If a sequence of weights arrays are provided the datasets will
+    be combined according to those weights. The result will be a composite
+    dataset where the data in each pixel is coming from the dataset having the
+    highest weight.
+
+    """
+    if weights:
+        return _stack_weighted(datasets, weights, combine_times)
+
     base = datasets[0].copy()
     for dataset in datasets[1:]:
-        base = base.where(dataset.isnull(), dataset)
+        try:
+            base = base.where(dataset == dataset.attrs["_FillValue"], dataset)
+        except KeyError:
+            base = base.where(dataset.isnull(), dataset)
+
     return base
 
 
+def _stack_weighted(datasets, weights, combine_times):
+    """Stack datasets using weights."""
+    weights = set_weights_to_zero_where_invalid(datasets, weights)
+
+    indices = da.argmax(da.dstack(weights), axis=-1)
+    attrs = combine_metadata(*[x.attrs for x in datasets])
+
+    if combine_times:
+        if 'start_time' in attrs and 'end_time' in attrs:
+            attrs['start_time'], attrs['end_time'] = _get_combined_start_end_times(*[x.attrs for x in datasets])
+
+    dims = datasets[0].dims
+    weighted_array = xr.DataArray(da.choose(indices, datasets), dims=dims, attrs=attrs)
+    return weighted_array
+
+
+def set_weights_to_zero_where_invalid(datasets, weights):
+    """Go through the weights and set to pixel values to zero where corresponding datasets are invalid."""
+    for i, dataset in enumerate(datasets):
+        try:
+            weights[i] = xr.where(dataset == dataset.attrs["_FillValue"], 0, weights[i])
+        except KeyError:
+            weights[i] = xr.where(dataset.isnull(), 0, weights[i])
+
+    return weights
+
+
+def _get_combined_start_end_times(*metadata_objects):
+    """Get the start and end times attributes valid for the entire dataset series."""
+    start_time = datetime.now()
+    end_time = datetime.fromtimestamp(0)
+    for md_obj in metadata_objects:
+        if md_obj['start_time'] < start_time:
+            start_time = md_obj['start_time']
+        if md_obj['end_time'] > end_time:
+            end_time = md_obj['end_time']
+
+    return start_time, end_time
+
+
 def timeseries(datasets):
     """Expand dataset with and concatenate by time dimension."""
     expanded_ds = []
@@ -508,7 +564,7 @@ def _get_single_frame(self, ds, enh_args, fill_value):
         enh_args = enh_args.copy()  # don't change caller's dict!
         if "decorate" in enh_args:
             enh_args["decorate"] = self._format_decoration(
-                    ds, enh_args["decorate"])
+                ds, enh_args["decorate"])
         img = get_enhanced_image(ds, **enh_args)
         data, mode = img.finalize(fill_value=fill_value)
         if data.ndim == 3:
@@ -632,7 +688,7 @@ def _get_writers_and_frames(
             info_datasets = [scn.get(dataset_id) for scn in info_scenes]
             this_fn, shape, this_fill = self._get_animation_info(info_datasets, filename, fill_value=fill_value)
             data_to_write = self._get_animation_frames(
-                    all_datasets, shape, this_fill, ignore_missing, enh_args)
+                all_datasets, shape, this_fill, ignore_missing, enh_args)
 
             writer = imageio.get_writer(this_fn, **imio_args)
             frames[dataset_id] = data_to_write
@@ -703,8 +759,8 @@ def save_animation(self, filename, datasets=None, fps=10, fill_value=None,
             raise ImportError("Missing required 'imageio' library")
 
         (writers, frames) = self._get_writers_and_frames(
-                filename, datasets, fill_value, ignore_missing,
-                enh_args, imio_args={"fps": fps, **kwargs})
+            filename, datasets, fill_value, ignore_missing,
+            enh_args, imio_args={"fps": fps, **kwargs})
 
         client = self._get_client(client=client)
         # get an ordered list of frames

satpy/tests/multiscene_tests/__init__.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2018-2023 Satpy developers
+#
+# This file is part of satpy.
+#
+# satpy is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# satpy is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# satpy.  If not, see <http://www.gnu.org/licenses/>.
+"""Unit tests for Multiscene."""