Download this notebook from github.

Ensembles

An important aspect of climate models is that they are run multiple times with some initial perturbations to see how they replicate the natural variability of the climate. Through xclim.ensembles, xclim provides an easy interface to compute ensemble statistics on different members. Most methods perform checks and conversion on top of simpler xarray methods, providing an easier interface to use.

create_ensemble

Our first step is to create an ensemble. This method takes a list of files defining the same variables over the same coordinates and concatenates them into one dataset with an added dimension realization.

Using xarray a very simple way of creating an ensemble dataset would be :

import xarray

xarray.open_mfdataset(files, concat_dim='realization')

However, this is only successful when the dimensions of all the files are identical AND only if the calendar type of each netcdf file is the same

xclim’s create_ensemble() method overcomes these constraints, selecting the common time period to all files and assigns a standard calendar type to the dataset.

Input netcdf files still require equal spatial dimension size (e.g. lon, lat dimensions).

Given files all named ens_tas_m[member number].nc, we use glob to get a list of all those files.

import glob

import xarray as xr

import xclim as xc

# Set display to HTML sytle (for fancy output)
xr.set_options(display_style="html", display_width=50)

import matplotlib.pyplot as plt

%matplotlib inline

from xclim import ensembles

ens = ensembles.create_ensemble(glob.glob("ens_tas_m*.nc")).load()
ens.close()

/home/docs/checkouts/readthedocs.org/user_builds/xclim/conda/v0.43.0/lib/python3.10/site-packages/xclim/indices/fire/_cffwis.py:207: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def _day_length(lat: int | float, mth: int):  # pragma: no cover
/home/docs/checkouts/readthedocs.org/user_builds/xclim/conda/v0.43.0/lib/python3.10/site-packages/xclim/indices/fire/_cffwis.py:227: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def _day_length_factor(lat: float, mth: int):  # pragma: no cover

plt.style.use("seaborn-v0_8-dark")
plt.rcParams["figure.figsize"] = (13, 5)
ens.tas.plot(hue="realization")
plt.show()

ens.tas  # Attributes of the first dataset to be opened are copied to the final output

<xarray.DataArray 'tas' (realization: 10,
                         time: 1461)>
array([[ -0.74661228,  -1.14834882,  -1.53530599, ..., -17.61971376,
        -17.56192347, -17.48840648],
       [ -1.60477585,  -1.53196209,  -1.45478034, ...,  -1.29130334,
         -1.32623196,  -1.35751854],
       [-14.54425643, -14.41094884, -14.26933018, ...,  -2.77476123,
         -2.3891895 ,  -1.98746361],
       ...,
       [-13.22850295, -12.24738429, -11.29336498, ..., -15.67208672,
        -15.77506349, -15.87533369],
       [-14.46163135, -14.49828035, -14.52805846, ..., -12.03841506,
        -11.51116215, -10.95658622],
       [ -6.31010317,  -5.97253002,  -5.64234933, ...,  -7.47462733,
         -7.38653016,  -7.29506706]])
Coordinates:
  * time         (time) datetime64[ns] 2000-01-...
  * realization  (realization) int64 0 1 2 ... 8 9
Attributes:
    units:          degC
    standard_name:  air_temperature
    long_name:      Mean air temperature at sur...
    title:          tas of member 01

xarray.DataArray

'tas'

• realization: 10
• time: 1461

• -0.7466 -1.148 -1.535 -1.907 -2.265 ... -7.559 -7.475 -7.387 -7.295

  array([[ -0.74661228,  -1.14834882,  -1.53530599, ..., -17.61971376,
          -17.56192347, -17.48840648],
         [ -1.60477585,  -1.53196209,  -1.45478034, ...,  -1.29130334,
           -1.32623196,  -1.35751854],
         [-14.54425643, -14.41094884, -14.26933018, ...,  -2.77476123,
           -2.3891895 ,  -1.98746361],
         ...,
         [-13.22850295, -12.24738429, -11.29336498, ..., -15.67208672,
          -15.77506349, -15.87533369],
         [-14.46163135, -14.49828035, -14.52805846, ..., -12.03841506,
          -11.51116215, -10.95658622],
         [ -6.31010317,  -5.97253002,  -5.64234933, ...,  -7.47462733,
           -7.38653016,  -7.29506706]])

