Add support for temporal composites

doc/source/composites.rst

@@ -304,6 +304,43 @@ categories representing clear sky (cat1/cat5), cloudy features (cat2-cat4) and m
     >>> compositor = CategoricalDataCompositor('binary_cloud_mask', lut=lut)
     >>> composite = compositor([cloud_type])  # 0 - cat1/cat5, 1 - cat2/cat3/cat4, nan - cat0
 
+Time-dependent compositors
+--------------------------
+
+.. versionadded:: 0.43
+
+.. warning::
+
+   Time-dependent compositors are experimental.
+
+The built-in compositor :class:`TemporalRGB` creates an RGB consisting of
+three timesteps from the same channel.  When creating the compositor, the
+user must define the times using a :class:`~satpy.DataQuery`::
+
+    >>> temporal = TemporalRGB(
+    ...     "temporal_rgb",
+    ...     [DataQuery(wavelength=0.6, time=0),
+    ...      DataQuery(wavelength=0.6, time="-10 min"),
+    ...      DataQuery(wavelength=0.6, time="-20 min")])
+
+When calling the compositor, it should be passed a dataset that has a
+time dimension.  This can be obtained by blending a :class:`satpy.multiscene.MultiScene`
+with the :func:`satpy.multiscene.timeseries` function::
+
+    >>> ms = MultiScene.from_files(
+    ...         glob("/media/nas/x21308/MTG_test_data/2022_05_MTG_Testdata/RC007[012]/*BODY*.nc"),
+    ...         reader="fci_l1c_nc",
+    ...         group_keys=["repeat_cycle_in_day"])
+    >>> ms.load(["vis_06"])
+    >>> sc = ms.blend(blend_function=satpy.multiscene.timeseries)
+    >>> comp = temporal([sc["vis_06"]]*3)  # pass three times: same source for each timestep
+
+When defining a time-dependent composite in a configuration file (see below)
+it should be loaded directly from the :class:`~satpy.multiscene.MultiScene`, but will be
+available only after timeseries blending::
+
+    >>> ms.load(["temporal_rgb_vis06"])
+    >>> sc = ms.blend(blend_function=timeseries)
 
 Creating composite configuration files
 ======================================

satpy/composites/__init__.py

@@ -1,4 +1,4 @@
 # Copyright (c) 2015-2023 Satpy developers
 #
 # This file is part of satpy.
 #
@@ -22,6 +22,7 @@
 
 import dask.array as da
 import numpy as np
+import pandas as pd
 import xarray as xr
 from trollimage.colormap import Colormap
 
@@ -245,6 +246,16 @@
                       "'{}'".format(self.attrs['name']))
             raise IncompatibleAreas("Areas are different")
 
+    def _concat_datasets(self, projectables, mode):
+        try:
+            data = xr.concat(projectables, 'bands', coords='minimal')
+            data['bands'] = list(mode)
+        except ValueError as e:
+            LOG.debug("Original exception for incompatible areas: {}".format(str(e)))
+            raise IncompatibleAreas
+
+        return data
+
 
 class DifferenceCompositor(CompositeBase):
     """Make the difference of two data arrays."""
@@ -397,16 +408,6 @@
             return ''.join(data_arr.coords['bands'].values)
         return cls.modes[data_arr.sizes['bands']]
 
-    def _concat_datasets(self, projectables, mode):
-        try:
-            data = xr.concat(projectables, 'bands', coords='minimal')
-            data['bands'] = list(mode)
-        except ValueError as e:
-            LOG.debug("Original exception for incompatible areas: {}".format(str(e)))
-            raise IncompatibleAreas
-
-        return data
-
     def _get_sensors(self, projectables):
         sensor = set()
         for projectable in projectables:
@@ -1701,3 +1702,85 @@
 
         masked_projectable = projectable.where(lon_min_max)
         return super().__call__([masked_projectable], **info)
