Add Hovmoeller MJO diagnostic (#4564)

Co-authored-by: Manuel Schlund <manuel.schlund@dlr.de>
Co-authored-by: Manuel Schlund <32543114+schlunma@users.noreply.github.com>
Co-authored-by: Lisa Bock <lisa.bock@dlr.de>

Author: 3 people

doc/sphinx/source/recipes/figures/mjo/era5_mjo_hovmoeller.png

@@ -51,6 +51,7 @@ Atmosphere
    recipe_iht_toa
    recipe_impact
    recipe_lifetime
+   recipe_mjo_hovmoeller
    recipe_modes_of_variability
    recipe_mpqb_xch4
    recipe_quantilebias

doc/sphinx/source/recipes/index.rst

@@ -0,0 +1,131 @@
+.. _recipes_mjo_hovmoeller:
+
+Madden-Julian Oscillation precipitation Hovmöller diagnostic
+==============================================================
+
+Overview
+--------
+
+This recipe computes lag-regression Hovmöller diagrams of the
+Madden-Julian Oscillation (MJO). For each dataset, daily tropical
+precipitation is averaged over a latitude band and its day-of-year
+climatology is removed. The resulting anomalies are Lanczos band-pass
+filtered to the 20-100 day MJO period range. A reference index is built
+by averaging the filtered field over a reference longitude sector, and
+the full filtered field is then regressed against this index at a range
+of lags.
+
+The result is a longitude-lag diagram: an eastward-propagating diagonal
+band is the signature of MJO convection.
+
+
+Available recipes and diagnostics
+----------------------------------
+
+Recipes are stored in recipes/
+
+* recipe_mjo_hovmoeller.yml
+
+Diagnostics are stored in diag_scripts/mjo/
+
+* mjo_hovmoeller.py: compute the lag regression and plot the Hovmöller
+  diagram.
+
+
+User settings in recipe
+------------------------
+
+#. Script mjo_hovmoeller.py
+
+   *Required settings for script*
+
+   * ``reference_longitudes``: longitude sector ``[lon0, lon1]`` (in
+     degrees East) used to build the MJO reference index that the
+     filtered field is regressed against.
+   * ``low_period``: lower period cutoff (in days) of the Lanczos
+     band-pass filter.
+   * ``high_period``: upper period cutoff (in days) of the Lanczos
+     band-pass filter.
+   * ``lanczos_weights``: number of weights of the Lanczos band-pass
+     filter. Must be an odd integer greater than 1.
+   * ``max_lag``: maximum lag (in days, in both directions) computed by
+     the lag regression.
+
+   *Optional settings for script*
+
+   * ``longitude_limits``: longitude axis limits of the Hovmöller plot
+     (default: ``[0.0, 360.0]``).
+   * ``contour_levels``: number of contour levels in the Hovmöller plot.
+     Must be at least 3 (default: ``21``).
+   * ``colormap``: matplotlib colormap used for the Hovmöller contour
+     plot (default: ``RdYlBu``).
+   * ``plot_title``: title of the Hovmöller plot (default: ``MJO
+     Hovmöller diagram``).
+   * ``colorbar_label``: label for the figure's colorbar (default:
+     ``Precipitation regression coefficient``).
+
+   *Required settings for variables*
+
+   * none beyond the standard ``short_name``, ``mip``, ``preprocessor``
+     and ``timerange``.
+
+   *Optional settings for variables*
+
+   * none
+
+   *Required settings for preprocessor*
+
+   * ``extract_region``: restrict the data to the tropical latitude
+     band used for the diagnostic.
+   * ``regrid``: regrid all datasets onto a common regular grid.
+   * ``meridional_statistics``: average over the extracted latitude
+     band (``operator: mean``).
+   * ``daily_statistics``: reduce the data to daily means
+     (``operator: mean``).
+   * ``anomalies``: remove the day-of-year climatology
+     (``period: day``).
+   * ``convert_units``: convert precipitation to ``kg m-2 day-1``.
+
+   *Optional settings for preprocessor*
+
+   * none
+
+   *Color tables*
+
+   * none
+
+
+Variables
+---------
+
+* pr (atmos, daily mean, longitude latitude time)
+
+
+Observations and reformat scripts
+----------------------------------
+
+*Note: ERA5 is read directly through ESMValCore's native6 support; no
+separate reformat script needs to be run beforehand.*
+
+* ERA5 (native6 project, tier 3, ``frequency: 1hr``)
+
+
+References
+----------
+
+* Hannah, W. M., Jones, C. R., Hillman, B. R., Norman, M. R., Bader, D. C., Taylor, M. A., et al. (2020).
+  Initial results from the super-parameterized E3SM. Journal of Advances in Modeling Earth Systems. 12,
+  e2019MS001863. https://doi.org/10.1029/2019MS001863
+
+
+Example plots
+-------------
+
+.. _fig_mjo_hovmoeller_1:
+.. figure:: /recipes/figures/mjo/era5_mjo_hovmoeller.png
+   :align: center
+
+   Lag regression of 20-100 day filtered ERA5 precipitation against a
+   precipitation index averaged over 80-100E, 1979-1983. Positive
+   longitude-lag slope through the reference sector shows the
+   eastward-propagating MJO precipitation signal.