• Coordinates: (2)
  • time
    (time)
    datetime64[ns]
    2000-01-01 ... 2003-12-31

    array(['2000-01-01T00:00:00.000000000',
           '2000-01-02T00:00:00.000000000',
           '2000-01-03T00:00:00.000000000', ...,
           '2003-12-29T00:00:00.000000000',
           '2003-12-30T00:00:00.000000000',
           '2003-12-31T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • realization
    (realization)
    int64
    0 1 2 3 4 5 6 7 8 9

    axis :

    E

    array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-02', '2000-01-03', '2000-01-04',
                   '2000-01-05', '2000-01-06', '2000-01-07', '2000-01-08',
                   '2000-01-09', '2000-01-10',
                   ...
                   '2003-12-22', '2003-12-23', '2003-12-24', '2003-12-25',
                   '2003-12-26', '2003-12-27', '2003-12-28', '2003-12-29',
                   '2003-12-30', '2003-12-31'],
                  dtype='datetime64[ns]', name='time', length=1461, freq=None))

  • realization
    PandasIndex

    PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='realization'))

• Attributes: (4)

  units :

  degC

  standard_name :

  air_temperature

  long_name :

  Mean air temperature at surface

  title :

  tas of member 01

Ensemble statistics

Beyond creating an ensemble dataset, the xclim.ensembles module contains functions for calculating statistics between realizations

Ensemble mean, standard-deviation, max & min

In the example below, we use xclim’s ensemble_mean_std_max_min() to calculate statistics across the 10 realizations in our test dataset. Output variables are created combining the original variable name tas with additional ending indicating the statistic calculated on the realization dimension : _mean, _stdev, _min, _max

The resulting output now contains 4 derived variables from the original single variable in our ensemble dataset.

ens_stats = ensembles.ensemble_mean_std_max_min(ens)
ens_stats

<xarray.Dataset>
Dimensions:    (time: 1461)
Coordinates:
  * time       (time) datetime64[ns] 2000-01-01...
Data variables:
    tas_mean   (time) float64 -8.165 ... -8.005
    tas_stdev  (time) float64 5.79 5.589 ... 5.398
    tas_max    (time) float64 0.633 ... -1.358
    tas_min    (time) float64 -15.26 ... -17.49
Attributes:
    history:  [2023-05-09 21:46:39] : Computati...

xarray.Dataset

• Dimensions:
  • time: 1461
• Coordinates: (1)
  • time
    (time)
    datetime64[ns]
    2000-01-01 ... 2003-12-31

    array(['2000-01-01T00:00:00.000000000',
           '2000-01-02T00:00:00.000000000',
           '2000-01-03T00:00:00.000000000', ...,
           '2003-12-29T00:00:00.000000000',
           '2003-12-30T00:00:00.000000000',
           '2003-12-31T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (4)
  • tas_mean
    (time)
    float64
    -8.165 -8.127 ... -8.078 -8.005

    units :

    degC

    standard_name :

    air_temperature

    long_name :

    Mean air temperature at surface

    title :

    tas of member 01

    array([-8.16530078, -8.12670256, -8.08432193, ..., -8.14376852,
           -8.07809341, -8.00524327])

  • tas_stdev
    (time)
    float64
    5.79 5.589 5.425 ... 5.486 5.398

    units :

    degC

    standard_name :

    air_temperature

    long_name :

    Mean air temperature at surface

    title :

    tas of member 01

    array([5.79020463, 5.58913219, 5.42523228, ..., 5.58371659, 5.48585542,
           5.39847414])

  • tas_max
    (time)
    float64
    0.633 -0.2453 ... -1.326 -1.358

    units :

    degC

    standard_name :

    air_temperature

    long_name :

    Mean air temperature at surface

    title :

    tas of member 01

    array([ 0.63295006, -0.24533423, -1.09541398, ..., -1.29130334,
           -1.32623196, -1.35751854])

  • tas_min
    (time)
    float64
    -15.26 -15.61 ... -17.56 -17.49

    units :

    degC

    standard_name :

    air_temperature

    long_name :

    Mean air temperature at surface

    title :

    tas of member 01

    array([-15.25890853, -15.61373987, -15.94539174, ..., -17.61971376,
           -17.56192347, -17.48840648])

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-02', '2000-01-03', '2000-01-04',
                   '2000-01-05', '2000-01-06', '2000-01-07', '2000-01-08',
                   '2000-01-09', '2000-01-10',
                   ...
                   '2003-12-22', '2003-12-23', '2003-12-24', '2003-12-25',
                   '2003-12-26', '2003-12-27', '2003-12-28', '2003-12-29',
                   '2003-12-30', '2003-12-31'],
                  dtype='datetime64[ns]', name='time', length=1461, freq=None))

