global_avg_temp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ba8ee18a-2ea7-4cc8-ab6b-bf8a9367f1de",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_2353/3090431912.py:6: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)\n",
      "  from tqdm.autonotebook import tqdm\n"
     ]
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import dask\n",
    "from dask.diagnostics import progress\n",
    "from tqdm.autonotebook import tqdm \n",
    "import intake\n",
    "import fsspec\n",
    "import seaborn as sns\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6d9c342d-f7d8-4c38-abaa-633a2d086a49",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>experiment_id</th>\n",
       "      <th>variable_id</th>\n",
       "      <th>table_id</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>source_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ACCESS-CM2</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AWI-CM-1-1-MR</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BCC-CSM2-MR</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CAMS-CSM1-0</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CESM2-WACCM</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CMCC-CM2-SR5</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CMCC-ESM2</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CanESM5</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>EC-Earth3</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>EC-Earth3-Veg</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>EC-Earth3-Veg-LR</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FGOALS-f3-L</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>FGOALS-g3</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GFDL-ESM4</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IITM-ESM</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>INM-CM4-8</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>INM-CM5-0</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IPSL-CM6A-LR</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>KACE-1-0-G</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MIROC6</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MPI-ESM1-2-HR</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MPI-ESM1-2-LR</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MRI-ESM2-0</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NorESM2-MM</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TaiESM1</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  experiment_id  variable_id  table_id\n",
       "source_id                                             \n",
       "ACCESS-CM2                    6            4         1\n",
       "AWI-CM-1-1-MR                 6            4         1\n",
       "BCC-CSM2-MR                   6            4         1\n",
       "CAMS-CSM1-0                   6            4         1\n",
       "CESM2-WACCM                   6            4         1\n",
       "CMCC-CM2-SR5                  6            4         1\n",
       "CMCC-ESM2                     6            4         1\n",
       "CanESM5                       6            4         1\n",
       "EC-Earth3                     6            4         1\n",
       "EC-Earth3-Veg                 6            4         1\n",
       "EC-Earth3-Veg-LR              6            4         1\n",
       "FGOALS-f3-L                   6            4         1\n",
       "FGOALS-g3                     6            4         1\n",
       "GFDL-ESM4                     6            4         1\n",
       "IITM-ESM                      6            4         1\n",
       "INM-CM4-8                     6            4         1\n",
       "INM-CM5-0                     6            4         1\n",
       "IPSL-CM6A-LR                  6            4         1\n",
       "KACE-1-0-G                    6            4         1\n",
       "MIROC6                        6            4         1\n",
       "MPI-ESM1-2-HR                 6            4         1\n",
       "MPI-ESM1-2-LR                 6            4         1\n",
       "MRI-ESM2-0                    6            4         1\n",
       "NorESM2-MM                    6            4         1\n",
       "TaiESM1                       6            4         1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col = intake.open_esm_datastore(\"https://storage.googleapis.com/cmip6/pangeo-cmip6.json\")\n",
    "col\n",
    "\n",
    "# there is currently a significant amount of data for these runs\n",
    "expts_full = ['historical', 'ssp126', 'ssp245','ssp370', 'ssp585', 'piControl']\n",
    "\n",
    "query = dict(\n",
    "    experiment_id=expts_full,\n",
    "    table_id='Amon',                           \n",
    "    variable_id=['tas', 'pr', 'ua', 'va'],\n",
    "    member_id = 'r1i1p1f1',                     \n",
    ")\n",
    "\n",
    "col_subset = col.search(require_all_on=[\"source_id\"], **query)\n",
    "col_subset_var = [col_subset.search(variable_id=var_name) for var_name in query['variable_id']]\n",
    "col_subset.df.groupby(\"source_id\")[\n",
    "    [\"experiment_id\", \"variable_id\", \"table_id\"]\n",
    "].nunique()\n",
    "#col_subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e38909b-f5e7-4193-9a64-8f233850ce2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def drop_all_bounds(ds):\n",
    "    drop_vars = [vname for vname in ds.coords\n",
    "                 if (('_bounds') in vname ) or ('_bnds') in vname]\n",
    "    return ds.drop(drop_vars)\n",
    "\n",
    "def open_dset(df):\n",
    "    #assert len(df) == 1\n",
    "    ds = xr.open_zarr(fsspec.get_mapper(df.zstore.values[0]), consolidated=True, decode_times=True, use_cftime=True)\n",
    "    return drop_all_bounds(ds)\n",
    "\n",
    "def open_delayed(df):\n",
    "    return dask.delayed(open_dset)(df)\n",
    "\n",
    "from collections import defaultdict\n",
    "dsets = []\n",
    "for col_subset in col_subset_var :\n",
    "    dset = defaultdict(dict)\n",
    "\n",
    "    for group, df in col_subset.df.groupby(by=['source_id', 'experiment_id']):\n",
    "        dset[group[0]][group[1]] = open_delayed(df)\n",
    "    dsets.append(dset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33b5aee3-f9ee-41a1-a7a7-8ffbf4b7d82c",
   "metadata": {},
   "outputs": [],
   "source": [
    "dsets_ = [dask.compute(dict(dset))[0]for dset in dsets]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "283f67b0-4533-4bdd-b767-95dce78992c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate global means\n",
    "\n",
    "def get_lat_name(ds):\n",
    "    for lat_name in ['lat', 'latitude']:\n",
    "        if lat_name in ds.coords:\n",
    "            return lat_name\n",
    "    raise RuntimeError(\"Couldn't find a latitude coordinate\")\n",
    "\n",
    "def global_mean(ds):\n",
    "    lat = ds[get_lat_name(ds)]\n",
    "    weight = np.cos(np.deg2rad(lat))\n",
    "    weight /= weight.mean()\n",
    "    other_dims = set(ds.dims) - {'time'}\n",
    "    return (ds * weight).mean(other_dims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17b9128f-1e0d-469b-a76c-86be8c12d87b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from toolz.functoolz import juxt\n",
    "expts = expts_full[:-1]\n",
    "expt_da = xr.DataArray(expts, dims='experiment_id', name='experiment_id',\n",
    "                       coords={'experiment_id': expts})\n",
    "\n",
    "dsets_aligned_list = []\n",
    "for dset_ in dsets_[:2]:\n",
    "    j=0\n",
    "    dsets_aligned = {}\n",
    "    for k, v in tqdm(dset_.items()):\n",
    "        \n",
    "        expt_dsets = v.values()\n",
    "        if any([d is None for d in expt_dsets]):\n",
    "            print(f\"Missing experiment for {k}\")\n",
    "            continue\n",
    "\n",
    "        for ds in expt_dsets:\n",
    "            ds.coords['year'] = ds.time.dt.year\n",
    "    \n",
    "        # workaround for\n",
    "        # https://github.com/pydata/xarray/issues/2237#issuecomment-620961663\n",
    "        dsets_ann_mean = [v[expt].pipe(global_mean)\n",
    "                                .swap_dims({'time': 'year'})\n",
    "                                .drop('time')\n",
    "                                .coarsen(year=12).mean()\n",
    "                                for expt in expts]\n",
    "\n",
    "        # align everything with the 4xCO2 experiment\n",
    "\n",
    "        dsets_aligned[k] = xr.concat(dsets_ann_mean, join='outer',\n",
    "                                    dim=expt_da)\n",
    "    dsets_aligned_list.append(dsets_aligned)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f724ed46-1d72-4d3a-84d4-476c4fe34bb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "with progress.ProgressBar():\n",
    "    dsets_aligned_list_1 = dask.compute(dsets_aligned_list[0])[0]\n",
    "        \n",
    "with progress.ProgressBar():\n",
    "    dsets_aligned_list_2 = dask.compute(dsets_aligned_list[1])[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f89cc47-df8a-45c6-8770-431c93ac0f9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "dsets_algned_list_ = [dsets_aligned_list_1, dsets_aligned_list_2]\n",
    "type(dsets_algned_list_[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1233280-e3d5-41f4-ab41-90daba2f5c15",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "source_ids = [list(dsets_aligned_.keys()) for dsets_aligned_ in dsets_algned_list_]\n",
    "#source_da = xr.DataArray(source_ids, dims='source_id', name='source_id',\n",
    "#                         coords={'source_id': source_ids})\n",
    "big_ds = []\n",
    "for idx, dsets_aligned_ in enumerate(dsets_algned_list_):\n",
    "    source_da = xr.DataArray(source_ids[idx], dims='source_id', name='source_id',\n",
    "                         coords={'source_id': source_ids[idx]})\n",
    "    big_ds.append(xr.concat([ds.reset_coords(drop=True)\n",
    "                        for ds in dsets_aligned_.values()],\n",
    "                        dim=source_da))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58c01d85-64a0-487c-ae36-a9b282c7166c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_all = big_ds[0].assign(pr=big_ds[1].pr)\n",
    "ds_all.to_netcdf('/home/jovyan/pangeo/data/tas_pr_global_timeseries.nc')\n",
    "ds_all"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60b48db9-334e-4736-b8c3-618cf0728a8e",
   "metadata": {},
   "source": [
    "# Processing and Analysing with Pandas - skipped\n",
    "\n",
    "This is a bit more dfficult since the dataset has multiple coordinates to deal with.\n",
    "\n",
    "### 1. Converting dataset to dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2877097d-2c1a-4045-b417-374ec2ecb1f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all = big_ds.sel(year=slice(1950, 2100)).to_dataframe().reset_index()\n",
    "df_all.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1569baf1-4183-47c3-827a-f4d4606fb90a",
   "metadata": {},
   "source": [
    "### Simple plotting using seaborn relplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbdb0c9f-12a1-42c8-8a4e-a4a78be48db4",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.relplot(data=df_all,\n",
    "            x=\"year\", y=\"pr\", hue='experiment_id',\n",
    "            kind=\"line\", ci=\"sd\", aspect=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7af80310-b641-4080-b731-cdfc9a7ea6dc",
   "metadata": {},
   "source": [
    "### 2. Finding climatic mean value for each models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03fa1e7e-4362-4ee7-afb6-df9861e86808",
   "metadata": {},
   "outputs": [],
   "source": [
    "clim_mean = df_all[(df_all['year']>1980) & (df_all['year']<2014)].groupby('source_id')['pr'].apply(np.nanmean).reset_index()\n",
    "clim_mean.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "319e5fe7-d734-43d0-95b1-7cf2fbe77469",
   "metadata": {},
   "source": [
    "### 3. Calculating the fractional deviation from climatic mean "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb1e7548-d0f9-43d6-b735-416cf772c5f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "def func(df):\n",
    "    df_m = df_all\n",
    "    for group in df.groups.keys():\n",
    "        df_m['pr'].iloc[df.groups[group]] = ((df_m['pr'].iloc[df.groups[group]].values - clim_mean['pr'].values) / clim_mean['pr'].values)*100\n",
    "    return df_m\n",
    "dd = df_all.groupby(['year', 'experiment_id'])[['pr']].pipe(func)\n",
    "dd.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88943a91-65ad-4d62-83f2-332c41aaf25e",
   "metadata": {},
   "source": [
    "### 4. Finding the ensemble mean and plotting with seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7afa561-c538-4593-9fa5-8647fad3286b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ens_mean = dd.groupby(['experiment_id', 'year'])['pr'].mean().reset_index()\n",
    "ax = sns.relplot(data=ens_mean, x='year', y='pr', hue='experiment_id', kind='line')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfabf75e-8e3b-4091-861b-daa712d5f695",
   "metadata": {},
   "source": [
    "# Processing and Analysing with Xarray - Global Mean\n",
    "\n",
    "Xarray is fluid and easy in terms of dealing with higher number of coordinates.\n",
    "\n",
    "## 1. Subsetting the data from big to small"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "89e7ec04-0d5c-40a3-8fe2-a96e5b6c2213",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_all = xr.open_dataset('/home/jovyan/pangeo/data/tas_pr_global_timeseries.nc')\n",
    "small_ds = ds_all.sel(year=slice(1950, 2100)).rolling(year=2).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "806303b5-7156-41bf-a767-1c9d8b98f0de",
   "metadata": {},
   "source": [
    "## 2. Finding Climatic Mean and calculating the fractional change\n",
    "This is because the values are so small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2854f486-e17c-4831-8c60-f45d69590385",
   "metadata": {},
   "outputs": [],
   "source": [
    "cl_t0, cl_tf = 1984, 2014\n",
    "clim_ds = small_ds.sel(experiment_id='historical', year=slice(cl_t0, cl_tf)).mean(dim=('year')) \n",
    "\n",
    "dev_ds_pr = ((small_ds.pr - clim_ds.pr)/clim_ds.pr)*100 \n",
    "dev_ds_tas = small_ds.tas - clim_ds.tas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5188e405-2acc-46ae-8fec-7b8cc9ffe1f2",
   "metadata": {},
   "source": [
    "## 3. Ensemble mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3d1e073-6b7a-45d0-9865-2354b7b9f66d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ens_ds_pr = dev_ds_pr.mean(dim='source_id')\n",
    "print(ens_ds_pr.shape)\n",
    "\n",
    "ens_ds_tas = dev_ds_tas.mean(dim='source_id')\n",
    "ens_ds_tas.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3efe8f5-ead6-4f79-b383-d041ec391e10",
   "metadata": {},
   "source": [
    "## 4. Finding the actual spread (not standard deviation) w.r.t ensemble mean for each models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8df049ba-0e73-493e-aebc-d5bfd361b20b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools as it\n",
    "spread_max_pr = np.empty(shape=(5, 151))\n",
    "spread_min_pr = np.empty(shape=(5, 151))\n",
    "\n",
    "spread_max_tas = np.empty(shape=(5, 151))\n",
    "spread_min_tas = np.empty(shape=(5, 151))\n",
    "\n",
    "for i,j in it.product(range(5), range(151)):  ##instead of nested for loops. works same.\n",
    "    spread_max_pr[i, j] = np.nanmax(dev_ds_pr[:, i, j]-ens_ds_pr[i, j])\n",
    "    spread_min_pr[i, j] = np.nanmin(dev_ds_pr[:, i, j]-ens_ds_pr[i, j])\n",
    "    \n",
    "    spread_max_tas[i, j] = np.nanmax(dev_ds_tas[:, i, j]-ens_ds_tas[i, j])\n",
    "    spread_min_tas[i, j] = np.nanmin(dev_ds_tas[:, i, j]-ens_ds_tas[i, j])\n",
    "\n",
    "spread_tas ={'max_vals' : spread_max_tas, \n",
    "             'min_vals' : spread_min_tas}\n",
    "\n",
    "spread_pr ={'max_vals' : spread_max_pr, \n",
    "             'min_vals' : spread_min_pr}\n",
    "\n",
    "spr = [spread_tas, spread_pr]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0f121dc-2a1f-49e8-b584-7ec9180e1f51",
   "metadata": {},
   "source": [
    "## 5. Visualizing the ensemble mean with spreads for the historical and future projection data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f32da74-958b-4a28-957e-6da526074b9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import cm\n",
    "import matplotlib.colors as mcl\n",
    "cmap = cm.turbo\n",
    "cl = ['k']+[mcl.rgb2hex(cmap(i)[:3]) for i in range(0,cmap.N,70)]\n",
    "cl = ['k'] + ['blue', 'orange', 'green','red']\n",
    "\n",
    "ncols, nrows=2,2\n",
    "fig, axs = plt.subplots(nrows, ncols, dpi=600, figsize = (14,9))\n",
    "ens_dss = [ens_ds_tas, ens_ds_pr]\n",
    "\n",
    "y_l = [f'Relative to {cl_t0}-{cl_tf} ($\\circ$C)', \n",
    "       f'Relative to {cl_t0}-{cl_tf} (%)']\n",
    "\n",
    "for i in range(nrows):\n",
    "    for j in range(ncols):\n",
    "        ax = axs[i, j]\n",
    "        ens_ds = ens_dss[j]\n",
    "        if i==0:\n",
    "            for idx,v in enumerate(ens_ds.experiment_id):\n",
    "                len(ens_ds.year)\n",
    "                ax.plot(ens_ds.year, ens_ds[ens_ds['experiment_id']==v.values].squeeze(), label = v.values, c=cl[idx])\n",
    "                ax.fill_between(ens_ds.year, ens_ds[ens_ds['experiment_id']==v.values].squeeze()+spr[j]['max_vals'][idx,:], \n",
    "                                 ens_ds[ens_ds['experiment_id']==v.values].squeeze()+spr[j]['min_vals'][idx, :], alpha=0.2, color=cl[idx] )\n",
    "                ax.set_xlim(1955,2100)\n",
    "        else:\n",
    "            for idx,v in enumerate(ens_ds.experiment_id[1:]):\n",
    "                ax.plot(ens_ds.year, ens_ds[ens_ds['experiment_id']==v.values].squeeze(), label=v.values, c=cl[idx+1])\n",
    "                ax.set_xlim(2020,2100)\n",
    "        ax.set_xlabel('Year')\n",
    "        ax.set_ylabel(y_l[j])\n",
    "        ax.legend(loc='upper left')\n",
    "        ax.grid(alpha=0.3)\n",
    "\n",
    "fig.suptitle('CMIP6 Global mean T$_s$ and P', x=0.5, y =0.95, fontsize=28, weight='bold')\n",
    "plt.savefig('/home/jovyan/pangeo/plot/tas_pr_global_timeseries.png', bbox_inches='tight', facecolor='white')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b492cede-c12e-4bfc-a61c-3abf5a8b4692",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}