Note

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

MESMER-X: calibrate scenario#

In this tutorial we use a conditional distribution to model the dependence of a climate variable on a global predictor using the approach from MESMER-X (Quilquaille et al., 2022). Specifically, we model annual maximum temperature (TXx) with a normal distribution (see note) with a linear dependence on global mean temperature. We calibrate the historical and one additional scenario but the approach to calibrate several scenarios is the same.

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 xarray as xr

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

Configuration#

scenario = "ssp585"
target_name = "tasmax"
option_2ndfit = False

THRESHOLD_LAND = 1 / 3
REFERENCE_PERIOD = slice("1850", "1900")

esm = "IPSL-CM6A-LR"

# path of the example data
# NOTE: this downloads the tutorial data and caches it locally if needed
cmip6_data_path = mesmer.example_data.cmip6_ng_path(relative=True)

Load and prepare predictor and target data#

MESMER-X expects the predictor and forcing data in the same format as MESMER and MESMER-M. First define helper functions to load data from the cmip6-ng example data archive.

Caution

In the original cmip6-ng archive annual tasmax is not the annual maximum temperature (TXx) but the annual mean daily maximum temperature. In MESMER’s example data archive the annual maximum is provided.

def find_files(variable, model, scenarios, member=None):

    CMIP_FILEFINDER = filefisher.FileFinder(
        path_pattern=str(cmip6_data_path / "{variable}/{time_res}/{resolution}"),
        file_pattern="{variable}_{time_res}_{model}_{scenario}_{member}_{resolution}.nc",
    )

    fc_scens_pred = CMIP_FILEFINDER.find_files(
        variable=variable,
        scenario=scenarios,
        model=model,
        member=member,
        resolution="g025",
        time_res="ann",
    )

    # get the historical members that are also in the future scenarios, but only once
    unique_scen_members = fc_scens_pred.df.member.unique()

    fc_hist_pred = CMIP_FILEFINDER.find_files(
        variable=variable,
        scenario="historical",
        model=model,
        resolution="g025",
        time_res="ann",
        member=unique_scen_members,
    )

    fc_pred = fc_hist_pred.concat(fc_scens_pred)

    return fc_pred

def load_data(fc):
    data = xr.DataTree()

    scenarios = fc.df.scenario.unique().tolist()

    for scen in scenarios:
        # load data for each scenario
        data_scen = []

        for var in fc.df.variable.unique():
            files = fc.search(variable=var, scenario=scen)

            # load all members for a scenario
            members = []
            for fN, meta in files.items():
                time_coder = xr.coders.CFDatetimeCoder(use_cftime=True)
                ds = xr.open_dataset(fN, decode_times=time_coder)
                # drop unnecessary variables
                ds = ds.drop_vars(
                    ["height", "time_bnds", "file_qf", "area"], errors="ignore"
                )
                # assign member-ID as coordinate
                ds = ds.assign_coords({"member": meta["member"]})
                members.append(ds)

            # create a Dataset that holds each member along the member dimension
            data_var = xr.concat(members, dim="member")
            data_scen.append(data_var)

        data_scen = xr.merge(data_scen)
        data[scen] = xr.DataTree(data_scen)
    return data

load predictor data: tas#

fc_pred = find_files(variable="tas", model=esm, scenarios=scenario, member="r1i1p1f1")
fc_pred.df

	variable	time_res	resolution	model	scenario	member
path
../data/cmip6-ng/tas/ann/g025/tas_ann_IPSL-CM6A-LR_historical_r1i1p1f1_g025.nc	tas	ann	g025	IPSL-CM6A-LR	historical	r1i1p1f1
../data/cmip6-ng/tas/ann/g025/tas_ann_IPSL-CM6A-LR_ssp585_r1i1p1f1_g025.nc	tas	ann	g025	IPSL-CM6A-LR	ssp585	r1i1p1f1

tas = load_data(fc_pred)
# calculate anomalies w.r.t. the reference period
tas_anom = mesmer.anomaly.calc_anomaly(tas, reference_period=REFERENCE_PERIOD)
# calculate global mean
tas_glob_mean = mesmer.weighted.global_mean(tas_anom)