+
+
+class MissingTime(Exception):
+    """Raised when temporal composite building lacks time dimension."""
+
+
+class BaseTemporalCompositor(CompositeBase):
+    """Compositors combining multiple time steps.
+
+    Base class for compositors that combine inputs from two or more time steps.
+
+    To use this, the user must start with a :class:`MultiScene` and load the
+    temporal composite from there.  Composite generation will fail due to missing
+    temporal information in the containing scenes, but it will add the temporal
+    composite to the scene wishlists.  Now we run ``MultiScene.blend(timeseries)``,
+    which will create DataArrays with a time dimension.  After blending, the
+    blend method tries again to generate the composites, and on this second run
+    the generation should work.
+    """
+
+    def __init__(self, name, prerequisites, optional_prerequisites=None, **kwargs):
+        """Initialise a temporal compositor."""
+        self._ensure_prerequisites_have_times(prerequisites)
+        super().__init__(name, prerequisites, optional_prerequisites, **kwargs)
+
+    def _check_time_dimension(self, projectables):
+        """Make sure all projectables have a time dimension."""
+        for projectable in projectables:
+            if "time" not in projectable.dims:
+                raise MissingTime(
+                    "Creating temporal composite needs time dimensions. "
+                    "Typically, this comes from starting with a MultiScene "
+                    "and then performing timeseries blending.")
+
+    def _ensure_prerequisites_have_times(self, prerequisites):
+        """Make sure all prerequisites have times defined."""
+        for preq in prerequisites:
+            if "time" not in preq.to_dict():
+                raise KeyError(
+                    "In a temporal composite, all prerequisites should have "
+                    "time information defined.  Time information is missing for "
+                    f"{preq!s}")
+
+    def _apply_temporal_prerequisites(self, projectables):
+        """Apply temporal prerequisites.
+
+        The generic Satpy composite loading logic does not understand time
+        dependencies.  When :meth:`MultiScene.blend` reruns the composite
+        loading for a composite combining the same channel for three time
+        steps, we will be passed the same 3-D data array three times.
+
+        This method matches the temporal prerequisites as defined in the
+        compositor with the time dimension associated with the dataarrays.
+        """
+        # reference time is the newest time
+        ref_time = projectables[0]["time"][-1]
+        new_projectables = []
+        for (dq, proj) in zip(self.attrs["prerequisites"], projectables):
+            new_da = proj.sel(time=ref_time+pd.Timedelta(dq["time"]), method='nearest')
+            new_projectables.append(new_da)
+        return new_projectables
+
+
+class TemporalRGB(BaseTemporalCompositor):
+    """Make an RGB where different timesteps of the same channel go in different bands.
+
+    See the note in the parent class :class:`BaseTemporalCompositor` for
+    usage instructions.
+    """
+
+    def __call__(self, projectables, optional_datasets=None, **info):
+        """Build the composite."""
+        projectables = self.match_data_arrays(projectables)
+        self._check_time_dimension(projectables)
+        new_projectables = self._apply_temporal_prerequisites(projectables)
+        # need to drop the time coordinate, otherwise we can't concatenate
+        # the bands
+        dataset = self._concat_datasets(
+            [proj.drop_vars("time") for proj in new_projectables], "RGB")
+        dataset.attrs.update(self.attrs)
+        dataset.attrs.update(info)
+        return dataset

satpy/etc/composites/visir.yaml

@@ -642,3 +642,20 @@ composites:
       - wavelength: 3.9
       - wavelength: 6.2
       - wavelength: 1.6
+
+  temporal_rgb_vis06:
+    description: >
+      An RGB to highlight changes at 0.64 µm. To use this, start with a
+      MultiScene containing three scenes that are 10 minutes apart, such
+      as three subsequent FCI or ABI FD scenes. Load the composite from
+      the MultiScene, then apply timeseries blending. The resulting Scene
+      should have the composites.
+    compositor: !!python/name:satpy.composites.TemporalRGB
+    prerequisites:
+      - wavelength: 0.64
+        time: 0
+      - wavelength: 0.64
+        time: -10 min
+      - wavelength: 0.64
+        time: -20 min
+    standard_name: temporal_rgb_vis06

