{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"[](https://colab.research.google.com/github/neuromatch/climate-course-content/blob/main/tutorials/W2D4_ClimateResponse-Extremes&Variability/W2D4_Tutorial8.ipynb) ![\"Open](\"https://kaggle.com/static/images/open-in-kaggle.svg\"/)
"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Bonus Tutorial 8: Thresholds\n",
"\n",
"**Week 2, Day 4, Extremes & Variability**\n",
"\n",
"**Content creators:** Matthias Aengenheyster, Joeri Reinders\n",
"\n",
"**Content reviewers:** Younkap Nina Duplex, Sloane Garelick, Paul Heubel, Zahra Khodakaramimaghsoud, Peter Ohue, Laura Paccini, Jenna Pearson, Derick Temfack, Peizhen Yang, Cheng Zhang, Chi Zhang, Ohad Zivan\n",
"\n",
"**Content editors:** Paul Heubel, Jenna Pearson, Chi Zhang, Ohad Zivan\n",
"\n",
"**Production editors:** Wesley Banfield, Paul Heubel, Jenna Pearson, Konstantine Tsafatinos, Chi Zhang, Ohad Zivan\n",
"\n",
"**Our 2024 Sponsors:** CMIP, NFDI4Earth"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Tutorial Objectives"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"*Estimated timing of tutorial:* 30 minutes\n",
"\n",
"The human body has physiological limits within which it can function. In hot conditions, the body cools itself through the process of sweating, where water evaporates from the skin, resulting in the loss of heat. The effectiveness of this cooling mechanism depends on the air's capacity to hold moisture. This is why sweating is more effective in dry heat, while humid heat feels \"hotter\" because it hampers the body's ability to cool down.\n",
"\n",
"As a result, the combination of temperature and humidity sets limits on the body's ability to regulate its temperature. One measure that captures this combined effect is the \"wet-bulb globe temperature,\" which combines information about ambient temperature, relative humidity, wind, and solar radiation to monitor heat stress risks while in direct sunlight. You can learn more about wet-bulb temperature on the following Wikipedia page: [Wet-bulb globe temperature](https://en.wikipedia.org/wiki/Wet-bulb_globe_temperature).\n",
"\n",
"In this tutorial, we will look at extreme levels of wet-bulb temperature spatially and consider the importance of thresholds. \n",
"\n",
"By the end of the tutorial, you will be able to:\n",
"1. Assess the risk of increasing wet-bulb globe temperatures.\n",
"2. Analyse how the probability of crossing threshold changes over time and between scenarios.\n",
"3. Assess the results of a spatial GEV analysis."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Setup"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"# installations ( uncomment and run this cell ONLY when using google colab or kaggle )\n",
"\n",
"# !pip install -q condacolab\n",
"# import condacolab\n",
"# condacolab.install()\n",
"# !mamba install numpy matplotlib seaborn pandas cartopy scipy texttable xarrayutils cf_xarray cmocean\n",
"\n",
"!wget https://raw.githubusercontent.com/neuromatch/climate-course-content/main/tutorials/W2D4_ClimateResponse-Extremes%26Variability/extremes_functions.py\n",
"!wget https://raw.githubusercontent.com/neuromatch/climate-course-content/main/tutorials/W2D4_ClimateResponse-Extremes%26Variability/gev_functions.py\n",
"!wget https://raw.githubusercontent.com/neuromatch/climate-course-content/main/tutorials/W2D4_ClimateResponse-Extremes%26Variability/sdfc_classes.py"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 2018,
"status": "ok",
"timestamp": 1682614821419,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"# imports\n",
"import xarray as xr\n",
"import matplotlib.pyplot as plt\n",
"import cartopy.crs as ccrs\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"import cmocean.cm as cmo\n",
"import os\n",
"import numpy as np\n",
"from scipy import stats\n",
"from scipy.stats import genextreme as gev\n",
"import pooch\n",
"import os\n",
"import tempfile\n",
"import gev_functions as gf\n",
"import extremes_functions as ef\n",
"import sdfc_classes as sd\n",
"\n",
"import warnings\n",
"\n",
"warnings.filterwarnings(\"ignore\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Note that `import gev_functions as gf` and `extremes_functions as ef` import the functions introduced in previous tutorials. Furthermore, `import sdfc_classes as sd` allows to use SDFC classes, that have been introduced in Tutorial 7."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Figure Settings\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"# @title Figure Settings\n",
"import ipywidgets as widgets # interactive display\n",
"\n",
"%config InlineBackend.figure_format = 'retina'\n",
"plt.style.use(\n",
" \"https://raw.githubusercontent.com/neuromatch/climate-course-content/main/cma.mplstyle\"\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Helper functions\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"# @title Helper functions\n",
"\n",
"\n",
"def pooch_load(filelocation=None, filename=None, processor=None):\n",
" shared_location = \"/home/jovyan/shared/Data/tutorials/W2D4_ClimateResponse-Extremes&Variability\" # this is different for each day\n",
" user_temp_cache = tempfile.gettempdir()\n",
"\n",
" if os.path.exists(os.path.join(shared_location, filename)):\n",
" file = os.path.join(shared_location, filename)\n",
" else:\n",
" file = pooch.retrieve(\n",
" filelocation,\n",
" known_hash=None,\n",
" fname=os.path.join(user_temp_cache, filename),\n",
" processor=processor,\n",
" )\n",
"\n",
" return file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Video 1: Thresholds\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"# @title Video 1: Thresholds\n",
"\n",
"from ipywidgets import widgets\n",
"from IPython.display import YouTubeVideo\n",
"from IPython.display import IFrame\n",
"from IPython.display import display\n",
"\n",
"\n",
"class PlayVideo(IFrame):\n",
" def __init__(self, id, source, page=1, width=400, height=300, **kwargs):\n",
" self.id = id\n",
" if source == 'Bilibili':\n",
" src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'\n",
" elif source == 'Osf':\n",
" src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'\n",
" super(PlayVideo, self).__init__(src, width, height, **kwargs)\n",
"\n",
"\n",
"def display_videos(video_ids, W=400, H=300, fs=1):\n",
" tab_contents = []\n",
" for i, video_id in enumerate(video_ids):\n",
" out = widgets.Output()\n",
" with out:\n",
" if video_ids[i][0] == 'Youtube':\n",
" video = YouTubeVideo(id=video_ids[i][1], width=W,\n",
" height=H, fs=fs, rel=0)\n",
" print(f'Video available at https://youtube.com/watch?v={video.id}')\n",
" else:\n",
" video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,\n",
" height=H, fs=fs, autoplay=False)\n",
" if video_ids[i][0] == 'Bilibili':\n",
" print(f'Video available at https://www.bilibili.com/video/{video.id}')\n",
" elif video_ids[i][0] == 'Osf':\n",
" print(f'Video available at https://osf.io/{video.id}')\n",
" display(video)\n",
" tab_contents.append(out)\n",
" return tab_contents\n",
"\n",
"\n",
"video_ids = [('Youtube', 'zC91CZCInGM'), ('Bilibili', 'BV1v8411D7jk')]\n",
"tab_contents = display_videos(video_ids, W=730, H=410)\n",
"tabs = widgets.Tab()\n",
"tabs.children = tab_contents\n",
"for i in range(len(tab_contents)):\n",
" tabs.set_title(i, video_ids[i][0])\n",
"display(tabs)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"pycharm": {
"name": "#%%\n"
},
"tags": [
"remove-input"
]
},
"outputs": [],
"source": [
"# @markdown\n",
"from ipywidgets import widgets\n",
"from IPython.display import IFrame\n",
"\n",
"link_id = \"fu8bv\"\n",
"\n",
"download_link = f\"https://osf.io/download/{link_id}/\"\n",
"render_link = f\"https://mfr.ca-1.osf.io/render?url=https://osf.io/{link_id}/?direct%26mode=render%26action=download%26mode=render\"\n",
"# @markdown\n",
"out = widgets.Output()\n",
"with out:\n",
" print(f\"If you want to download the slides: {download_link}\")\n",
" display(IFrame(src=f\"{render_link}\", width=730, height=410))\n",
"display(out)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Section 1: Downloading the Data\n",
"\n",
"In this tutorial, we will utilize wet-bulb globe temperature data derived from the MPI-ESM1-2-HR climate model, developed by the [Max Planck Institute for Meteorology](https://mpimet.mpg.de/en/homepage) in Hamburg, Germany. The data covers the historical period (hist) as well as three future [CMIP6](https://wcrp-cmip.org/) climate scenarios (SSP1-2.6, SSP2-4.5, and SSP5-8.5), which were introduced on W2D1 and in the previous Tutorial 6.\n",
"\n",
"You can learn more about CMIP, including the different scenarios, through our [CMIP Resource Bank](https://github.com/neuromatch/climate-course-content/blob/main/tutorials/CMIP/CMIP_resource_bank.md) and the [CMIP website](https://wcrp-cmip.org/).\n",
"\n",
"During the pre-processing phase, the data was subjected to a 7-day averaging process, followed by the computation of the annual maximum. As a result, the data for each grid point represents the wet-bulb temperature during the most extreme 7-day period within each year."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"# download file: 'WBGT_day_MPI-ESM1-2-HR_historical_r1i1p1f1_raw_runmean7_yearmax.nc'\n",
"\n",
"filename_WBGT_day = \"WBGT_day_MPI-ESM1-2-HR_historical_r1i1p1f1_raw_runmean7_yearmax.nc\"\n",
"url_WBGT_day = \"https://osf.io/69ms8/download\"\n",
"\n",
"wetbulb_hist = xr.open_dataset(pooch_load(url_WBGT_day, filename_WBGT_day)).WBGT\n",
"wetbulb_hist.attrs[\"units\"] = \"degC\""
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"The dataset consists of one entry per year. However, due to the inclusion of leap years, the data processing step resulted in different days assigned to each year. This discrepancy is deemed undesirable for analysis purposes. To address this, we resample the data by grouping all the data points belonging to each year and taking their average. Since there is only one data point per year, this resampling process does not alter the data itself but rather adjusts the time coordinate. This serves as a reminder to thoroughly inspect datasets before analysis, as overlooking such issues can lead to workflow failures."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"wetbulb_hist = wetbulb_hist.resample(time=\"1Y\").mean()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Let's load the data for the remaining scenarios:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 6,
"status": "ok",
"timestamp": 1682614821421,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"# SSP1-2.6 - 'WBGT_day_MPI-ESM1-2-HR_ssp126_r1i1p1f1_raw_runmean7_yearmax.nc'\n",
"filename_SSP126 = \"WBGT_day_MPI-ESM1-2-HR_ssp126_r1i1p1f1_raw_runmean7_yearmax.nc\"\n",
"url_SSP126 = \"https://osf.io/67b8m/download\"\n",
"wetbulb_126 = xr.open_dataset(pooch_load(url_SSP126, filename_SSP126)).WBGT\n",
"wetbulb_126.attrs[\"units\"] = \"degC\"\n",
"wetbulb_126 = wetbulb_126.resample(time=\"1Y\").mean()\n",
"\n",
"# SSP2-4.5 - WBGT_day_MPI-ESM1-2-HR_ssp245_r1i1p1f1_raw_runmean7_yearmax.nc\n",
"filename_SSP245 = \"WBGT_day_MPI-ESM1-2-HR_ssp245_r1i1p1f1_raw_runmean7_yearmax.nc\"\n",
"url_SSP245 = \"https://osf.io/fsx5y/download\"\n",
"wetbulb_245 = xr.open_dataset(pooch_load(url_SSP245, filename_SSP245)).WBGT\n",
"wetbulb_245.attrs[\"units\"] = \"degC\"\n",
"wetbulb_245 = wetbulb_245.resample(time=\"1Y\").mean()\n",
"\n",
"# SSP5-8.5 - WBGT_day_MPI-ESM1-2-HR_ssp585_r1i1p1f1_raw_runmean7_yearmax.nc\n",
"filename_SSP585 = \"WBGT_day_MPI-ESM1-2-HR_ssp585_r1i1p1f1_raw_runmean7_yearmax.nc\"\n",
"url_SSP585 = \"https://osf.io/pr456/download\"\n",
"wetbulb_585 = xr.open_dataset(pooch_load(url_SSP585, filename_SSP585)).WBGT\n",
"wetbulb_585.attrs[\"units\"] = \"degC\"\n",
"wetbulb_585 = wetbulb_585.resample(time=\"1Y\").mean()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Let's look at how the data is structured:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 2032,
"status": "ok",
"timestamp": 1682614897119,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"wetbulb_hist"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"There is one data point per year on a latitude-longitude grid. Let's compute the grid spacing in the longitude and latitude directions:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 3,
"status": "ok",
"timestamp": 1682614897449,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"wetbulb_hist.lon.diff(\"lon\").values.mean()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 2,
"status": "ok",
"timestamp": 1682614897910,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"wetbulb_hist.lat.diff(\"lat\").values.mean()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Each grid box in the dataset has an approximate size of 1 degree by 1 degree, which translates to about 100 km by 100 km at the equator. However, this size decreases as we move towards the poles due to the convergence of the meridians. It is important to consider the limitations imposed by this grid resolution.\n",
"\n",
"As a result, at the equator, the grid boxes cover an area of approximately 100 km by 100 km, while their size decreases in the mid-latitudes.\n",
"\n",
"Considering these grid box limitations, can you identify any potential limitations or challenges they may introduce in the analysis?"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Section 1.1: Focus on New Delhi, India"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Before doing spatial analysis, in the following, we focus on the capital of India: New Delhi."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 3,
"status": "ok",
"timestamp": 1682614597150,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"# find the nearest model point to the latitude and longitude of New Delhi\n",
"wetbulb_hist_delhi = wetbulb_hist.sel(lon=77.21, lat=28.61, method=\"nearest\")\n",
"wetbulb_126_delhi = wetbulb_126.sel(lon=77.21, lat=28.61, method=\"nearest\")\n",
"wetbulb_245_delhi = wetbulb_245.sel(lon=77.21, lat=28.61, method=\"nearest\")\n",
"wetbulb_585_delhi = wetbulb_585.sel(lon=77.21, lat=28.61, method=\"nearest\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 1155,
"status": "ok",
"timestamp": 1682614598302,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"# plot time series of all scenarios for New Delhi\n",
"fig, ax = plt.subplots()\n",
"wetbulb_hist_delhi.plot(linestyle=\"-\", marker=\".\", label=\"hist\", ax=ax)\n",
"wetbulb_126_delhi.plot(linestyle=\"-\", marker=\".\", label=\"ssp126\", ax=ax)\n",
"wetbulb_245_delhi.plot(linestyle=\"-\", marker=\".\", label=\"ssp245\", ax=ax)\n",
"wetbulb_585_delhi.plot(linestyle=\"-\", marker=\".\", label=\"ssp585\", ax=ax)\n",
"\n",
"ax.legend()\n",
"ax.set_title(\"\")\n",
"ax.set_xlabel(\"Time (years)\")\n",
"ax.set_ylabel(\"Maximum 7-day Mean\\nWet-Bulb Globe Temperature ($\\degree$C)\")\n",
"ax.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Note:\n",
"1. Trends are visible in the historical period\n",
"2. Distinct differences between climate scenarios are apparent\n",
"3. Strong variability - each years maximum is not necessarily warmer than the previous one\n",
"\n",
"Let's fit the data to a GEV distribution and get the associated return levels."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"# GEV fit for historical period in New Delhi with bootstrapping for uncertainty estimation\n",
"shape_hist, loc_hist, scale_hist = gev.fit(wetbulb_hist_delhi.values, 0)\n",
"return_levels_hist = gf.fit_return_levels(\n",
" wetbulb_hist_delhi.values, years=np.arange(1.1, 1000), N_boot=100\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Now we can plot the probability density functions, the return levels, and assess the fit using the QQ plot from previous tutorials."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"# create figures\n",
"fig, axs = plt.subplots(2, 2, constrained_layout=True)\n",
"# reduce dimension of axes to have only one index\n",
"ax = axs.flatten()\n",
"\n",
"# create x-data for QQ-plot (top left)\n",
"x = np.linspace(0, 1, 100)\n",
"\n",
"# calculate quantiles by sorting from min to max\n",
"# and plot vs. return levels calculated with ppf (percent point function) and GEV fit parameters,\n",
"# cf. e.g. Tut4)\n",
"ax[0].plot(\n",
" gev.ppf(x, shape_hist, loc=loc_hist, scale=scale_hist),\n",
" np.quantile(wetbulb_hist_delhi, x),\n",
" \"o\",\n",
")\n",
"# plot identity line\n",
"xlim = ax[0].get_xlim()\n",
"ylim = ax[0].get_ylim()\n",
"ax[0].plot(\n",
" [min(xlim[0], ylim[0]), max(xlim[1], ylim[1])],\n",
" [min(xlim[0], ylim[0]), max(xlim[1], ylim[1])],\n",
" \"k\",\n",
")\n",
"\n",
"#ax[0].set_xlim(xlim)\n",
"#ax[0].set_ylim(ylim)\n",
"\n",
"# create time data array for PDF comparison (bottom left)\n",
"x = np.linspace(wetbulb_hist_delhi.min() - 1, wetbulb_hist_delhi.max() + 1, 1000)\n",
"wetbulb_hist_delhi.plot.hist(\n",
" bins=np.arange(29, 33, 0.25),\n",
" histtype=\"step\",\n",
" density=True,\n",
" lw=1,\n",
" color=\"k\",\n",
" ax=ax[2],\n",
" label=\"histogram\",\n",
")\n",
"# methodology comparison (bottom left)\n",
"# plot PDF of modeled GEV and empirical PDF estimate (kernel density estimation)\n",
"ax[2].plot(x, gev.pdf(x, shape_hist, loc=loc_hist, scale=scale_hist), label=\"Modeled\")\n",
"sns.kdeplot(wetbulb_hist_delhi, ax=ax[2], label=\"Empirical\")\n",
"ax[2].legend()\n",
"\n",
"# plot return levels (bottom right) using functions also used in Tutorial 6 and 7\n",
"gf.plot_return_levels(return_levels_hist, ax=ax[3])\n",
"ax[3].set_xlim(1.5, 1000)\n",
"# ax[3].set_ylim(0,None)\n",
"\n",
"# aesthetics\n",
"ax[0].set_title(\"QQ-plot\")\n",
"ax[2].set_title(\"PDF\")\n",
"ax[3].set_title(\"Return levels\")\n",
"\n",
"ax[0].set_xlabel(\"Wet-Bulb Globe Temperature ($\\degree$C)\")\n",
"ax[0].set_ylabel(\"Return Levels ($\\degree$C)\")\n",
"\n",
"ax[2].set_xlabel(\"Wet-Bulb Globe Temperature ($\\degree$C)\")\n",
"ax[3].set_xlabel(\"Return Periods (years)\")\n",
"ax[3].set_ylabel(\"Return Levels ($\\degree$C)\")\n",
"\n",
"ax[1].remove()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"\n",
" Click here for a detailed description of the plot.
\n",
"\n",
"*Top left*: \n",
"The QQ-plot shows the quantiles versus the return levels both in $\\degree$C, calculated from the wet-bulb globe temperature data and emphasized as blue dots. For small quantiles, the return levels follow the identity line (black), in contrast to a few larger quantiles, where return levels increase with a smaller slope than the id line. This indicates a slight difference in the distributions of the two properties.\n",
" \n",
"*Bottom left*: \n",
"A comparison between modeled and empirical probability density functions (PDFs) shows slightly different distributions. The empirical PDF weights the density of larger anomalies more, so it appears that the GEV PDF underestimates large wet-bulb globe temperatures. In addition, the modeled scale parameter yields a higher maximum than the empirical estimate. The tails are comparable, however the empirical PDF sees more events in the extremes.\n",
"\n",
"*Bottom right*: \n",
"The fitted return levels of the wet-bulb globe temperatures in $\\degree$C versus return periods log-scaled in years. Blue dots highlight the empirical calculations (using the algorithm introduced in Tutorial 2), while the blue line emphasizes the fitted GEV corresponding relationship. Shaded areas represent the uncertainty ranges that were calculated with a bootstrap approach before. The larger the return periods, the less data on such events is given, hence the uncertainty increases. \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"We use the function `estimate_return_level_period()` to compute the X-year return level:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 4,
"status": "ok",
"timestamp": 1682614629210,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"# X-year of interest\n",
"x = 100\n",
"print(\n",
" \"{}-year return level: {:.2f} C\".format(x,\n",
" gf.estimate_return_level_period(x, loc_hist, scale_hist, shape_hist).mean()\n",
" )\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Now let's compare with the last 30 years of the SSP-245 scenario, the middle scenario we looked at before. 2050-2100 are approximately stationary here, so we can leave out that utility."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 22216,
"status": "ok",
"timestamp": 1682614651424,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"shape_245, loc_245, scale_245 = gev.fit(\n",
" wetbulb_245_delhi.sel(time=slice(\"2070\", \"2100\")).values, 0\n",
")\n",
"return_levels_245 = gf.fit_return_levels(\n",
" wetbulb_245_delhi.sel(time=slice(\"2070\", \"2100\")).values,\n",
" years=np.arange(1.1, 1000),\n",
" N_boot=100,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 2162,
"status": "ok",
"timestamp": 1682614653581,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"# create figures\n",
"fig, axs = plt.subplots(2, 2, constrained_layout=True)\n",
"# reduce dimension of axes to have only one index\n",
"ax = axs.flatten()\n",
"\n",
"# create x-data for QQ-plot (top left)\n",
"x = np.linspace(0, 1, 100)\n",
"\n",
"# calculate quantiles by sorting from min to max\n",
"# and plot vs. return levels calculated with ppf (percent point function) and GEV fit parameters,\n",
"# cf. e.g. Tut4)\n",
"ax[0].plot(\n",
" gev.ppf(x, shape_245, loc=loc_245, scale=scale_245),\n",
" np.quantile(wetbulb_245_delhi.sel(time=slice(\"2051\", \"2100\")), x),\n",
" \"o\",\n",
")\n",
"# plot identity line\n",
"xlim = ax[0].get_xlim()\n",
"ylim = ax[0].get_ylim()\n",
"ax[0].plot(\n",
" [min(xlim[0], ylim[0]), max(xlim[1], ylim[1])],\n",
" [min(xlim[0], ylim[0]), max(xlim[1], ylim[1])],\n",
" \"k\",\n",
")\n",
"\n",
"#ax[0].set_xlim(xlim)\n",
"#ax[0].set_ylim(ylim)\n",
"\n",
"# create time data array for PDF comparison (bottom left)\n",
"x = np.linspace(\n",
" wetbulb_245_delhi.sel(time=slice(\"2051\", \"2100\")).min() - 1,\n",
" wetbulb_245_delhi.sel(time=slice(\"2051\", \"2100\")).max() + 1,\n",
" 1000,\n",
")\n",
"# plot histogram\n",
"wetbulb_245_delhi.sel(time=slice(\"2051\", \"2100\")).plot.hist(\n",
" bins=np.arange(29, 33, 0.25),\n",
" histtype=\"step\",\n",
" density=True,\n",
" lw=1,\n",
" color=\"k\",\n",
" ax=ax[2],\n",
" label=\"histogram\",\n",
")\n",
"# add GEV PDF defined by fit parameters\n",
"ax[2].plot(x, gev.pdf(x, shape_245, loc=loc_245, scale=scale_245), label=\"Modeled\")\n",
"# add kernel density estimation (kde) of PDF\n",
"sns.kdeplot(\n",
" wetbulb_245_delhi.sel(time=slice(\"2051\", \"2100\")), ax=ax[2], label=\"Empirical\"\n",
")\n",
"ax[2].legend()\n",
"\n",
"# plot return levels (bottom right) using functions also used in Tutorial 6 and 7\n",
"gf.plot_return_levels(return_levels_245, ax=ax[3])\n",
"ax[3].set_xlim(1.5, 1000)\n",
"# ax[3].set_ylim(0,None)\n",
"\n",
"# aesthetics\n",
"ax[0].set_title(\"QQ-plot\")\n",
"ax[2].set_title(\"PDF\")\n",
"ax[3].set_title(\"Return levels\")\n",
"\n",
"ax[0].set_xlabel(\"Wet-Bulb Globe Temperature ($\\degree$C)\")\n",
"ax[0].set_ylabel(\"Return Levels ($\\degree$C)\")\n",
"\n",
"ax[2].set_xlabel(\"Wet-Bulb Globe Temperature ($\\degree$C)\")\n",
"ax[3].set_xlabel(\"Return Periods (years)\")\n",
"ax[3].set_ylabel(\"Return Levels ($\\degree$C)\")\n",
"\n",
"ax[1].remove()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"\n",
" Click here for a detailed description of the plot.
\n",
"\n",
"*Top left*: \n",
"The QQ-plot shows the quantiles versus the return levels both in $\\degree$C, calculated from the wet-bulb globe temperature data and emphasized as blue dots. For small quantiles, the return levels follow the identity line (black) less, additionally, in a few larger quantiles corresponding return levels increase with a smaller slope than the id line. This indicates a slight difference in the distributions of the two properties.\n",
" \n",
"*Bottom left*: \n",
"A comparison between modeled and empirical probability density functions (PDFs) shows slightly different distributions. The empirical PDF weights the density of larger anomalies more, so it appears that the GEV PDF underestimates large wet-bulb globe temperatures. In addition, the modeled scale parameter yields a higher maximum than the empirical estimate. The tails are comparable, however the empirical PDF sees more events in the extremes, especially in the maximum range.\n",
"\n",
"*Bottom right*: \n",
"The fitted return levels of the wet-bulb globe temperatures in $\\degree$C versus return periods log-scaled in years. Blue dots highlight the empirical calculations (using the algorithm introduced in Tutorial 2), while the blue line emphasizes the fitted GEV corresponding relationship. Shaded areas represent the uncertainty ranges that were calculated with a bootstrap approach before. The larger the return periods, the less data on such events is given, hence the uncertainty increases. \n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"# X-year of interest\n",
"x = 100\n",
"print(\n",
" \"{}-year return level: {:.2f} C\".format(x,\n",
" gf.estimate_return_level_period(x, loc_245, scale_245, shape_245).mean()\n",
" )\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Compute as well the fit and return levels for the remaining two scenarios (SSP-126 and SSP-585). Save the QQ plot etc. for your own testing later on.\n",
"\n",
"You can then plot all return level curves together to compare."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"return_levels_126 = gf.fit_return_levels(\n",
" wetbulb_126_delhi.sel(time=slice(\"2070\", \"2100\")).values,\n",
" years=np.arange(1.1, 1000),\n",
" N_boot=100,\n",
")\n",
"return_levels_585 = gf.fit_return_levels(\n",
" wetbulb_585_delhi.sel(time=slice(\"2070\", \"2100\")).values,\n",
" years=np.arange(1.1, 1000),\n",
" N_boot=100,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"gf.plot_return_levels(return_levels_hist, c=\"k\", label=\"historical\", ax=ax)\n",
"gf.plot_return_levels(return_levels_126, c=\"C0\", label=\"ssp126\", ax=ax)\n",
"gf.plot_return_levels(return_levels_245, c=\"C1\", label=\"ssp245\", ax=ax)\n",
"gf.plot_return_levels(return_levels_585, c=\"C2\", label=\"ssp585\", ax=ax)\n",
"\n",
"ax.set_xlim(1, 100)\n",
"ax.set_ylim(29.5, 37)\n",
"ax.legend()\n",
"ax.grid(True, which=\"both\")\n",
"ax.set_xlabel(\"Return Period (years)\")\n",
"ax.set_ylabel(\"Return Level ($\\degree$ C)\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"\n",
" Click here for a description of the plot.
\n",
"\n",
"The fitted return levels of the wet-bulb globe temperatures in $\\degree$C versus return periods log-scaled in years for all four scenarios. Dots highlight the empirical calculations (using the algorithm introduced in Tutorial 2), while lines emphasize the fitted GEV corresponding relationship, respectively. Shaded areas represent the uncertainty ranges that were calculated with a bootstrap approach before. The larger the return periods, the less data on such events is given, hence the uncertainty increases. \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Questions 1.1\n",
"\n",
"Compare the common event (return period << 10 years) to the very rare event (return period ~ 100 years) events under the different scenarios.\n",
"\n",
"1. What is the return level of a 3-year event under the SSP5-8.5 scenario? Note down the level. What would be the return period of such an event in the other scenarios?\n",
"\n",
"2. What is the return level of a 100-year event in the historical scenario. How often would such an event occur under the other scenarios?"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"execution": {}
},
"source": [
"[*Click for solution*](https://github.com/neuromatch/climate-course-content/tree/main/tutorials/W2D4_ClimateResponse-Extremes&Variability/solutions/W2D4_Tutorial8_Solution_a557abd9.py)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Section 1.2: Return Levels Over Different Intervals"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Besides the late period (2070-2011), we compute the return levels over the near future (2015-2050). Then let's plot the time series, and overlay the 100-year return level, as computed over 2015-2050, 2070-2100, and the historical period:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"return_levels_126_2015_2050 = gf.fit_return_levels(\n",
" wetbulb_126_delhi.sel(time=slice(\"2015\", \"2050\")).values,\n",
" years=np.arange(1.1, 1000),\n",
" N_boot=100,\n",
")\n",
"return_levels_245_2015_2050 = gf.fit_return_levels(\n",
" wetbulb_245_delhi.sel(time=slice(\"2015\", \"2050\")).values,\n",
" years=np.arange(1.1, 1000),\n",
" N_boot=100,\n",
")\n",
"return_levels_585_2015_2050 = gf.fit_return_levels(\n",
" wetbulb_585_delhi.sel(time=slice(\"2015\", \"2050\")).values,\n",
" years=np.arange(1.1, 1000),\n",
" N_boot=100,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"# calculate annual mean for all scenarios and plot time series\n",
"wetbulb_hist_delhi.groupby(\"time.year\").mean().plot.line(\n",
" alpha=0.5, c=\"k\", label=\"hist\", ax=ax\n",
")\n",
"wetbulb_126_delhi.groupby(\"time.year\").mean().plot.line(\n",
" alpha=0.5, c=\"C0\", label=\"ssp126\", ax=ax\n",
")\n",
"wetbulb_245_delhi.groupby(\"time.year\").mean().plot.line(\n",
" alpha=0.5, c=\"C1\", label=\"ssp245\", ax=ax\n",
")\n",
"wetbulb_585_delhi.groupby(\"time.year\").mean().plot.line(\n",
" alpha=0.5, c=\"C2\", label=\"ssp585\", ax=ax\n",
")\n",
"# aesthetics\n",
"ax.set_title(\"\")\n",
"ax.legend()\n",
"\n",
"# add horizontal lines of 100 year event return levels for certain periods\n",
"# historical 1950-2014\n",
"ax.hlines(\n",
" return_levels_hist.GEV.sel(period=100, method=\"nearest\").values,\n",
" 1950,\n",
" 2014,\n",
" \"k\",\n",
" linestyle=\"-.\",\n",
" lw=2,\n",
")\n",
"# SSP1-2.6 2015-2050\n",
"ax.hlines(\n",
" return_levels_126_2015_2050.GEV.sel(period=100, method=\"nearest\").values,\n",
" 2015,\n",
" 2050,\n",
" \"C0\",\n",
" linestyle=\"--\",\n",
" lw=2,\n",
")\n",
"# SSP2-4.5 2015-2050\n",
"ax.hlines(\n",
" return_levels_245_2015_2050.GEV.sel(period=100, method=\"nearest\").values,\n",
" 2015,\n",
" 2050,\n",
" \"C1\",\n",
" linestyle=\"--\",\n",
" lw=2,\n",
")\n",
"# SSP5-8.5 2015-2050\n",
"ax.hlines(\n",
" return_levels_585_2015_2050.GEV.sel(period=100, method=\"nearest\").values,\n",
" 2015,\n",
" 2050,\n",
" \"C2\",\n",
" linestyle=\"--\",\n",
" lw=2,\n",
")\n",
"# SSP1-2.6 2070-2100\n",
"ax.hlines(\n",
" return_levels_126.GEV.sel(period=100, method=\"nearest\").values,\n",
" 2070,\n",
" 2100,\n",
" \"C0\",\n",
" linestyle=\":\",\n",
" lw=2,\n",
")\n",
"# SSP2-4.5 2070-2100\n",
"ax.hlines(\n",
" return_levels_245.GEV.sel(period=100, method=\"nearest\").values,\n",
" 2070,\n",
" 2100,\n",
" \"C1\",\n",
" linestyle=\":\",\n",
" lw=2,\n",
")\n",
"# SSP5-8.5 2070-2100\n",
"ax.hlines(\n",
" return_levels_585.GEV.sel(period=100, method=\"nearest\").values,\n",
" 2070,\n",
" 2100,\n",
" \"C2\",\n",
" linestyle=\":\",\n",
" lw=2,\n",
")\n",
"# aesthetics\n",
"plt.title(\"\")\n",
"plt.xlabel(\"Time (years)\")\n",
"plt.ylabel(\"Maximum 7-day Mean\\nWet-Bulb Globe Temperature ($\\degree$C)\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"\n",
" Click here for a detailed description of the plot.
\n",
"Time series of the wet-bulb globe temperature in $\\degree$C in India, New Delhi for three scenarios and the historical period, overlaid by the 100-year return level of certain periods, as computed over 1950-2014 (dot-dashed), 2015-2050 (dashed), and 2070-2100 (dotted). In summary, the return levels increase over time, only the blue scenario (SSP1-2.6) keeps it at a temporally constant level at the end of the century, while the more severe scenarios still show rising return levels. The SSP5-8.5, often called the 'business-as-usual' scenario due to omitted implementation of policies, shows a return level of 36 $\\degree$C for 100-year events, 3 years more than the SSP2-4.5 scenario.\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Section 1.3: Time-Dependent Return Levels"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Looking at the previous plot we see trends present in our datasets, with the 100-year event return levels varying across the time periods we have choosen. This suggests that our location parameter is changing with time.\n",
"\n",
"Now, similar to the previous tutorial, we assume that the location parameter is a function of time and proceed to estimate the GEV distribution for the four scenarios:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"def estimate_return_level_model(quantile, model):\n",
" loc, scale, shape = model.loc, model.scale, model.shape\n",
" level = loc - scale / shape * (1 - (-np.log(quantile)) ** (-shape))\n",
" return level"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 20,
"status": "ok",
"timestamp": 1682614653581,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"# fit the GEV to the data, while specifying that the location parameter ('loc') is meant\n",
"# to be a covariate ('c_') of the time axis (data.index)\n",
"\n",
"law_ns_hist = sd.GEV()\n",
"for i in range(250):\n",
" law_ns_hist.fit(wetbulb_hist_delhi.values, c_loc=np.arange(wetbulb_hist_delhi.time.size))\n",
" #print(law_ns_hist.coef_)\n",
" # if the first coefficient is not zero, we stop fitting\n",
" if law_ns_hist.coef_[0] != 0:\n",
" print(f'Found non-trivial solution for law_ns_hist after {i} fitting iterations.')\n",
" print(law_ns_hist.coef_)\n",
" break\n",
"\n",
"\n",
"law_ns_126 = sd.GEV()\n",
"for i in range(250):\n",
" law_ns_126.fit(wetbulb_126_delhi.values, c_loc=np.arange(wetbulb_126_delhi.time.size))\n",
" # if the first coefficient is not zero, we stop fitting\n",
" if law_ns_126.coef_[0] != 0:\n",
" print(f'Found non-trivial solution for law_ns_126 after {i} fitting iterations.')\n",
" print(law_ns_126.coef_)\n",
" break\n",
"\n",
"law_ns_245 = sd.GEV()\n",
"for i in range(250):\n",
" law_ns_245.fit(wetbulb_245_delhi.values, c_loc=np.arange(wetbulb_245_delhi.time.size))\n",
" # if the first coefficient is not zero, we stop fitting\n",
" if law_ns_245.coef_[0] != 0:\n",
" print(f'Found non-trivial solution for law_ns_245 after {i} fitting iterations.')\n",
" print(law_ns_245.coef_)\n",
" break\n",
"\n",
"law_ns_585 = sd.GEV()\n",
"for i in range(250):\n",
" law_ns_585.fit(wetbulb_585_delhi.values, c_loc=np.arange(wetbulb_585_delhi.time.size))\n",
" # if the first coefficient is not zero, we stop fitting\n",
" if law_ns_585.coef_[0] != 0:\n",
" print(f'Found non-trivial solution for law_ns_585 after {i} fitting iterations.')\n",
" print(law_ns_585.coef_)\n",
" break"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 20,
"status": "ok",
"timestamp": 1682614653581,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"wetbulb_hist_delhi.plot.line(c=\"k\", ax=ax)\n",
"wetbulb_126_delhi.plot.line(c=\"C0\", ax=ax)\n",
"wetbulb_245_delhi.plot.line(c=\"C1\", ax=ax)\n",
"wetbulb_585_delhi.plot.line(c=\"C2\", ax=ax)\n",
"\n",
"ax.plot(\n",
" wetbulb_hist_delhi.time,\n",
" estimate_return_level_model(1 - 1 / 100, law_ns_hist),\n",
" \"k--\",\n",
" label=\"100-year return level: hist\",\n",
")\n",
"ax.plot(\n",
" wetbulb_126_delhi.time,\n",
" estimate_return_level_model(1 - 1 / 100, law_ns_126),\n",
" \"C0--\",\n",
" label=\"100-year return level: ssp126\",\n",
")\n",
"ax.plot(\n",
" wetbulb_245_delhi.time,\n",
" estimate_return_level_model(1 - 1 / 100, law_ns_245),\n",
" \"C1--\",\n",
" label=\"100-year return level: ssp245\",\n",
")\n",
"ax.plot(\n",
" wetbulb_585_delhi.time,\n",
" estimate_return_level_model(1 - 1 / 100, law_ns_585),\n",
" \"C2--\",\n",
" label=\"100-year return level: ssp585\",\n",
")\n",
"\n",
"ax.legend()\n",
"ax.set_title(\"\")\n",
"ax.set_xlabel(\"Time (years)\")\n",
"ax.set_ylabel(\"Maximum 7-day Mean\\nWet-Bulb Globe Temperature ($\\degree$C)\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Now we again compute the AIC for the constant and time-dependent models, and compare their performance:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"executionInfo": {
"elapsed": 13,
"status": "ok",
"timestamp": 1682614653581,
"user": {
"displayName": "Matthias Aengenheyster",
"userId": "16322208118439170907"
},
"user_tz": -120
},
"tags": []
},
"outputs": [],
"source": [
"def compute_aic(model):\n",
" return 2 * len(model.coef_) + 2 * model.info_.mle_optim_result.fun"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"# compute stationary models:\n",
"law_ss_hist = sd.GEV()\n",
"law_ss_hist.fit(wetbulb_hist_delhi.values)\n",
"\n",
"law_ss_126 = sd.GEV()\n",
"law_ss_126.fit(wetbulb_126_delhi.values)\n",
"\n",
"law_ss_126 = sd.GEV()\n",
"law_ss_126.fit(wetbulb_126_delhi.values)\n",
"\n",
"law_ss_245 = sd.GEV()\n",
"law_ss_245.fit(wetbulb_245_delhi.values)\n",
"\n",
"law_ss_585 = sd.GEV()\n",
"law_ss_585.fit(wetbulb_585_delhi.values)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"aics = pd.DataFrame(\n",
" columns=[\"hist\", \"ssp126\", \"ssp245\", \"ssp585\"], index=[\"constant\", \"covariate\"]\n",
")\n",
"\n",
"aics[\"hist\"] = compute_aic(law_ss_hist), compute_aic(law_ns_hist)\n",
"aics[\"ssp126\"] = compute_aic(law_ss_126), compute_aic(law_ns_126)\n",
"aics[\"ssp245\"] = compute_aic(law_ss_245), compute_aic(law_ns_245)\n",
"aics[\"ssp585\"] = compute_aic(law_ss_585), compute_aic(law_ns_585)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"aics.round(2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"The AIC is lower when using a covariate, suggesting that including the time dependence in the location parameter improves the quality of the model. The exception is the SSP1-2.6 scenario, which does not perform as well. This is because, unlike the other scenarios and the historical period, the wet-bulb globe temperatures stabilize, hence this location parameter is less dependent on time. In this instance, making other parameters depend on time could potentially improve the performance."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Section 2: Spatial Analysis"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"After looking at New Delhi, India, now we can make use of the spatial information:"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"\n",
"The code provided below is commented and is used to fit the GEV distribution for each grid point. For the historical scenario, the entire time range is used, while for the selected scenarios, the period from 2071 to 2100 (the last 30 years of the data) is used.\n",
"\n",
"Please note that the computation for this code takes some time (approximately 9 minutes per dataset). To save time, we have already precomputed the data, so there is no need to run the commented code. However, you are free to uncomment and run the code, make modifications, or include time-dependent parameters (as shown above) at your convenience. If desired, you can also focus on specific regions."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Expensive code that fits a GEV distribution to each grid point:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"# this code requires the authors' extremes_functions.py file and SDFC library from github: https://github.com/yrobink/SDFC\n",
"# The code takes roughly 30 minutes to execute, in the next cell we load in the precomputed data. Uncomment the following lines if you want to rerun.\n",
"# fit_sp_hist = ef.fit_return_levels_sdfc_2d(wetbulb_hist.rename({'lon':'longitude','lat':'latitude'}),times=np.arange(1.1,1000),periods_per_year=1,kind='GEV',N_boot=0,full=True)\n",
"# fit_sp_126 = ef.fit_return_levels_sdfc_2d(wetbulb_126.sel(time=slice('2071','2100')).rename({'lon':'longitude','lat':'latitude'}),times=np.arange(1.1,1000),periods_per_year=1,kind='GEV',N_boot=0,full=True)\n",
"# fit_sp_245 = ef.fit_return_levels_sdfc_2d(wetbulb_245.sel(time=slice('2071','2100')).rename({'lon':'longitude','lat':'latitude'}),times=np.arange(1.1,1000),periods_per_year=1,kind='GEV',N_boot=0,full=True)\n",
"# fit_sp_585 = ef.fit_return_levels_sdfc_2d(wetbulb_585.sel(time=slice('2071','2100')).rename({'lon':'longitude','lat':'latitude'}),times=np.arange(1.1,1000),periods_per_year=1,kind='GEV',N_boot=0,full=True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Section 2.1: Load Pre-Computed Data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"# historical - wbgt_hist_raw_runmean7_gev.nc\n",
"fname_wbgt_hist = \"wbgt_hist_raw_runmean7_gev.nc\"\n",
"url_wbgt_hist = \"https://osf.io/dakv3/download\"\n",
"fit_sp_hist = xr.open_dataset(pooch_load(url_wbgt_hist, fname_wbgt_hist))\n",
"\n",
"# SSP-126 - wbgt_126_raw_runmean7_gev_2071-2100.nc\n",
"fname_wbgt_126 = \"wbgt_126_raw_runmean7_gev_2071-2100.nc\"\n",
"url_wbgt_126 = \"https://osf.io/ef9pv/download\"\n",
"fit_sp_126 = xr.open_dataset(pooch_load(url_wbgt_126, fname_wbgt_126))\n",
"\n",
"# SSP-245 - wbgt_245_raw_runmean7_gev_2071-2100.nc\n",
"fname_wbgt_245 = \"wbgt_245_raw_runmean7_gev_2071-2100.nc\"\n",
"url_wbgt_245 = \"https://osf.io/j4hfc/download\"\n",
"fit_sp_245 = xr.open_dataset(pooch_load(url_wbgt_245, fname_wbgt_245))\n",
"\n",
"# SSP-585 - wbgt_585_raw_runmean7_gev_2071-2100.nc\n",
"fname_bgt_58 = \"wbgt_585_raw_runmean7_gev_2071-2100.nc\"\n",
"url_bgt_585 = \"https://osf.io/y6edw/download\"\n",
"fit_sp_585 = xr.open_dataset(pooch_load(url_bgt_585, fname_bgt_58))"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Also load the area for each grid box - we will use this later to compute global averages:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"# filename - area_mpi.nc\n",
"filename_area_mpi = \"area_mpi.nc\"\n",
"url_area_mpi = \"https://osf.io/zqd86/download\"\n",
"area = xr.open_dataarray(pooch_load(url_area_mpi, filename_area_mpi))\n",
"\n",
"# filename - area_land_mpi.nc\n",
"filename_area_mpi = \"area_land_mpi.nc\"\n",
"url_area_land_mpi = \"https://osf.io/dxq98/download\"\n",
"area_land = xr.open_dataarray(pooch_load(url_area_land_mpi, filename_area_mpi)).fillna(\n",
" 0.0\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Now, let's examine the 100-year return level in the historical run and compare it to the period from 2071-2100 in the three scenarios. The colorbar has been set to start at 28 $\\degree$C, which is approximately the temperature reached during the severe heatwaves in Europe in 2003 and Russia in 2010."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"# setup plots\n",
"fig, axs = plt.subplots(\n",
" 2,\n",
" 2,\n",
" constrained_layout=True,\n",
" figsize=(12, 8),\n",
" subplot_kw=dict(projection=ccrs.Robinson()),\n",
")\n",
"# reduce axis dimension to have one index only\n",
"ax = axs.flatten()\n",
"\n",
"# summarize all kwargs used for al plots\n",
"kwargs = dict(\n",
" vmin=28, vmax=38, cmap=cmo.amp, transform=ccrs.PlateCarree(), add_colorbar=False\n",
")\n",
"\n",
"# select 100-year event for every grid cell in every scenario dataset\n",
"p = (\n",
" fit_sp_hist[\"return level\"]\n",
" .sel({\"return period\": 100}, method=\"nearest\")\n",
" .plot(ax=ax[0], **kwargs)\n",
")\n",
"fit_sp_126[\"return level\"].sel({\"return period\": 100}, method=\"nearest\").plot(\n",
" ax=ax[1], **kwargs\n",
")\n",
"fit_sp_245[\"return level\"].sel({\"return period\": 100}, method=\"nearest\").plot(\n",
" ax=ax[2], **kwargs\n",
")\n",
"fit_sp_585[\"return level\"].sel({\"return period\": 100}, method=\"nearest\").plot(\n",
" ax=ax[3], **kwargs\n",
")\n",
"\n",
"# define colorbar\n",
"cbar = fig.colorbar(\n",
" p,\n",
" ax=ax,\n",
" pad=0.025,\n",
" orientation=\"horizontal\",\n",
" shrink=0.75,\n",
" label=\"100-year return level (degree)\",\n",
" extend=\"max\",\n",
")\n",
"\n",
"# aesthetics for all subplots\n",
"ax[0].set_title(\"Historical\")\n",
"ax[1].set_title(\"SSP1-2.6, 2091-2100\")\n",
"ax[2].set_title(\"SSP2-4.5, 2091-2100\")\n",
"ax[3].set_title(\"SSP5-8.5, 2091-2100\")\n",
"\n",
"[axi.set_facecolor(\"grey\") for axi in ax]\n",
"[axi.coastlines(lw=0.5) for axi in ax]"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"In the following regions where the hottest heatwave is above 31 degrees wet-bulb globe temperature are given by the red shading, which is considered a \"critical temperature\" above which a human will die within a few hours without non-evaporative cooling like air conditioning:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"fig, axs = plt.subplots(\n",
" 2,\n",
" 2,\n",
" constrained_layout=True,\n",
" figsize=(12, 8),\n",
" subplot_kw=dict(projection=ccrs.Robinson()),\n",
")\n",
"ax = axs.flatten()\n",
"\n",
"kwargs = dict(vmin=0, cmap=cmo.amp, transform=ccrs.PlateCarree(), add_colorbar=False)\n",
"\n",
"p = (wetbulb_hist.sel(time=slice(\"2005\", \"2014\")).max(\"time\") > 31).plot(\n",
" ax=ax[0], **kwargs\n",
")\n",
"(wetbulb_126.sel(time=slice(\"2091\", \"2100\")).max(\"time\") > 31).plot(ax=ax[1], **kwargs)\n",
"(wetbulb_245.sel(time=slice(\"2091\", \"2100\")).max(\"time\") > 31).plot(ax=ax[2], **kwargs)\n",
"(wetbulb_585.sel(time=slice(\"2091\", \"2100\")).max(\"time\") > 31).plot(ax=ax[3], **kwargs)\n",
"\n",
"# cbar = fig.colorbar(p,ax=ax,pad=0.025,orientation='horizontal',shrink=0.75,label='Most extreme 7-day mean WBGT')\n",
"\n",
"ax[0].set_title(\"Historical, 2005-2014\")\n",
"ax[1].set_title(\"SSP1-2.6, 2091-2100\")\n",
"ax[2].set_title(\"SSP2-4.5, 2091-2100\")\n",
"ax[3].set_title(\"SSP5-8.5, 2091-2100\")\n",
"\n",
"[axi.set_facecolor(\"grey\") for axi in ax]\n",
"[axi.coastlines(lw=0.5) for axi in ax]\n",
"\n",
"fig.suptitle(\"Shaded regions for most extreme heatwave is > 31 WBGT\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Now we will examine the changes over time in the portion of the Earth's land surface affected by extreme heatwaves. To accomplish this, we utilize the previously loaded grid box area data.\n",
"\n",
"Next, we assign a value of \"1\" to the temporal-spatial data if it surpasses the defined threshold, and a value of \"0\" if it does not. By conducting an area-weighted average across the entire land area of the world, we determine the fraction of land area experiencing a heatwave above the threshold for each year."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"((wetbulb_hist > 31) * 1).weighted(area_land).mean([\"lon\", \"lat\"]).plot.line(\n",
" \"k.-\", label=\"historical\", ax=ax\n",
")\n",
"((wetbulb_126 > 31) * 1).weighted(area_land).mean([\"lon\", \"lat\"]).plot.line(\n",
" \".-\", label=\"ssp126\", ax=ax\n",
")\n",
"((wetbulb_245 > 31) * 1).weighted(area_land).mean([\"lon\", \"lat\"]).plot.line(\n",
" \".-\", label=\"ssp245\", ax=ax\n",
")\n",
"((wetbulb_585 > 31) * 1).weighted(area_land).mean([\"lon\", \"lat\"]).plot.line(\n",
" \".-\", label=\"ssp585\", ax=ax\n",
")\n",
"\n",
"ax.grid(True)\n",
"\n",
"ax.legend()\n",
"ax.set_xlabel(\"Time (years)\")\n",
"ax.set_ylabel(\"Land Area Fraction\")\n",
"ax.set_title(\"Fraction of land area\\nwith 7 days of wet-bulb temperature > 31 degrees\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"Let's average the fraction of the land area of the world that experiences a heatwave above wet-bulb temperature of 31 in the last 10 years of each run."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"tags": []
},
"outputs": [],
"source": [
"print(\n",
" \"Fraction of the land area of the world that experiences a heatwave above wet-bulb temperature of 31 in the last 10 years of each run in percent:\"\n",
")\n",
"(\n",
" pd.Series(\n",
" index=[\"historical\", \"SSP-126\", \"SSP-245\", \"SSP-585\"],\n",
" data=[\n",
" ((wetbulb_hist > 31) * 1)\n",
" .weighted(area_land)\n",
" .mean([\"lon\", \"lat\"])\n",
" .isel(time=slice(-10, None))\n",
" .mean()\n",
" .values,\n",
" ((wetbulb_126 > 31) * 1)\n",
" .weighted(area_land)\n",
" .mean([\"lon\", \"lat\"])\n",
" .isel(time=slice(-10, None))\n",
" .mean()\n",
" .values,\n",
" ((wetbulb_245 > 31) * 1)\n",
" .weighted(area_land)\n",
" .mean([\"lon\", \"lat\"])\n",
" .isel(time=slice(-10, None))\n",
" .mean()\n",
" .values,\n",
" ((wetbulb_585 > 31) * 1)\n",
" .weighted(area_land)\n",
" .mean([\"lon\", \"lat\"])\n",
" .isel(time=slice(-10, None))\n",
" .mean()\n",
" .values,\n",
" ],\n",
" ).astype(float)\n",
" * 100\n",
").round(2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"### Questions 2.1\n",
"1. What observations can you make when examining the time evolution and comparing the different scenarios?\n",
"2. Do you think these results would change if you used a different threshold such as 28 or 33 degrees? Why would it be important to look at this?"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"execution": {}
},
"source": [
"[*Click for solution*](https://github.com/neuromatch/climate-course-content/tree/main/tutorials/W2D4_ClimateResponse-Extremes&Variability/solutions/W2D4_Tutorial8_Solution_43e64815.py)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Summary\n",
"In this tutorial, you learned what the **wet-bulb globe temperature** is and how it affects human health. You analyzed the likelihood of crossing critical thresholds under historical and future climate scenarios, using data from a specific climate model. You learned how to conduct a spatial GEV analysis and evaluated the potential human impact under extreme heat waves."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Resources\n",
"\n",
"The data for this tutorial was accessed through the [Pangeo Cloud platform](https://pangeo.io/cloud.html). \n",
" \n",
"This tutorial uses data from the simulations conducted as part of the [CMIP6](https://pcmdi.llnl.gov/CMIP6/) multi-model ensemble, in particular the models MPI-ESM1-2-HR. \n",
"\n",
"MPI-ESM1-2-HR was developed and the runs were conducted by the [Max Planck Institute for Meteorology](https://mpimet.mpg.de/en/homepage) in Hamburg, Germany. \n",
"\n",
"For references on particular model experiments see this [database](https://www.wdc-climate.de/ords/f?p=127:2).\n",
"\n",
"For more information on what CMIP is and how to access the data, please see this [page](https://github.com/neuromatch/climate-course-content/blob/main/tutorials/CMIP/CMIP_resource_bank.md)."
]
}
],
"metadata": {
"colab": {
"collapsed_sections": [],
"include_colab_link": true,
"name": "W2D4_Tutorial8",
"provenance": [],
"toc_visible": true
},
"kernel": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.19"
}
},
"nbformat": 4,
"nbformat_minor": 4
}