{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../../thumbnail.png\" width=250 alt=\"CESM LENS image\"></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reproducing Key Figures from Kay et al. (2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "This notebook demonstrates how one might use the NCAR Community Earth System Model (CESM) Large Ensemble (LENS) data hosted on AWS S3. The notebook shows how to reproduce figures 2 and 4 from the [Kay et al. (2015) paper](https://doi.org/10.1175/BAMS-D-13-00255.1) describing the CESM LENS dataset.\n",
    "\n",
    "This resource is intended to be helpful for people not familiar with elements of the [Pangeo](https://pangeo.io/) framework including Jupyter Notebooks, [Xarray](https://docs.xarray.dev/en/stable/), and [Zarr](https://zarr.readthedocs.io/en/stable/) data format, or with the original paper, so it includes additional explanation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "\n",
    "| Concepts | Importance | Notes |\n",
    "| --- | --- | --- |\n",
    "| [Intro to Xarray](https://foundations.projectpythia.org/core/xarray/xarray-intro.html) | Necessary | |\n",
    "| Dask | Helpful | |\n",
    "\n",
    "- **Time to learn**: 30 minutes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\"><b>NOTE: </b>In this notebook, we access very large cloud-served datasets and use <a href=\"https://dask.org\">Dask</a> to parallelize our workflow. The end-to-end execution time may be on the order of an hour or more, depending on your computing resources.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import intake\n",
    "import matplotlib.pyplot as plt\n",
    "from dask.distributed import Client\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import cmaps  # for NCL colormaps\n",
    "import cartopy.crs as ccrs\n",
    "import dask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dask.config.set({\"distributed.scheduler.worker-saturation\": 1.0})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Create and Connect to Dask Distributed Cluster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we'll use a dask cluster to parallelize our analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "platform = sys.platform\n",
    "\n",
    "if (platform == 'win32'):\n",
    "    import multiprocessing.popen_spawn_win32\n",
    "else:\n",
    "    import multiprocessing.popen_spawn_posix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client = Client()\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Load and Prepare Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_url = 'https://ncar-cesm-lens.s3-us-west-2.amazonaws.com/catalogs/aws-cesm1-le.json'\n",
    "col = intake.open_esm_datastore(catalog_url)\n",
    "col"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the first few lines of the catalog:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col.df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show expanded version of collection structure with details:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col.keys_info().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract data needed to construct Figure 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Search the catalog to find the desired data, in this case the reference height temperature of the atmosphere, at daily time resolution, for the Historical, 20th Century, and RCP8.5 (IPCC Representative Concentration Pathway 8.5) experiments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_subset = col.search(frequency=[\"daily\", \"monthly\"], component=\"atm\", variable=\"TREFHT\",\n",
    "                        experiment=[\"20C\", \"RCP85\", \"HIST\"])\n",
    "\n",
    "col_subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_subset.df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load catalog entries for subset into a dictionary of Xarray Datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsets = col_subset.to_dataset_dict(zarr_kwargs={\"consolidated\": True}, storage_options={\"anon\": True})\n",
    "print(f\"\\nDataset dictionary keys:\\n {dsets.keys()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define Xarray Datasets corresponding to the three experiments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_HIST = dsets['atm.HIST.monthly']\n",
    "ds_20C = dsets['atm.20C.daily']\n",
    "ds_RCP85 = dsets['atm.RCP85.daily']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `dask.distributed` utility function to display size of each dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.utils import format_bytes\n",
    "print(f\"Historical: {format_bytes(ds_HIST.nbytes)}\\n\"\n",
    "      f\"20th Century: {format_bytes(ds_20C.nbytes)}\\n\"\n",
    "      f\"RCP8.5: {format_bytes(ds_RCP85.nbytes)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, extract the Reference Height Temperature data variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_hist = ds_HIST[\"TREFHT\"]\n",
    "t_20c = ds_20C[\"TREFHT\"]\n",
    "t_rcp = ds_RCP85[\"TREFHT\"]\n",
    "t_20c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The global surface temperature anomaly was computed relative to the 1961-90 base period in the Kay et al. paper, so extract that time slice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_ref = t_20c.sel(time=slice(\"1961\", \"1990\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Figure 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read grid cell areas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cell size varies with latitude, so this must be accounted for when computing the global mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat = col.search(frequency=\"static\", component=\"atm\", experiment=[\"20C\"])\n",
    "_, grid = cat.to_dataset_dict(aggregate=False, storage_options={'anon':True}, zarr_kwargs={\"consolidated\": True}).popitem()\n",
    "grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_area = grid.area.load()\n",
    "total_area = cell_area.sum()\n",
    "cell_area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define weighted means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: `resample(time=\"AS\")` does an annual resampling based on start of calendar year. See documentation for [Pandas resampling options](https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#dateoffset-objects)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_ref_ts = (\n",
    "    (t_ref.resample(time=\"AS\").mean(\"time\") * cell_area).sum(dim=(\"lat\", \"lon\"))\n",
    "    / total_area\n",
    ").mean(dim=(\"time\", \"member_id\"))\n",
    "\n",
    "t_hist_ts = (\n",
    "    (t_hist.resample(time=\"AS\").mean(\"time\") * cell_area).sum(dim=(\"lat\", \"lon\"))\n",
    ") / total_area\n",
    "\n",
    "t_20c_ts = (\n",
    "    (t_20c.resample(time=\"AS\").mean(\"time\") * cell_area).sum(dim=(\"lat\", \"lon\"))\n",
    ") / total_area\n",
    "\n",
    "t_rcp_ts = (\n",
    "    (t_rcp.resample(time=\"AS\").mean(\"time\") * cell_area).sum(dim=(\"lat\", \"lon\"))\n",
    ") / total_area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read data and compute means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dask’s \"lazy execution\" philosophy means that until this point we have not actually read the bulk of the data. Steps 1, 3, and 4 take a while to complete, so we include the Notebook \"cell magic\" directive `%%time` to display elapsed and CPU times after computation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 1 (takes a while)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# this cell takes a while, be patient\n",
    "t_ref_mean = t_ref_ts.load()\n",
    "t_ref_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 2 (executes quickly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time \n",
    "t_hist_ts_df = t_hist_ts.to_series().T\n",
    "#t_hist_ts_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 3 (takes even longer than Step 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "t_20c_ts_df = t_20c_ts.to_series().unstack().T\n",
    "t_20c_ts_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 4 (similar to Step 3 in its execution time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# This also takes a while\n",
    "t_rcp_ts_df = t_rcp_ts.to_series().unstack().T\n",
    "t_rcp_ts_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get observations for Figure 2 (HadCRUT4; Morice et al. 2012)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observational time series data for comparison with ensemble average:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obsDataURL = \"https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/cru/hadcrut4/air.mon.anom.median.nc\"\n",
    "ds = xr.open_dataset(obsDataURL).load()\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def weighted_temporal_mean(ds):\n",
    "    \"\"\"\n",
    "    weight by days in each month\n",
    "    \"\"\"\n",
    "    time_bound_diff = ds.time_bnds.diff(dim=\"nbnds\")[:, 0]\n",
    "    wgts = time_bound_diff.groupby(\"time.year\") / time_bound_diff.groupby(\n",
    "        \"time.year\"\n",
    "    ).sum(xr.ALL_DIMS)\n",
    "    obs = ds[\"air\"]\n",
    "    cond = obs.isnull()\n",
    "    ones = xr.where(cond, 0.0, 1.0)\n",
    "    obs_sum = (obs * wgts).resample(time=\"AS\").sum(dim=\"time\")\n",
    "    ones_out = (ones * wgts).resample(time=\"AS\").sum(dim=\"time\")\n",
    "    obs_s = (obs_sum / ones_out).mean((\"lat\", \"lon\")).to_series()\n",
    "    return obs_s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Limit observations to 20th century:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_s = weighted_temporal_mean(ds)\n",
    "obs_s = obs_s['1920':]\n",
    "obs_s.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_ts_anom = pd.concat([t_20c_ts_df, t_rcp_ts_df]) - t_ref_mean.data\n",
    "years = [val.year for val in all_ts_anom.index]\n",
    "obs_years = [val.year for val in obs_s.index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine ensemble member 1 data from historical and 20th century experiments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_anom = t_hist_ts_df - t_ref_mean.data\n",
    "member1 = pd.concat([hist_anom.iloc[:-2], all_ts_anom.iloc[:,0]], verify_integrity=True)\n",
    "member1_years = [val.year for val in member1.index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Figure 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Global surface temperature anomaly (1961-90 base period) for individual ensemble members, and observations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.axes()\n",
    "\n",
    "ax.tick_params(right=True, top=True, direction=\"out\", length=6, width=2, grid_alpha=0.5)\n",
    "ax.plot(years, all_ts_anom.iloc[:,1:], color=\"grey\")\n",
    "ax.plot(obs_years, obs_s['1920':], color=\"red\")\n",
    "ax.plot(member1_years, member1, color=\"black\")\n",
    "\n",
    "\n",
    "ax.text(\n",
    "    0.35,\n",
    "    0.4,\n",
    "    \"observations\",\n",
    "    verticalalignment=\"bottom\",\n",
    "    horizontalalignment=\"left\",\n",
    "    transform=ax.transAxes,\n",
    "    color=\"red\",\n",
    "    fontsize=10,\n",
    ")\n",
    "ax.text(\n",
    "    0.35,\n",
    "    0.33,\n",
    "    \"members 2-40\",\n",
    "    verticalalignment=\"bottom\",\n",
    "    horizontalalignment=\"left\",\n",
    "    transform=ax.transAxes,\n",
    "    color=\"grey\",\n",
    "    fontsize=10,\n",
    ")\n",
    "ax.text(\n",
    "    0.05,\n",
    "    0.2,\n",
    "    \"member 1\",\n",
    "    verticalalignment=\"bottom\",\n",
    "    horizontalalignment=\"left\",\n",
    "    transform=ax.transAxes,\n",
    "    color=\"black\",\n",
    "    fontsize=10,\n",
    ")\n",
    "\n",
    "ax.set_xticks([1850, 1920, 1950, 2000, 2050, 2100])\n",
    "plt.ylim(-1, 5)\n",
    "plt.xlim(1850, 2100)\n",
    "plt.ylabel(\"Global Surface\\nTemperature Anomaly (K)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Figure 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute linear trend for winter seasons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_trend(da, dim=\"time\"):\n",
    "    da_chunk = da.chunk({dim: -1})\n",
    "    trend = xr.apply_ufunc(\n",
    "        calc_slope,\n",
    "        da_chunk,\n",
    "        vectorize=True,\n",
    "        input_core_dims=[[dim]],\n",
    "        output_core_dims=[[]],\n",
    "        output_dtypes=[np.float64],\n",
    "        dask=\"parallelized\",\n",
    "    )\n",
    "    return trend\n",
    "\n",
    "\n",
    "def calc_slope(y):\n",
    "    \"\"\"ufunc to be used by linear_trend\"\"\"\n",
    "    x = np.arange(len(y))\n",
    "\n",
    "    # drop missing values (NaNs) from x and y\n",
    "    finite_indexes = ~np.isnan(y)\n",
    "    slope = np.nan if (np.sum(finite_indexes) < 2) else np.polyfit(x[finite_indexes], y[finite_indexes], 1)[0]\n",
    "    return slope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute ensemble trends"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time \n",
    "# Takes several minutes\n",
    "t = xr.concat([t_20c, t_rcp], dim=\"time\")\n",
    "seasons = t.sel(time=slice(\"1979\", \"2012\")).resample(time=\"QS-DEC\").mean(\"time\")\n",
    "# Include only full seasons from 1979 and 2012\n",
    "seasons = seasons.sel(time=slice(\"1979\", \"2012\")).load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "winter_seasons = seasons.sel(\n",
    "    time=seasons.time.where(seasons.time.dt.month == 12, drop=True)\n",
    ")\n",
    "winter_trends = linear_trend(\n",
    "    winter_seasons.chunk({\"lat\": 20, \"lon\": 20, \"time\": -1})\n",
    ").load() * len(winter_seasons.time)\n",
    "\n",
    "# Compute ensemble mean from the first 30 members\n",
    "winter_trends_mean = winter_trends.isel(member_id=range(30)).mean(dim='member_id')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure that we have 34 seasons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(winter_seasons.time) == 34"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get observations for Figure 4 (NASA GISS GisTemp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is observational time series data for comparison with ensemble average (NASA GISS Surface Temperature Analysis, https://data.giss.nasa.gov/gistemp/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obsDataURL = \"https://data.giss.nasa.gov/pub/gistemp/gistemp1200_GHCNv4_ERSSTv5.nc.gz\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download, unzip, and load file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.system(\"wget \" + obsDataURL)\n",
    "\n",
    "obsDataFileName = obsDataURL.split('/')[-1]\n",
    "os.system(\"gunzip \" + obsDataFileName)\n",
    "\n",
    "obsDataFileName = obsDataFileName[:-3]\n",
    "ds = xr.open_dataset(obsDataFileName).load()\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remap longitude range from [-180, 180] to [0, 360] for plotting purposes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = ds.assign_coords(lon=((ds.lon + 360) % 360)).sortby('lon')\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute observed trends"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Include only full seasons from 1979 through 2012:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_seasons = ds.sel(time=slice(\"1979\", \"2012\")).resample(time=\"QS-DEC\").mean(\"time\")\n",
    "obs_seasons = obs_seasons.sel(time=slice(\"1979\", \"2012\")).load()\n",
    "obs_winter_seasons = obs_seasons.sel(\n",
    "    time=obs_seasons.time.where(obs_seasons.time.dt.month == 12, drop=True)\n",
    ")\n",
    "obs_winter_seasons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And compute observed winter trends:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_winter_trends = linear_trend(\n",
    "    obs_winter_seasons.chunk({\"lat\": 20, \"lon\": 20, \"time\": -1})\n",
    ").load() * len(obs_winter_seasons.time)\n",
    "obs_winter_trends"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Figure 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Global maps of historical (1979 - 2012) boreal winter (DJF) surface air trends:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contour_levels = [-6, -5, -4, -3, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 3, 4, 5, 6]\n",
    "color_map = cmaps.ncl_default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_map_plot(nplot_rows, nplot_cols, plot_index, data, plot_label):\n",
    "    \"\"\" Create a single map subplot. \"\"\"\n",
    "    ax = plt.subplot(nplot_rows, nplot_cols, plot_index, projection = ccrs.Robinson(central_longitude = 180))\n",
    "    cplot = plt.contourf(lons, lats, data,\n",
    "                         levels = contour_levels,\n",
    "                         cmap = color_map,\n",
    "                         extend = 'both',\n",
    "                         transform = ccrs.PlateCarree())\n",
    "    ax.coastlines(color = 'grey')\n",
    "    ax.text(0.01, 0.01, plot_label, fontsize = 14, transform = ax.transAxes)\n",
    "    return cplot, ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Generate plot (may take a while as many individual maps are generated)\n",
    "numPlotRows = 8\n",
    "numPlotCols = 4\n",
    "figWidth = 20\n",
    "figHeight = 30\n",
    "\n",
    "fig, axs = plt.subplots(numPlotRows, numPlotCols, figsize=(figWidth,figHeight))\n",
    "\n",
    "lats = winter_trends.lat\n",
    "lons = winter_trends.lon\n",
    "\n",
    "# Create ensemble member plots\n",
    "for ensemble_index in range(30):\n",
    "    plot_data = winter_trends.isel(member_id = ensemble_index)\n",
    "    plot_index = ensemble_index + 1\n",
    "    plot_label = str(plot_index)\n",
    "    plotRow = ensemble_index // numPlotCols\n",
    "    plotCol = ensemble_index % numPlotCols\n",
    "    # Retain axes objects for figure colorbar\n",
    "    cplot, axs[plotRow, plotCol] = make_map_plot(numPlotRows, numPlotCols, plot_index, plot_data, plot_label)\n",
    "\n",
    "# Create plots for the ensemble mean, observations, and a figure color bar.\n",
    "cplot, axs[7,2] = make_map_plot(numPlotRows, numPlotCols, 31, winter_trends_mean, 'EM')\n",
    "\n",
    "lats = obs_winter_trends.lat\n",
    "lons = obs_winter_trends.lon\n",
    "cplot, axs[7,3] = make_map_plot(numPlotRows, numPlotCols, 32, obs_winter_trends.tempanomaly, 'OBS')\n",
    "\n",
    "cbar = fig.colorbar(cplot, ax=axs, orientation='horizontal', shrink = 0.7, pad = 0.02)\n",
    "cbar.ax.set_title('1979-2012 DJF surface air temperature trends (K/34 years)', fontsize = 16)\n",
    "cbar.set_ticks(contour_levels)\n",
    "cbar.set_ticklabels(contour_levels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Close our client:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "In this notebook, we used CESM LENS data hosted on AWS to recreate two key figures in the paper that describes the project.\n",
    "\n",
    "### What's next?\n",
    "More example workflows using these datasets may be added in the future."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources and references\n",
    "[Original notebook in the Pangeo Gallery](https://gallery.pangeo.io/repos/NCAR/cesm-lens-aws/notebooks/kay-et-al-2015.v3.html)"
   ]
  }
 ],
 "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.11.0"
  },
  "nbdime-conflicts": {
   "local_diff": [
    {
     "diff": [
      {
       "diff": [
        {
         "key": 0,
         "op": "addrange",
         "valuelist": [
          "Python 3"
         ]
        },
        {
         "key": 0,
         "length": 1,
         "op": "removerange"
        }
       ],
       "key": "display_name",
       "op": "patch"
      }
     ],
     "key": "kernelspec",
     "op": "patch"
    }
   ],
   "remote_diff": [
    {
     "diff": [
      {
       "diff": [
        {
         "key": 0,
         "op": "addrange",
         "valuelist": [
          "Python3"
         ]
        },
        {
         "key": 0,
         "length": 1,
         "op": "removerange"
        }
       ],
       "key": "display_name",
       "op": "patch"
      }
     ],
     "key": "kernelspec",
     "op": "patch"
    }
   ]
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}