satpy/multiscene/_multiscene.py

@@ -359,8 +359,19 @@ def blend(
             datasets = [scn[ds_id] for scn in self.scenes if ds_id in scn]
             new_scn[ds_id] = blend_function(datasets)
 
+        new_scn._wishlist = self._shared_wishlist()
+        # without copying the dependency tree, there is a KeyError in dependency_tree.trunk
+        new_scn._dependency_tree = self.first_scene._dependency_tree.copy()
+        new_scn.generate_possible_composites(True)
         return new_scn
 
+    def _shared_wishlist(self):
+        """Get shared wishlist."""
+        shared_wishlist = self.scenes[0]._wishlist
+        for scene in self.scenes[1:]:
+            shared_wishlist &= scene._wishlist
+        return shared_wishlist
+
     def group(self, groups):
         """Group datasets from the multiple scenes.
 

satpy/scene.py

@@ -1,4 +1,4 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # Copyright (c) 2010-2022 Satpy developers
 #
@@ -28,7 +28,7 @@
 from pyresample.geometry import AreaDefinition, BaseDefinition, SwathDefinition
 from xarray import DataArray
 
-from satpy.composites import IncompatibleAreas
+from satpy.composites import IncompatibleAreas, MissingTime
 from satpy.composites.config_loader import load_compositor_configs_for_sensors
 from satpy.dataset import DataID, DataQuery, DatasetDict, combine_metadata, dataset_walker, replace_anc
 from satpy.dependency_tree import DependencyTree
@@ -1285,7 +1285,7 @@
 
         missing_str = ", ".join(str(x) for x in missing)
         LOG.warning("The following datasets were not created and may require "
-                    "resampling to be generated: {}".format(missing_str))
+                    "resampling or temporal blending to be generated: {}".format(missing_str))
 
     def unload(self, keepables=None):
         """Unload all unneeded datasets.
@@ -1521,8 +1521,9 @@
                 self._wishlist.remove(comp_node.name)
                 self._wishlist.add(cid)
             self._dependency_tree.update_node_name(comp_node, cid)
-        except IncompatibleAreas:
-            LOG.debug("Delaying generation of %s because of incompatible areas", str(compositor.id))
+        except (IncompatibleAreas, MissingTime):
+            LOG.debug("Delaying generation of %s because of incompatible areas "
+                      "or missing time dimension", str(compositor.id))
             preservable_datasets = set(self._datasets.keys())
             prereq_ids = set(p.name for p in prereqs)
             opt_prereq_ids = set(p.name for p in optional_prereqs)

satpy/tests/multiscene_tests/test_temporal_composites.py

@@ -0,0 +1,146 @@
+# Copyright (c) 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/>.
+
+"""Tests for loading temporal composites."""
+
+import datetime
+
+import numpy as np
+import pytest
+import xarray as xr
+
+composite_definition = """sensor_name: visir
+
+composites:
+  temporal:
+    compositor: !!python/name:satpy.composites.TemporalRGB
+    prerequisites:
+      - name: ir
+        time: 0
+      - name: ir
+        time: -10 min
+      - name: ir
+        time: -20 min
+    standard_name: temporal
+"""
+
+
+@pytest.fixture
+def fake_config(tmp_path):
+    """Make a configuration path with a temporal composite definition."""
+    confdir = tmp_path / "etc"
+    conffile = tmp_path / "etc" / "composites" / "visir.yaml"
+    conffile.parent.mkdir(parents=True)
+    with conffile.open(mode="wt", encoding="ascii") as fp:
+        fp.write(composite_definition)
+    return confdir
+
+
+@pytest.fixture
+def fake_dataset():
+    """Create minimal fake Satpy CF NC dataset."""
+    ds = xr.Dataset()
+    nx = ny = 4
+    ds["ir"] = xr.DataArray(
+            np.zeros((nx, ny)),
+            dims=("y", "x"),
+            attrs={"sensor": "visir"})
+    return ds
+
+
+@pytest.fixture
+def fake_dataset_3d():
+    """Create fake 3D dataset."""
+    ds = xr.Dataset()
+    nt = 3
+    nx = ny = 4
+    times = np.array(
+            ["2017-09-01T10:00:00"], dtype="M8[s]") + np.array(
+                    [0, 600, 1200], dtype="m8[s]")
+    ds["ir"] = xr.DataArray(
+            np.zeros((nt, nx, ny)),
+            dims=("time", "y", "x"),
+            attrs={"sensor": "visir"},
+            coords={"time": times})
+    return ds
+
+
+@pytest.fixture
+def fake_files(tmp_path, fake_dataset):
+    """Make fake files for the Satpy CF reader."""
+    start_time = datetime.datetime(2050, 5, 3, 12, 0, 0, tzinfo=None)
+    delta = datetime.timedelta(minutes=10)
+    n_timesteps = 5
+    offsets = [-2.2, 1.3, 0.5, -0.9, 0.7]
+    ofs = []
+    for i in range(n_timesteps):
+        begin = start_time + i*delta + datetime.timedelta(seconds=offsets[i])
+        end = start_time + (i+1)*delta + datetime.timedelta(seconds=offsets[i])
+        of = tmp_path / f"Meteosat99-imager-{begin:%Y%m%d%H%M%S}-{end:%Y%m%d%H%M%S}.nc"
+        fd = fake_dataset.copy()
+        fd["ir"][...] = i
+        fd["ir"].attrs["start_time"] = f"{start_time:%Y-%m-%dT%H:%M:%S}"
+        fd.to_netcdf(of)
+        ofs.append(of)
+    return ofs
+
+
+@pytest.fixture
+def time_compositor():
+    """Construct a time compositor."""
+    from satpy import DataQuery
+    from satpy.composites import TemporalRGB
+    return TemporalRGB(
+            "temporal_rgb",
+            [DataQuery(wavelength=0.6, time=0),
+             DataQuery(wavelength=0.6, time="-10 min"),
+             DataQuery(wavelength=0.6, time="-20 min")])
+
+
+def test_load_temporal_composite(fake_files, fake_config):
+    """Test loading a temporal composite."""
+    from satpy import config
+    from satpy.multiscene import MultiScene, timeseries
+    with config.set(config_path=[fake_config]):
+        ms = MultiScene.from_files(fake_files, reader="satpy_cf_nc")
+        ms.load(["temporal"])
+        sc = ms.blend(blend_function=timeseries)
+        assert sc["temporal"].shape == (3, 4, 4)
+        np.testing.assert_array_equal(sc["temporal"][0, :, :], np.full((4, 4), 4))
+        np.testing.assert_array_equal(sc["temporal"][1, :, :], np.full((4, 4), 3))
+        np.testing.assert_array_equal(sc["temporal"][2, :, :], np.full((4, 4), 2))
+
+
+def test_init_compositor():
+    """Initialise the temporal compositor."""
+    from satpy import DataQuery
+    from satpy.composites import TemporalRGB
+    TemporalRGB(
+            "temporal_rgb",
+            [DataQuery(wavelength=0.6, time=0),
+             DataQuery(wavelength=0.6, time="-10 min"),
+             DataQuery(wavelength=0.6, time="-20 min")])
+    with pytest.raises(KeyError):
+        TemporalRGB(
+            "temporal_rgb",
+            [DataQuery(wavelength=0.6, time=0),
+             DataQuery(wavelength=0.6),
+             DataQuery(wavelength=0.6, time="-20 min")])
+
+
+def test_call_compositor(time_compositor, fake_dataset_3d):
+    """Test calling the temporal compositor."""
+    time_compositor([fake_dataset_3d]*3)