diff --git a/doc/source/examples/index.rst b/doc/source/examples/index.rst
index d86e12431e..90cde0681f 100644
--- a/doc/source/examples/index.rst
+++ b/doc/source/examples/index.rst
@@ -15,6 +15,7 @@ as explanations in the various sections of this documentation.
 
     fci_l1c_natural_color
     vii_l1b_nc
+    mean_time
 
 .. list-table::
     :header-rows: 1
@@ -45,3 +46,5 @@ as explanations in the various sections of this documentation.
       - Generate Natural Color RGB from Meteosat Third Generation (MTG) FCI Level 1c data
     * - :doc:`Reading EPS-SG Visible and Infrared Imager (VII) with Pytroll <vii_l1b_nc>`
       - Read and visualize EPS-SG VII L1B test data and save it to an image
+    * - :doc:`Storing approximate measurement time <mean_time>`
+      - Read, resample, and store the approximate measurement time for supported readers and writers.
diff --git a/doc/source/examples/mean_time.rst b/doc/source/examples/mean_time.rst
new file mode 100644
index 0000000000..64bfc266f4
--- /dev/null
+++ b/doc/source/examples/mean_time.rst
@@ -0,0 +1,44 @@
+Tracking measurement time
+=========================
+
+.. versionadded:: v0.58
+
+For some readers and writers, it is possible to keep track of pixel-level
+measurement times and store the average measurement time in the metadata
+for the resampled image.  This can be stored in the filename or the file
+headers (or both).  It is only supported for selected readers and writers.
+
+By default, Satpy does not keep track of measurement times.  To keep track
+of measurement times, we must first tell the reader to add such times to
+the metadata of each dataset.  With supported readers, this can be done
+by passing ``reader_kwargs={"track_time": True}`` to :meth:`~satpy.scene.Scene`:
+
+.. code-block:: python
+
+   sc = Scene(filenames={"seviri_l1b_hrit": seviri_files}, reader_kwargs={"track_time": True})
+   sc.load(["IR_108"])
+
+The time is stored as a coordinate:
+
+.. code-block:: python
+
+   sc["IR_108"].coords["time"]
+
+To retain it upon resampling, pass ``resample_coords=True`` to :meth:`~satpy.scene.Scene.resample`:
+
+.. code-block:: python
+
+   ls = sc.resample("eurol", resample_coords=True)
+
+For supported writers, it can be stored in the filename by passing ``dynamic_fields={"mean_time"}``
+to :meth:`~satpy.scene.Scene.save_datasets`:
+
+.. code-block:: python
+
+   ls.save_datasets(
+       writer="geotiff",
+       filename="{platform_name}-{sensor}-{name}-{area.area_id}-{start_time:%Y%m%d%H%M}-{mean_time:%Y%m%d%H%M%S}.tif",
+       dynamic_fields={"mean_time"})
+
+For supported writers, valid time may also be written to the headers.  Consult
+the documentation of your writer for details.
diff --git a/doc/source/writing.rst b/doc/source/writing.rst
index 9cd6145bbe..bc84b48c61 100644
--- a/doc/source/writing.rst
+++ b/doc/source/writing.rst
@@ -17,10 +17,15 @@ One common parameter across almost all Writers is ``filename`` and
     ...     filename="{name}_{start_time:%Y%m%d_%H%M%S}.tif",
     ...     base_dir="/tmp/my_ouput_dir")
 
-.. versionchanged:: 0.10
-
-    The `file_pattern` keyword argument was renamed to `filename` to match
-    the `save_dataset` method"s keyword argument.
+The ``filename`` argument can specify Python string formatting fields.
+Those fields are mostly filled by attributes available or the individual
+datasets.  Following Python string formatting rules, attributes of
+attributes can be referenced as well, for example ``area.name``.  In
+addition to dataset attributes, some reader/writer combinations support
+dynamically calculated field values.  This currently exists only for
+``{mean_time}`` if the keyword argument ``dynamic_fields={"mean_time"}``
+is passed to :meth:`~satpy.scene.Scene.save_datasets`.  See the
+:doc:`example on storing valid time </examples/mean_time>` for details.
 
 .. _writer_table:
 
@@ -53,7 +58,7 @@ One common parameter across almost all Writers is ``filename`` and
       -
     * - GeoTIFF with NinJo tags (from NinJo 7)
       - :class:`ninjogeotiff <satpy.writers.ninjogeotiff.NinJoGeoTIFFWriter>`
-      - Beta
+      - Nominal
       -
 
 Available Writers
diff --git a/satpy/cf/decoding.py b/satpy/cf/decoding.py
index 2515f6bd38..040d74bd23 100644
--- a/satpy/cf/decoding.py
+++ b/satpy/cf/decoding.py
@@ -22,6 +22,8 @@
 import datetime as dt
 import json
 
+import numpy as np
+
 
 def decode_attrs(attrs):
     """Decode CF-encoded attributes to Python object.
@@ -74,3 +76,39 @@ def _str2datetime(string):
         return dt.datetime.fromisoformat(string)
     except (TypeError, ValueError):
         return None
+
+
+def lazy_decode_cf_time(xrda_encoded):
+    """Lazily decode CF-encoded time in limited situations.
+
+    This is a restricted alternative to xarray.coding.times.decode_cf_datetime
+    that avoids dask computations on the values to be decoded.  It is
+    restricted to standard calendars (proleptic gregorian).
+
+    Args:
+        xrda_encoded (array-like):
+            Xarray data array with units attribute.  The units attribute should
+            be a string describing the time units using "x since timestamp",
+            UDUNITS-style.
+    """
+    # An early iteration of this function was written with the
+    # assistance of GPT-5.4.
+
+    (unit, ref_str) = xrda_encoded.attrs["units"].split(" since ")
+    ref = np.datetime64(ref_str)
+
+    unit_map = {
+        "days": "timedelta64[D]",
+        "day": "timedelta64[D]",
+        "hours": "timedelta64[h]",
+        "hour": "timedelta64[h]",
+        "minutes": "timedelta64[m]",
+        "minute": "timedelta64[m]",
+        "seconds": "timedelta64[s]",
+        "second": "timedelta64[s]",
+        "milliseconds": "timedelta64[ms]",
+        "microseconds": "timedelta64[us]",
+        "nanoseconds": "timedelta64[ns]",
+    }
+
+    return ref + xrda_encoded.astype(unit_map[unit.lower()])
diff --git a/satpy/composites/arithmetic.py b/satpy/composites/arithmetic.py
index 668b9a7972..0353504e46 100644
--- a/satpy/composites/arithmetic.py
+++ b/satpy/composites/arithmetic.py
@@ -34,7 +34,8 @@ def __call__(self, projectables, nonprojectables=None, **attrs):
         """Generate the composite."""
         if len(projectables) != 2:
             raise ValueError("Expected 2 datasets, got %d" % (len(projectables),))
-        projectables = self.match_data_arrays(projectables)
+        projectables = self.match_data_arrays(projectables,
+                                              drop_coordinates=False)
         info = combine_metadata(*projectables)
         info["name"] = self.attrs["name"]
         info.update(self.attrs)  # attrs from YAML/__init__
@@ -52,7 +53,8 @@ def __call__(self, projectables, nonprojectables=None, **info):
         """Generate the composite."""
         if len(projectables) != 2:
             raise ValueError("Expected 2 datasets, got %d" % (len(projectables),))
-        projectables = self.match_data_arrays(projectables)
+        projectables = self.match_data_arrays(projectables,
+                                              drop_coordinates=False)
         info = combine_metadata(*projectables)
         info.update(self.attrs)
 
@@ -68,7 +70,8 @@ def __call__(self, projectables, nonprojectables=None, **info):
         """Generate the composite."""
         if len(projectables) != 2:
             raise ValueError("Expected 2 datasets, got %d" % (len(projectables),))
-        projectables = self.match_data_arrays(projectables)
+        projectables = self.match_data_arrays(projectables,
+                                              drop_coordinates=False)
         info = combine_metadata(*projectables)
         info["name"] = self.attrs["name"]
 
diff --git a/satpy/composites/core.py b/satpy/composites/core.py
index 1b124eb88a..eb6b64cbe8 100644
--- a/satpy/composites/core.py
+++ b/satpy/composites/core.py
@@ -23,7 +23,6 @@
 from typing import Optional, Sequence
 
 import dask.array as da
-import numpy as np
 import xarray as xr
 
 from satpy.dataset import DataID, combine_metadata
@@ -35,8 +34,6 @@
 NEGLIGIBLE_COORDS = ["time"]
 """Keywords identifying non-dimensional coordinates to be ignored during composite generation."""
 
-TIME_COMPATIBILITY_TOLERANCE = np.timedelta64(1, "s")
-
 
 class IncompatibleAreas(Exception):
     """Error raised upon compositing things of different shapes."""
@@ -124,7 +121,8 @@ def _add_missing_modifier(self, key, destination, origin):
         elif origin.get(key) is not None:
             destination[key] = origin[key]
 
-    def match_data_arrays(self, data_arrays: Sequence[xr.DataArray]) -> list[xr.DataArray]:
+    def match_data_arrays(self, data_arrays: Sequence[xr.DataArray],
+                          drop_coordinates: bool=True) -> list[xr.DataArray]:
         """Match data arrays so that they can be used together in a composite.
 
         For the purpose of this method, "can be used together" means:
