New preprocessor: local_solar_time (#2258)

Co-authored-by: Bouwe Andela <[email protected]>

doc/recipe/preprocessor.rst

@@ -1219,6 +1219,7 @@ The ``_time.py`` module contains the following preprocessor functions:
 * regrid_time_: Aligns the time axis of each dataset to have common time
  points and calendars.
 * timeseries_filter_: Allows application of a filter to the time-series data.
+* local_solar_time_: Convert cube with UTC time to local solar time.
 
 Statistics functions are applied by default in the order they appear in the
 list. For example, the following example applied to hourly data will retrieve
@@ -1653,6 +1654,55 @@ Examples:
 
 See also :func:`esmvalcore.preprocessor.timeseries_filter`.
 
+.. _local_solar_time:
+
+``local_solar_time``
+--------------------
+
+Many variables in the Earth system show a strong diurnal cycle.
+The reason for that is of course Earth's rotation around its own axis, which
+leads to a diurnal cycle of the incoming solar radiation.
+While UTC time is a very good absolute time measure, it is not really suited to
+analyze diurnal cycles over larger regions.
+For example, diurnal cycles over Russia and the USA are phase-shifted by ~180°
+= 12 hr in UTC time.
+
+This is where the `local solar time (LST)
+<https://en.wikipedia.org/wiki/Solar_time>`__ comes into play:
+For a given location, 12:00 noon LST is defined as the moment when the sun
+reaches its highest point in the sky.
+By using this definition based on the origin of the diurnal cycle (the sun), we
+can directly compare diurnal cycles across the globe.
+LST is mainly determined by the longitude of a location, but due to the
+eccentricity of Earth's orbit, it also depends on the day of year (see
+`equation of time <https://en.wikipedia.org/wiki/Equation_of_time>`__).
+However, this correction is at most ~15 min, which is usually smaller than the
+highest frequency output of CMIP6 models (1 hr) and smaller than the time scale
+for diurnal evolution of meteorological phenomena (which is in the order of
+hours, not minutes).
+Thus, instead, we use the **mean** LST, which solely depends on longitude:
+
+.. math::
+
+ LST = UTC + 12 \cdot \frac{lon}{180°}
+
+where the times are given in hours and `lon` in degrees in the interval [-180,
+180].
+To transform data from UTC to LST, this preprocessor shifts data along the time
+axis based on the longitude.
+
+This preprocessor does not need any additional parameters.
+
+Example:
+
+.. code-block:: yaml
+
+ calculate_local_solar_time:
+ local_solar_time:
+
+See also :func:`esmvalcore.preprocessor.local_solar_time`.
+
+
 .. _area operations:
 
 Area manipulation

esmvalcore/iris_helpers.py

@@ -1,11 +1,14 @@
 """Auxiliary functions for :mod:`iris`."""
-from typing import Dict, List, Sequence
+from __future__ import annotations
+
+from typing import Dict, Iterable, List, Literal, Sequence
 
 import dask.array as da
 import iris
 import iris.cube
 import iris.util
 import numpy as np
+from iris.coords import Coord
 from iris.cube import Cube
 from iris.exceptions import CoordinateMultiDimError
 
