Setting up your own output

This demo demonstrates how you can setup your raw model output with climpred.preprocessing to match climpred’s expectations.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import climpred

from climpred.preprocessing.shared import load_hindcast, set_integer_time_axis
from climpred.preprocessing.mpi import get_path

Assuming your raw model output is stored in multiple files per member and initialization, load_hindcast is a nice wrapper function based on get_path designed for the output format of MPI-ESM to aggregated all hindcast output into one file as expected by climpred.

The basic idea is to look over the output of all members and concatinate, then loop over all initializations and concatinate. Before concatination, it is important to make the time dimension identical in all input datasets for concatination.

To reduce the data size, use the preprocess function provided to xr.open_mfdataset wisely in combination with set_integer_axis, e.g. additionally extracting only a certain region, time-step, time-aggregation or only few variables for a multi-variable input file as in MPI-ESM standard output.

# check the code of load_hindcast
# load_hindcast??

v = "global_primary_production"

def preprocess_1var(ds, v=v):
    """Only leave one variable `v` in dataset """
    return ds[v].to_dataset(name=v).squeeze()

# lead_offset because yearmean output
%time ds = load_hindcast(inits=range(1961, 1965), members=range(1, 3), preprocess=preprocess_1var, get_path=get_path)

Processing init 1961 ...
Processing init 1962 ...
Processing init 1963 ...
Processing init 1964 ...
CPU times: user 5.07 s, sys: 2.06 s, total: 7.13 s
Wall time: 5.19 s

# what we need for climpred
ds.coords

Coordinates:
    depth    float64 0.0
    lat      float64 0.0
    lon      float64 0.0
  * lead     (lead) int64 1 2 3 4 5 6 7 8 9 10
  * member   (member) int64 1 2
  * init     (init) int64 1961 1962 1963 1964

ds[v].data

Array Chunk Bytes 320 B 4 B Shape (4, 2, 10) (1, 1, 1) Count 720 Tasks 80 Chunks Type float32 numpy.ndarray

# loading the data into memory
# if not rechunk
%time ds = ds.load()

CPU times: user 216 ms, sys: 44 ms, total: 260 ms
Wall time: 220 ms

# go on with creation of PredictionEnsemble
# climpred.HindcastEnsemble(ds).add_observations(obs).verify(metric='acc', comparison='e2o', dim='init', alignment='maximize')
# climpred.PerfectModelEnsemble(ds).add_control(control).verify(metric='acc', comparison='m2e', dim=['init','member'])

intake-esm for cmorized output

In case you have access to cmorized output of CMIP experiments, consider using intake-esm. With the preprocess function you can align the time dimension of all input files. Finally, rename_to_climpred_dims only renames.

from climpred.preprocessing.shared import rename_to_climpred_dims, set_integer_time_axis

# make to have intake-esm installed
import intake # this is enough for intake-esm to work

col_url = "/home/mpim/m300524/intake-esm-datastore/catalogs/mistral-cmip6.json"
col_url = "https://raw.githubusercontent.com/NCAR/intake-esm-datastore/master/catalogs/pangeo-cmip6.json"
col = intake.open_esm_datastore(col_url)

col.df.columns

Index(['activity_id', 'institution_id', 'source_id', 'experiment_id',
       'member_id', 'table_id', 'variable_id', 'grid_label', 'dcpp_init_year',
       'version', 'time_range', 'path'],
      dtype='object')

# load 2 members for 2 inits for one variable from one model
query = dict(experiment_id=[
    'dcppA-hindcast'], table_id='Amon', member_id=['r1i1p1f1', 'r2i1p1f1'], dcpp_init_year=[1970, 1971],
    variable_id='tas', source_id='MPI-ESM1-2-HR')
cat = col.search(**query)
cdf_kwargs = {'chunks': {'time': 12}, 'decode_times': False}

cat.df.head()

	activity_id	institution_id	source_id	experiment_id	member_id	table_id	variable_id	grid_label	dcpp_init_year	version	time_range	path
0	DCPP	MPI-M	MPI-ESM1-2-HR	dcppA-hindcast	r1i1p1f1	Amon	tas	gn	1971.0	v20190906	197111-198112	/work/ik1017/CMIP6/data/CMIP6/DCPP/MPI-M/MPI-E...
1	DCPP	MPI-M	MPI-ESM1-2-HR	dcppA-hindcast	r1i1p1f1	Amon	tas	gn	1970.0	v20190906	197011-198012	/work/ik1017/CMIP6/data/CMIP6/DCPP/MPI-M/MPI-E...
2	DCPP	MPI-M	MPI-ESM1-2-HR	dcppA-hindcast	r2i1p1f1	Amon	tas	gn	1971.0	v20190906	197111-198112	/work/ik1017/CMIP6/data/CMIP6/DCPP/MPI-M/MPI-E...
3	DCPP	MPI-M	MPI-ESM1-2-HR	dcppA-hindcast	r2i1p1f1	Amon	tas	gn	1970.0	v20190906	197011-198012	/work/ik1017/CMIP6/data/CMIP6/DCPP/MPI-M/MPI-E...

def preprocess(ds):
    # extract tiny spatial and temporal subset to make this fast
    ds = ds.isel(lon=[50, 51, 52], lat=[50, 51, 52],
                 time=np.arange(12 * 2))
    # make time dim identical
    ds = set_integer_time_axis(ds,time_dim='time')
    return ds

dset_dict = cat.to_dataset_dict(
    cdf_kwargs=cdf_kwargs, preprocess=preprocess)

Progress: |███████████████████████████████████████████████████████████████████████████████| 100.0%

--> The keys in the returned dictionary of datasets are constructed as follows:
        'activity_id.institution_id.source_id.experiment_id.table_id.grid_label'

--> There are 1 group(s)

# get first dict value
_, ds = dset_dict.popitem()
ds.coords

Coordinates:
    height          float64 ...
  * dcpp_init_year  (dcpp_init_year) float64 1.97e+03 1.971e+03
  * lat             (lat) float64 -42.55 -41.61 -40.68
  * time            (time) int64 1 2 3 4 5 6 7 8 9 ... 17 18 19 20 21 22 23 24
  * lon             (lon) float64 46.88 47.81 48.75
  * member_id       (member_id) <U8 'r1i1p1f1' 'r2i1p1f1'

# rename to comply with climpred's required dimension names
ds = rename_to_climpred_dims(ds)

# what we need for climpred
ds.coords

Coordinates:
    height   float64 ...
  * init     (init) float64 1.97e+03 1.971e+03
  * lat      (lat) float64 -42.55 -41.61 -40.68
  * lead     (lead) int64 1 2 3 4 5 6 7 8 9 10 ... 15 16 17 18 19 20 21 22 23 24
  * lon      (lon) float64 46.88 47.81 48.75
  * member   (member) <U8 'r1i1p1f1' 'r2i1p1f1'

ds['tas'].data

Array Chunk Bytes 3.46 kB 432 B Shape (2, 2, 24, 3, 3) (1, 1, 12, 3, 3) Count 176 Tasks 8 Chunks Type float32 numpy.ndarray

# loading the data into memory
# if not rechunk
# this is here quite fast before we only select 9 grid cells
%time ds = ds.load()

CPU times: user 218 ms, sys: 74 ms, total: 292 ms
Wall time: 237 ms

# go on with creation of PredictionEnsemble
# climred.HindcastEnsemble(ds).add_observations(obs).verify(metric='acc', comparison='e2o', dim='init', alignment='maximize')
# climred.PerfectModelEnsemble(ds).add_control(control).verify(metric='acc', comparison='m2e', dim=['init','member'])