@@ -133,12 +131,14 @@ def match_data_arrays(self, data_arrays: Sequence[xr.DataArray]) -> list[xr.Data
         - Either all arrays should have an area, or none should.
         - If all have an area, the areas should be all the same.
 
-        In addition, negligible non-dimensional coordinates are dropped (see
+        In addition, negligible non-dimensional coordinates can be dropped (see
         :meth:`drop_coordinates`) and dask chunks are unified (see
         :func:`satpy.utils.unify_chunks`).
 
         Args:
             data_arrays: Arrays to be checked
+            drop_coordinates: If true, drop non-dimensional coordinates.
+                If false, unify them too.
 
         Returns:
             Arrays with negligible non-dimensional coordinates removed.
@@ -150,7 +150,7 @@ def match_data_arrays(self, data_arrays: Sequence[xr.DataArray]) -> list[xr.Data
                 If some, but not all data arrays lack an area attribute.
         """
         self.check_geolocation(data_arrays)
-        new_arrays = self.drop_coordinates(data_arrays)
+        new_arrays = self.drop_coordinates(data_arrays) if drop_coordinates else self.combine_coordinates(data_arrays)
         new_arrays = self.align_geo_coordinates(new_arrays)
         new_arrays = list(unify_chunks(*new_arrays))
         return new_arrays
@@ -210,6 +210,8 @@ def drop_coordinates(data_arrays: Sequence[xr.DataArray]) -> list[xr.DataArray]:
         dimension.  Negligible coordinates are defined in the
         :attr:`NEGLIGIBLE_COORDS` module attribute.
 
+        See also: :meth:`combine_coordinates`.
+
         Args:
             data_arrays: Arrays to be checked
         """
@@ -225,6 +227,47 @@ def drop_coordinates(data_arrays: Sequence[xr.DataArray]) -> list[xr.DataArray]:
 
         return new_arrays
 
+    def combine_coordinates(self, data_arrays: Sequence[xr.DataArray]) -> list[xr.DataArray]:
+        """Combine time coordinates.
+
+        Combine coordinates if they do not correspond to any
+        dimension.  Currently only supports time coordinates.
+
+        See also: :meth:`drop_coordinates`.
+
+        Args:
+            data_arrays: Arrays to be checked
+        """
+        new_time = self.combine_times(data_arrays)
+        return [data.assign_coords({"time": new_time}) if "time" in data.coords else data
+                    for data in data_arrays]
+
+    @staticmethod
+    def combine_times(projectables):
+        """Get a combined time coordinate between projectables.
+
+        Returns the arithmetic mean between the available time coordinates,
+        provided they have a common dimensionality and units.  If no projectables
+        have time coordinates, return None.  Time coordinates should be CF-encoded,
+        i.e. have a numeric dtype and a units attribute.
+        """
+        _verify_times(projectables)
+        timed_projectables = [proj for proj in projectables
+                              if hasattr(proj, "coords")
+                              and "time" in proj.coords]
+        if len(timed_projectables) > 0:
+            with xr.set_options(keep_attrs=True):
+                # when the time coordinates are different, can't use xarray to
+                # sum them without this leading to expensive computations!
+                da_new_time = sum([x.coords["time"].data for x in timed_projectables]) / len(timed_projectables)
+                first = timed_projectables[0].coords["time"]
+                xr_new_time = xr.DataArray(
+                        da_new_time,
+                        dims=first.dims,
+                        attrs=first.attrs,
+                        coords={"y": first.coords["y"], "x": first.coords["x"]})
+            return xr_new_time.assign_coords(time=(first.dims, da_new_time))
+
     @staticmethod
     def align_geo_coordinates(data_arrays: Sequence[xr.DataArray]) -> list[xr.DataArray]:
         """Align DataArrays along geolocation coordinates.
@@ -428,7 +471,7 @@ def _concat_datasets(self, projectables, mode):
             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)))
+            LOG.exception("Original exception for incompatible areas: {}".format(str(e)))
             raise IncompatibleAreas("Areas do not match.")
 
         return data
@@ -483,17 +526,14 @@ def _get_mode(self, attrs, num):
         return mode
 
     def _check_datasets_and_data(self, datasets, mode):
+        time = self.combine_times(datasets)
         datasets = self.match_data_arrays(datasets)
         data = self._concat_datasets(datasets, mode)
         # Skip masking if user wants it or a specific alpha channel is given.
         if self.common_channel_mask and mode[-1] != "A":
             data = data.where(data.notnull().all(dim="bands"))
-        # if inputs have a time coordinate that may differ slightly between
-        # themselves then find the mid time and use that as the single
-        # time coordinate value
-        time = check_times(datasets)
-        if time is not None and "time" in data.dims:
-            data["time"] = [time]
+        if time is not None:
+            data = data.assign_coords({"time": time})
 
         return datasets, data
 
@@ -518,39 +558,58 @@ def _get_updated_attrs(self, datasets, attrs, mode):
         return new_attrs
 
 
-def check_times(projectables):
-    """Check that *projectables* have compatible times."""
-    times = []
-    for proj in projectables:
-        status = _collect_time_from_proj(times, proj)
-        if not status:
-            break
-    else:
-        return _get_average_time(times)
-
-
-def _collect_time_from_proj(times, proj):
-    status = False
-    try:
-        if proj["time"].size and proj["time"][0] != 0:
-            times.append(proj["time"][0].values)
-            status = True
-    except KeyError:
-        # the datasets don't have times
-        pass
-    except IndexError:
-        # time is a scalar
-        if proj["time"].values != 0:
-            times.append(proj["time"].values)
-            status = True
-    return status
-
-
-def _get_average_time(times):
-    # Is there a more gracious way to handle this ?
-    if np.max(times) - np.min(times) > TIME_COMPATIBILITY_TOLERANCE:
-        raise IncompatibleTimes("Times do not match.")
-    return (np.max(times) - np.min(times)) / 2 + np.min(times)
+    @staticmethod
+    def combine_times(projectables):
+        """Get a combined time coordinate between projectables.
+
+        Returns the arithmetic mean between the available time coordinates,
+        provided they have a common dimensionality and units.  If no projectables
+        have time coordinates, return None.  Time coordinates should be CF-encoded,
+        i.e. have a numeric dtype and a units attribute.
+        """
+        _verify_times(projectables)
+        timed_projectables = [proj for proj in projectables
+                              if hasattr(proj, "coords")
+                              and "time" in proj.coords]
+        if len(timed_projectables) > 0:
+            with xr.set_options(keep_attrs=True):
+                # when the time coordinates are different, can't use xarray to
+                # sum them without this leading to expensive computations.  At
+                # this point, we should make sure NO time coordinates are
+                # assigned to the time coordinate or things will get messed up
+                # later.
+                first = timed_projectables[0].coords["time"]
+                new_coords = {}
+                if "y" in first.coords:
+                    new_coords["y"] = first.coords["y"].copy()
+                if "x" in first.coords:
+                    new_coords["x"] = first.coords["x"].copy()
+                da_new_time = sum([x.coords["time"].data for x in timed_projectables]) / len(timed_projectables)
+                xr_new_time = xr.DataArray(
+                        da_new_time,
+                        dims=first.dims,
+                        attrs=first.attrs.copy(),
+                        coords=new_coords)
+            return xr_new_time.assign_coords(time=(first.dims, da_new_time))
+
+
+def _verify_times(projectables):
+    """Verify that the times can be combined.
+
+    Times can be combined if they have consistent units and dimensions.
+    """
+    times = [p.coords["time"] for p in projectables
+             if hasattr(p, "coords") and "time" in p.coords]
+    for time in times:
+        if "units" not in time.attrs:
+            raise ValueError("Time coordinate lacks units attribute")
+        if time.attrs["units"] != times[0].attrs["units"]:
+            raise IncompatibleTimes("Time coordinates have inconsistent units.  "
+                                    "Conversions not implemented.")
+        if time.dims != times[0].dims:
+            raise IncompatibleTimes(
+                "Time coordinates have inconsistent dimensionality.  Only "
+                "consistent dimensionality is implemented.")
 
 
 class RGBCompositor(GenericCompositor):
diff --git a/satpy/composites/fill.py b/satpy/composites/fill.py
index bd5b51577e..76fed83511 100644
--- a/satpy/composites/fill.py
+++ b/satpy/composites/fill.py
@@ -42,7 +42,8 @@ class FillingCompositor(GenericCompositor):
 
     def __call__(self, projectables, nonprojectables=None, **info):
         """Generate the composite."""
-        projectables = self.match_data_arrays(projectables)
+        projectables = self.match_data_arrays(projectables,
+                                              drop_coordinates=False)
         projectables[1] = projectables[1].fillna(projectables[0])
         projectables[2] = projectables[2].fillna(projectables[0])
         projectables[3] = projectables[3].fillna(projectables[0])
@@ -54,7 +55,8 @@ class Filler(GenericCompositor):
 
     def __call__(self, projectables, nonprojectables=None, **info):
         """Generate the composite."""
-        projectables = self.match_data_arrays(projectables)
+        projectables = self.match_data_arrays(projectables,
+                                              drop_coordinates=False)
         filled_projectable = projectables[0].fillna(projectables[1])
         return super(Filler, self).__call__([filled_projectable], **info)
 
@@ -114,7 +116,7 @@ def __call__(
             **attrs
     ) -> xr.DataArray:
         """Generate the composite."""
-        datasets = self.match_data_arrays(datasets)
+        datasets = self.match_data_arrays(datasets, drop_coordinates=False)
         # At least one composite is requested.
         foreground_data = datasets[0]
         weights = self._get_coszen_blending_weights(datasets)
@@ -250,7 +252,14 @@ def _merge_bands_with_weights(self, day_data, night_data, weights, attrs):
                 night_band = night_band * (1 - weights)
             band = day_band + night_band
             band.attrs = attrs
+            # in theory, combining the times like the data values would be
+            # better, but averaging is good enough for any use case I can think
+            # of
+            time = self.combine_times([day_band, night_band])
+            if time is not None:
+                band.coords["time"] = time
             data.append(band)
+
         return data
 
     def _is_weightable(self, band):
diff --git a/satpy/etc/readers/seviri_l1b_hrit.yaml b/satpy/etc/readers/seviri_l1b_hrit.yaml
index d83135092e..b0754c3713 100644
--- a/satpy/etc/readers/seviri_l1b_hrit.yaml
+++ b/satpy/etc/readers/seviri_l1b_hrit.yaml
@@ -11,6 +11,7 @@ reader:
     HRIT reader for EUMETSAT MSG (Meteosat 8 to 11) SEVIRI Level 1b files.
   status: Nominal
   supports_fsspec: true
+  supports_time_tracking: true
   sensors: [seviri]
   reader: !!python/name:satpy.readers.core.yaml_reader.GEOSegmentYAMLReader
 
diff --git a/satpy/modifiers/atmosphere.py b/satpy/modifiers/atmosphere.py
index c7144c27ca..e73919d888 100644
--- a/satpy/modifiers/atmosphere.py
+++ b/satpy/modifiers/atmosphere.py
@@ -188,7 +188,8 @@ class CO2Corrector(ModifierBase):
 
     def __call__(self, projectables, optional_datasets=None, **info):
         """Apply correction."""
-        ir_039, ir_108, ir_134 = projectables
+        (ir_039, ir_108, ir_134) = self.match_data_arrays(projectables,
+                                                          drop_coordinates=False)
         logger.info("Applying CO2 correction")
         dt_co2 = (ir_108 - ir_134) / 4.0
         rcorr = ir_108 ** 4 - (ir_108 - dt_co2) ** 4
diff --git a/satpy/modifiers/geometry.py b/satpy/modifiers/geometry.py
index f2ae4a50c0..abaff9c223 100644
--- a/satpy/modifiers/geometry.py
+++ b/satpy/modifiers/geometry.py
@@ -47,7 +47,8 @@ def __init__(self, max_sza=95.0, **kwargs):  # noqa: D417
 
     def __call__(self, projectables, **info):
         """Generate the composite."""
-        projectables = self.match_data_arrays(list(projectables) + list(info.get("optional_datasets", [])))
+        projectables = self.match_data_arrays(list(projectables) + list(info.get("optional_datasets", [])),
+                                              drop_coordinates=False)
         vis = projectables[0]
         if vis.attrs.get("sunz_corrected"):
             logger.debug("Sun zenith correction already applied")
diff --git a/satpy/readers/core/seviri.py b/satpy/readers/core/seviri.py
index aaaa766117..a371877cf9 100644
--- a/satpy/readers/core/seviri.py
+++ b/satpy/readers/core/seviri.py
@@ -458,6 +458,41 @@ def add_scanline_acq_time(dataset, acq_time):
         "long_name"] = "Mean scanline acquisition time"
 
 
+def add_pixel_acq_time(dataset, label="time"):
+    """Add estimated pixel acquisition time for the given dataset.
+
+    Pixel acquisition time is added as a coordinate with the indicated name.
+    The coordinate will have dtype float32 and units of seconds since the
+    nominal start time.
+
+    For simplicity, takes scanline time as measurement time for an entire scanline.
+
+    Args:
+        dataset (xarray.DataArray): dataset with acq_time coordinate
+        label (str): Label for time coordinate.  Defaults to "time".  Choosing
+            a different label might impact the ability of writers trying to
+            process time coordinate information.
+    """
+    st_time = dataset.attrs["nominal_start_time"]
+    delta_per_line = (dataset["acq_time"].chunk({"y": dataset.sizes["y"]}) - np.datetime64(st_time)).astype("m8[s]")
+    # convert to float, because datetime64 cannot be resampled and ints cannot
+    # cover missing data
+    #
+    # NaT does not cast to nan but to -9.223372036854776e+18, hence detour via
+    # where
+    delta_per_line = xr.where(
+            delta_per_line.notnull(),
+            delta_per_line.astype(np.float32),
+            np.nan)
+    # this constant assumption could be improved
+    delta_per_pixel = delta_per_line.expand_dims({"x": dataset.sizes["x"]})
+    dataset.coords[label] = (
+            ("y", "x"),
+            delta_per_pixel.rename("time").transpose("y", "x").data)
+    dataset.coords[label].attrs["long_name"] = "Approximate pixel acquisition time"
+    dataset.coords[label].attrs["units"] = f"Seconds since {st_time:%Y-%m-%d %H:%M:%SZ}"
+
+
 def dec10216(inbuf):
     """Decode 10 bits data into 16 bits words.
 
diff --git a/satpy/readers/seviri_l1b_hrit.py b/satpy/readers/seviri_l1b_hrit.py
index b3811d6f6c..229abce3be 100644
--- a/satpy/readers/seviri_l1b_hrit.py
+++ b/satpy/readers/seviri_l1b_hrit.py
@@ -71,7 +71,7 @@
 ``nominal_start_time`` and ``nominal_end_time`` are also available directly
 via ``start_time`` and ``end_time`` respectively.
 
-Here is an exmaple of the content of the start/end time and ``time_parameters`` attibutes
+Here is an example of the content of the start/end time and ``time_parameters`` attibutes
 
 .. code-block:: python
 
@@ -246,6 +246,7 @@
     OrbitPolynomialFinder,
     ScanParams,
     SEVIRICalibrationHandler,
+    add_pixel_acq_time,
     add_scanline_acq_time,
     create_coef_dict,
     get_cds_time,
@@ -336,7 +337,8 @@ class HRITMSGPrologueFileHandler(HRITMSGPrologueEpilogueBase):
 
     def __init__(self, filename, filename_info, filetype_info, calib_mode="nominal",
                  ext_calib_coefs=None, include_raw_metadata=False,
-                 mda_max_array_size=None, fill_hrv=None, mask_bad_quality_scan_lines=None):
+                 mda_max_array_size=None, fill_hrv=None,
+                 mask_bad_quality_scan_lines=None, track_time=False):
         """Initialize the reader."""
         super().__init__(filename, filename_info,
                          filetype_info,
@@ -407,7 +409,8 @@ class HRITMSGEpilogueFileHandler(HRITMSGPrologueEpilogueBase):
 
     def __init__(self, filename, filename_info, filetype_info, calib_mode="nominal",
                  ext_calib_coefs=None, include_raw_metadata=False,
-                 mda_max_array_size=None, fill_hrv=None, mask_bad_quality_scan_lines=None):
+                 mda_max_array_size=None, fill_hrv=None,
+                 mask_bad_quality_scan_lines=None, track_time=False):
         """Initialize the reader."""
         super(HRITMSGEpilogueFileHandler, self).__init__(filename, filename_info,
                                                          filetype_info,
@@ -453,6 +456,19 @@ class HRITMSGFileHandler(HRITFileHandler):
                             reader='seviri_l1b_hrit',
                             reader_kwargs={'fill_hrv': False})
 
+    **Time tracking**
+
+    The reader supports adding per-pixel time estimates as coordinates to
+    loaded variables::
+
+        scene = satpy.Scene(filenames,
+                            reader='seviri_l1b_hrit',
+                            reader_kwargs={'track_time': True})
+
+    Those coordinates are kept in resampling if passing the argument
+    ``resample_coords=True`` to :meth:`~satpy.scene.Scene.resample`.
+    See the :doc:`example on storing valid time </examples/mean_time>` for details.
+
     **Metadata**
 
     See :mod:`satpy.readers.core.seviri`.
@@ -463,7 +479,8 @@ def __init__(self, filename, filename_info, filetype_info,
                  prologue, epilogue, calib_mode="nominal",
                  ext_calib_coefs=None, include_raw_metadata=False,
                  mda_max_array_size=100, fill_hrv=True,
-                 mask_bad_quality_scan_lines=True):
+                 mask_bad_quality_scan_lines=True,
+                 track_time=False):
         """Initialize the reader."""
         super(HRITMSGFileHandler, self).__init__(filename, filename_info,
                                                  filetype_info,
@@ -483,6 +500,7 @@ def __init__(self, filename, filename_info, filetype_info,
         self.ext_calib_coefs = ext_calib_coefs or {}
         self.mask_bad_quality_scan_lines = mask_bad_quality_scan_lines
         self._get_header()
+        self.track_time = track_time
 
     def _get_header(self):
         """Read the header info, and fill the metadata dictionary."""
@@ -690,6 +708,8 @@ def get_dataset(self, key, info):
             res = self.pad_hrv_data(res)
         self._update_attrs(res, info)
         self._add_scanline_acq_time(res)
+        if self.track_time:
+            self._add_pixel_acq_time(res)
         return res
 
     def pad_hrv_data(self, res):
@@ -769,6 +789,10 @@ def _add_scanline_acq_time(self, dataset):
         acq_time = get_cds_time(days=tline["days"], msecs=tline["milliseconds"])
         add_scanline_acq_time(dataset, acq_time)
 
+    def _add_pixel_acq_time(self, dataset):
+        """Estimate pixel acquisition time to the given dataset."""
+        add_pixel_acq_time(dataset)
+
     def _update_attrs(self, res, info):
         """Update dataset attributes."""
         res.attrs["units"] = info["units"]
diff --git a/satpy/scene.py b/satpy/scene.py
index 8693cc77b4..0cb5b9e5d7 100644
--- a/satpy/scene.py
+++ b/satpy/scene.py
@@ -862,14 +862,13 @@ def _slice_data(self, source_area, slices, dataset):
         return dataset
 
     def _resampled_scene(self, new_scn, destination_area, reduce_data=True,
+                         resample_coords=False,
                          **resample_kwargs):
         """Resample `datasets` to the `destination` area.
 
-        If data reduction is enabled, some local caching is perfomed in order to
+        If data reduction is enabled, some local caching is performed in order to
         avoid recomputation of area intersections.
         """
-        from satpy.resample.base import resample_dataset
-
         new_datasets = {}
         datasets = list(new_scn._datasets.values())
         destination_area = self._get_finalized_destination_area(destination_area, new_scn)
@@ -878,33 +877,86 @@ def _resampled_scene(self, new_scn, destination_area, reduce_data=True,
         reductions = {}
         for dataset, parent_dataset in dataset_walker(datasets):
             ds_id = DataID.from_dataarray(dataset)
-            pres = None
-            if parent_dataset is not None:
-                pres = new_datasets[DataID.from_dataarray(parent_dataset)]
-            if ds_id in new_datasets:
-                replace_anc(new_datasets[ds_id], pres)
-                if ds_id in new_scn._datasets:
-                    new_scn._datasets[ds_id] = new_datasets[ds_id]
+            pres = self._get_new_datasets_from_parent(new_datasets, parent_dataset)
+            if self._replace_anc_for_new_datasets(new_datasets, ds_id, pres, new_scn):
                 continue
-            if dataset.attrs.get("area") is None:
-                if parent_dataset is None:
-                    new_scn._datasets[ds_id] = dataset
-                else:
-                    replace_anc(dataset, pres)
+            if self._update_area(dataset, parent_dataset, new_scn, ds_id, pres):
                 continue
             LOG.debug("Resampling %s", ds_id)
             source_area = dataset.attrs["area"]
-            dataset, source_area = self._reduce_data(dataset, source_area, destination_area,
-                                                     reduce_data, reductions, resample_kwargs)
-            self._prepare_resampler(source_area, destination_area, resamplers, resample_kwargs)
-            kwargs = resample_kwargs.copy()
-            kwargs["resampler"] = resamplers[source_area]
-            res = resample_dataset(dataset, destination_area, **kwargs)
+            res = self._reduce_and_resample(dataset, source_area,
+                                            destination_area, reduce_data,
+                                            resamplers, reductions, resample_kwargs)
             new_datasets[ds_id] = res
             if ds_id in new_scn._datasets:
                 new_scn._datasets[ds_id] = res
             if parent_dataset is not None:
                 replace_anc(res, pres)
+            if resample_coords:
+                self._resample_coords(dataset, res,
+                                      destination_area, reduce_data,
+                                      resamplers, reductions, resample_kwargs)
+
+    def _reduce_and_resample(self, dataset, source_area, destination_area, reduce_data,
+                             resamplers, reductions, resample_kwargs):
+        from satpy.resample.base import resample_dataset
+
+        dataset, source_area = self._reduce_data(dataset, source_area, destination_area,
+                                                 reduce_data, reductions, resample_kwargs)
+        self._prepare_resampler(source_area, destination_area, resamplers, resample_kwargs)
+        kwargs = resample_kwargs.copy()
+        kwargs["resampler"] = resamplers[source_area]
+        return resample_dataset(dataset, destination_area, **kwargs)
+
+    def _resample_coords(self, orig_dataset, dataset,
+                         destination_area, reduce_data, resamplers, reductions,
+                         resample_kwargs):
+        """Resample matching coordinates on dataset."""
+        for cv in orig_dataset.coords:
+            if orig_dataset.coords[cv].dims == orig_dataset.dims:
+                LOG.debug(f"resampling coordinate {cv:s}")
+                orig_dataset.coords[cv].attrs["area"] = orig_dataset.attrs["area"]
+                res = self._reduce_and_resample(orig_dataset.coords[cv],
+                                                orig_dataset.attrs["area"],
+                                                destination_area, reduce_data,
+                                                resamplers, reductions, resample_kwargs)
+                dataset.coords[cv] = res
+
+
+    @classmethod
+    def _get_new_datasets_from_parent(self, new_datasets, parent_dataset):
+        if parent_dataset is not None:
+            return new_datasets[DataID.from_dataarray(parent_dataset)]
+        return None
+
+    @classmethod
+    def _replace_anc_for_new_datasets(self, new_datasets, ds_id, pres, new_scn):
+        if ds_id in new_datasets:
+            replace_anc(new_datasets[ds_id], pres)
+            if ds_id in new_scn._datasets:
+                new_scn._datasets[ds_id] = new_datasets[ds_id]
+            return True
+        return False
+
+    @classmethod
+    def _update_area(self, dataset, parent_dataset, new_scn, ds_id, pres):
+        if dataset.attrs.get("area") is None:
+            if parent_dataset is None:
+                new_scn._datasets[ds_id] = dataset
+            else:
+                replace_anc(dataset, pres)
+            return True
+        return False
+
+    def _resample_dataset(self, source_area, destination_area, dataset, resamplers, resample_kwargs):
+        from satpy.resample.base import resample_dataset
+
+        self._prepare_resampler(source_area, destination_area, resamplers, resample_kwargs)
+        kwargs = resample_kwargs.copy()
+        kwargs["resampler"] = resamplers[source_area]
+        res = resample_dataset(dataset, destination_area, **kwargs)
+
+        return res
 
     def _get_finalized_destination_area(self, destination_area, new_scn):
         if isinstance(destination_area, str):
@@ -927,32 +979,39 @@ def _prepare_resampler(self, source_area, destination_area, resamplers, resample
             resamplers[source_area] = resampler
             self._resamplers[key] = resampler
 
-    def _reduce_data(self, dataset, source_area, destination_area, reduce_data, reductions, resample_kwargs):
+    def _reduce_data(self, dataset, source_area, destination_area, reduce_data,
+                     reductions, resample_kwargs):
+        if not reduce_data:
+            LOG.debug("Data reduction disabled by the user")
+            return dataset, source_area
+
         try:
-            if reduce_data:
-                key = source_area
-                try:
-                    (slice_x, slice_y), source_area = reductions[key]
-                except KeyError:
-                    if resample_kwargs.get("resampler") == "gradient_search":
-                        factor = resample_kwargs.get("shape_divisible_by", 2)
-                    else:
-                        factor = None
-                    try:
-                        slice_x, slice_y = source_area.get_area_slices(
-                            destination_area, shape_divisible_by=factor)
-                    except TypeError:
-                        slice_x, slice_y = source_area.get_area_slices(
-                            destination_area)
-                    source_area = source_area[slice_y, slice_x]
-                    reductions[key] = (slice_x, slice_y), source_area
-                dataset = self._slice_data(source_area, (slice_x, slice_y), dataset)
-            else:
-                LOG.debug("Data reduction disabled by the user")
+            slice_x, slice_y = self._get_source_dest_slices(source_area, destination_area, reductions, resample_kwargs)
+            source_area = source_area[slice_y, slice_x]
+            reductions[source_area] = (slice_x, slice_y), source_area
+            dataset = self._slice_data(source_area, (slice_x, slice_y), dataset)
         except NotImplementedError:
             LOG.info("Not reducing data before resampling.")
+
         return dataset, source_area
 
+    @classmethod
+    def _get_source_dest_slices(self, source_area, destination_area, reductions, resample_kwargs):
+        try:
+            (slice_x, slice_y), source_area = reductions[source_area]
+        except KeyError:
+            if resample_kwargs.get("resampler") == "gradient_search":
+                factor = resample_kwargs.get("shape_divisible_by", 2)
+            else:
+                factor = None
+            try:
+                slice_x, slice_y = source_area.get_area_slices(
+                    destination_area, shape_divisible_by=factor)
+            except TypeError:
+                slice_x, slice_y = source_area.get_area_slices(
+                    destination_area)
+        return slice_x, slice_y
+
     def resample(
             self,
             destination: AreaDefinition | CoordinateDefinition | str | None = None,
@@ -961,6 +1020,7 @@ def resample(
             unload: bool = True,
             resampler: str | None = None,
             reduce_data: bool = True,
+            resample_coords: bool = False,
             **resample_kwargs,
     ) -> Scene:
         """Resample datasets and return a new scene.
@@ -989,6 +1049,11 @@ def resample(
                 information.
             reduce_data: Reduce data by matching the input and output
                 areas and slicing the data arrays (default: True)
+            resample_coords: If true, resample coordinates with (y, x)
+                dimensions.  If false (default), drop those coordinates.  Such
+                coordinates might be time coordinates if a scene was created
+                while passing ``track_time=True`` to the readers and those
+                readers support doing so.
             resample_kwargs: Remaining keyword arguments to pass to individual
                 resampler classes. See the individual resampler class
                 documentation :mod:`here <satpy.resample>` for available
@@ -999,7 +1064,8 @@ def resample(
             destination = self.finest_area(datasets)
         new_scn = self.copy(datasets=datasets)
         self._resampled_scene(new_scn, destination, resampler=resampler,
-                              reduce_data=reduce_data, **resample_kwargs)
+                              reduce_data=reduce_data,
+                              resample_coords=resample_coords, **resample_kwargs)
 
         # regenerate anything from the wishlist that needs it (combining
         # multiple resolutions, etc.)
diff --git a/satpy/tests/cf_tests/test_decoding.py b/satpy/tests/cf_tests/test_decoding.py
index 6309b62263..9b159b1ecf 100644
--- a/satpy/tests/cf_tests/test_decoding.py
+++ b/satpy/tests/cf_tests/test_decoding.py
@@ -20,9 +20,14 @@
 
 import datetime as dt
 
+import dask
+import dask.array as da
+import numpy as np
 import pytest
+import xarray as xr
 
 import satpy.cf.decoding
+from satpy.tests.utils import CustomScheduler
 
 
 class TestDecodeAttrs:
@@ -64,3 +69,17 @@ def test_decoding_doesnt_modify_original(self, attrs):
         """Test that decoding doesn't modify the original attributes."""
         satpy.cf.decoding.decode_attrs(attrs)
         assert isinstance(attrs["my_dict"], str)
+
+    def test_lazy_decode(self):
+        """Test that lazy decoding is lazy."""
+        xrda = xr.DataArray(
+                da.array([1, 2, 3]), dims=("y",),
+                attrs={"units": "seconds since 1950-05-01 01:00:00"})
+
+        with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+            res = satpy.cf.decoding.lazy_decode_cf_time(xrda)
+        resc = res.compute()
+        np.testing.assert_array_equal(
+            resc.values,
+            np.array(["1950-05-01T01:00:01", "1950-05-01T01:00:02",
+                      "1950-05-01T01:00:03"], dtype="M8[s]"))
diff --git a/satpy/tests/compositor_tests/test_arithmetic.py b/satpy/tests/compositor_tests/test_arithmetic.py
index f353cdfd37..e974cfa107 100644
--- a/satpy/tests/compositor_tests/test_arithmetic.py
+++ b/satpy/tests/compositor_tests/test_arithmetic.py
@@ -20,11 +20,14 @@
 import datetime as dt
 import unittest
 
+import dask
 import dask.array as da
 import numpy as np
 import pytest
 import xarray as xr
 
+from satpy.tests.utils import CustomScheduler
+
 
 class TestDifferenceCompositor(unittest.TestCase):
     """Test case for the difference compositor."""
@@ -100,6 +103,28 @@ def fake_dataset_pair(fake_area):
     return (ds1, ds2)
 
 
+@pytest.fixture
+def fake_dataset_pair_with_times(fake_area):
+    """Return a fake pair of 2×2 datasets with daskified time coordinates."""
+    ds1 = xr.DataArray(da.full((2, 2), 8, chunks=2, dtype=np.float32),
+                       dims=("y", "x"),
+                       coords={"time": xr.DataArray(
+                           da.array([[1, 2], [3, 4]]),
+                           dims=("y", "x"),
+                           attrs={"area": fake_area,
+                                  "units": (u:="seconds since 1999-08-07 06:45:00")})},
+                       attrs={"area": fake_area, "standard_name": "toa_bidirectional_reflectance"})
+    ds2 = xr.DataArray(da.full((2, 2), 4, chunks=2, dtype=np.float32),
+                       dims=("y", "x"),
+                       coords={"time": xr.DataArray(
+                           da.array([[1.1, 2], [3.1, 4]]),
+                           dims=("y", "x"),
+                           attrs={"area": fake_area,
+                                  "units": u})},
+                       attrs={"area": fake_area, "standard_name": "toa_bidirectional_reflectance"})
+    return (ds1, ds2)
+
+
 @pytest.mark.parametrize("kwargs", [{}, {"standard_name": "channel_ratio", "foo": "bar"}])
 def test_ratio_compositor(fake_dataset_pair, kwargs):
     """Test the ratio compositor."""
@@ -124,3 +149,30 @@ def test_sum_compositor(fake_dataset_pair):
     comp = SumCompositor(name="sum", standard_name="channel_sum")
     res = comp(fake_dataset_pair)
     np.testing.assert_allclose(res.values, 12)
+
+
+def test_difference_compositor_time_dask(fake_dataset_pair_with_times):
+    """Test difference compositor does not compute time coordinates."""
+    from satpy.composites.arithmetic import DifferenceCompositor
+    comp = DifferenceCompositor("diff")
+    with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+        res = comp(fake_dataset_pair_with_times)
+    assert "time" in res.coords
+
+
+def test_sum_compositor_time_dask(fake_dataset_pair_with_times):
+    """Test sum compositor does not compute time coordinates."""
+    from satpy.composites.arithmetic import SumCompositor
+    comp = SumCompositor("sum")
+    with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+        res = comp(fake_dataset_pair_with_times)
+    assert "time" in res.coords
+
+
+def test_ratio_compositor_time_dask(fake_dataset_pair_with_times):
+    """Test ratio compositor does not compute time coordinates."""
+    from satpy.composites.arithmetic import RatioCompositor
+    comp = RatioCompositor("ratio")
+    with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+        res = comp(fake_dataset_pair_with_times)
+    assert "time" in res.coords
diff --git a/satpy/tests/compositor_tests/test_core.py b/satpy/tests/compositor_tests/test_core.py
index 315431aa4c..caa024fed1 100644
--- a/satpy/tests/compositor_tests/test_core.py
+++ b/satpy/tests/compositor_tests/test_core.py
@@ -184,41 +184,149 @@ def test_inline_composites(self):
         assert comps["seviri"][fog_dep_ids[1]].attrs["prerequisites"] == ["IR_108", "IR_087"]
 
 
-class TestSingleBandCompositor(unittest.TestCase):
-    """Test the single-band compositor."""
-
-    def setUp(self):
-        """Create test data."""
-        from satpy.composites.core import SingleBandCompositor
-        self.comp = SingleBandCompositor(name="test")
-
-        all_valid = np.ones((2, 2))
-        self.all_valid = xr.DataArray(all_valid, dims=["y", "x"])
-
-    def test_call(self):
-        """Test calling the compositor."""
-        # Dataset with extra attributes
-        all_valid = self.all_valid
-        all_valid.attrs["sensor"] = "foo"
-        attrs = {
-            "foo": "bar",
-            "resolution": 333,
-            "units": "K",
-            "sensor": {"fake_sensor1", "fake_sensor2"},
-            "calibration": "BT",
-            "wavelength": 10.8
-        }
-        self.comp.attrs["resolution"] = None
-        res = self.comp([all_valid], **attrs)
-        # Verify attributes
-        assert res.attrs.get("sensor") == "foo"
-        assert "foo" in res.attrs
-        assert res.attrs.get("foo") == "bar"
-        assert "units" in res.attrs
-        assert "calibration" in res.attrs
-        assert "modifiers" not in res.attrs
-        assert res.attrs["wavelength"] == 10.8
-        assert res.attrs["resolution"] == 333
+@pytest.fixture
+def simple_dataset():
+    """A simple xarray DataArray with dimensions (y, x)."""
+    return xr.DataArray(
+            np.ones((2, 2)),
+            dims=["y", "x"])
+
+
+@pytest.fixture
+def dataset_with_1d_time_coordinate(simple_dataset):
+    """An xarray DataArray with time dimensions."""
+    return simple_dataset.assign_coords(
+            {"time": xr.DataArray(
+                np.array([1, 1], dtype=np.float32),
+                dims=("y",),
+                attrs={"units": "Seconds since 2100-03-04 05:30:00Z"})})
+
+@pytest.fixture
+def dataset_with_2d_time_coordinate(simple_dataset):
+    """An xarray DataArray with 2D time dimensions."""
+    return simple_dataset.assign_coords(
+            {"time": xr.DataArray(
+                np.array([[1, 1], [1, 1]], dtype=np.float32),
+                dims=("y", "x"),
+                attrs={"units": "Seconds since 2100-03-04 05:30:00Z"})})
+
+@pytest.fixture
+def datasets_rgb_with_identical_1d_time_coordinate(dataset_with_1d_time_coordinate):
+    """List of three xarray DataArray with same 1d time dimensions."""
+    return [dataset_with_1d_time_coordinate]*3
+
+@pytest.fixture
+def datasets_rgb_with_identical_2d_time_coordinate(dataset_with_2d_time_coordinate):
+    """List of three xarray DataArray with same 1d time dimensions."""
+    return [dataset_with_2d_time_coordinate]*3
+
+@pytest.fixture
+def datasets_rgb_with_non_identical_1d_time_coordinate(simple_dataset):
+    """List of three xarray DataArray with differing 1d time dimensions."""
+    r = xr.DataArray(
+        da.array([[0, 1], [0, 1]]),
+        dims=("y", "x"),
+        coords={
+            "time": xr.DataArray(
+                da.array([100, 200], dtype=np.float32),
+                dims=("y", ),
+                attrs={"units": "Seconds since 2100-03-04 05:30:00Z"})})
+    g = xr.DataArray(
+        da.array([[1, 0], [1, 0]]),
+        dims=("y", "x"),
+        coords={
+            "time": xr.DataArray(
+                da.array([100.1, 200.1], dtype=np.float32),
+                dims=("y", ),
+                attrs=r.coords["time"].attrs.copy())})
+    b = xr.DataArray(
+        da.array([[0, 1], [1, 1]]),
+        dims=("y", "x"),
+        coords={
+            "time": xr.DataArray(
+                da.array([100.2, 200.2], dtype=np.float32),
+                dims=("y", ),
+                attrs=r.coords["time"].attrs.copy())})
+    return [r, g, b]
+
+@pytest.fixture
+def datasets_rgb_with_non_identical_2d_time_coordinate(simple_dataset):
+    """List of three xarray DataArray with differing 2d time dimensions."""
+    r = simple_dataset.assign_coords(
+            {"time": xr.DataArray(
+                da.array([[100, 200], [101, 201]], dtype=np.float32),
+                dims=("y", "x"),
+                attrs={"units": "Seconds since 2100-03-04 05:30:00Z"})})
+    g = r.copy()
+    g[:] = da.array([[100.1, 200.1], [101.1, 201.1]])
+    b = r.copy()
+    b[:] = da.array([[100.2, 200.2], [101.2, 201.2]])
+    return [r, g, b]
+
+
+@pytest.fixture(params=["datasets_rgb_with_identical_1d_time_coordinate",
+                        "datasets_rgb_with_non_identical_1d_time_coordinate",
+                        "datasets_rgb_with_identical_2d_time_coordinate",
+                        "datasets_rgb_with_non_identical_2d_time_coordinate"])
+def datasets_rgb_with_time(request):
+    """Return one of several RGB datasets with time coordinates."""
+    return request.getfixturevalue(request.param)
+
+
+@pytest.fixture
+def single_band_compositor():
+    """A minimal SingleBandCompositor."""
+    from satpy.composites.core import SingleBandCompositor
+    return SingleBandCompositor(name="test")
+
+
+@pytest.fixture
+def generic_compositor():
+    """A minimal GenericCompositor."""
+    from satpy.composites.core import GenericCompositor
+    return GenericCompositor(name="test")
+
+
+def test_call_single_band_compositor(single_band_compositor, simple_dataset):
+    """Test calling the compositor."""
+    # Dataset with extra attributes
+    all_valid = simple_dataset
+    all_valid.attrs["sensor"] = "foo"
+    attrs = {
+        "foo": "bar",
+        "resolution": 333,
+        "units": "K",
+        "sensor": {"fake_sensor1", "fake_sensor2"},
+        "calibration": "BT",
+        "wavelength": 10.8
+    }
+    single_band_compositor.attrs["resolution"] = None
+    res = single_band_compositor([all_valid], **attrs)
+    # Verify attributes
+    assert res.attrs.get("sensor") == "foo"
+    assert "foo" in res.attrs
+    assert res.attrs.get("foo") == "bar"
+    assert "units" in res.attrs
+    assert "calibration" in res.attrs
+    assert "modifiers" not in res.attrs
+    assert res.attrs["wavelength"] == 10.8
+    assert res.attrs["resolution"] == 333
+
+
+def test_single_band_compositor_retain_time_coordinate(
+        single_band_compositor, dataset_with_1d_time_coordinate):
+    """Test that a SingleBandCompositor retains the time coordinate."""
+    res = single_band_compositor([dataset_with_1d_time_coordinate])
+    assert "time" in res.coords
+
+
+def test_generic_compositor_retain_time_coordinate(
+        generic_compositor, datasets_rgb_with_time):
+    """Test that a GenericCompositor retains the time coordinate."""
+    with assert_maximum_dask_computes(0):
+        res = generic_compositor(datasets_rgb_with_time)
+    assert "time" in res.coords
+    assert res.coords["time"].attrs["units"] == datasets_rgb_with_time[0].coords["time"].attrs["units"]
 
 
 class TestGenericCompositor(unittest.TestCase):
@@ -292,9 +400,8 @@ def test_get_sensors(self):
 
     @mock.patch("satpy.composites.core.GenericCompositor._get_sensors")
     @mock.patch("satpy.composites.core.combine_metadata")
-    @mock.patch("satpy.composites.core.check_times")
     @mock.patch("satpy.composites.core.GenericCompositor.match_data_arrays")
-    def test_call_with_mock(self, match_data_arrays, check_times, combine_metadata, get_sensors):
+    def test_call_with_mock(self, match_data_arrays, combine_metadata, get_sensors):
         """Test calling generic compositor."""
         from satpy.composites.core import IncompatibleAreas
         combine_metadata.return_value = dict()
@@ -343,7 +450,7 @@ def test_call(self):
         assert res.attrs["resolution"] == 333
 
     def test_deprecation_warning(self):
-        """Test deprecation warning for dcprecated composite recipes."""
+        """Test deprecation warning for deprecated composite recipes."""
         warning_message = "foo is a deprecated composite. Use composite bar instead."
         self.comp.attrs["deprecation_warning"] = warning_message
         with pytest.warns(UserWarning, match=warning_message):
diff --git a/satpy/tests/compositor_tests/test_fill.py b/satpy/tests/compositor_tests/test_fill.py
index 1cbc2c1dc5..153cc27555 100644
--- a/satpy/tests/compositor_tests/test_fill.py
+++ b/satpy/tests/compositor_tests/test_fill.py
@@ -31,6 +31,30 @@
 from satpy.tests.utils import CustomScheduler
 
 
+@pytest.fixture
+def data_with_time():
+    """Return fake data array with time coordinate."""
+    return xr.DataArray(
+            da.array([[1, 2, 3], [4, 3, 2]]),
+            dims=("y", "x"),
+            coords={"time": xr.DataArray(
+                da.array([[0, 1, 2], [3, 4, 5]]),
+                dims=("y", "x"),
+                attrs={"units": "seconds since 1900-01-01 00:00:00"})})
+
+
+@pytest.fixture
+def data_with_different_time():
+    """Return fake data array with a different time coordinate."""
+    return xr.DataArray(
+            da.array([[1, 2, 3], [4, 3, 2]]),
+            dims=("y", "x"),
+            coords={"time": xr.DataArray(
+                da.array([[0, 1, 2.1], [3, 4, 5]]),
+                dims=("y", "x"),
+                attrs={"units": "seconds since 1900-01-01 00:00:00"})})
+
+
 class TestDayNightCompositor(unittest.TestCase):
     """Test DayNightCompositor."""
 
@@ -222,8 +246,26 @@ def test_day_only_area_without_alpha(self):
         np.testing.assert_allclose(res.values[0], expected)
         assert "A" not in res.bands
 
+    def test_day_night_retains_time_coordinate(self):
+        """Test that time coordinate is retained."""
+        from satpy.composites.fill import DayNightCompositor
 
-class TestFillingCompositor(unittest.TestCase):
+        comp = DayNightCompositor(name="dn_test", day_night="day_night")
+        data_a = self.data_a.copy()
+        data_b = self.data_b.copy()
+        data_a.coords["time"] = xr.DataArray(
+                np.zeros_like(data_a.values),
+                dims=data_a.dims,
+                attrs={"area": data_a.attrs["area"],
+                       "units": "seconds since 2123-12-12 12:12:30"})
+        with xr.set_options(keep_attrs=True):
+            data_b.coords["time"] = data_a.coords["time"].copy()+0.1
+        res = comp((data_a, data_b, self.sza))
+        assert "time" in res.coords
+        np.testing.assert_array_equal(res.coords["bands"], ["R", "G", "B"])
+
+
+class TestFillingCompositor:
     """Test case for the filling compositor."""
 
     def test_fill(self):
@@ -239,6 +281,15 @@ def test_fill(self):
         np.testing.assert_allclose(res.sel(bands="G").data, filler.data)
         np.testing.assert_allclose(res.sel(bands="B").data, blue.data)
 
+    def test_fill_retains_time_coordinate(self, data_with_time,
+                                          data_with_different_time):
+        """Test that time coordinate is retained."""
+        from satpy.composites.fill import FillingCompositor
+        with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+            comp = FillingCompositor(name="fill_test")
+            res = comp([data_with_time]*2 + [data_with_different_time]*2)
+        assert "time" in res.coords
+
 
 class TestMultiFiller(unittest.TestCase):
     """Test case for the MultiFiller compositor."""
@@ -260,6 +311,16 @@ def test_fill(self):
         assert res.attrs["units"] == "K"
 
 
+def test_filler_retains_time_coordinate(data_with_time,
+                                        data_with_different_time):
+    """Test that Filler retains the time coordinate."""
+    from satpy.composites.fill import Filler
+    with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+        comp = Filler(name="filler_test")
+        res = comp([data_with_time, data_with_different_time])
+    assert "time" in res.coords
+
+
 def _enhance2dataset(dataset, convert_p=False):
     """Mock the enhance2dataset to return the original data."""
     return dataset
diff --git a/satpy/tests/compositor_tests/test_moved_compositors.py b/satpy/tests/compositor_tests/test_moved_compositors.py
index f74d13727e..a080f530fb 100644
--- a/satpy/tests/compositor_tests/test_moved_compositors.py
+++ b/satpy/tests/compositor_tests/test_moved_compositors.py
@@ -27,7 +27,6 @@
         "add_bands",
         "BackgroundCompositor",
         "CategoricalDataCompositor",
-        "check_times",
         "CloudCompositor",
         "ColorizeCompositor",
         "ColormapCompositor",
diff --git a/satpy/tests/modifier_tests/test_atmosphere.py b/satpy/tests/modifier_tests/test_atmosphere.py
new file mode 100644
index 0000000000..8e27556a6a
--- /dev/null
+++ b/satpy/tests/modifier_tests/test_atmosphere.py
@@ -0,0 +1,320 @@
+# Copyright (c) 2026 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 modifiers in modifiers/atmosphere.py."""
+import datetime as dt
+import unittest
+from unittest import mock
+
+import dask
+import dask.array as da
+import numpy as np
+import pytest
+import xarray as xr
+from pyresample.geometry import AreaDefinition
+
+from satpy.tests.utils import CustomScheduler
+
+
+class TestPSPRayleighReflectance:
+    """Test the pyspectral-based Rayleigh correction modifier."""
+
+    def _make_data_area(self):
+        """Create test area definition and data."""
+        rows = 3
+        cols = 5
+        area = AreaDefinition(
+            "some_area_name", "On-the-fly area", "geosabii",
+            {"a": "6378137.0", "b": "6356752.31414", "h": "35786023.0", "lon_0": "-89.5", "proj": "geos", "sweep": "x",
+             "units": "m"},
+            cols, rows,
+            (-5434894.954752679, -5434894.964451744, 5434894.964451744, 5434894.954752679))
+
+        data = np.zeros((rows, cols)) + 25
+        data[1, :] += 25
+        data[2, :] += 50
+        data = da.from_array(data, chunks=2)
+        return area, data
+
+    def _create_test_data(self, name, wavelength, resolution):
+        area, dnb = self._make_data_area()
+        input_band = xr.DataArray(dnb,
+                                  dims=("y", "x"),
+                                  attrs={
+                                      "platform_name": "Himawari-8",
+                                      "calibration": "reflectance", "units": "%", "wavelength": wavelength,
+                                      "name": name, "resolution": resolution, "sensor": "ahi",
+                                      "start_time": "2017-09-20 17:30:40.800000",
+                                      "end_time": "2017-09-20 17:41:17.500000",
+                                      "area": area, "ancillary_variables": [],
+                                      "orbital_parameters": {
+                                          "satellite_nominal_longitude": -89.5,
+                                          "satellite_nominal_latitude": 0.0,
+                                          "satellite_nominal_altitude": 35786023.4375,
+                                      },
+                                  })
+
+        red_band = xr.DataArray(dnb,
+                                dims=("y", "x"),
+                                attrs={
+                                    "platform_name": "Himawari-8",
+                                    "calibration": "reflectance", "units": "%", "wavelength": (0.62, 0.64, 0.66),
+                                    "name": "B03", "resolution": 500, "sensor": "ahi",
+                                    "start_time": "2017-09-20 17:30:40.800000",
+                                    "end_time": "2017-09-20 17:41:17.500000",
+                                    "area": area, "ancillary_variables": [],
+                                    "orbital_parameters": {
+                                        "satellite_nominal_longitude": -89.5,
+                                        "satellite_nominal_latitude": 0.0,
+                                        "satellite_nominal_altitude": 35786023.4375,
+                                    },
+                                })
+        fake_angle_data = da.ones_like(dnb, dtype=np.float32) * 90.0
+        angle1 = xr.DataArray(fake_angle_data,
+                              dims=("y", "x"),
+                              attrs={
+                                  "platform_name": "Himawari-8",
+                                  "calibration": "reflectance", "units": "%", "wavelength": wavelength,
+                                  "name": "satellite_azimuth_angle", "resolution": resolution, "sensor": "ahi",
+                                  "start_time": "2017-09-20 17:30:40.800000",
+                                  "end_time": "2017-09-20 17:41:17.500000",
+                                  "area": area, "ancillary_variables": [],
+                              })
+        return input_band, red_band, angle1, angle1, angle1, angle1
+
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
+    @pytest.mark.parametrize(
+        ("name", "wavelength", "resolution", "aerosol_type", "reduce_lim_low", "reduce_lim_high", "reduce_strength"),
+        [
+            ("B01", (0.45, 0.47, 0.49), 1000, "rayleigh_only", 70, 95, 1),
+            ("B02", (0.49, 0.51, 0.53), 1000, "rayleigh_only", 70, 95, 1),
+            ("B03", (0.62, 0.64, 0.66), 500, "rayleigh_only", 70, 95, 1),
+            ("B01", (0.45, 0.47, 0.49), 1000, "rayleigh_only", -95, -70, -1),
+        ]
+    )
+    def test_rayleigh_corrector(
+            self, tmp_path, name, wavelength, resolution, aerosol_type,
+            reduce_lim_low, reduce_lim_high, reduce_strength, dtype):
+        """Test PSPRayleighReflectance with fake data."""
+        from pyspectral.testing import mock_rayleigh
+
+        from satpy.modifiers.atmosphere import PSPRayleighReflectance
+
+        ray_cor = PSPRayleighReflectance(name=name, atmosphere="us-standard", aerosol_types=aerosol_type,
+                                         reduce_lim_low=reduce_lim_low, reduce_lim_high=reduce_lim_high,
+                                         reduce_strength=reduce_strength)
+        assert ray_cor.attrs["name"] == name
+        assert ray_cor.attrs["atmosphere"] == "us-standard"
+        assert ray_cor.attrs["aerosol_types"] == aerosol_type
+        assert ray_cor.attrs["reduce_lim_low"] == reduce_lim_low
+        assert ray_cor.attrs["reduce_lim_high"] == reduce_lim_high
+        assert ray_cor.attrs["reduce_strength"] == reduce_strength
+
+        input_band, red_band, *_ = self._create_test_data(name, wavelength, resolution)
+        with mock_rayleigh(rayleigh_dir=tmp_path):
+            res = ray_cor([input_band.astype(dtype), red_band.astype(dtype)])
+
+        assert isinstance(res, xr.DataArray)
+        assert isinstance(res.data, da.Array)
+        assert res.dtype == dtype
+        data = res.values
+        assert data.shape == (3, 5)
+        assert data.dtype == dtype
+
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
+    @pytest.mark.parametrize("as_optionals", [False, True])
+    def test_rayleigh_with_angles(self, tmp_path, as_optionals, dtype):
+        """Test PSPRayleighReflectance with angles provided."""
+        from pyspectral.testing import mock_rayleigh
+
+        from satpy.modifiers.atmosphere import PSPRayleighReflectance
+
+        aerosol_type = "rayleigh_only"
+        ray_cor = PSPRayleighReflectance(name="B01", atmosphere="us-standard", aerosol_types=aerosol_type)
+        prereqs, opt_prereqs = self._get_angles_prereqs_and_opts(as_optionals, dtype)
+        with mock.patch("satpy.modifiers.atmosphere.get_angles") as get_angles, mock_rayleigh(rayleigh_dir=tmp_path):
+            res = ray_cor(prereqs, opt_prereqs)
+        get_angles.assert_not_called()
+
+        assert isinstance(res, xr.DataArray)
+        assert isinstance(res.data, da.Array)
+        assert res.dtype == dtype
+        data = res.values
+        assert data.shape == (3, 5)
+        assert data.dtype == dtype
+
+    def _get_angles_prereqs_and_opts(self, as_optionals, dtype):
+        wavelength = (0.45, 0.47, 0.49)
+        resolution = 1000
+        input_band, red_band, *angles = self._create_test_data("B01", wavelength, resolution)
+        prereqs = [input_band.astype(dtype), red_band.astype(dtype)]
+        opt_prereqs = []
+        angles = [a.astype(dtype) for a in angles]
+        if as_optionals:
+            opt_prereqs = angles
+        else:
+            prereqs += angles
+        return prereqs, opt_prereqs
+
+
+class TestPSPAtmosphericalCorrection(unittest.TestCase):
+    """Test the pyspectral-based atmospheric correction modifier."""
+
+    def test_call(self):
+        """Test atmospherical correction."""
+        from pyresample.geometry import SwathDefinition
+
+        from satpy.modifiers.atmosphere import PSPAtmosphericalCorrection
+
+        # Patch methods
+        lons = np.zeros((5, 5))
+        lons[1, 1] = np.inf
+        lons = da.from_array(lons, chunks=5)
+        lats = np.zeros((5, 5))
+        lats[1, 1] = np.inf
+        lats = da.from_array(lats, chunks=5)
+        area = SwathDefinition(lons, lats)
+        stime = dt.datetime(2020, 1, 1, 12, 0, 0)
+        orb_params = {
+            "satellite_actual_altitude": 12345678,
+            "nadir_longitude": 0.0,
+            "nadir_latitude": 0.0,
+        }
+        band = xr.DataArray(da.zeros((5, 5)),
+                            attrs={"area": area,
+                                   "start_time": stime,
+                                   "name": "name",
+                                   "platform_name": "platform",
+                                   "sensor": "sensor",
+                                   "orbital_parameters": orb_params},
+                            dims=("y", "x"))
+
+        # Perform atmospherical correction
+        psp = PSPAtmosphericalCorrection(name="dummy")
+        res = psp(projectables=[band])
+        res.compute()
+
+
+@pytest.fixture
+def fake_ir_039():
+    """Return fake IR 3.9 µm data."""
+    return xr.DataArray(
+            np.array([
+                [258.85126, 265.5641 , 279.3434 , 288.12598],
+                [280.36703, 286.23367, 289.83792, 279.60886]]),
+            dims=("y", "x"))
+
+
+@pytest.fixture
+def fake_ir_039_dask_with_time():
+    """Return fake IR 3.9 µm data, daskified with time coordinate."""
+    return xr.DataArray(
+            da.array([
+                [258.85126, 265.5641 , 279.3434 , 288.12598],
+                [280.36703, 286.23367, 289.83792, 279.60886]]),
+            dims=("y", "x"),
+            coords={"time":
+                    xr.DataArray(
+                        da.array([
+                            [0, 1, 2, 3],
+                            [0, 1, 2, 3]]),
+                        dims=("y", "x"),
+                        attrs={"units": "seconds since 1950-05-05 05:00:00"})})
+
+
+
+@pytest.fixture
+def fake_ir_108():
+    """Return fake IR 10.8 µm data."""
+    return xr.DataArray(
+            np.array([
+                [251.42664, 244.1881 , 281.82275, 276.40485],
+                [280.53705, 286.7915 , 290.83807, 265.59827]]),
+            dims=("y", "x"))
+
+
+@pytest.fixture
+def fake_ir_108_dask_with_time():
+    """Return fake IR 10.8 µm data, daskified with time coordinate."""
+    return xr.DataArray(
+            da.array([
+                [251.42664, 244.1881 , 281.82275, 276.40485],
+                [280.53705, 286.7915 , 290.83807, 265.59827]]),
+            dims=("y", "x"),
+            coords={"time":
+                    xr.DataArray(
+                        da.array([
+                            [0, 1, 2, 3.1],
+                            [0, 1, 2.1, 3]]),
+                        dims=("y", "x"),
+                        attrs={"units": "seconds since 1950-05-05 05:00:00"})})
+
+
+@pytest.fixture
+def fake_ir_134():
+    """Return fake IR 13.4 µm data."""
+    return xr.DataArray(
+            np.array([
+                [238.38925, 231.85939, 254.27702, 251.15706],
+                [256.5422 , 259.16937, 260.9417 , 246.25409]]),
+            dims=("y", "x"))
+
+
+@pytest.fixture
+def fake_ir_134_dask_with_time():
+    """Return fake IR 13.4 µm data, daskified with time coordinate."""
+    return xr.DataArray(
+            da.array([
+                [238.38925, 231.85939, 254.27702, 251.15706],
+                [256.5422 , 259.16937, 260.9417 , 246.25409]]),
+            dims=("y", "x"),
+            coords={"time":
+                    xr.DataArray(
+                        da.array([
+                            [0, 1, 2.1, 3],
+                            [0, 1, 2.1, 3]]),
+                        dims=("y", "x"),
+                        attrs={"units": "seconds since 1950-05-05 05:00:00"})})
+
+
+class TestCO2Corrector:
+    """Test the 3.9 µm CO2 corrector."""
+
+    def test_co2_corrector(self, fake_ir_039, fake_ir_108, fake_ir_134):
+        """Test the CO2 corrector."""
+        from satpy.modifiers.atmosphere import CO2Corrector
+
+        corrector = CO2Corrector(name="co2_corrector")
+        res = corrector([fake_ir_039, fake_ir_108, fake_ir_134])
+        np.testing.assert_array_almost_equal(
+                res,
+                np.array([
+                    [261.732083, 267.884717, 285.923655, 293.365848],
+                    [286.013747, 292.709672, 296.845263, 283.557481]]))
+
+    def test_co2_corrector_time_dask(self, fake_ir_039_dask_with_time,
+                                     fake_ir_108_dask_with_time,
+                                     fake_ir_134_dask_with_time):
+        """Test CO2 corrector does not compute with time coordinates."""
+        from satpy.modifiers.atmosphere import CO2Corrector
+
+        corrector = CO2Corrector(name="co2_corrector")
+
+        with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+            res = corrector([fake_ir_039_dask_with_time,
+                             fake_ir_108_dask_with_time,
+                             fake_ir_134_dask_with_time])
+        assert "time" in res.coords
diff --git a/satpy/tests/reader_tests/test_seviri_l1b_hrit.py b/satpy/tests/reader_tests/test_seviri_l1b_hrit.py
index 15167af040..84d746fd24 100644
--- a/satpy/tests/reader_tests/test_seviri_l1b_hrit.py
+++ b/satpy/tests/reader_tests/test_seviri_l1b_hrit.py
@@ -681,3 +681,26 @@ def test_read_real_segment_zipped_with_upath(compressed_seviri_hrit_files):
     res = filehandler.get_dataset(dict(name="VIS008", calibration="counts"),
                                   dict(units="", wavelength=0.8, standard_name="counts"))
     res.compute()
+
+
+def test_track_time(prologue_file, segment_file, epilogue_file):
+    """Check tracking time with a virtual aux dataset."""
+    info = dict(start_time=dt.datetime(2222, 2, 22, 22, 0), service="")
+    prologue_fh = HRITMSGPrologueFileHandler(prologue_file, info, dict())
+    epilogue_fh = HRITMSGEpilogueFileHandler(epilogue_file, info, dict())
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", category=UserWarning,
+            message="No orbit polynomial valid for")
+        filehandler = HRITMSGFileHandler(segment_file, info, dict(),
+                                     prologue_fh, epilogue_fh,
+                                     track_time=True)
+    fake_acq_time = (np.datetime64("2022-02-22T22:00:00") +
+                     np.linspace(0, 900, 464).astype("m8[s]"))
+    fake_acq_time[:2] = np.datetime64("NaT")
+    fake_acq_time[-2:] = np.datetime64("NaT")
+    with mock.patch("satpy.readers.seviri_l1b_hrit.get_cds_time") as srsg:
+        srsg.return_value = fake_acq_time
+        res = filehandler.get_dataset(dict(name="VIS008", calibration="counts"),
+                                      dict(units="", wavelength=0.8, standard_name="counts"))
+    assert "time" in res.coords
+    assert res.dims == res.coords["time"].dims
diff --git a/satpy/tests/test_modifiers.py b/satpy/tests/test_modifiers.py
index ffea32e29e..0e955f1fbb 100644
--- a/satpy/tests/test_modifiers.py
+++ b/satpy/tests/test_modifiers.py
@@ -175,7 +175,7 @@ def test_basic_lims_provided(self, data_arr, sunz_sza, dtype):
         assert res.dtype == dtype
         assert values.dtype == dtype
 
-    def test_imcompatible_areas(self, sunz_ds2, sunz_sza):
+    def test_incompatible_areas(self, sunz_ds2, sunz_sza):
         """Test sunz correction on incompatible areas."""
         from satpy.composites.core import IncompatibleAreas
         from satpy.modifiers.geometry import SunZenithCorrector
@@ -183,6 +183,18 @@ def test_imcompatible_areas(self, sunz_ds2, sunz_sza):
         with pytest.raises(IncompatibleAreas):
             comp((sunz_ds2, sunz_sza), test_attr="test")
 
+    def test_retain_time_coordinate(self, sunz_ds1):
+        """Test sunz correction retains time coordinate."""
+        from satpy.modifiers.geometry import SunZenithCorrector
+        sunz_ds1.coords["time"] = xr.DataArray(
+                np.zeros_like(sunz_ds1.values),
+                dims=sunz_ds1.dims,
+                attrs={"area": sunz_ds1.attrs["area"],
+                       "units": "seconds since 2022-02-22 22:22:00"})
+        comp = SunZenithCorrector(name="sza_test", modifiers=tuple())
+        res = comp([sunz_ds1])
+        assert "time" in res.coords
+        np.testing.assert_array_equal(sunz_ds1.coords["time"], res.coords["time"])
 
 class TestSunZenithReducer:
     """Test case for the sun zenith reducer."""
@@ -423,182 +435,3 @@ def test_compositor(self, calculator, apply_modifier_info, sza):
         assert res.attrs["name"] == chan_name
         calculator.assert_called_with("NOAA-20", "viirs", "M12", sunz_threshold=86.0,
                                       masking_limit=NIRReflectance.MASKING_LIMIT)
-
-
-class TestPSPRayleighReflectance:
-    """Test the pyspectral-based Rayleigh correction modifier."""
-
-    def _make_data_area(self):
-        """Create test area definition and data."""
-        rows = 3
-        cols = 5
-        area = AreaDefinition(
-            "some_area_name", "On-the-fly area", "geosabii",
-            {"a": "6378137.0", "b": "6356752.31414", "h": "35786023.0", "lon_0": "-89.5", "proj": "geos", "sweep": "x",
-             "units": "m"},
-            cols, rows,
-            (-5434894.954752679, -5434894.964451744, 5434894.964451744, 5434894.954752679))
-
-        data = np.zeros((rows, cols)) + 25
-        data[1, :] += 25
-        data[2, :] += 50
-        data = da.from_array(data, chunks=2)
-        return area, data
-
-    def _create_test_data(self, name, wavelength, resolution):
-        area, dnb = self._make_data_area()
-        input_band = xr.DataArray(dnb,
-                                  dims=("y", "x"),
-                                  attrs={
-                                      "platform_name": "Himawari-8",
-                                      "calibration": "reflectance", "units": "%", "wavelength": wavelength,
-                                      "name": name, "resolution": resolution, "sensor": "ahi",
-                                      "start_time": "2017-09-20 17:30:40.800000",
-                                      "end_time": "2017-09-20 17:41:17.500000",
-                                      "area": area, "ancillary_variables": [],
-                                      "orbital_parameters": {
-                                          "satellite_nominal_longitude": -89.5,
-                                          "satellite_nominal_latitude": 0.0,
-                                          "satellite_nominal_altitude": 35786023.4375,
-                                      },
-                                  })
-
-        red_band = xr.DataArray(dnb,
-                                dims=("y", "x"),
-                                attrs={
-                                    "platform_name": "Himawari-8",
-                                    "calibration": "reflectance", "units": "%", "wavelength": (0.62, 0.64, 0.66),
-                                    "name": "B03", "resolution": 500, "sensor": "ahi",
-                                    "start_time": "2017-09-20 17:30:40.800000",
-                                    "end_time": "2017-09-20 17:41:17.500000",
-                                    "area": area, "ancillary_variables": [],
-                                    "orbital_parameters": {
-                                        "satellite_nominal_longitude": -89.5,
-                                        "satellite_nominal_latitude": 0.0,
-                                        "satellite_nominal_altitude": 35786023.4375,
-                                    },
-                                })
-        fake_angle_data = da.ones_like(dnb, dtype=np.float32) * 90.0
-        angle1 = xr.DataArray(fake_angle_data,
-                              dims=("y", "x"),
-                              attrs={
-                                  "platform_name": "Himawari-8",
-                                  "calibration": "reflectance", "units": "%", "wavelength": wavelength,
-                                  "name": "satellite_azimuth_angle", "resolution": resolution, "sensor": "ahi",
-                                  "start_time": "2017-09-20 17:30:40.800000",
-                                  "end_time": "2017-09-20 17:41:17.500000",
-                                  "area": area, "ancillary_variables": [],
-                              })
-        return input_band, red_band, angle1, angle1, angle1, angle1
-
-    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
-    @pytest.mark.parametrize(
-        ("name", "wavelength", "resolution", "aerosol_type", "reduce_lim_low", "reduce_lim_high", "reduce_strength"),
-        [
-            ("B01", (0.45, 0.47, 0.49), 1000, "rayleigh_only", 70, 95, 1),
-            ("B02", (0.49, 0.51, 0.53), 1000, "rayleigh_only", 70, 95, 1),
-            ("B03", (0.62, 0.64, 0.66), 500, "rayleigh_only", 70, 95, 1),
-            ("B01", (0.45, 0.47, 0.49), 1000, "rayleigh_only", -95, -70, -1),
-        ]
-    )
-    def test_rayleigh_corrector(
-            self, tmp_path, name, wavelength, resolution, aerosol_type,
-            reduce_lim_low, reduce_lim_high, reduce_strength, dtype):
-        """Test PSPRayleighReflectance with fake data."""
-        from pyspectral.testing import mock_rayleigh
-
-        from satpy.modifiers.atmosphere import PSPRayleighReflectance
-
-        ray_cor = PSPRayleighReflectance(name=name, atmosphere="us-standard", aerosol_types=aerosol_type,
-                                         reduce_lim_low=reduce_lim_low, reduce_lim_high=reduce_lim_high,
-                                         reduce_strength=reduce_strength)
-        assert ray_cor.attrs["name"] == name
-        assert ray_cor.attrs["atmosphere"] == "us-standard"
-        assert ray_cor.attrs["aerosol_types"] == aerosol_type
-        assert ray_cor.attrs["reduce_lim_low"] == reduce_lim_low
-        assert ray_cor.attrs["reduce_lim_high"] == reduce_lim_high
-        assert ray_cor.attrs["reduce_strength"] == reduce_strength
-
-        input_band, red_band, *_ = self._create_test_data(name, wavelength, resolution)
-        with mock_rayleigh(rayleigh_dir=tmp_path):
-            res = ray_cor([input_band.astype(dtype), red_band.astype(dtype)])
-
-        assert isinstance(res, xr.DataArray)
-        assert isinstance(res.data, da.Array)
-        assert res.dtype == dtype
-        data = res.values
-        assert data.shape == (3, 5)
-        assert data.dtype == dtype
-
-    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
-    @pytest.mark.parametrize("as_optionals", [False, True])
-    def test_rayleigh_with_angles(self, tmp_path, as_optionals, dtype):
-        """Test PSPRayleighReflectance with angles provided."""
-        from pyspectral.testing import mock_rayleigh
-
-        from satpy.modifiers.atmosphere import PSPRayleighReflectance
-
-        aerosol_type = "rayleigh_only"
-        ray_cor = PSPRayleighReflectance(name="B01", atmosphere="us-standard", aerosol_types=aerosol_type)
-        prereqs, opt_prereqs = self._get_angles_prereqs_and_opts(as_optionals, dtype)
-        with mock.patch("satpy.modifiers.atmosphere.get_angles") as get_angles, mock_rayleigh(rayleigh_dir=tmp_path):
-            res = ray_cor(prereqs, opt_prereqs)
-        get_angles.assert_not_called()
-
-        assert isinstance(res, xr.DataArray)
-        assert isinstance(res.data, da.Array)
-        assert res.dtype == dtype
-        data = res.values
-        assert data.shape == (3, 5)
-        assert data.dtype == dtype
-
-    def _get_angles_prereqs_and_opts(self, as_optionals, dtype):
-        wavelength = (0.45, 0.47, 0.49)
-        resolution = 1000
-        input_band, red_band, *angles = self._create_test_data("B01", wavelength, resolution)
-        prereqs = [input_band.astype(dtype), red_band.astype(dtype)]
-        opt_prereqs = []
-        angles = [a.astype(dtype) for a in angles]
-        if as_optionals:
-            opt_prereqs = angles
-        else:
-            prereqs += angles
-        return prereqs, opt_prereqs
-
-
-class TestPSPAtmosphericalCorrection(unittest.TestCase):
-    """Test the pyspectral-based atmospheric correction modifier."""
-
-    def test_call(self):
-        """Test atmospherical correction."""
-        from pyresample.geometry import SwathDefinition
-
-        from satpy.modifiers import PSPAtmosphericalCorrection
-
-        # Patch methods
-        lons = np.zeros((5, 5))
-        lons[1, 1] = np.inf
-        lons = da.from_array(lons, chunks=5)
-        lats = np.zeros((5, 5))
-        lats[1, 1] = np.inf
-        lats = da.from_array(lats, chunks=5)
-        area = SwathDefinition(lons, lats)
-        stime = dt.datetime(2020, 1, 1, 12, 0, 0)
-        orb_params = {
-            "satellite_actual_altitude": 12345678,
-            "nadir_longitude": 0.0,
-            "nadir_latitude": 0.0,
-        }
-        band = xr.DataArray(da.zeros((5, 5)),
-                            attrs={"area": area,
-                                   "start_time": stime,
-                                   "name": "name",
-                                   "platform_name": "platform",
-                                   "sensor": "sensor",
-                                   "orbital_parameters": orb_params},
-                            dims=("y", "x"))
-
-        # Perform atmospherical correction
-        psp = PSPAtmosphericalCorrection(name="dummy")
-        res = psp(projectables=[band])
-        res.compute()
diff --git a/satpy/tests/test_resample.py b/satpy/tests/test_resample.py
index 32784e9696..1bdcba43a7 100644
--- a/satpy/tests/test_resample.py
+++ b/satpy/tests/test_resample.py
@@ -716,3 +716,52 @@ def test_moved_import_warns(name):
     import satpy.resample
     with pytest.warns(UserWarning, match=".*has been moved.*"):
         _ = getattr(satpy.resample, name)
+
+
+@pytest.fixture
+def scene_with_time_coords():
+    """Return a scene with time coordinates."""
+    from pyresample import create_area_def
+
+    from satpy.tests.utils import make_fake_scene
+
+    ar1 = create_area_def("test", 4087, shape=(5, 5), resolution=1000, center=(0, 0))
+
+    sc = make_fake_scene(
+            {"ir": np.arange(25, dtype="f4").reshape(5, 5)},
+            area=ar1)
+    sc["ir"].coords["time"] = (
+            ("y", "x"),
+            np.linspace(0, 900, 25, dtype="uint16").reshape(5, 5))
+    sc["ir"].coords["time"].attrs["units"] = "seconds since 2222-02-22T22:22:22"
+    return sc
+
+@pytest.mark.parametrize("reduce_data", [True, False])
+def test_resample_time_coordinate(scene_with_time_coords, reduce_data):
+    """Test that resampling retains the time coordinate."""
+    from pyresample import create_area_def
+
+    ar2 = create_area_def("test", 4087, shape=(4, 4), resolution=1200, center=(100, 100))
+    ls = scene_with_time_coords.resample(ar2, resampler="nearest",
+                                         reduce_data=reduce_data)
+    assert "time" not in ls["ir"].coords  # drop by default
+    ls = scene_with_time_coords.resample(ar2, resampler="nearest",
+                                         resample_coords=False,
+                                         reduce_data=reduce_data)
+    assert "time" not in ls["ir"].coords
+    ls = scene_with_time_coords.resample(ar2, resampler="nearest",
+                                         resample_coords=True,
+                                         reduce_data=reduce_data)
+    assert "time" in ls["ir"].coords
+    assert ls["ir"].coords["time"].sizes == ls["ir"].sizes
+    assert ls["ir"].coords["time"].dtype == scene_with_time_coords["ir"].coords["time"].dtype
+    np.testing.assert_allclose(ls["ir"].coords["time"].mean(), 449.75)
+
+
+def test_slice_scene_time_coordinate(scene_with_time_coords):
+    """Test that slicing retains the time coordinate."""
+    # this test function may fit better elsewhere, but shares a fixture with
+    # test_resample_time_coordinate
+    sc2 = scene_with_time_coords[::2, ::2]
+    assert "time" in sc2["ir"].coords
+    assert sc2["ir"].coords["time"].sizes == sc2["ir"].sizes
diff --git a/satpy/tests/writer_tests/test_ninjogeotiff.py b/satpy/tests/writer_tests/test_ninjogeotiff.py
index a200aa73d9..e393ad8394 100644
--- a/satpy/tests/writer_tests/test_ninjogeotiff.py
+++ b/satpy/tests/writer_tests/test_ninjogeotiff.py
@@ -20,8 +20,10 @@
 import datetime
 import logging
 import os
+from math import prod
 from unittest.mock import Mock
 
+import dask
 import dask.array as da
 import numpy as np
 import pytest
@@ -31,16 +33,8 @@
 
 from satpy import Scene
 from satpy.enhancements.enhancer import get_enhanced_image
-
-try:
-    from math import prod
-except ImportError:  # Remove when dropping Python < 3.8
-    from functools import reduce
-    from operator import mul
-
-    def prod(iterable):  # type: ignore
-        """Drop-in replacement for math.prod."""
-        return reduce(mul, iterable, 1)
+from satpy.tests.utils import CustomScheduler
+from satpy.writers.core.compute import compute_writer_results
 
 # NOTE:
 # The following fixtures are not defined in this file, but are used and injected by Pytest:
@@ -338,6 +332,28 @@ def test_image_latlon(test_area_epsg4326):
     return get_enhanced_image(arr)
 
 
+@pytest.fixture(scope="module")
+def test_image_with_time_coords_no_dask(test_image_small_mid_atlantic_L):
+    """Get test image with time coordinates."""
+    time_coor = xr.DataArray(
+            np.ones_like(test_image_small_mid_atlantic_L.data.sel(bands="L")),
+            dims=("y", "x"),
+            attrs={"units": "seconds since 1985-08-13 13:00:00"})
+    test_image_small_mid_atlantic_L.data.coords["time"] = time_coor
+    return test_image_small_mid_atlantic_L
+
+
+@pytest.fixture(scope="module")
+def test_image_with_time_coords_dask(test_image_small_mid_atlantic_L):
+    """Get test image with time coordinates."""
+    time_coor = xr.DataArray(
+            da.ones_like(test_image_small_mid_atlantic_L.data.sel(bands="L")),
+            dims=("y", "x"),
+            attrs={"units": "seconds since 1985-08-13 13:00:00"})
+    test_image_small_mid_atlantic_L.data.coords["time"] = time_coor
+    return test_image_small_mid_atlantic_L
+
+
 @pytest.fixture(scope="module")
 def ntg1(test_image_small_mid_atlantic_L):
     """Create instance of NinJoTagGenerator class."""
@@ -514,8 +530,8 @@ def test_write_and_read_file(test_image_small_mid_atlantic_L, tmp_path):
         ChannelID=900015,
         DataType="GORN",
         DataSource="dowsing rod")
-    src = rasterio.open(fn)
-    tgs = src.tags()
+    with rasterio.open(fn) as src:
+        tgs = src.tags()
     assert tgs["ninjo_FileName"] == fn
     assert tgs["ninjo_DataSource"] == "dowsing rod"
     np.testing.assert_allclose(float(tgs["ninjo_Gradient"]),
@@ -541,8 +557,8 @@ def test_write_and_read_file_RGB(test_image_large_asia_RGB, tmp_path):
         ChannelID=900015,
         DataType="GORN",
         DataSource="dowsing rod")
-    src = rasterio.open(fn)
-    tgs = src.tags()
+    with rasterio.open(fn) as src:
+        tgs = src.tags()
     assert tgs["ninjo_FileName"] == fn
     assert tgs["ninjo_DataSource"] == "dowsing rod"
     assert "ninjo_Gradient" not in tgs.keys()
@@ -567,9 +583,9 @@ def test_write_and_read_file_LA(test_image_latlon, tmp_path):
         ChannelID=900015,
         DataType="GORN",
         DataSource="dowsing rod")
-    src = rasterio.open(fn)
-    assert len(src.indexes) == 2  # mode LA
-    tgs = src.tags()
+    with rasterio.open(fn) as src:
+        assert len(src.indexes) == 2  # mode LA
+        tgs = src.tags()
     assert tgs["ninjo_FileName"] == fn
     assert tgs["ninjo_DataSource"] == "dowsing rod"
     np.testing.assert_allclose(float(tgs["ninjo_Gradient"]), 0.31058823679007746)
@@ -598,10 +614,10 @@ def test_write_and_read_file_P(test_image_small_arctic_P, tmp_path):
         DataSource="dowsing rod",
         keep_palette=True,
         cmap=Colormap(*enumerate(zip(*([np.linspace(0, 1, 256)]*3)))))
-    src = rasterio.open(fn)
-    assert len(src.indexes) == 1  # mode P
-    assert src.colorinterp[0] == rasterio.enums.ColorInterp.palette
-    tgs = src.tags()
+    with rasterio.open(fn) as src:
+        assert len(src.indexes) == 1  # mode P
+        assert src.colorinterp[0] == rasterio.enums.ColorInterp.palette
+        tgs = src.tags()
     assert tgs["ninjo_FileName"] == fn
     assert tgs["ninjo_DataSource"] == "dowsing rod"
     assert tgs["ninjo_Gradient"] == "1.0"
@@ -636,8 +652,8 @@ def test_write_and_read_file_units(
     # all, but that currently fails due to
     # https://github.com/pytroll/satpy/issues/2022
     assert test_image_small_mid_atlantic_K_L.data.attrs["enhancement_history"][0] != {"scale": 1, "offset": 273.15}
-    src = rasterio.open(fn)
-    tgs = src.tags()
+    with rasterio.open(fn) as src:
+        tgs = src.tags()
     assert tgs["ninjo_FileName"] == fn
     assert tgs["ninjo_DataSource"] == "dowsing rod"
     np.testing.assert_allclose(float(tgs["ninjo_Gradient"]),
@@ -687,8 +703,8 @@ def test_write_and_read_no_quantity(
         ChannelID=900015,
         DataType="GORN",
         DataSource="dowsing rod")
-    src = rasterio.open(fn)
-    tgs = src.tags()
+    with rasterio.open(fn) as src:
+        tgs = src.tags()
     assert "ninjo_Gradient" not in tgs.keys()
     assert "ninjo_AxisIntercept" not in tgs.keys()
 
@@ -714,28 +730,24 @@ def test_write_and_read_via_scene(test_image_small_mid_atlantic_L, tmp_path):
         SatelliteNameID=6400014,
         ChannelID=900015,
         DataType="GORN")
-    src = rasterio.open(tmp_path / "test-montanha-do-pico.tif")
-    tgs = src.tags()
+    with rasterio.open(tmp_path / "test-montanha-do-pico.tif") as src:
+        tgs = src.tags()
     assert tgs["ninjo_FileName"] == os.fspath(tmp_path / "test-montanha-do-pico.tif")
 
 
 def test_get_all_tags(ntg1, ntg3, ntg_latlon, ntg_northpole, caplog):
     """Test getting all tags from dataset."""
-    # test that passed, dynamic, and mandatory tags are all included, and
-    # nothing more
+    # test that passed, dynamic, and mandatory tags are all included
     t1 = ntg1.get_all_tags()
-    assert set(t1.keys()) == (
+    delta = set(t1.keys()) - (
             ntg1.fixed_tags.keys() |
             ntg1.passed_tags |
             ntg1.dynamic_tags.keys() |
             {"DataSource"})
+    assert delta - ntg1.optional_tags == set()
     # test that when extra tag is passed this is also included
     t3 = ntg3.get_all_tags()
-    assert t3.keys() == (
-            ntg3.fixed_tags.keys() |
-            ntg3.passed_tags |
-            ntg3.dynamic_tags.keys() |
-            {"OverFlightTime"})
+    assert "OverFlightTime" in t3.keys()
     assert t3["OverFlightTime"] == 42
     # test that CentralMeridian skipped and warning logged
     with caplog.at_level(logging.DEBUG):
@@ -975,7 +987,7 @@ def test_create_unknown_tags(test_image_small_arctic_P):
 
 
 def test_str_ids(test_image_small_arctic_P):
-    """Test that channel and satellit IDs can be str."""
+    """Test that channel and satellite IDs can be str."""
     from satpy.writers.ninjogeotiff import NinJoTagGenerator
     NinJoTagGenerator(
         test_image_small_arctic_P,
@@ -986,3 +998,61 @@ def test_str_ids(test_image_small_arctic_P):
         PhysicUnit="N/A",
         PhysicValue="N/A",
         SatelliteNameID="trollsat")
+
+
+def test_write_mean_time(test_image_with_time_coords_no_dask, tmp_path):
+    """Test that mean time is written to GeoTIFF."""
+    import rasterio
+
+    from satpy.writers.ninjogeotiff import NinJoGeoTIFFWriter
+    fn = os.fspath(tmp_path / "test-{mean_time:%Y%m%d%H%M%S}.tif")
+    ngtw = NinJoGeoTIFFWriter(filename=fn)
+    ngtw.save_dataset(
+        test_image_with_time_coords_no_dask.data,
+        blockxsize=128,
+        blockysize=128,
+        compress="lzw",
+        predictor=2,
+        dynamic_fields={"mean_time"},
+        PhysicUnit="N/A",
+        PhysicValue="N/A",
+        SatelliteNameID="trollsat",
+        ChannelID="trollchannel",
+        DataType="GORN",
+        DataSource="dowsing rod")
+    with rasterio.open(fn.replace("{mean_time:%Y%m%d%H%M%S}",
+                                  "19850813130001")) as src:
+        tgs = src.tags()
+    assert tgs["ninjo_DateID"] == "492786000"
+    assert tgs["ninjo_ValidDateID"] == "492786001"
+
+
+def test_write_mean_time_dask(test_image_with_time_coords_dask, tmp_path):
+    """Test that mean time is written to GeoTIFF."""
+    import rasterio
+
+    from satpy.writers.ninjogeotiff import NinJoGeoTIFFWriter
+
+    # actual mean time in filename not supported
+    fn = os.fspath(tmp_path / "test-no-mean-time-in-filename.tif")
+    with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+        ngtw = NinJoGeoTIFFWriter(filename=fn)
+        res = ngtw.save_dataset(
+            test_image_with_time_coords_dask.data,
+            blockxsize=128,
+            blockysize=128,
+            compress="lzw",
+            predictor=2,
+            dynamic_fields={"mean_time"},
+            PhysicUnit="N/A",
+            PhysicValue="N/A",
+            SatelliteNameID="trollsat",
+            ChannelID="trollchannel",
+            DataType="GORN",
+            DataSource="dowsing rod",
+            compute=False)
+    compute_writer_results([res])
+    with rasterio.open(fn) as src:
+        tgs = src.tags()
+    assert tgs["ninjo_DateID"] == "492786000"
+    assert tgs["ninjo_ValidDateID"] == "492786001"
diff --git a/satpy/tests/writer_tests/test_utils.py b/satpy/tests/writer_tests/test_utils.py
new file mode 100644
index 0000000000..62e4d55ec4
--- /dev/null
+++ b/satpy/tests/writer_tests/test_utils.py
@@ -0,0 +1,65 @@
+# Copyright (c) 2025 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/>.
+
+"""Testing of writer helper functions."""
+
+import datetime
+
+import numpy as np
+import pytest
+import xarray as xr
+
+
+def test_get_mean_time():
+    """Test extracting valid time from a dataset."""
+    from satpy.writers.core.utils import get_mean_time
+
+    xrda = xr.DataArray(
+            np.zeros((3, 3)),
+            name="Čuonajóhka",
+            dims=("y", "x"),
+            coords={"time":
+                    xr.DataArray(
+                        np.full((3, 3), 2),
+                        dims=("y", "x"),
+                        attrs={"units": "seconds since 2222-02-02T22:22:20"})})
+
+    assert get_mean_time(xrda) == datetime.datetime(2222, 2, 2, 22, 22, 22)
+
+
+def test_get_mean_time_no_time_coordinate_dataset_name():
+    """Test raising of ValueError in the absence of a time coordinate.
+
+    Takes name from dataset.
+    """
+    from satpy.writers.core.utils import get_mean_time
+
+    xrda = xr.DataArray(np.zeros((3, 3)), dims=("y", "x"), name="Liedik")
+    with pytest.raises(ValueError,
+                       match="Dataset Liedik has no time coordinate."):
+        get_mean_time(xrda)
+
+
+def test_get_mean_time_no_time_coordinate_attribute_name():
+    """Test raising of ValueError in the absence of a time coordinate.
+
+    Takes name from dataset attribute (more common in satpy).
+    """
+    from satpy.writers.core.utils import get_mean_time
+    xrda = xr.DataArray(np.zeros((3, 3)), dims=("y", "x"), name="Liedik", attrs={"name": "Gámasčearru"})
+    with pytest.raises(ValueError,
+                       match="Dataset Gámasčearru has no time coordinate."):
+        get_mean_time(xrda)
diff --git a/satpy/writers/cf_writer.py b/satpy/writers/cf_writer.py
index d98df05b2c..7ebe085745 100644
--- a/satpy/writers/cf_writer.py
+++ b/satpy/writers/cf_writer.py
@@ -232,7 +232,8 @@ def save_dataset(self, dataset, filename=None, fill_value=None, **kwargs):
 
     def save_datasets(self, datasets, filename=None, groups=None, header_attrs=None, engine=None, epoch=None,  # noqa: D417
                       flatten_attrs=False, exclude_attrs=None, include_lonlats=True, pretty=False,
-                      include_orig_name=True, numeric_name_prefix="CHANNEL_", **to_netcdf_kwargs):
+                      include_orig_name=True, numeric_name_prefix="CHANNEL_",
+                      dynamic_fields=set(), **to_netcdf_kwargs):
         """Save the given datasets in one netCDF file.
 
         Note that all datasets (if grouping: in one group) must have the same projection coordinates.
@@ -268,7 +269,9 @@ def save_datasets(self, datasets, filename=None, groups=None, header_attrs=None,
 
         # Define netCDF filename if not provided
         # - It infers the name from the first DataArray
-        filename = filename or self.get_filename(**datasets[0].attrs)
+        filename = filename or self.get_filename(**datasets[0].attrs,
+                                                 **self._get_dynamic_fields(datasets[0],
+                                                                            dynamic_fields))
 
         # Collect xr.Dataset for each group
         grouped_datasets, header_attrs = collect_cf_datasets(list_dataarrays=datasets,  # list of xr.DataArray
diff --git a/satpy/writers/core/base.py b/satpy/writers/core/base.py
index e0a53d0f7e..de88ea74d4 100644
--- a/satpy/writers/core/base.py
+++ b/satpy/writers/core/base.py
@@ -25,6 +25,8 @@
 from satpy.plugin_base import Plugin
 from satpy.writers.core.compute import compute_writer_results, split_results
 
+from . import utils
+
 if typing.TYPE_CHECKING:
     from collections.abc import Iterable
     from os import PathLike
@@ -256,3 +258,19 @@ def save_dataset(
         """
         raise NotImplementedError(
             "Writer '%s' has not implemented dataset saving" % (self.name, ))
+
+
+    def _get_dynamic_fields(self, dataset, dynamic_fields):
+        """Report dynamic fields.
+
+        Apart from attributes stored in the dataset, we can retrieve attributes
+        derived from the data.  This may be useful in case we want to include
+        such information in the filename.
+        """
+        return {
+            field: getattr(self, f"get_{field}")(dataset)
+            for field in dynamic_fields}
+
+    def get_mean_time(self, dataset):
+        """Get mean time for dataset."""
+        return utils.get_mean_time(dataset)
diff --git a/satpy/writers/core/image.py b/satpy/writers/core/image.py
index b33014e7c5..1c99337e0a 100644
--- a/satpy/writers/core/image.py
+++ b/satpy/writers/core/image.py
@@ -66,13 +66,6 @@ def __init__(self, name=None, filename=None, base_dir=None, enhance=None, **kwar
 
             kwargs (dict): Additional keyword arguments to pass to the
                 :class:`~satpy.writers.core.base.Writer` base class.
-
-        .. versionchanged:: 0.10
-
-            Deprecated `enhancement_config_file` and 'enhancer' in favor of
-            `enhance`. Pass an instance of the `Enhancer` class to `enhance`
-            instead.
-
         """
         super().__init__(name, filename, base_dir, **kwargs)
         if enhance is False:
diff --git a/satpy/writers/core/utils.py b/satpy/writers/core/utils.py
new file mode 100644
index 0000000000..f9e54f8b5e
--- /dev/null
+++ b/satpy/writers/core/utils.py
@@ -0,0 +1,46 @@
+# Copyright (c) 2025 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/>.
+"""Utilities shared between writers."""
+
+import xarray as xr
+
+from satpy.cf.decoding import lazy_decode_cf_time
+
+
+def get_mean_time(dataset):
+    """Get the arithmetic mean time for a dataset.
+
+    If a dataset has time coordinates, get the arithmetic mean time.
+    This can be used as a representative time for the dataset.
+
+    Args:
+        dataset: xarray dataarray with time coordinates
+
+    Returns:
+        datetime.datetime object with arithmetic mean time for dataset
+    """
+    if "time" not in dataset.coords:
+        nm = dataset.attrs.get("name", dataset.name)
+        raise ValueError(
+           f"Dataset {nm:s} has no time coordinate. "
+            "No mean time can be calculated.  To track time, "
+            "pass `reader_kwargs = {'track_time': True}` to `Scene.__init__` "
+            "for a supported reader and, if applicable, `resample_coords=True` "
+            "to `Scene.resample`.")
+
+    with xr.set_options(keep_attrs=True):
+        tm = dataset.coords["time"].mean()
+    return lazy_decode_cf_time(tm)
diff --git a/satpy/writers/geotiff.py b/satpy/writers/geotiff.py
index 6b80bcadf9..3601f3ed95 100644
--- a/satpy/writers/geotiff.py
+++ b/satpy/writers/geotiff.py
@@ -162,6 +162,7 @@ def save_image(  # noqa: D417
             colormap_tag: str | None = None,
             driver: str | None = None,
             tiled: bool = True,
+            dynamic_fields: set[str] = set(),
             **kwargs
     ) -> list[da.Array | Delayed] | tuple[list[da.Array], list[Any]] | list[str | PathLike | None]:
         """Save the image to the given ``filename`` in geotiff_ format.
@@ -244,11 +245,15 @@ def save_image(  # noqa: D417
             include_scale_offset: Deprecated.
                 Use ``scale_offset_tags=("scale", "offset")`` to include scale
                 and offset tags.
+            dynamic_fields: set of strings of fields that are calculated
+                dynamically to enter the filename.
 
         .. _geotiff: http://trac.osgeo.org/geotiff/
 
         """
-        filename = filename or self.get_filename(**img.data.attrs)
+        filename = filename or self.get_filename(**img.data.attrs,
+                                                 **self._get_dynamic_fields(img.data,
+                                                                            dynamic_fields))
 
         gdal_options = self._get_gdal_options(kwargs)
         if fill_value is None:
diff --git a/satpy/writers/ninjogeotiff.py b/satpy/writers/ninjogeotiff.py
index 7eca65d38a..29df5ff3e6 100644
--- a/satpy/writers/ninjogeotiff.py
+++ b/satpy/writers/ninjogeotiff.py
@@ -89,6 +89,7 @@
 
 import numpy as np
 
+from .core.utils import get_mean_time
 from .geotiff import GeoTIFFWriter
 
 logger = logging.getLogger(__name__)
@@ -110,6 +111,7 @@ def save_image(  # noqa: D417
             compute=True, keep_palette=False, cmap=None, overviews=None,
             overviews_minsize=256, overviews_resampling=None,
             tags=None, config_files=None,
+            dynamic_fields=set(),
             *, ChannelID, DataType, PhysicUnit, PhysicValue,
             SatelliteNameID, **kwargs):
         """Save image along with NinJo tags.
@@ -138,6 +140,8 @@ def save_image(  # noqa: D417
                 As for :meth:`~satpy.writers.geotiff.GeoTIFFWriter.save_image`.
             overviews_resampling (str):
                 As for :meth:`~satpy.writers.geotiff.GeoTIFFWriter.save_image`.
+            dynamic_fields (set[str]):
+                As for :meth:`~satpy.writers.geotiff.GeoTIFFWriter.save_image`.
             tags (dict): Extra (not NinJo) tags to add to GDAL MetaData
             config_files (Any): Not directly used by this writer, supported
                 for compatibility with other writers.
@@ -167,12 +171,16 @@ def save_image(  # noqa: D417
             PhysicValue (str)
                 NinJo label for quantity (example: "temperature")
 
+        The tag ``ValidDateID`` is included if the dataset has a ``time``
+        coordinate.  See the :doc:`example on storing mean time </examples/mean_time>`
+        for details.  All other tags are filled automatically without user action.
         """
         dataset = image.data
 
         # filename not passed on to writer by Scene.save_dataset, but I need
         # it!
-        filename = filename or self.get_filename(**dataset.attrs)
+        filename = filename or self.get_filename(**dataset.attrs,
+                                                 **self._get_dynamic_fields(dataset, dynamic_fields))
 
         gdal_opts = {}
         ntg_opts = {}
@@ -263,13 +271,16 @@ class NinJoTagGenerator:
     Tags are gathered from three sources:
 
     - Fixed tags, contained in the attribute ``fixed_tags``.  The value of
-      those tags is hardcoded and never changes.
+      those tags is hardcoded and does not change unless there is a major
+      change in the fileformat definition.
     - Tags passed by the user, contained in the attribute ``passed_tags``.
       Those tags must be passed by the user as arguments to the writer, which
       will pass them on when instantiating this class.
     - Tags calculated from data and metadata.  Those tags are defined in the
       attribute ``dynamic_tags``.  They are either calculated from image data,
       from image metadata, or from arguments passed by the user to the writer.
+      Some dynamic tags, such as ``ValidDateID``, are only calculated if
+      necessary metadata are available in the input data.
 
     Some tags are mandatory (defined in ``mandatory_tags``).  All tags that are
     not mandatory are optional.  By default, optional tags are generated if and
@@ -293,6 +304,7 @@ class NinJoTagGenerator:
         "ColorDepth": "color_depth",
         "CreationDateID": "creation_date_id",
         "DateID": "date_id",
+        "ValidDateID": "valid_date_id",
         "EarthRadiusLarge": "earth_radius_large",
         "EarthRadiusSmall": "earth_radius_small",
         "FileName": "filename",
@@ -323,7 +335,7 @@ class NinJoTagGenerator:
                      "OverFlightTime", "IsBlackLinesCorrection",
                      "IsAtmosphereCorrected", "IsCalibrated", "IsNormalized",
                      "OriginalHeader", "IsValueTableAvailable",
-                     "ValueTableFloatField"}
+                     "ValueTableFloatField", "ValidDateID"}
 
     # tags that are added later in other ways
     postponed_tags = {"AxisIntercept", "Gradient"}
@@ -357,7 +369,7 @@ def get_all_tags(self):
         for tag in self.tag_names:
             try:
                 tags[tag] = self.get_tag(tag)
-            except (AttributeError, KeyError) as e:
+            except (AttributeError, KeyError, ValueError) as e:
                 if tag in self.mandatory_tags:
                     raise
                 logger.debug(
@@ -418,12 +430,21 @@ def get_date_id(self):
         """Calculate the date ID.
 
         That's seconds since UNIX Epoch for the time corresponding to the
-        satellite image start of measurement time.
+        satellite image reference time or start of measurement time.
         """
         tm = self.dataset.attrs["start_time"]
         delta = tm.replace(tzinfo=datetime.timezone.utc) - self._epoch
         return int(delta.total_seconds())
 
+    def get_valid_date_id(self):
+        """Calculate the valid date ID.
+
+        That's seconds since UNIX Epoch for a representative time for the
+        image.
+        """
+        mean_time = get_mean_time(self.dataset)
+        return (mean_time - np.datetime64(self._epoch)).astype("m8[s]").astype("int64")
+
     def get_earth_radius_large(self):
         """Return the Earth semi-major axis in metre."""
         return self.dataset.attrs["area"].crs.ellipsoid.semi_major_metre
diff --git a/satpy/writers/simple_image.py b/satpy/writers/simple_image.py
index b24294599b..a0afa4d9d0 100644
--- a/satpy/writers/simple_image.py
+++ b/satpy/writers/simple_image.py
@@ -49,6 +49,7 @@ def save_image(
         img: XRImage,
         filename: str | None = None,
         compute: bool = True,
+        dynamic_fields: set[str] = set(),
         **kwargs,
     ) -> list[da.Array | Delayed] | tuple[list[da.Array], list[Any]] | list[str | PathLike | None]:
         """Save Image object to a given ``filename``.
@@ -63,6 +64,8 @@ def save_image(
                 If `False` return either a `dask.delayed.Delayed`
                 object or tuple of (source, target). See the
                 return values below for more information.
+            dynamic_fields: Fields calculated dynamically from the
+                data for the purposes of filling in the filename.
             **kwargs: Keyword arguments to pass to the images `save` method.
 
         Returns:
@@ -77,7 +80,9 @@ def save_image(
             this method.
 
         """
-        filename = filename or self.get_filename(**img.data.attrs)
+        filename = filename or self.get_filename(**img.data.attrs,
+                                                 **self._get_dynamic_fields(img.data,
+                                                                            dynamic_fields))
 
         LOG.debug("Saving to image: %s", filename)
         res = img.save(filename, compute=compute, **kwargs)