{ "cells": [ { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/neuromatch/climate-course-content/blob/main/tutorials/W2D3_ExtremesandVariability/instructor/W2D3_Tutorial4.ipynb)   \"Open" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Tutorial 4: Return Levels Using Normal and GEV Distributions\n", "\n", "**Week 2, Day 3, 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, Agustina Pesce, 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", "Now that we have learned how to create a probability density function (PDF) for the generalized extreme value (GEV) distribution, we can utilize it to calculate return values based on this distribution.\n", "\n", "Return levels are computed from a PDF using the [cumulative distribution function (CDF)](https://en.wikipedia.org/wiki/Cumulative_distribution_function), which provides the probability that a randomly drawn variable from our distribution will be less than a certain value X (for our specific variable). Extreme events are rarer, and thus the probability of a precipitation event being at or lower than it is high. For example, if there is a 99% chance of observing a storm with 80mm of rainfall or lower, it means that there is a 1% chance of observing a storm with at least 80mm of rainfall. This implies that the 100-year storm would bring 80mm of rain. In simple terms, the return level is the inverse of the CDF.\n", "\n", "By the end of this tutorial, you will be able to:\n", "\n", "- Estimate return levels using a given quantile and the parameters of a GEV distribution.\n", "- Compare return level plots for GEV and normal distributions." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Setup" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 2204, "status": "ok", "timestamp": 1681924180496, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "# imports\n", "import xarray as xr\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import os\n", "import pooch\n", "import tempfile\n", "from scipy import stats\n", "from scipy.stats import genextreme as gev" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Install and import feedback gadget\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Install and import feedback gadget\n", "\n", "!pip3 install vibecheck datatops --quiet\n", "\n", "from vibecheck import DatatopsContentReviewContainer\n", "def content_review(notebook_section: str):\n", " return DatatopsContentReviewContainer(\n", " \"\", # No text prompt\n", " notebook_section,\n", " {\n", " \"url\": \"https://pmyvdlilci.execute-api.us-east-1.amazonaws.com/klab\",\n", " \"name\": \"comptools_4clim\",\n", " \"user_key\": \"l5jpxuee\",\n", " },\n", " ).render()\n", "\n", "\n", "feedback_prefix = \"W2D3_T4\"" ] }, { "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/W2D3_ExtremesandVariability\" # 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: Return Levels Using Distributions\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "remove-input" ] }, "outputs": [], "source": [ "# @title Video 1: Return Levels Using Distributions\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', 'DZwh4H2TqDs'), ('Bilibili', 'BV1Bh4y1Z7gu')]\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": "markdown", "metadata": {}, "source": [ "## Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Return_Levels_Using_Distributions_Video\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n" ] }, { "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 = \"sdezg\"\n", "\n", "print(f\"If you want to download the slides: https://osf.io/download/{link_id}/\")\n", "IFrame(src=f\"https://mfr.ca-1.osf.io/render?url=https://osf.io/{link_id}/?direct%26mode=render%26action=download%26mode=render\", width=854, height=480)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Return_Levels_Using_Distributions_Slides\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Section 1: Return Levels and Return Period" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "As before let's load the annual maximum precipitation data from Germany:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 530, "status": "ok", "timestamp": 1681924181024, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "# download file: 'precipitationGermany_1920-2022.csv'\n", "filename_precipitationGermany = \"precipitationGermany_1920-2022.csv\"\n", "url_precipitationGermany = \"https://osf.io/xs7h6/download\"\n", "data = pd.read_csv(\n", " pooch_load(url_precipitationGermany, filename_precipitationGermany), index_col=0\n", ").set_index(\"years\")\n", "data.columns = [\"precipitation\"]\n", "precipitation = data.precipitation" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "And fit the Generalized Extreme Value (GEV) distribution:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 267, "status": "ok", "timestamp": 1681924238489, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "shape, loc, scale = gev.fit(precipitation.values, 0)\n", "print(f\"Fitted parameters:\\nShape: {shape:.5f}, Location: {loc:.5f}, Scale: {scale:.5f}\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "To compute return levels, you can use the function ```estimate_return_level(x, location, scale, shape)``` defined below. In this case, x represents the _quantile_, the probability of a random value from our distribution being lower. For example, for the 100-year storm, x would be 0.99, and for the 1000-year storm, it would be 0.999.\n", "\n", "This utility function takes a given *quantile* and the GEV parameters (loc, scale, shape) and computes the corresponding return level." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "def estimate_return_level(quantile, loc, scale, shape):\n", " level = loc + scale / shape * (1 - (-np.log(quantile)) ** (shape))\n", " return level" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now, let's utilize this function to calculate the 2-year return level, which represents the precipitation level we anticipate with a 50% chance each year." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 7, "status": "ok", "timestamp": 1681924241704, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "estimate_return_level(0.5, loc, scale, shape)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "This function is also available as part of the scipy GEV implementation - here called the \"Percent point function\" or ```ppf``` which you used previously to make the quantile-quantile plots in the last tutorial." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "gev.ppf(0.5, shape, loc=loc, scale=scale)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "To make things easier, let's replace x with the formula $1-(1/\\text{return period})$. Note this is different than the exceedance probability discussed in Tutorial 2. So for a 100-year storm, we would use $1-(1/100) = 0.99$.\n", "\n", "Let us compute return levels for the 100- and 1000-year storms. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 188, "status": "ok", "timestamp": 1681924244959, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "estimate_return_level(1 - 1 / 100, loc, scale, shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 2, "status": "ok", "timestamp": 1681924245377, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "estimate_return_level(1 - 1 / 1000, loc, scale, shape)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now we will make a figure in which we plot return levels against return periods: \n", "1. Create a vector “periods” which is a sequence of 2 to 1000 years (with a step size of 2).\n", "2. Calculate the associated return levels for those return periods.\n", "3. Plot the return levels against the return periods. Typically, the return periods go on the x-axis at a log scale.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 784, "status": "ok", "timestamp": 1681924247138, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "periods = np.arange(2, 1000, 2)\n", "quantiles = 1 - 1 / periods\n", "levels = estimate_return_level(quantiles, loc, scale, shape)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(periods, levels, \".-\")\n", "ax.set_xlabel(\"Return Period (years)\")\n", "ax.set_ylabel(\"Return Level (mm/day)\")\n", "ax.set_xscale(\"log\")\n", "ax.grid(True)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "
\n", " Click to see a description of the plot \n", "Return levels versus return periods calculated from fitted GEV parameters and equal-spaced periods, log scaled on the x-axis. Hence, we can easily access the return level for an X-year event, e.g. 20 year-event corresponds to an annual maximum precipitation of 70 millimeters per day.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Coding Exercise 1\n", "\n", "Compute and plot the 50-, 100-, 500- and 1000-year return levels based on a normal distribution and based on the Normal and GEV distributions.\n", "\n", "Plot these three things:\n", "1. the empirical return level\n", "2. the estimate based on a normal distribution\n", "3. the estimate based on a GEV distribution\n", "\n", "*Note that you can get the empirical return levels from tutorial 2, we will provide this code below.*" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "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" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "# setup plots\n", "fig, ax = plt.subplots()\n", "\n", "# get empirical return levels and plot them\n", "_ = ...\n", "\n", "# create vector of years\n", "years = np.arange(1.1, 100, 0.1)\n", "\n", "# calculate and plot the normal return levels\n", "_ = ...\n", "\n", "# calculate and plot the GEV distribution, note the negtive shape parameter\n", "_ = ...\n", "\n", "# set x axis to log scale\n", "ax.set_xscale(\"log\")\n", "\n", "# show legend\n", "ax.legend([\"empirical\", \"normal\", \"GEV\"])\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 357, "status": "ok", "timestamp": 1681924249195, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# setup plots\n", "fig, ax = plt.subplots()\n", "\n", "# get empirical return levels and plot them\n", "_ = empirical_return_level(precipitation).plot(ax=ax, marker=\".\", linestyle=\"None\")\n", "\n", "# create vector of years\n", "years = np.arange(1.1, 100, 0.1)\n", "\n", "# calculate and plot the normal return levels\n", "ax.plot(\n", " years,\n", " stats.norm.ppf(1 - 1 / years, loc=precipitation.mean(), scale=precipitation.std()),\n", ")\n", "\n", "# calculate and plot the GEV distribution, note the negtive shape parameter\n", "ax.plot(years, gev.ppf(1 - 1 / years, shape, loc=loc, scale=scale))\n", "# set x axis to log scale\n", "ax.set_xscale(\"log\")\n", "\n", "# show legend\n", "ax.legend([\"empirical\", \"normal\", \"GEV\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Coding_Exercise_1\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Question 1\n", "\n", "1. What can you say about the plot and how the distributions describe the data? How do the distributions differ? At short/long return periods? What will happen at even longer return periods?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\"\"\"\n", "Put your observations here. GEV distribution describes the data better than the normal distribution. However, they are similar for shorter return periods. This suggests that for common, less extreme events (short return periods), both the Normal and GEV distributions provide a similar representation of the data. This similarity is likely due to the fact that for less extreme events, the tail of the distribution (which captures extreme events) is less relevant. When we move toards longer return periods (which correspond to more extreme events), the two curves start to diverge.\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Questions_1\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Bonus: Find Confidence Intervals\n", "\n", "When computing return levels of extreme events - or really when doing any quantitative assessment in science - it is useful and essential to provide uncertainty estimates, that is how much the result may be expected to fluctuate around the central estimate. \n", "\n", "In Tutorial 2 the exercise discussed one approach to construct such uncertainty ranges - bootstrapping. You can now extend this method by resampling the data and re-fitting the GEV, and in doing so generate uncertainty estimates of the GEV parameters and GEV-based return levels.\n", "\n", "One possible approach could involve the following steps:\n", "\n", "1. Begin with the `gev` class which can be employed to fit GEV parameters using the `gev.fit(data)` method. This function returns the three parameters, taking note of a sign change in the shape parameter (as different conventions exist).\n", "2. Determine the number of observations in the data, denoted as N.\n", "3. Perform resampling by randomly drawing N samples from the data with replacement. This process generates an artificial ensemble that differs from the true data due to resampling.\n", "4. Estimate the parameters for each resampling.\n", "5. Utilize the `gev.ppf()` function for each parameter set to compute the return level.\n", "6. Visualize the resulting uncertainty bands by plotting multiple lines or calculating confidence intervals using `np.quantile()`." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "# initalize list to store parameters from samples\n", "params = []\n", "\n", "# generate 1000 samples by resampling data with replacement\n", "for i in range(1000):\n", " ...\n", "\n", "# print the estimate of the mean of each parameter and it's confidence intervals\n", "print(\n", " \"Mean estimate: \",\n", " np.mean(np.array(params), axis=0),\n", " \" and 95% confidence intervals: \",\n", " ...,\n", ")\n", "\n", "# generate years vector\n", "years = np.arange(1.1, 1000, 0.1)\n", "\n", "# intialize list for return levels\n", "levels = []\n", "\n", "# calculate return levels for each of the 1000 samples\n", "for i in range(1000):\n", " levels.append(...)\n", "levels = np.array(levels)\n", "\n", "# setup plots\n", "fig, ax = plt.subplots()\n", "\n", "# find empirical return levels\n", "_ = ...\n", "\n", "# plot return mean levels\n", "_ = ...\n", "\n", "# plot confidence intervals\n", "_ = ...\n", "\n", "# aesthetics\n", "ax.set_xlim(1.5, 1000)\n", "ax.set_ylim(20, 110)\n", "ax.set_xscale(\"log\")\n", "ax.set_xlabel(\"Return Period (years)\")\n", "ax.set_ylabel(\"Return Level (mm/day)\")\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# initalize list to store parameters from samples\n", "params = []\n", "\n", "# generate 1000 samples by resampling data with replacement\n", "for i in range(1000):\n", " params.append(\n", " gev.fit(np.random.choice(precipitation, size=precipitation.size, replace=True))\n", " )\n", "\n", "# print the estimate of the mean of each parameter and it's confidence intervals\n", "print(\n", " \"Mean estimate: \",\n", " np.mean(np.array(params), axis=0),\n", " \" and 95% confidence intervals: \",\n", " np.quantile(np.array(params), [0.025, 0.975], axis=0),\n", ")\n", "\n", "# generate years vector\n", "years = np.arange(1.1, 1000, 0.1)\n", "\n", "# intialize list for return levels\n", "levels = []\n", "\n", "# calculate return levels for each of the 1000 samples\n", "for i in range(1000):\n", " levels.append(gev.ppf(1 - 1 / years, *params[i]))\n", "levels = np.array(levels)\n", "\n", "# setup plots\n", "fig, ax = plt.subplots()\n", "\n", "# find empirical return levels\n", "_ = empirical_return_level(precipitation).plot(ax=ax, marker=\".\", linestyle=\"None\")\n", "\n", "# plot return mean levels\n", "_ = ax.plot(years, levels.mean(axis=0))\n", "\n", "# plot confidence intervals\n", "_ = ax.plot(years, np.quantile(levels, [0.025, 0.975], axis=0).T, \"k--\")\n", "\n", "# aesthetics\n", "ax.set_xlim(1.5, 1000)\n", "ax.set_ylim(20, 110)\n", "ax.set_xscale(\"log\")\n", "ax.set_xlabel(\"Return Period (years)\")\n", "ax.set_ylabel(\"Return Level (mm/day)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Bonus_Coding_Exercise_1\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Summary\n", "In this tutorial, you learned how to calculate return values based on the Generalized Extreme Value (GEV) distribution. You learned how they are derived from the cumulative distribution function (CDF), and compared them to those from the normal distribution. Finally, you learned how to create confidence intervals for the parameters that define the GEV distribution for a given dataset." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Resources\n", "\n", "Data from this tutorial uses the 0.25 degree precipitation dataset E-OBS. It combines precipitation observations to generate a gridded (i.e. no \"holes\") precipitation over Europe. We used the precipitation data from the gridpoint at 51 N, 6 E. \n", "\n", "The dataset can be accessed using the KNMI Climate Explorer [here](https://climexp.knmi.nl/select.cgi?id=someone@somewhere&field=ensembles_025_rr). The Climate Explorer is a great resource to access, manipulate and visualize climate data, including observations and climate model simulations. It is freely accessible - feel free to explore!" ] } ], "metadata": { "colab": { "collapsed_sections": [], "include_colab_link": true, "name": "W2D3_Tutorial4", "toc_visible": true, "version": "" }, "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 }