• Attributes: (1)

  history :

  [2023-05-09 21:46:39] : Computation of statistics on 10 ensemble members. - xclim version: 0.43.0

Ensemble percentiles

Here, we use xclim’s ensemble_percentiles() to calculate percentile values across the 10 realizations. The output has now a percentiles dimension instead of realization. Split variables can be created instead, by specifying split=True (the variable name tas will be appended with _p{x}). Compared to NumPy’s percentile() and xarray’s quantile(), this method handles more efficiently dataset with invalid values and the chunking along the realization dimension (which is automatic when dask arrays are used).

ens_perc = ensembles.ensemble_percentiles(ens, values=[15, 50, 85], split=False)
ens_perc

<xarray.Dataset>
Dimensions:      (time: 1461, percentiles: 3)
Coordinates:
  * time         (time) datetime64[ns] 2000-01-...
  * percentiles  (percentiles) int64 15 50 85
Data variables:
    tas          (time, percentiles) float64 -1...
Attributes:
    units:          degC
    standard_name:  air_temperature
    long_name:      Mean air temperature at sur...
    title:          tas of member 01
    history:        [2023-05-09 21:46:39] : Com...

xarray.Dataset

• Dimensions:
  • time: 1461
  • percentiles: 3
• Coordinates: (2)
  • time
    (time)
    datetime64[ns]
    2000-01-01 ... 2003-12-31

    array(['2000-01-01T00:00:00.000000000',
           '2000-01-02T00:00:00.000000000',
           '2000-01-03T00:00:00.000000000', ...,
           '2003-12-29T00:00:00.000000000',
           '2003-12-30T00:00:00.000000000',
           '2003-12-31T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • percentiles
    (percentiles)
    int64
    15 50 85

    array([15, 50, 85])

• Data variables: (1)
  • tas
    (time, percentiles)
    float64
    -14.52 -8.066 ... -7.064 -2.533

    units :

    degC

    standard_name :

    air_temperature

    long_name :

    Mean air temperature at surface

    title :

    tas of member 01

    history :

    [2023-05-09 21:46:39] : Computation of the percentiles on 10 ensemble members. - xclim version: 0.43.0

    array([[-14.51533765,  -8.06558363,  -1.04696953],
           [-14.46771432,  -7.79924856,  -1.28261346],
           [-14.43750356,  -7.53961213,  -1.48296432],
           ...,
           [-14.48397046,  -6.91629283,  -2.81763514],
           [-14.34807023,  -6.98922344,  -2.67832622],
           [-14.19967223,  -7.06413705,  -2.53270064]])

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2000-01-01', '2000-01-02', '2000-01-03', '2000-01-04',
                   '2000-01-05', '2000-01-06', '2000-01-07', '2000-01-08',
                   '2000-01-09', '2000-01-10',
                   ...
                   '2003-12-22', '2003-12-23', '2003-12-24', '2003-12-25',
                   '2003-12-26', '2003-12-27', '2003-12-28', '2003-12-29',
                   '2003-12-30', '2003-12-31'],
                  dtype='datetime64[ns]', name='time', length=1461, freq=None))

  • percentiles
    PandasIndex

    PandasIndex(Index([15, 50, 85], dtype='int64', name='percentiles'))

• Attributes: (5)

  units :

  degC

  standard_name :

  air_temperature

  long_name :

  Mean air temperature at surface

  title :

  tas of member 01

  history :

  [2023-05-09 21:46:39] : Computation of the percentiles on 10 ensemble members. - xclim version: 0.43.0

fig, ax = plt.subplots()
ax.fill_between(
    ens_stats.time.values,
    ens_stats.tas_min,
    ens_stats.tas_max,
    alpha=0.3,
    label="Min-Max",
)
ax.fill_between(
    ens_perc.time.values,
    ens_perc.tas.sel(percentiles=15),
    ens_perc.tas.sel(percentiles=85),
    alpha=0.5,
    label="Perc. 15-85",
)
ax._get_lines.get_next_color()  # Hack to get different line
ax._get_lines.get_next_color()
ax.plot(ens_stats.time.values, ens_stats.tas_mean, linewidth=2, label="Mean")
ax.plot(
    ens_perc.time.values, ens_perc.tas.sel(percentiles=50), linewidth=2, label="Median"
)
ax.legend()
plt.show()