{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/ClimateMatchAcademy/course-content/blob/main/tutorials/W2D4_ClimateResponse-Extremes&Variability/W2D4_Tutorial6.ipynb)   \"Open" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Tutorial 6: Scenario-dependence of Future Changes in Extremes**\n", "\n", "**Week 2, Day 4, Extremes & Vulnerability**\n", "\n", "**Content creators:** Matthias Aengenheyster, Joeri Reinders\n", "\n", "**Content reviewers:** Younkap Nina Duplex, Sloane Garelick, Zahra Khodakaramimaghsoud, Peter Ohue, Laura Paccini, Jenna Pearson, Agustina Pesce, Derick Temfack, Peizhen Yang, Cheng Zhang, Chi Zhang, Ohad Zivan\n", "\n", "**Content editors:** Jenna Pearson, Chi Zhang, Ohad Zivan\n", "\n", "**Production editors:** Wesley Banfield, Jenna Pearson, Chi Zhang, Ohad Zivan\n", "\n", "**Our 2023 Sponsors:** NASA TOPS and Google DeepMind" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Tutorial Objectives**" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "In this tutorial, we will analyze climate model output for various cities worldwide to investigate the changes in extreme temperature and precipitation patterns over time under different emission scenarios.\n", "\n", "The data we will be using consists of climate model simulations for the historical period and three future climate scenarios:\n", "\n", "1. **Historical (hist)**: This scenario covers the time range from 1851 to 2014 and incorporates information about volcanic eruptions, greenhouse gas emissions, and other factors relevant to historical climate conditions.\n", "\n", "The Socio-economic pathway (SSP) scenarios represent potential climate futures beyond 2014. It's important to note that these scenarios are predictions (\"this is what we think will happen\") but are not certainties (\"this is plausible given the assumptions\"). Each scenario is based on different assumptions, primarily concerning the speed and effectiveness of global efforts to address global warming and reduce greenhouse gas emissions and other pollutants. Here are the scenarios we will be using today, along with desciptions [(see Section 1.6 of the recent IPCC AR6 WG1 report for more detail)](https://www.ipcc.ch/report/ar6/wg1/chapter/chapter-1/#1.6)\n", "\n", "2. **SSP1-2.6** *(ambitious climate scenario)*: This scenario stays below 2.0°C warming relative to 1850–1900 (median) with implied net zero CO2 emissions in the second half of the century.\n", "3. **SSP2-4.5** *(middle-ground climate scenario)*: This scenario is in line with the upper end of aggregate NDC (nationally-defined contribution) emissions levels by 2030. CO2 emissions remain around current levels until the middle of the century. The SSP2-4.5 scenario deviates mildly from a ‘no-additional-climate-policy’ reference scenario, resulting in a best-estimate warming around 2.7°C by the end of the 21st century relative to 1850–1900 (Chapter 4).\n", "4. **SSP5-8.5** *(pessimistic climate scenario)*: The highest emissions scenario with no additional climate policy and where CO2 emissions roughly double from current levels by 2050. Emissions levels as high as SSP5-8.5 are not obtained by integrated assessment models (IAMs) under any of the SSPs other than the fossil-fuelled SSP5 socio-economic development pathway. It exhibits the highest level of warming among all scenarios and is often used as a \"worst-case scenario\" (though it may not be the most likely outcome). It is worth noting that many experts today consider this scenario to be unlikely due to significant improvements in mitigation policies over the past decade or so.\n", "\n", "By the end of this tutorial, you will be able to:\n", "\n", "- Utilize climate model output from scenario runs to assess changes during the historical period.\n", "- Compare potential future climate scenarios, focusing on their impact on extreme events." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Setup**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 4521, "status": "ok", "timestamp": 1682954196975, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 } }, "outputs": [], "source": [ "# !pip install -q condacolab\n", "# import condacolab\n", "# condacolab.install()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 2969, "status": "ok", "timestamp": 1682954246617, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "# imports\n", "import xarray as xr\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import pandas as pd\n", "import cartopy.crs as ccrs\n", "from scipy import stats\n", "from scipy.stats import genextreme as gev\n", "from datetime import datetime\n", "import os\n", "import pooch\n", "import tempfile" ] }, { "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/ClimateMatchAcademy/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 estimate_return_level_period(period, loc, scale, shape):\n", " \"\"\"\n", " Compute GEV-based return level for a given return period, and GEV parameters\n", " \"\"\"\n", " return stats.genextreme.ppf(1 - 1 / period, shape, loc=loc, scale=scale)\n", "\n", "\n", "def empirical_return_level(data):\n", " \"\"\"\n", " Compute empirical return level using the algorithm introduced in Tutorial 2\n", " \"\"\"\n", " df = pd.DataFrame(index=np.arange(data.size))\n", " # sort the data\n", " df[\"sorted\"] = np.sort(data)[::-1]\n", " # rank via scipy instead to deal with duplicate values\n", " df[\"ranks_sp\"] = np.sort(stats.rankdata(-data))\n", " # find exceedence probability\n", " n = data.size\n", " df[\"exceedance\"] = df[\"ranks_sp\"] / (n + 1)\n", " # find return period\n", " df[\"period\"] = 1 / df[\"exceedance\"]\n", "\n", " df = df[::-1]\n", "\n", " out = xr.DataArray(\n", " dims=[\"period\"],\n", " coords={\"period\": df[\"period\"]},\n", " data=df[\"sorted\"],\n", " name=\"level\",\n", " )\n", " return out\n", "\n", "\n", "def fit_return_levels(data, years, N_boot=None, alpha=0.05):\n", " \"\"\"\n", " Fit GEV to data, compute return levels and confidence intervals\n", " \"\"\"\n", " empirical = (\n", " empirical_return_level(data)\n", " .rename({\"period\": \"period_emp\"})\n", " .rename(\"empirical\")\n", " )\n", " shape, loc, scale = gev.fit(data, 0)\n", " print(\"Location: %.1e, scale: %.1e, shape: %.1e\" % (loc, scale, shape))\n", " central = estimate_return_level_period(years, loc, scale, shape)\n", "\n", " out = xr.Dataset(\n", " # dims = ['period'],\n", " coords={\"period\": years, \"period_emp\": empirical[\"period_emp\"]},\n", " data_vars={\n", " \"empirical\": ([\"period_emp\"], empirical.data),\n", " \"GEV\": ([\"period\"], central),\n", " },\n", " )\n", "\n", " if N_boot:\n", " levels = []\n", " shapes, locs, scales = [], [], []\n", " for i in range(N_boot):\n", " datai = np.random.choice(data, size=data.size, replace=True)\n", " # print(datai.mean())\n", " shapei, loci, scalei = gev.fit(datai, 0)\n", " shapes.append(shapei)\n", " locs.append(loci)\n", " scales.append(scalei)\n", " leveli = estimate_return_level_period(years, loci, scalei, shapei)\n", " levels.append(leveli)\n", "\n", " levels = np.array(levels)\n", " quant = alpha / 2, 1 - alpha / 2\n", " quantiles = np.quantile(levels, quant, axis=0)\n", "\n", " print(\"Ranges with alpha=%.3f :\" % alpha)\n", " print(\"Location: [%.2f , %.2f]\" % tuple(np.quantile(locs, quant).tolist()))\n", " print(\"Scale: [%.2f , %.2f]\" % tuple(np.quantile(scales, quant).tolist()))\n", " print(\"Shape: [%.2f , %.2f]\" % tuple(np.quantile(shapes, quant).tolist()))\n", "\n", " quantiles = xr.DataArray(\n", " dims=[\"period\", \"quantiles\"],\n", " coords={\"period\": out.period, \"quantiles\": np.array(quant)},\n", " data=quantiles.T,\n", " )\n", " out[\"range\"] = quantiles\n", " return out\n", "\n", "\n", "def plot_return_levels(obj, c=\"C0\", label=\"\", ax=None):\n", " \"\"\"\n", " Plot fitted data:\n", " - empirical return level\n", " - GEV-fitted return level\n", " - alpha-confidence ranges with bootstrapping (if N_boot is given)\n", " \"\"\"\n", " if not ax:\n", " ax = plt.gca()\n", " obj[\"GEV\"].plot.line(\"%s-\" % c, lw=3, _labels=False, label=label, ax=ax)\n", " obj[\"empirical\"].plot.line(\"%so\" % c, mec=\"k\", markersize=5, _labels=False, ax=ax)\n", " if \"range\" in obj:\n", " # obj['range'].plot.line('k--',hue='quantiles',label=obj['quantiles'].values)\n", " ax.fill_between(obj[\"period\"], *obj[\"range\"].T, alpha=0.3, lw=0, color=c)\n", " ax.semilogx()\n", " ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Video 1: Speaker Introduction\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Video 1: Speaker Introduction\n", "# Tech team will add code to format and display the video" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# 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" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Section 1: Load CMIP6 Data**\n", "\n", "\n", "As in IPCC Physical Basis Day, you will be loading (Coupled Model Intercomparison Project 6) CMIP6 data from [Pangeo](https://pangeo.io/). In this way you can access large amounts of climate model output that has been stored in the cloud. Here, we have already accessed the data of interest and collected it into a .nc file for you. However, the information on how to access this data directly is provided in the Resources section at the end of this notebook.\n", "\n", "Here you will load precipitation and maximum daily near-surface air temperature." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# download file: 'cmip6_data_city_daily_scenarios_tasmax_pr_models.nc'\n", "filename_cmip6_data = \"cmip6_data_city_daily_scenarios_tasmax_pr_models.nc\"\n", "url_cmip6_data = \"https://osf.io/ngafk/download\"\n", "\n", "data = xr.open_dataset(pooch_load(url_cmip6_data, filename_cmip6_data))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Section 2: Inspect Data**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "data" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "data.pr" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "data.tasmax" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "You can take a little time to explore the different variables within this dataset." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Section 2: Processing**" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Let's look at the data for one selected city, for one climate model. In this case here, we choose Madrid and the `MPI-ESM1-2-HR` model." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 441, "status": "ok", "timestamp": 1682954270647, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "city = \"Madrid\"\n", "data_city = data.sel(city=city, model=\"MPI-ESM1-2-HR\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 1456, "status": "ok", "timestamp": 1682954272100, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "data_city" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "The data has daily resolution, for three climate scenarios. Until 2014 the scenarios are identical (the 'historical' scenario). After 2014 they diverge given different climate change trajectories. Let's plot these two variables over time to get a sense of this." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 13213, "status": "ok", "timestamp": 1682954288444, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "# setup plot\n", "fig, ax = plt.subplots(2, sharex=True, figsize=(10, 5), constrained_layout=True)\n", "data_city[\"tasmax\"].plot(hue=\"scenario\", ax=ax[0], lw=0.5)\n", "data_city[\"pr\"].plot(hue=\"scenario\", ax=ax[1], lw=0.5, add_legend=False)\n", "\n", "ax[0].set_title(\"Maximum Daily Near-Surface Air Temperature\")\n", "ax[1].set_title(\"Precipitation\")\n", "\n", "ax[0].set_xlabel(\"\")\n", "\n", "ax[0].set_ylabel(\"K\")\n", "ax[1].set_ylabel(\"mm/day\")\n", "\n", "ax[0].set_xlim(datetime(1850, 1, 1), datetime(2100, 12, 31))\n", "\n", "ax[1].set_ylim(0, None);" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "In the previous tutorials we have been operating on annual maxima data - looking at the most extreme event observed in each year. We will do the same here: take the day in each year with the highest temperature or the largest amount of rainfall" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "fig, ax = plt.subplots(2, sharex=True, figsize=(10, 5), constrained_layout=True)\n", "data_city[\"tasmax\"].resample(time=\"1Y\").max().plot(hue=\"scenario\", ax=ax[0], lw=2)\n", "data_city[\"pr\"].resample(time=\"1Y\").max().plot(\n", " hue=\"scenario\", ax=ax[1], lw=2, add_legend=False\n", ")\n", "\n", "ax[0].set_title(\"Annual Maximum of Daily Maximum Near-Surface Air Temperature\")\n", "ax[1].set_title(\"Annual Maximum of Daily Precipitation\")\n", "\n", "ax[0].set_xlabel(\"\")\n", "\n", "ax[0].set_ylabel(\"K\")\n", "ax[1].set_ylabel(\"mm/day\")\n", "\n", "ax[0].set_xlim(datetime(1850, 1, 1), datetime(2100, 12, 31))\n", "\n", "ax[1].set_ylim(0, None);" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## **Questions 2** \n", "1. Describe the plot - what do you see for the two variables, over time, between scenarios?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# to_remove explanation\n", "\"\"\"\n", "1. Possible answer: For near-surface air tempereature, there tends to be an increase from ssp 126 to ssp 585. For rainfall, the variability rather than the trend seems to be changing the most.\n", "\"\"\";" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Section 3: Differences Between Historical Periods**" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "As in the previous tutorial we want to compare consecutive 30-year periods in the past: therefore take the historical run (1850-2014), and look at the annual maximum daily precipitation for the last three 30-year periods. We only need to look at one scenario because they all use the historical run until 2014." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 971, "status": "ok", "timestamp": 1682954330425, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "# take max daily precip values from Madrid for three climate normal periods\n", "pr_city = data_city[\"pr\"]\n", "pr_city_max = pr_city.resample(time=\"1Y\").max()\n", "\n", "data_period1 = (\n", " pr_city_max.sel(scenario=\"ssp245\", time=slice(\"2014\"))\n", " .sel(time=slice(\"1925\", \"1954\"))\n", " .to_dataframe()[\"pr\"]\n", ")\n", "data_period2 = (\n", " pr_city_max.sel(scenario=\"ssp245\", time=slice(\"2014\"))\n", " .sel(time=slice(\"1955\", \"1984\"))\n", " .to_dataframe()[\"pr\"]\n", ")\n", "data_period3 = (\n", " pr_city_max.sel(scenario=\"ssp245\", time=slice(\"2014\"))\n", " .sel(time=slice(\"1985\", \"2014\"))\n", " .to_dataframe()[\"pr\"]\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Plot the histograms of annual maximum daily precipitation for the three climate normals covering the historical period. What do you see? Compare to the analysis in the previous tutorial where we analyzed sea level height. Any similarities or differences? Why do you think that is?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 444, "status": "ok", "timestamp": 1682954338685, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "# plot histograms for climate normals during historical period\n", "fig, ax = plt.subplots()\n", "sns.histplot(\n", " data_period1,\n", " bins=np.arange(20, 90, 5),\n", " color=\"C0\",\n", " element=\"step\",\n", " alpha=0.5,\n", " kde=True,\n", " label=\"1925-1954\",\n", " ax=ax,\n", ")\n", "sns.histplot(\n", " data_period2,\n", " bins=np.arange(20, 90, 5),\n", " color=\"C1\",\n", " element=\"step\",\n", " alpha=0.5,\n", " kde=True,\n", " label=\"1955-1984\",\n", " ax=ax,\n", ")\n", "sns.histplot(\n", " data_period3,\n", " bins=np.arange(20, 90, 5),\n", " color=\"C2\",\n", " element=\"step\",\n", " alpha=0.5,\n", " kde=True,\n", " label=\"1985-2014\",\n", " ax=ax,\n", ")\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Daily Precipitation (mm/day)\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 7, "status": "ok", "timestamp": 1682954338884, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "# calculate moments of the data\n", "periods_stats = pd.DataFrame(index=[\"Mean\", \"Standard Deviation\", \"Skew\"])\n", "periods_stats[\"1925-1954\"] = [\n", " data_period1.mean(),\n", " data_period1.std(),\n", " data_period1.skew(),\n", "]\n", "periods_stats[\"1955-1984\"] = [\n", " data_period2.mean(),\n", " data_period2.std(),\n", " data_period2.skew(),\n", "]\n", "periods_stats[\"1985-2014\"] = [\n", " data_period3.mean(),\n", " data_period3.std(),\n", " data_period3.skew(),\n", "]\n", "\n", "periods_stats = periods_stats.T\n", "periods_stats" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now, we fit a GEV to the three time periods, and plot the distributions using the gev.pdf function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# reminder of how to fit function works and the output shape parameters\n", "gev = stats.genextreme\n", "gev.fit(data_period1.values)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# fit GEV distribution for three climate normals during historical period\n", "params_period1 = gev.fit(data_period1, 0)\n", "shape_period1, loc_period1, scale_period1 = params_period1\n", "params_period2 = gev.fit(data_period2, 0)\n", "shape_period2, loc_period2, scale_period2 = params_period2\n", "params_period3 = gev.fit(data_period3, 0)\n", "shape_period3, loc_period3, scale_period3 = params_period3" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# plot corresponding PDFs of GEV distribution fits\n", "fig, ax = plt.subplots()\n", "x = np.linspace(20, 90, 1000)\n", "ax.plot(\n", " x,\n", " gev.pdf(x, shape_period1, loc=loc_period1, scale=scale_period1),\n", " c=\"C0\",\n", " lw=3,\n", " label=\"1925-1954\",\n", ")\n", "ax.plot(\n", " x,\n", " gev.pdf(x, shape_period2, loc=loc_period2, scale=scale_period2),\n", " c=\"C1\",\n", " lw=3,\n", " label=\"1955-1984\",\n", ")\n", "ax.plot(\n", " x,\n", " gev.pdf(x, shape_period3, loc=loc_period3, scale=scale_period3),\n", " c=\"C2\",\n", " lw=3,\n", " label=\"1985-2014\",\n", ")\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Daily Precipitation (mm/day)\")\n", "ax.set_ylabel(\"Density\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# show the parameters of the GEV fit\n", "parameters = pd.DataFrame(index=[\"Location\", \"Scale\", \"Shape\"])\n", "parameters[\"1925-1954\"] = [loc_period1, scale_period1, shape_period1]\n", "parameters[\"1955-1984\"] = [loc_period2, scale_period2, shape_period2]\n", "parameters[\"1985-2014\"] = [loc_period3, scale_period3, shape_period3]\n", "\n", "parameters = parameters.T\n", "parameters.round(4) # .astype('%.2f')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now we will create a return level plot for the three periods. To do so we will be using some helper functions defined at the beginning of the tutorial, most of which you have seen before. In particular we will use `fit_return_levels` to generate an xr.Dataset that contains empirical and GEV fits, as well as confidence intervals, and `plot_return_levels` to generate a plot from this xr.Dataset with calculated confidence intervals shaded (alpha printed below).\n", "\n", "These functions can also be found in gev_functions (`import gev_functions as gf`)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "fit_period1 = fit_return_levels(\n", " data_period1, np.arange(1.1, 1000, 0.1), N_boot=100, alpha=0.05\n", ")\n", "fit_period2 = fit_return_levels(\n", " data_period2, np.arange(1.1, 1000, 0.1), N_boot=100, alpha=0.05\n", ")\n", "fit_period3 = fit_return_levels(\n", " data_period3, np.arange(1.1, 1000, 0.1), N_boot=100, alpha=0.05\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "plot_return_levels(fit_period1, c=\"C0\", label=\"1925-1954\", ax=ax)\n", "plot_return_levels(fit_period2, c=\"C1\", label=\"1955-1984\", ax=ax)\n", "plot_return_levels(fit_period3, c=\"C2\", label=\"1985-2014\", ax=ax)\n", "ax.set_xlim(1.5, 1000)\n", "ax.set_ylim(40, 140)\n", "\n", "ax.legend()\n", "ax.set_ylabel(\"Return Level (mm/day)\")\n", "ax.set_xlabel(\"Return Period (years)\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## **Questions 3**\n", "\n", "1. What do you conclude for the historical change in extreme precipitation in this city? What possible limitations could this analysis have? (How) could we address this?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# to_remove explanation\n", "\"\"\"\n", "1. The changes over the historical period are relatively small, however there are larger return levels noted for the most recent climate normal from 1985-2014. One limitation is the length of time used, as the shape parameter is typically not well defined. Using a longer time interval and likely a better shape parameters could give us more confidence in our return level estimates.\n", "\"\"\";" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Section 4: Climate Futures**" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now let's look at maximum precipitation in possible climate futures: the years 2071-2100 (the last 30 years). For comparison we use the historical period, 1850-2014." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "In the next box we select the data, then you will plot it as a coding exercise:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 238, "status": "ok", "timestamp": 1682954382837, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "data_city = data.sel(city=city, model=\"MPI-ESM1-2-HR\")\n", "# data_city = data.sel(city=city,model='MIROC6')\n", "data_city_fut = (\n", " data_city[\"pr\"].sel(time=slice(\"2071\", \"2100\")).resample(time=\"1Y\").max()\n", ")\n", "\n", "# select the different time periods\n", "data_hist = (\n", " data_city[\"pr\"]\n", " .sel(scenario=\"ssp126\", time=slice(\"1850\", \"2014\"))\n", " .resample(time=\"1Y\")\n", " .max()\n", " .to_dataframe()[\"pr\"]\n", ")\n", "data_ssp126 = data_city_fut.sel(scenario=\"ssp126\").to_dataframe()[\"pr\"]\n", "data_ssp245 = data_city_fut.sel(scenario=\"ssp245\").to_dataframe()[\"pr\"]\n", "data_ssp585 = data_city_fut.sel(scenario=\"ssp585\").to_dataframe()[\"pr\"]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## **Coding Exercises 4**\n", "Differences between climate scenarios:\n", "\n", "Repeat the analysis that we did above for three different time periods, but now for three different climate scenarios, and the historical period for comparison.\n", "1. Create a figure that displays the histograms of the four records. Find a useful number and spacing of bins (via the bins= keyword to sns.histplot). Calculate the moments.\n", "2. Fit GEV distributions to the four records using the same commands as above. Use the `gev.pdf` function to plot the fitted distributions.\n", "3. Inspect location, scale and shape parameters\n", "4. Create a return-level plot using the `ef.plot_levels_from_obj` function." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "# setup plot\n", "fig, ax = plt.subplots()\n", "\n", "# create histograms for each scenario and historical\n", "# Uncomment and complete\n", "# sns.histplot(\n", "# ...,\n", "# bins=np.arange(30, 100, 5),\n", "# color=\"k\",\n", "# element=\"step\",\n", "# stat=\"density\",\n", "# alpha=0.3,\n", "# lw=0.5,\n", "# line_kws=dict(lw=3),\n", "# kde=True,\n", "# label=\"Historical, 1850-2014\",\n", "# ax=ax,\n", "# )\n", "# sns.histplot(\n", "# ...,\n", "# bins=np.arange(30, 100, 5),\n", "# color=\"C0\",\n", "# element=\"step\",\n", "# stat=\"density\",\n", "# alpha=0.3,\n", "# lw=0.5,\n", "# line_kws=dict(lw=3),\n", "# kde=True,\n", "# label=\"SSP-126, 2071-2100\",\n", "# ax=ax,\n", "# )\n", "# sns.histplot(\n", "# ...,\n", "# bins=np.arange(30, 100, 5),\n", "# color=\"C0\",\n", "# element=\"step\",\n", "# stat=\"density\",\n", "# alpha=0.3,\n", "# lw=0.5,\n", "# line_kws=dict(lw=3),\n", "# kde=True,\n", "# label=\"SSP-126, 2071-2100\",\n", "# ax=ax,\n", "# )\n", "# sns.histplot(\n", "# ...,\n", "# bins=np.arange(30, 100, 5),\n", "# color=\"C1\",\n", "# element=\"step\",\n", "# stat=\"density\",\n", "# alpha=0.3,\n", "# lw=0.5,\n", "# line_kws=dict(lw=3),\n", "# kde=True,\n", "# label=\"SSP-245, 2071-2100\",\n", "# ax=ax,\n", "# )\n", "# sns.histplot(\n", "# ...,\n", "# bins=np.arange(30, 100, 5),\n", "# color=\"C2\",\n", "# element=\"step\",\n", "# stat=\"density\",\n", "# alpha=0.3,\n", "# lw=0.5,\n", "# line_kws=dict(lw=3),\n", "# kde=True,\n", "# label=\"SSP-585, 2071-2100\",\n", "# ax=ax,\n", "# )\n", "\n", "# aesthetics\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Daily Precipitation (mm/day)\")\n", "\n", "# calculate moments\n", "periods_stats = pd.DataFrame(index=[\"Mean\", \"Standard Deviation\", \"Skew\"])\n", "periods_stats[...] = [..., ..., ...]\n", "periods_stats[...] = [..., ..., ...]\n", "periods_stats[...] = [..., ..., ...]\n", "periods_stats[...] = [..., ..., ...]\n", "periods_stats = periods_stats.T\n", "periods_stats\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# setup plot\n", "fig, ax = plt.subplots()\n", "\n", "# create histograms for each scenario and historical\n", "# Uncomment and complete\n", "sns.histplot(\n", " data_hist,\n", " bins=np.arange(30, 100, 5),\n", " color=\"k\",\n", " element=\"step\",\n", " stat=\"density\",\n", " alpha=0.3,\n", " lw=0.5,\n", " line_kws=dict(lw=3),\n", " kde=True,\n", " label=\"Historical, 1850-2014\",\n", " ax=ax,\n", ")\n", "sns.histplot(\n", " data_ssp126,\n", " bins=np.arange(30, 100, 5),\n", " color=\"C0\",\n", " element=\"step\",\n", " stat=\"density\",\n", " alpha=0.3,\n", " lw=0.5,\n", " line_kws=dict(lw=3),\n", " kde=True,\n", " label=\"SSP-126, 2071-2100\",\n", " ax=ax,\n", ")\n", "sns.histplot(\n", " data_ssp126,\n", " bins=np.arange(30, 100, 5),\n", " color=\"C0\",\n", " element=\"step\",\n", " stat=\"density\",\n", " alpha=0.3,\n", " lw=0.5,\n", " line_kws=dict(lw=3),\n", " kde=True,\n", " label=\"SSP-126, 2071-2100\",\n", " ax=ax,\n", ")\n", "sns.histplot(\n", " data_ssp245,\n", " bins=np.arange(30, 100, 5),\n", " color=\"C1\",\n", " element=\"step\",\n", " stat=\"density\",\n", " alpha=0.3,\n", " lw=0.5,\n", " line_kws=dict(lw=3),\n", " kde=True,\n", " label=\"SSP-245, 2071-2100\",\n", " ax=ax,\n", ")\n", "sns.histplot(\n", " data_ssp585,\n", " bins=np.arange(30, 100, 5),\n", " color=\"C2\",\n", " element=\"step\",\n", " stat=\"density\",\n", " alpha=0.3,\n", " lw=0.5,\n", " line_kws=dict(lw=3),\n", " kde=True,\n", " label=\"SSP-585, 2071-2100\",\n", " ax=ax,\n", ")\n", "\n", "# aesthetics\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Daily Precipitation (mm/day)\")\n", "\n", "# calculate moments\n", "periods_stats = pd.DataFrame(index=[\"Mean\", \"Standard Deviation\", \"Skew\"])\n", "periods_stats[\"hist\"] = [data_hist.mean(), data_hist.std(), data_hist.skew()]\n", "periods_stats[\"ssp126\"] = [data_ssp126.mean(), data_ssp126.std(), data_ssp126.skew()]\n", "periods_stats[\"ssp245\"] = [data_ssp245.mean(), data_ssp245.std(), data_ssp245.skew()]\n", "periods_stats[\"ssp585\"] = [data_ssp585.mean(), data_ssp585.std(), data_ssp585.skew()]\n", "periods_stats = periods_stats.T\n", "periods_stats" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "# fit GEV distribution\n", "# Uncomment and complete\n", "# shape_hist, loc_hist, scale_hist = ...\n", "# shape_ssp126, loc_ssp126, scale_ssp126 = ...\n", "# shape_ssp245, loc_ssp245, scale_ssp245 = ...\n", "# shape_ssp585, loc_ssp585, scale_ssp585 = ...\n", "\n", "# make plots\n", "fig, ax = plt.subplots()\n", "x = np.linspace(20, 120, 1000)\n", "# Uncomment and complete\n", "# ax.plot(\n", "# x,\n", "# gev.pdf(x, ..., loc=..., scale=...),\n", "# c=\"k\",\n", "# lw=3,\n", "# label=\"historical, 1850-2014\",\n", "# )\n", "# ax.plot(\n", "# x,\n", "# gev.pdf(x, ..., loc=..., scale=...),\n", "# c=\"C0\",\n", "# lw=3,\n", "# label=\"SSP-126, 2071-2100\",\n", "# )\n", "# ax.plot(\n", "# x,\n", "# gev.pdf(x, ..., loc=..., scale=...),\n", "# c=\"C1\",\n", "# lw=3,\n", "# label=\"SSP-245, 2071-2100\",\n", "# )\n", "# ax.plot(\n", "# x,\n", "# gev.pdf(x, ..., loc=..., scale=...),\n", "# c=\"C2\",\n", "# lw=3,\n", "# label=\"SSP-585, 2071-2100\",\n", "# )\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Daily Precipitation (mm/day)\")\n", "ax.set_ylabel(\"Density\");\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 42287, "status": "ok", "timestamp": 1682954432623, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -120 }, "tags": [] }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# fit GEV distribution\n", "shape_hist, loc_hist, scale_hist = gev.fit(data_hist)\n", "shape_ssp126, loc_ssp126, scale_ssp126 = gev.fit(data_ssp126)\n", "shape_ssp245, loc_ssp245, scale_ssp245 = gev.fit(data_ssp245)\n", "shape_ssp585, loc_ssp585, scale_ssp585 = gev.fit(data_ssp585)\n", "\n", "# make plots\n", "fig, ax = plt.subplots()\n", "x = np.linspace(20, 120, 1000)\n", "ax.plot(\n", " x,\n", " gev.pdf(x, shape_hist, loc=loc_hist, scale=scale_hist),\n", " c=\"k\",\n", " lw=3,\n", " label=\"historical, 1850-2014\",\n", ")\n", "ax.plot(\n", " x,\n", " gev.pdf(x, shape_ssp126, loc=loc_ssp126, scale=scale_ssp126),\n", " c=\"C0\",\n", " lw=3,\n", " label=\"SSP-126, 2071-2100\",\n", ")\n", "ax.plot(\n", " x,\n", " gev.pdf(x, shape_ssp245, loc=loc_ssp245, scale=scale_ssp245),\n", " c=\"C1\",\n", " lw=3,\n", " label=\"SSP-245, 2071-2100\",\n", ")\n", "ax.plot(\n", " x,\n", " gev.pdf(x, shape_ssp585, loc=loc_ssp585, scale=scale_ssp585),\n", " c=\"C2\",\n", " lw=3,\n", " label=\"SSP-585, 2071-2100\",\n", ")\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Daily Precipitation (mm/day)\")\n", "ax.set_ylabel(\"Density\");" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "# find return levels\n", "# Uncomment and complete\n", "# fit_hist = fit_return_levels(\n", "# ..., np.arange(1.1, 200, 0.1), N_boot=100, alpha=0.05\n", "# )\n", "fit_ssp126 = ...\n", "fit_ssp245 = ...\n", "fit_ssp585 = ...\n", "\n", "\n", "# plot\n", "fig, ax = plt.subplots()\n", "# Uncomment and complete\n", "# plot_return_levels(..., c=\"k\", label=\"historical, 1850-2014\", ax=ax)\n", "# plot_return_levels(..., c=\"C0\", label=\"SSP-126, 2071-2100\", ax=ax)\n", "# plot_return_levels(..., c=\"C1\", label=\"SSP-245, 2071-2100\", ax=ax)\n", "# plot_return_levels(..., c=\"C2\", label=\"SSP-585, 2071-2100\", ax=ax)\n", "ax.set_xlim(1, 200)\n", "ax.set_ylim(30, 110)\n", "ax.set_ylabel(\"Return Level (mm/day)\")\n", "ax.set_xlabel(\"Return Period (years)\")\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# find return levels\n", "fit_hist = fit_return_levels(\n", " data_hist, np.arange(1.1, 200, 0.1), N_boot=100, alpha=0.05\n", ")\n", "fit_ssp126 = fit_return_levels(\n", " data_ssp126, np.arange(1.1, 200, 0.1), N_boot=100, alpha=0.05\n", ")\n", "fit_ssp245 = fit_return_levels(\n", " data_ssp245, np.arange(1.1, 200, 0.1), N_boot=100, alpha=0.05\n", ")\n", "fit_ssp585 = fit_return_levels(\n", " data_ssp585, np.arange(1.1, 200, 0.1), N_boot=100, alpha=0.05\n", ")\n", "\n", "\n", "# plot\n", "fig, ax = plt.subplots()\n", "plot_return_levels(fit_hist, c=\"k\", label=\"historical, 1850-2014\", ax=ax)\n", "plot_return_levels(fit_ssp126, c=\"C0\", label=\"SSP-126, 2071-2100\", ax=ax)\n", "plot_return_levels(fit_ssp245, c=\"C1\", label=\"SSP-245, 2071-2100\", ax=ax)\n", "plot_return_levels(fit_ssp585, c=\"C2\", label=\"SSP-585, 2071-2100\", ax=ax)\n", "ax.set_xlim(1, 200)\n", "ax.set_ylim(30, 110)\n", "ax.set_ylabel(\"Return Level (mm/day)\")\n", "ax.set_xlabel(\"Return Period (years)\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## **Questions 4**\n", "\n", "1. What can you say about how extreme precipitation differs between the climate scenarios? Are the differences large or small compared to periods in the historical records? What are the limitations? Consider the x-axis in the return-level plot compared to the space covered by the data (only 30 years). How could we get more information for longer return periods?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\"\"\"\n", "The differences between the scenarios and the historical records are large for SSP 245 and SSP 558, while such difference is small for SSP 126. In general, we might expect larger differences as we move towards more extreme emission scenarios and further into the future. One limitation is the inherent uncertainty in any climate model. These models make assumptions about future emissions and other factors, which might not pan out exactly as predicted. Another limitation is the relatively short period covered by the data (only 30 years). To gain more information about longer return periods, you could use longer time series of data if available, for example from paleo-simulations, proxies, or reanalysis products.\n", "\"\"\";" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# **Summary**\n", "In this tutorial, you've learned how to analyze climate model output to investigate the changes in extreme precipitation patterns over time under various emission scenarios and the historical period. Specifically, we've focused on three future Shared Socioeconomic Pathways (SSPs) scenarios, which represent potential futures based on different assumptions about greenhouse gas emissions.\n", "\n", "You've explored how to fit Generalized Extreme Value (GEV) distributions to the data and used these fitted distributions to create return-level plots. These plots allow us to visualize the probability of extreme events under different climate scenarios." ] }, { "attachments": {}, "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). Additionally, a [notebook](https://github.com/ClimateMatchAcademy/course-content/blob/main/tutorials/W2D4_ClimateResponse-Extremes%26Variability/get_CMIP6_data_from_pangeo.ipynb) is available here that downloads the specific datasets we used in this tutorial.\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 and ACCESS-CM2. \n", "\n", "MPI-ESM1-2-HR was developed and the runs conducted by the [Max Planck Institute for Meteorology](https://mpimet.mpg.de/en/homepage) in Hamburg, Germany. \n", "ACCESS-CM2 was developed and the runs conducted by [CSIRO](https://www.csiro.au) (Commonwealth Scientific and Industrial Research Organisation) ARCCSS (Australian Research Council Centre of Excellence for Climate System Science), see also climateextremes.org.au. \n", "\n", "For references on particular model experiments see this [database](https://www.wdc-climate.de/ords/f?p=127:2)." ] } ], "metadata": { "colab": { "collapsed_sections": [], "include_colab_link": true, "name": "W2D4_Tutorial6", "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.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }