Note

This page was generated from an Jupyter notebook that can be accessed from github.

MESMER-X: create emulations

Using calibrated parameters, we can emulate a large ensemble of realisations of climate variables with MESMER-X, see tutorials on calibrating MESMER-X in the tutorial section. Here, we continue the tutorial and emulate the annual maximum temperature (annual maximum tasmax or TXx) with a normal distribution with a linear dependence on global mean temperature. As before the following warnings apply:

Note

• We fit annual maximum data (i.e. a block maxima) using a normal distribution. This is not the appropriate distribution for this type of data but is used for illustrative purposes and speed.

• A normal distribution with a linear dependence of the mean on a variable is equivalent to a linear regression with is much faster to use directly. Again, this is only for illustrative purposes.

import pathlib

import filefisher
import matplotlib.pyplot as plt
import xarray as xr

import mesmer
from mesmer.distrib import (
    ConditionalDistribution,
    Expression,
    ProbabilityIntegralTransform,
)

To generate emulations the workflow of the calibration is reversed, using the estimated parameters. Here, we use the same local annual mean temperatures to force the emulations, but temperatures from other models, scenarios, ensemble members or emulated annual local temperatures can be used as well.

Configuration

# set some configuration parameters
REFERENCE_PERIOD = slice("1850", "1900")


n_realizations = 3
seed = 121638161038369005437100875176885634488
buffer = 10
esm = "IPSL-CM6A-LR"
scenario_train = "ssp585"

targ_var = "tasmax"

Where the seed is random but constant (for reporducibility) and was created using

import secrets

secrets.randbits(64)

Load data

In contrast to the calibration we concatenate the historical simulations and the projection.

cmip6_data_path = mesmer.example_data.cmip6_ng_path()


def load_cmip6_ng(variable):

    # load predictor data
    path = cmip6_data_path / variable / "ann" / "g025"

    fN_hist = path / f"{variable}_ann_{esm}_historical_r1i1p1f1_g025.nc"
    fN_ssp585 = path / f"{variable}_ann_{esm}_ssp585_r1i1p1f1_g025.nc"

    time_coder = xr.coders.CFDatetimeCoder(use_cftime=True)
    data_hist = xr.open_dataset(fN_hist, decode_times=time_coder).drop_vars(
        ["height", "file_qf", "time_bnds"], errors="ignore"
    )
    data_ssp585 = xr.open_dataset(fN_ssp585, decode_times=time_coder).drop_vars(
        ["height", "file_qf", "time_bnds"], errors="ignore"
    )

    # concat
    data = xr.concat([data_hist, data_ssp585], dim="time")

    data = data - data.sel(time=REFERENCE_PERIOD).mean("time")

    return data

tas_gridded = load_cmip6_ng("tas")

# make global mean
predictor = mesmer.weighted.global_mean(tas_gridded)

time = predictor.time

Load tasmax data - for illustration purposes only

tasmax_gridded = load_cmip6_ng("tasmax")

THRESHOLD_LAND = 1 / 3


# mask ocean, Antarctica and stack the gridpoints
def mask_and_stack(ds, threshold_land):
    ds = mesmer.mask.mask_ocean_fraction(ds, threshold_land)
    ds = mesmer.mask.mask_antarctica(ds)
    ds = mesmer.grid.stack_lat_lon(ds)
    return ds


tasmax_stacked = mask_and_stack(tasmax_gridded, threshold_land=THRESHOLD_LAND)

Load parameters

param_path = pathlib.Path("./output/calibrated_parameters/mesmer-x/")
variable = "txx"

# load the parameters
PARAM_FILEFINDER = filefisher.FileFinder(
    path_pattern=param_path / "{variable}" / "{esm}_{scen}",
    file_pattern="params_{module}_{esm}_{scen}.nc",
)


distrib_file = PARAM_FILEFINDER.create_full_name(
    variable=variable, module="distrib", esm=esm, scen=scenario_train
)
local_ar_file = PARAM_FILEFINDER.create_full_name(
    variable=variable, module="local-variability", esm=esm, scen=scenario_train
)
localized_ecov_file = PARAM_FILEFINDER.create_full_name(
    variable=variable, module="covariance", esm=esm, scen=scenario_train
)

transform_coeffs = xr.open_dataset(distrib_file)
local_ar_params = xr.open_dataset(local_ar_file)
localized_ecov = xr.open_dataset(localized_ecov_file)

Create temporally and spatially correlated samples

# generate realizations based on the auto-regression with spatially correlated innovations
transf_emus = mesmer.stats.draw_auto_regression_correlated(
    local_ar_params,
    localized_ecov.localized_covariance_adjusted,
    time=time,
    realisation=n_realizations,
    seed=seed,
    buffer=buffer,
)
transf_emus = transf_emus.rename({"samples": targ_var})

Transform samples back to the conditional distribution

Initialize the conditional distribution with the calibrated values and the expression:

distrib_orig = ConditionalDistribution.from_dataset(transform_coeffs)

define a standard normal distribution and the probability integral transform:

# back-transform the realizations
expr_tranf = Expression("norm(loc=0, scale=1)", "standard_normal")
distrib_transf = ConditionalDistribution(expr_tranf)

back_pit = ProbabilityIntegralTransform(
    distrib_orig=distrib_transf,
    distrib_targ=distrib_orig,
)

and finally back-transform the realizations

emus = back_pit.transform(
    data=transf_emus, target_name=targ_var, preds_orig=None, preds_targ=predictor
)

Illustration

gridcell = 0


f, ax = plt.subplots()

emus.tasmax.isel(gridcell=gridcell).plot.line(
    ax=ax, x="time", color="#6baed6", add_legend=False
)

ax.plot([], [], label="realisations", color="#6baed6")

tasmax_stacked.tasmax.isel(gridcell=gridcell).plot(
    ax=ax, x="time", color="0.1", label="training data", lw=1
)


ax.legend()

ax.set_title("Training data and emulations at random grid point")

Text(0.5, 1.0, 'Training data and emulations at random grid point')