Note

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

MESMER-M: calibrate multiple scenarios#

Training monthly local temperature from yearly local temperature for multiple scenarios and ensemble members. We use an example data set on a coarse grid. This roughly follows the approach outlined in Nath et al. (2022).

MESMER-M trains the local monthly temperature using the local annual temperature (i.e. the temperature from the same grid point) as forcing. This is different from MESMER which uses global mean values as predictors for where local annual mean temperatures. Training MESMER-M consists of 4 steps:

harmonic model: fit the seasonal cycle with a harmonic model
power transformer: make the resulting residuals more normal by using a Yeo-Johnson transformation
cyclo-stationary AR(1) process: the monthly residuals are assumed to follow a cyclo-stationary AR(1) process, where one months value depends on the previous one
local variability: estimate parameters needed to generate local variability

This example can be extended to more scenarios, ensemble members and higher resolution data. See also the MESMER-M calibration and emulation tests in tests/integration/.

import pathlib

import filefisher
import matplotlib.pyplot as plt
import pandas as pd
import scipy as sp
import xarray as xr

import mesmer

Configuration#

# target model
model = "IPSL-CM6A-LR"

# scenarios used for training
scenarios = ["ssp126", "ssp585"]

# gridcells considered as land
THRESHOLD_LAND = 1 / 3

# reference period
REFERENCE_PERIOD = slice("1850", "1900")

# localisation radius for covariance matrix
LOCALISATION_RADII = list(range(7_500, 12_501, 500))

# 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)

# mesmer-m downscales annual temperature to
# monthly resolution

# variable, shared for input and output
variable = "tas"

# spatial resolution, shared for input and output
spatial_resolution = "g025"

# targeted temporal resolution
temporal_resolution_target = "mon"

# temporal resolution of the forcing data
temporal_resolution_forcing = "ann"

Load Data#

We search monthly and annual mean temperatures using the library filefisher

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

We search for the initial (annual) data for ssp126 and ssp585 as well as the corresponding historical data:

fc_scens_y = CMIP_FILEFINDER.find_files(
    variable=variable,
    scenario=scenarios,
    model=model,
    resolution=spatial_resolution,
    time_res=temporal_resolution_forcing,
)

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

fc_hist_y = CMIP_FILEFINDER.find_files(
    variable=variable,
    scenario="historical",
    model=model,
    resolution=spatial_resolution,
    time_res=temporal_resolution_forcing,
    member=unique_scen_members_y,
)

fc_all_y = fc_hist_y.concat(fc_scens_y)
fc_all_y.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_historical_r2i1p1f1_g025.nc	tas	ann	g025	IPSL-CM6A-LR	historical	r2i1p1f1
../data/cmip6-ng/tas/ann/g025/tas_ann_IPSL-CM6A-LR_ssp126_r1i1p1f1_g025.nc	tas	ann	g025	IPSL-CM6A-LR	ssp126	r1i1p1f1
../data/cmip6-ng/tas/ann/g025/tas_ann_IPSL-CM6A-LR_ssp585_r1i1p1f1_g025.nc	tas	ann	g025	IPSL-CM6A-LR	ssp585	r1i1p1f1
../data/cmip6-ng/tas/ann/g025/tas_ann_IPSL-CM6A-LR_ssp585_r2i1p1f1_g025.nc	tas	ann	g025	IPSL-CM6A-LR	ssp585	r2i1p1f1

Similarly, we search for the target (monthly) data for ssp126 and ssp585 as well as the corresponding historical data:

fc_scens_m = CMIP_FILEFINDER.find_files(
    variable=variable,
    scenario=scenarios,
    model=model,
    resolution=spatial_resolution,
    time_res=temporal_resolution_target,
)

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

fc_hist_m = CMIP_FILEFINDER.find_files(
    variable="tas",
    scenario="historical",
    model=model,
    resolution="g025",
    time_res="mon",
    member=unique_scen_members_m,
)

fc_all_m = fc_hist_m.concat(fc_scens_m)
fc_all_m.df

	variable	time_res	resolution	model	scenario	member
