notebooks/00_calculate_simulated_global_warming.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating historical global warming in climate model simulations\n",
    "\n",
    "#### Python Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import xesmf as xe\n",
    "import cartopy\n",
    "from tqdm.autonotebook import tqdm  # Fancy progress bars for our loops!\n",
    "import intake\n",
    "# util.py is in the local directory\n",
    "# it contains code that is common across project notebooks\n",
    "# or routines that are too extensive and might otherwise clutter\n",
    "# the notebook design\n",
    "import util\n",
    "\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = 12, 6\n",
    "%config InlineBackend.figure_format = 'retina' "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data organization and processing workflow\n",
    "\n",
    "This notebook uses [`intake-esm`](https://intake-esm.readthedocs.io/en/latest/) to ingest and organize climate model output from the fresh-off-the-supercomputers Phase 6 of the Coupled Model Intercomparison Project (CMIP6).\n",
    "\n",
    "Learn more about the Pangeo mirror of CMIP6 data hosted in Google Storage in [this recent blog post](https://cloud.google.com/blog/products/data-analytics/new-climate-model-data-now-google-public-datasets) or the [official dataset reference page](https://console.cloud.google.com/marketplace/details/noaa-public/cmip6).\n",
    "\n",
    "\n",
    "#### Data catalogs\n",
    "\n",
    "The file `https://storage.googleapis.com/cmip6/cmip6-zarr-consolidated-stores.csv` in google cloud storage contains thousands of lines of metadata, each describing an individual climate model experiment's simulated data.\n",
    "\n",
    "For example, the first line in the csv file contains the precipitation rate (`variable_id = 'pr'`), as a function of latitude, longitude, and time, in an individual climate model experiment with the BCC-ESM1 model (`source_id = 'BCC-ESM1'`) developed by the Beijing Climate Center (`institution_id = 'BCC'`). The model is *forced* by the forcing experiment SSP370 (`experiment_id = 'ssp370'`), which stands for the Shared Socio-Economic Pathway 3 that results in a change in radiative forcing of $\\Delta F = 7.0$ W/m$^{2}$ from pre-industrial to 2100. This simulation was run as part of the `AerChemMIP` activity, which is a spin-off of the CMIP activity that focuses specifically on how aerosol chemistry affects climate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('https://storage.googleapis.com/cmip6/cmip6-zarr-consolidated-stores.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file `pangeo-cmip6.json` describes the structure of the CMIP6 metadata and is formatted so as to be read in by the `intake.open_esm_datastore` method, which categorizes all of the data pointers into a tiered collection. For example, this collection contains the simulated data from 28691 individual experiments, representing 48 different models from 23 different scientific institutions. There are 190 different climate variables (e.g. sea surface temperature, sea ice concentration, atmospheric winds, dissolved organic carbon in the ocean, etc.) available for 29 different forcing experiments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col = intake.open_esm_datastore(\"../catalogs/pangeo-cmip6.json\")\n",
    "col"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we show the various forcing experiments that climate modellers ran in these simulations. A few examples are:\n",
    " - `piControl` which fixes CO2 levels at pre-industrial concentrations of 300 ppm\n",
    " - `historical` which includes the historical evolution of greenhouse concentrations as well as historical volcanic eruptions, changes in solar luminosity, and changes in atmospheric aerosol concentrations (and some other, less impactful forcings).\n",
    " - `abrupt-4xCO2` in which the CO2 concentrations in an pre-industrial control simulation are abrupted quadrupled from 300 ppm to 1200 ppm.\n",
    " - `ssp585`, a `worst-case scenario` in which fossil-fueled development leads to a disastrous increase of $\\Delta F = 8.5$ W/m$^{2}$ in radiative forcing (see description of the various SSP narratives below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['experiment_id'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Narratives of the CMIP6 Shared Socio-economic Pathways\n",
    "\n",
    "See this excellent [in-depth explainer](https://www.carbonbrief.org/explainer-how-shared-socioeconomic-pathways-explore-future-climate-change) of the SSP's by Carbonbrief for more info.\n",
    "\n",
    "#### SSP1: Sustainability – Taking the Green Road (Low challenges to mitigation and adaptation)\n",
    "The world shifts gradually, but pervasively, toward a more sustainable path, emphasizing more inclusive development that respects perceived environmental boundaries. Management of the global commons slowly improves, educational and health investments accelerate the demographic transition, and the emphasis on economic growth shifts toward a broader emphasis on human well-being. Driven by an increasing commitment to achieving development goals, inequality is reduced both across and within countries. Consumption is oriented toward low material growth and lower resource and energy intensity.\n",
    "#### SSP2: Middle of the Road (Medium challenges to mitigation and adaptation)\n",
    "The world follows a path in which social, economic, and technological trends do not shift markedly from historical patterns. Development and income growth proceeds unevenly, with some countries making relatively good progress while others fall short of expectations. Global and national institutions work toward but make slow progress in achieving sustainable development goals. Environmental systems experience degradation, although there are some improvements and overall the intensity of resource and energy use declines. Global population growth is moderate and levels off in the second half of the century. Income inequality persists or improves only slowly and challenges to reducing vulnerability to societal and environmental changes remain.\n",
    "#### SSP3: Regional Rivalry – A Rocky Road (High challenges to mitigation and adaptation)\n",
    "A resurgent nationalism, concerns about competitiveness and security, and regional conflicts push countries to increasingly focus on domestic or, at most, regional issues. Policies shift over time to become increasingly oriented toward national and regional security issues. Countries focus on achieving energy and food security goals within their own regions at the expense of broader-based development. Investments in education and technological development decline. Economic development is slow, consumption is material-intensive, and inequalities persist or worsen over time. Population growth is low in industrialized and high in developing countries. A low international priority for addressing environmental concerns leads to strong environmental degradation in some regions.\n",
    "#### SSP4: Inequality – A Road Divided (Low challenges to mitigation, high challenges to adaptation)\n",
    "Highly unequal investments in human capital, combined with increasing disparities in economic opportunity and political power, lead to increasing inequalities and stratification both across and within countries. Over time, a gap widens between an internationally-connected society that contributes to knowledge- and capital-intensive sectors of the global economy, and a fragmented collection of lower-income, poorly educated societies that work in a labor intensive, low-tech economy. Social cohesion degrades and conflict and unrest become increasingly common. Technology development is high in the high-tech economy and sectors. The globally connected energy sector diversifies, with investments in both carbon-intensive fuels like coal and unconventional oil, but also low-carbon energy sources. Environmental policies focus on local issues around middle and high income areas.\n",
    "#### SSP5: Fossil-fueled Development – Taking the Highway (High challenges to mitigation, low challenges to adaptation)\n",
    "This world places increasing faith in competitive markets, innovation and participatory societies to produce rapid technological progress and development of human capital as the path to sustainable development. Global markets are increasingly integrated. There are also strong investments in health, education, and institutions to enhance human and social capital. At the same time, the push for economic and social development is coupled with the exploitation of abundant fossil fuel resources and the adoption of resource and energy intensive lifestyles around the world. All these factors lead to rapid growth of the global economy, while global population peaks and declines in the 21st century. Local environmental problems like air pollution are successfully managed. There is faith in the ability to effectively manage social and ecological systems, including by geo-engineering if necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of Climate Model Output Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading data\n",
    "\n",
    "`intake-esm` enables loading data directly into an [xarray.DataArray](http://xarray.pydata.org/en/stable/api.html#dataset), a metadata-aware extension of numpy arrays. `xarray` objects leverage [dask](https://dask.org/) to only read data into memory as needed for any specific operation (i.e. lazy evaluation). Think of `xarray` Datasets as ways of conveniently organizing large arrays of floating point numbers (e.g. climate model data) on an n-dimensional discrete grid, with important metadata such as units, variable, names, etc.\n",
    "\n",
    "Note that data on the cloud are in [zarr](https://zarr.readthedocs.io/en/stable/) format, an extension of the metadata-aware format [netcdf](https://www.unidata.ucar.edu/software/netcdf/) commonly used in geosciences.\n",
    "\n",
    "`intake-esm` has rules for aggegating datasets; these rules are defined in the collection-specification file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat = col.search(experiment_id='historical',  # pick the `historical` forcing experiment\n",
    "                 table_id='Amon',             # choose to look at atmospheric variables (A) saved at monthly resolution (mon)\n",
    "                 variable_id='tas',           # choose to look at near-surface air temperature (tas) as our variable\n",
    "                 member_id = 'r1i1p1f1')      # arbitrarily pick one realization for each model (i.e. just one set of initial conditions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_slice = slice('1850','2019') # specific years that bracket our period of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert data catalog into a dictionary of xarray datasets\n",
    "dset_dict = cat.to_dataset_dict(zarr_kwargs={'consolidated': True, 'decode_times': False})\n",
    "\n",
    "ds_dict = {}\n",
    "gmst_dict = {}\n",
    "for name, ds in tqdm(dset_dict.items()):\n",
    "    # rename spatial dimensions if necessary\n",
    "    if ('longitude' in ds.dims) and ('latitude' in ds.dims):\n",
    "        ds = ds.rename({'longitude':'lon', 'latitude': 'lat'}) # some models labelled dimensions differently...\n",
    "        \n",
    "    ds = xr.decode_cf(ds) # temporary hack, not sure why I need this but has to do with calendar-aware metadata on the time variable\n",
    "    ds = ds.sel(time=time_slice) # subset the data for the time period of interest\n",
    "    \n",
    "    # drop redundant variables (like \"height: 2m\")\n",
    "    for coord in ds.coords:\n",
    "        if coord not in ['lat','lon','time']:\n",
    "            ds = ds.drop(coord)\n",
    "\n",
    "    ## Calculate global-mean surface temperature (GMST)\n",
    "    cos_lat_2d = np.cos(np.deg2rad(ds['lat'])) * xr.ones_like(ds['lon']) # effective area weights\n",
    "    gmst = (\n",
    "        (ds['tas'] * cos_lat_2d).sum(dim=['lat','lon']) /\n",
    "        cos_lat_2d.sum(dim=['lat','lon'])\n",
    "    )\n",
    "    \n",
    "    # Add GMST to dictionary\n",
    "    gmst_dict[name] = gmst.squeeze()\n",
    "    \n",
    "    # Add near-surface air temperature to dictionary\n",
    "    ds_dict[name] = ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example `xarray.DataArray`\n",
    "\n",
    "The following is an example of an xarray.DataArray for the near-surface air temperature (`tas`) output from the MCM-UA climate model.\n",
    "\n",
    "The dimensions `Dimensions:   (lat: YY, lon: XX, member_id: 1, time: 1980)` signify that the 165 years of output is discretized into 1980 monthly-averages on a globe discretized into `YY` latitdue bins and `XX` longitude bins. The attributes of the DataArray describes various relevant metadata, such as a url for more info about the model, a reference publication which describes the simulation in detail, and the name of the model.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds['tas'].data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualizing warming patterns from 1890-1910 to 1990-2010"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature_change = (\n",
    "    ds['tas'].sel(time=slice('1990','2010')).mean(dim='time') -\n",
    "    ds['tas'].sel(time=slice('1890','1910')).mean(dim='time')\n",
    ").compute()\n",
    "temperature_change.attrs.update(ds.attrs)\n",
    "temperature_change = temperature_change.rename(\n",
    "    r'temperature change ($^{\\circ}$C) from 1890-1910 to 1990-2010'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "ortho = ccrs.Orthographic(-90, 20) # define target coordinate frame\n",
    "geo = ccrs.PlateCarree() # define origin coordinate frame\n",
    "\n",
    "plt.figure(figsize=(9,7))\n",
    "ax = plt.subplot(1, 1, 1, projection=ortho)\n",
    "\n",
    "q = temperature_change.plot(ax=ax, transform = geo) # plot a colormap in transformed coordinates\n",
    "\n",
    "ax.add_feature(cartopy.feature.COASTLINE)\n",
    "ax.add_feature(cartopy.feature.BORDERS, linestyle='-')\n",
    "plt.title('Patterns of global warming over the Americas',fontsize=16, ha='center');\n",
    "plt.savefig('../figures/historical_warming_patterns.png',dpi=100,bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Eagerly compute GMST for each model, in preparation for plotting\n",
    "\n",
    "The operations we have done up to this point to calculate the global-mean surface temperature were evaluated lazily with [`xarray`](http://xarray.pydata.org/en/stable/). In other worse, we have created a blueprint for how we want to evaluate the calculations, but have not yet computing them. This lets us do things like multiply two 1 Tb arrays together even though they are each individually larger-than-memory.\n",
    "\n",
    "Now we call xarray's `compute()` method to carry out the computations we defined in the for loop above for calculation the global-mean surface temperature anomaly relative to the 1850-1900 mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, gmst in tqdm(gmst_dict.items()): # loop through each model\n",
    "    \n",
    "    # remove the climatological seasonal cycle from the signal, for visual clarity\n",
    "    gmst_climatology = gmst.groupby('time.month').mean('time')\n",
    "    gmst_anomalies = gmst.groupby('time.month') - gmst_climatology\n",
    "    \n",
    "    # reference the GMST relative to the 1850 to 1900 pre-industrial period\n",
    "    gmst_anomalies = gmst_anomalies - gmst_anomalies.sel(time=slice('1850','1900')).mean(dim='time')\n",
    "    \n",
    "    # calculate annual-averages, to smooth the time-series a bit\n",
    "    gmst_ann_anomalies = gmst_anomalies.groupby('time.year').mean('time')\n",
    "    gmst_ann_anomalies = gmst_ann_anomalies.rename(\n",
    "        r'global-mean surface temperature ($^{\\circ}$C) (relative to 1850-1900)'\n",
    "    )\n",
    "    \n",
    "    gmst_dict[name] = gmst_ann_anomalies.compute() # eager calculation of global-mean surface temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmst_ann_anomalies.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting a timeseries of global warming from 1850-2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "\n",
    "for name, gmst in gmst_dict.items():\n",
    "    gmst.plot()\n",
    "    \n",
    "plt.xlabel('time')\n",
    "plt.ylabel(r'global-mean temperature ($^{\\circ}$C),'\n",
    "           '\\n'\n",
    "           'relative to 1850-1900 mean')\n",
    "plt.xlim([1850,2015]);\n",
    "plt.savefig('../figures/historical_global_warming.png',dpi=100,bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise\n",
    "\n",
    "#### Part 1 (Beginner)\n",
    "Edit the code above to plot (for an SSP of your choice):\n",
    " - a map of the temperature change from 2015-2025 to 2090-2100\n",
    " - a time-series of the global-mean surface temperature from 2015-2100\n",
    " \n",
    "#### Part 2 (Advanced)\n",
    "Edit the code above to plot (for an SSP of your choice):\n",
    " - a map of the temperature change from 1980-2000 to 2080-2100\n",
    " - a time-series of the global-mean surface temperature from 1850-2100\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The trickiest part of part 2 of the exercise is splicing together the historical simulation (which ends in 2014) and the SSP simulation (which starts in 2015). You may find the following code helpful:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "gmst_ssp_dict = {} # dictionary that will hold spliced DataArrays\n",
    "for name, gmst in gmst_dict.items(): # Loop through dictionary\n",
    "    if 'ssp' not in name: continue # keep only SSP simulations\n",
    "            \n",
    "    add_hist_gmst = None\n",
    "    \n",
    "    # Loop through dictionary to find matching historical simulation\n",
    "    for hist_name, hist_gmst in gmst_dict.items():\n",
    "        model_name = name.split(\".\")[2] # extract model name from name of simulation\n",
    "        \n",
    "        # find historical simulation that matches SSP simulation (same model)\n",
    "        if (model_name in hist_name) and ('historical' in hist_name):\n",
    "            add_hist_gmst = hist_gmst\n",
    "    \n",
    "    # if we found one, slice it together with the SSP and add it to the dictionary\n",
    "    if add_hist_gmst is not None:\n",
    "        gmst_ssp = xr.concat([add_hist_gmst, gmst],dim='year')\n",
    "        gmst_ssp_dict[name] = gmst_ssp\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cmip6-temperature-demo",
   "language": "python",
   "name": "cmip6-temperature-demo"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}