@@ -157,3 +160,111 @@ def merge_cube_attributes(
  # Step 3: modify the cubes in-place
  for cube in cubes:
  cube.attributes = final_attributes
+
+
+def _rechunk(
+ array: da.core.Array,
+ complete_dims: list[int],
+ remaining_dims: int | Literal['auto'],
+) -> da.core.Array:
+ """Rechunk a given array so that it is not chunked along given dims."""
+ new_chunks: list[str | int] = [remaining_dims] * array.ndim
+ for dim in complete_dims:
+ new_chunks[dim] = -1
+ return array.rechunk(new_chunks)
+
+
+def _rechunk_dim_metadata(
+ cube: Cube,
+ complete_dims: Iterable[int],
+ remaining_dims: int | Literal['auto'] = 'auto',
+) -> None:
+ """Rechunk dimensional metadata of a cube (in-place)."""
+ # Non-dimensional coords that span complete_dims
+ # Note: dimensional coords are always realized (i.e., numpy arrays), so no
+ # chunking is necessary
+ for coord in cube.coords(dim_coords=False):
+ dims = cube.coord_dims(coord)
+ complete_dims_ = [dims.index(d) for d in complete_dims if d in dims]
+ if complete_dims_:
+ if coord.has_lazy_points():
+ coord.points = _rechunk(
+ coord.lazy_points(), complete_dims_, remaining_dims
+ )
+ if coord.has_bounds() and coord.has_lazy_bounds():
+ coord.bounds = _rechunk(
+ coord.lazy_bounds(), complete_dims_, remaining_dims
+ )
+
+ # Rechunk cell measures that span complete_dims
+ for measure in cube.cell_measures():
+ dims = cube.cell_measure_dims(measure)
+ complete_dims_ = [dims.index(d) for d in complete_dims if d in dims]
+ if complete_dims_ and measure.has_lazy_data():
+ measure.data = _rechunk(
+ measure.lazy_data(), complete_dims_, remaining_dims
+ )
+
+ # Rechunk ancillary variables that span complete_dims
+ for anc_var in cube.ancillary_variables():
+ dims = cube.ancillary_variable_dims(anc_var)
+ complete_dims_ = [dims.index(d) for d in complete_dims if d in dims]
+ if complete_dims_ and anc_var.has_lazy_data():
+ anc_var.data = _rechunk(
+ anc_var.lazy_data(), complete_dims_, remaining_dims
+ )
+
+
+def rechunk_cube(
+ cube: Cube,
+ complete_coords: Iterable[Coord | str],
+ remaining_dims: int | Literal['auto'] = 'auto',
+) -> Cube:
+ """Rechunk cube so that it is not chunked along given dimensions.
+
+ This will rechunk the cube's data, but also all non-dimensional
+ coordinates, cell measures, and ancillary variables that span at least one
+ of the given dimensions.
+
+ Note
+ ----
+ This will only rechunk `dask` arrays. `numpy` arrays are not changed.
+
+ Parameters
+ ----------
+ cube:
+ Input cube.
+ complete_coords:
+ (Names of) coordinates along which the output cubes should not be
+ chunked. The given coordinates must span exactly 1 dimension.
+ remaining_dims:
+ Chunksize of the remaining dimensions.
+
+ Returns
+ -------
+ Cube
+ Rechunked cube. This will always be a copy of the input cube.
+
+ """
+ cube = cube.copy() # do not modify input cube
+
+ # Make sure that complete_coords span exactly 1 dimension
+ complete_dims = []
+ for coord in complete_coords:
+ coord = cube.coord(coord)
+ dims = cube.coord_dims(coord)
+ if len(dims) != 1:
+ raise CoordinateMultiDimError(
+ f"Complete coordinates must be 1D coordinates, got "
+ f"{len(dims):d}D coordinate '{coord.name()}'"
+ )
+ complete_dims.append(dims[0])
+
+ # Rechunk data
+ if cube.has_lazy_data():
+ cube.data = _rechunk(cube.lazy_data(), complete_dims, remaining_dims)
+
+ # Rechunk dimensional metadata
+ _rechunk_dim_metadata(cube, complete_dims, remaining_dims=remaining_dims)
+
+ return cube

esmvalcore/preprocessor/__init__.py

@@ -70,6 +70,7 @@
  extract_season,
  extract_time,
  hourly_statistics,
+ local_solar_time,
  monthly_statistics,
  regrid_time,
  resample_hours,
@@ -148,17 +149,14 @@
  'extract_volume',
  'extract_trajectory',
  'extract_transect',
- # 'average_zone': average_zone,
- # 'cross_section': cross_section,
  'detrend',
  'extract_named_regions',
  'axis_statistics',
  'depth_integration',
  'area_statistics',
  'volume_statistics',
  # Time operations
- # 'annual_cycle': annual_cycle,
- # 'diurnal_cycle': diurnal_cycle,
+ 'local_solar_time',
  'amplitude',
  'zonal_statistics',
  'meridional_statistics',