path
../data/cmip6-ng/tas/mon/g025/tas_mon_IPSL-CM6A-LR_historical_r1i1p1f1_g025.nc	tas	mon	g025	IPSL-CM6A-LR	historical	r1i1p1f1
../data/cmip6-ng/tas/mon/g025/tas_mon_IPSL-CM6A-LR_historical_r2i1p1f1_g025.nc	tas	mon	g025	IPSL-CM6A-LR	historical	r2i1p1f1
../data/cmip6-ng/tas/mon/g025/tas_mon_IPSL-CM6A-LR_ssp126_r1i1p1f1_g025.nc	tas	mon	g025	IPSL-CM6A-LR	ssp126	r1i1p1f1
../data/cmip6-ng/tas/mon/g025/tas_mon_IPSL-CM6A-LR_ssp585_r1i1p1f1_g025.nc	tas	mon	g025	IPSL-CM6A-LR	ssp585	r1i1p1f1
../data/cmip6-ng/tas/mon/g025/tas_mon_IPSL-CM6A-LR_ssp585_r2i1p1f1_g025.nc	tas	mon	g025	IPSL-CM6A-LR	ssp585	r2i1p1f1

This found 1 ensemble member for SSP1-2.6 and two for SSP5-8.5 and the corresponding ones in the historical scenario.

Now we load all the files we found into a DataTree, a data structure provided by xarray. It is a container to hold xarray Dataset objects that are not alignable. This is useful for us since we have historical and future data, which have different time coordinates. Moreover, the scenarios may also have different numbers of members (as e.g., SSP1-2.6, which only has one). Thus, we store the data of each scenario in a Dataset with all its ensemble members along a member dimension. Then we store all the scenario datasets in one DataTree node. The DataTree allows us to perform computations on each of the scenarios separately.

We define a helper function to load the data from the cmip6_ng example data repository:

def load_data(filecontainer):

    out = xr.DataTree()

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

    # load data for each scenario
    for scen in scenarios:
        files = filecontainer.search(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"], 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
        scen_data = xr.concat(members, dim="member")
        # put the scenario dataset into the DataTree
        out[scen] = xr.DataTree(scen_data)

    return out

Load annual (yearly - y) and monthly (m) data:

tas_y_orig = load_data(fc_all_y)
tas_m_orig = load_data(fc_all_m)

This results in two DataTree objects, with 3 nodes, one for each scenario (click on Groups to see the individual Datasets for the three scenarios):

Saving the parameters#

Finally, we have calibrated all needed parameters and can save them. We can use filefisher to nicely create file names and save the parameters.

# define path relative to this notebook & create folder
param_path = pathlib.Path("./output/calibrated_parameters/mesmer-m/")

Create time coordinate#

We overwrote the time coordinate with ‘sample’. Thus, we need to get the original time coordinate to be able to validate our results later on. If it is not needed to align the final emulations with the original data, this can be omitted, the time coordinates can later be generated for example with

monthly_time = xr.cftime_range("1850-01-01", "2100-12-31", freq="MS", calendar="gregorian")
monthly_time = xr.DataArray(monthly_time, dims="time", coords={"time": monthly_time})

hist_time = tas_m.historical.time
scen_time = tas_m.ssp585.time
m_time = xr.concat([hist_time, scen_time], dim="time")

Store data for all modules#

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

scen_str = "-".join(scenarios)

folder = PARAM_FILEFINDER.create_path_name(esm=model, scen=scen_str)
pathlib.Path(folder).mkdir(exist_ok=True, parents=True)


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

params = {
    "harmonic-model": harmonic_model.params,
    "power-transformer": pt_coefficients,
    "local-variability": ar1_fit,
    "covariance": localized_ecov,
    "grid-orig": grid_orig,
    "monthly-time": m_time,
}

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(
            module=module,
            esm=model,
            scen=scen_str,
            variable=variable,
        )

        param.to_netcdf(filename)

MESMER-M: calibrate multiple scenarios

Contents

MESMER-M: calibrate multiple scenarios#

Configuration#

Load Data#

Preprocessing#

Calculate anomalies#

Extract land gridcells & stack#

Calibration#

Fit the harmonic model#

Approximate the residuals#

Train the power transformer#

Fit cyclo-stationary AR(1) process#

AR parameters#

Find localized empirical covariance#

Saving the parameters#

Create time coordinate#

Store data for all modules#