doc/sphinx/source/recipes/recipe_mjo_hovmoeller.rst

@@ -83,7 +83,9 @@
 
 
 def _convert_units(
-    cubes: iris.cube.CubeList, short_name: str, var: dict
+    cubes: iris.cube.CubeList,
+    short_name: str,
+    var: dict,
 ) -> iris.cube.Cube:
     """Perform variable-specific conversion of units.
 
@@ -187,7 +189,7 @@ def _extract_variable(in_files, var, cfg, out_dir, year, month):
             [
                 timecoord.units.date2num(start_date),
                 timecoord.units.date2num(end_date),
-            ]
+            ],
         )
 
     # Add longitude coordinate to cube only for o3_sage_omps.

esmvaltool/cmorizers/data/formatters/datasets/esacci_ozone.py

@@ -550,6 +550,10 @@ authors:
     institute: DLR, Germany
     orcid:
     github: ellensarauer
+  schoenfeld_jurij:
+    name: Schönfeld, Jurij
+    institute: DLR, Germany
+    orcid: https://orcid.org/0009-0000-7453-9517
   schulze_kirsten:
     name: Schulze, Kirsten
     institute: Uni Bremen, Germany

esmvaltool/config-references.yml

@@ -0,0 +1 @@
+"""Diagnostics for evaluating the Madden-Julian Oscillation."""

esmvaltool/diag_scripts/mjo/__init__.py

@@ -0,0 +1,100 @@
+"""Lanczos filtering utilities for MJO diagnostics."""
+
+import numpy as np
+
+
+def lanczos_weights(low_period, high_period, weights):
+    """Return Lanczos band-pass weights for daily data."""
+    if weights < 3 or weights % 2 == 0:
+        msg = "weights must be an odd integer greater than 1"
+        raise ValueError(msg)
+    if low_period <= 0 or high_period <= low_period:
+        msg = "Expected 0 < low_period < high_period, e.g. 20 and 100 days"
+        raise ValueError(msg)
+
+    half_width = weights // 2
+    n = np.arange(-half_width, half_width + 1, dtype=float)
+
+    # Cutoff frequencies in cycles per day.
+    low_frequency = 1.0 / high_period
+    high_frequency = 1.0 / low_period
+
+    coefficients = np.empty(weights, dtype=float)
+    coefficients[half_width] = 2.0 * (high_frequency - low_frequency)
+
+    nonzero = n != 0
+    coefficients[nonzero] = (
+        np.sin(2.0 * np.pi * high_frequency * n[nonzero])
+        - np.sin(2.0 * np.pi * low_frequency * n[nonzero])
+    ) / (np.pi * n[nonzero])
+
+    # Lanczos sigma window.
+    coefficients *= np.sinc(n / (half_width + 1.0))
+
+    # Do not normalize by coefficients.sum(): a band-pass filter has
+    # approximately zero response at zero frequency.
+    return coefficients
+
+
+def _filter_series(series, coefficients):
+    """Filter one time series and mask invalid endpoints."""
+    series = np.ma.asarray(series, dtype=float)
+    values = series.filled(np.nan)
+
+    valid = np.isfinite(values)
+    numerator = np.convolve(
+        np.where(valid, values, 0.0),
+        coefficients,
+        mode="same",
+    )
+
+    # Require every value in the filter window to be valid.
+    valid_count = np.convolve(
+        valid.astype(int),
+        np.ones(coefficients.size, dtype=int),
+        mode="same",
+    )
+    output = np.ma.masked_where(valid_count < coefficients.size, numerator)
+
+    half_width = coefficients.size // 2
+    output[:half_width] = np.ma.masked
+    output[-half_width:] = np.ma.masked
+    return output
+
+
+def lanczos_bandpass(
+    cube,
+    low_period=20,
+    high_period=100,
+    weights=91,
+):
+    """Apply a Lanczos band-pass filter along the cube time dimension."""
+    if not cube.coords("time"):
+        msg = "Input cube has no time coordinate"
+        raise ValueError(msg)
+
+    time_axis = cube.coord_dims("time")[0]
+    coefficients = lanczos_weights(low_period, high_period, weights)
+
+    data = np.moveaxis(np.ma.asarray(cube.data), time_axis, 0)
+    flattened = data.reshape(data.shape[0], -1)
+
+    filtered = np.ma.empty(flattened.shape, dtype=float)
+    for column in range(flattened.shape[1]):
+        filtered[:, column] = _filter_series(
+            flattened[:, column],
+            coefficients,
+        )
+
+    filtered = filtered.reshape(data.shape)
+    filtered = np.moveaxis(filtered, 0, time_axis)
+
+    result = cube.copy(data=filtered)
+    result.long_name = (
+        f"{cube.name()} {low_period}-{high_period} day filtered anomalies"
+    )
+    result.attributes["temporal_filter"] = (
+        f"Lanczos band-pass; periods={low_period}-{high_period} days; "
+        f"weights={weights}"
+    )
+    return result