Define the Expression and ConditionalDistribution#

The Expression defines the used distribution and the covariate structure and needs to be created as a string and has the following elements:

distribution: name of the desired distribution as written in scipy.stats
parameters: the parameters of the distribution as written in scipy (e.g. loc, scale, …)
predictors: predictors that the distribution will be conditional on, must be written as variable in python and surrounded by __ (e.g. __tas__)
coefficients: of the predictors and parameters, must be named c1, c2, …, c9.
mathematical operations: mathematical operations for the covariates are written as python code. Only operators (e.g. *) functions from numpy (np.*) are supported.

Thus for the normal distribution with a linear dependence on global mean temperature (tas) this looks as follows:

expr = "norm(loc=c1 + c2 * __tas__, scale=c3)"

expr_name = "expr1"
expression = Expression(expr, expr_name, boundaries_params={}, boundaries_coeffs={})

We then pass the Expression to the ConditionalDistribution class, which finds the best estimate of the coefficients for given climate model data.

distrib = ConditionalDistribution(expression=expression)

Transform to standard normal distribution#

After fitting the coefficients the original data is transformed to a standard normal distribution, also removing the fitted linear trend:

# probability integral transform on non-stacked data for AR(1) process
target_expression = Expression("norm(loc=0, scale=1)", "standard_normal")

pit = ProbabilityIntegralTransform(
    distrib_orig=distrib,
    distrib_targ=ConditionalDistribution(target_expression),
)
transf_target = pit.transform(
    data=targ_data, target_name=target_name, preds_orig=pred_data, preds_targ=None
)

Local variability#

The resulting residuals constitute the local variability and are now treated the same way as in MESMER by fitting a local AR(1) process with a spatially correlated noise term.

# training of auto-regression with spatially correlated innovations
local_ar_params = mesmer.stats.fit_auto_regression_scen_ens(
    transf_target,
    ens_dim="member",
    dim="time",
    lags=1,
)

# estimate covariance matrix
# prep distance matrix
geodist = mesmer.geospatial.geodist_exact(
    lon=targ_data["historical"].lon, lat=targ_data["historical"].lat
)
# prep localizer
localisation_radii = range(3_500, 15_001, 250)
phi_gc_localizer = mesmer.stats.gaspari_cohn_correlation_matrices(
    geodist=geodist, localisation_radii=localisation_radii
)

pooled_transf_target = mesmer.datatree.pool_scen_ens(transf_target)


localized_ecov = mesmer.stats.find_localized_empirical_covariance(
    data=pooled_transf_target[target_name],
    weights=pooled_weights.weights.squeeze(),
    localizer=phi_gc_localizer,
    dim="sample",
    k_folds=30,
)

localized_ecov["localized_covariance_adjusted"] = mesmer.stats.adjust_covariance_ar1(
    localized_ecov.localized_covariance, local_ar_params.coeffs
)

Potentially save the params:

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

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

# NOTE: to save disk space the "covariance" and "localized_covariance" can be
# removed from localized_ecov
localized_ecov = localized_ecov.drop_vars(
    ["covariance", "localized_covariance"], errors="ignore"
)

params = {
    "distrib": transform_coeffs,
    "local-variability": local_ar_params,
    "covariance": localized_ecov,
}


save_files = False  # we don't save them here in the example
if save_files:

    for module, param in params.items():
        filename = PARAM_FILEFINDER.create_full_name(
            variable=variable, module=module, esm=esm, scen=scenario
        )

        filename = pathlib.Path(filename)
        filename.parent.mkdir(parents=True, exist_ok=True)

        param.to_netcdf(filename)

MESMER-X: calibrate scenario

Contents

MESMER-X: calibrate scenario#

Configuration#

Load and prepare predictor and target data#

load predictor data: tas#

load target data: annual maximum tasmax (txx)#

Define the Expression and ConditionalDistribution#

Fit conditional distribution#

Transform to standard normal distribution#

Local variability#