{ "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_Tutorial5.ipynb)   \"Open" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Tutorial 5: Non-stationarity in Historical Records\n", "\n", "**Week 2, Day 3, Extremes & Variability**\n", "\n", "**Content creators:** Matthias Aengenheyster, Joeri Reinders\n", "\n", "**Content reviewers:** Yosemley Bermúdez, 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:* 25 minutes\n", "\n", "In this tutorial, we will analyze the annual maximum sea level heights from a measurement station near Washington DC. Coastal storms, particularly when combined with high tides, can result in exceptionally high sea levels, posing significant challenges for coastal cities. Understanding the magnitude of extreme events, such as the X-year storm, is crucial for implementing effective safety measures.\n", "\n", "By the end of this tutorial, you will gain the following abilities:\n", "\n", "- Analyze time series data during different climate normal periods.\n", "- Evaluate changes in statistical moments and parameter values over time to identify non-stationarity." ] }, { "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 cartopy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 1430, "status": "ok", "timestamp": 1681921157378, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "# imports\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import pandas as pd\n", "from scipy import stats\n", "from scipy.stats import genextreme as gev\n", "import os\n", "import pooch\n", "import tempfile" ] }, { "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_T5\"" ] }, { "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='d1725e1e5c6ccd7561644771c9ad17ab07e553c05b4156a38cb6780aaca36edb',\n", " fname=os.path.join(user_temp_cache, filename),\n", " processor=processor,\n", " )\n", "\n", " return file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Video 1: Non-Stationarity\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "remove-input" ] }, "outputs": [], "source": [ "# @title Video 1: Non-Stationarity\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', 'JAzC8_lpCxw'), ('Bilibili', 'BV1pV4y1h7bB')]\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}_Non_Stationarity_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 = \"e6vgz\"\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": { "execution": {} }, "source": [ "# Section 1: Washington DC's Maximum Sea Surface Height\n", "Let's inspect the annual maximum sea surface height data and create a plot over time. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 472, "status": "ok", "timestamp": 1681921184150, "user": { "displayName": "Matthias Aengenheyster", "userId": "16322208118439170907" }, "user_tz": -60 }, "tags": [] }, "outputs": [], "source": [ "# download file: 'WashingtonDCSSH1930-2022.csv'\n", "\n", "filename_WashingtonDCSSH1 = \"WashingtonDCSSH1930-2022.csv\"\n", "url_WashingtonDCSSH1 = \"https://osf.io/4zynp/download\"\n", "\n", "data = pd.read_csv(\n", " pooch_load(url_WashingtonDCSSH1, filename_WashingtonDCSSH1), index_col=0\n", ").set_index(\"years\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "data" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 1097, "status": "ok", "timestamp": 1680721204016, "user": { "displayName": "Yosmely Tamira Bermudez Gutierrez", "userId": "07776907551108334395" }, "user_tz": 180 }, "tags": [] }, "outputs": [], "source": [ "data.ssh.plot(\n", " linestyle=\"-\", marker=\".\", xlabel=\"Time (years)\", ylabel=\"Annual Maximum Sea Surface Height (mm)\", grid=True\n", ")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "
\n", " Click here for a description of the plot \n", "Time series plot of the annual maximum sea surface height in millimeters versus time in years from a measurement station near Washington DC. The lowest annual maximum was registered in 1931, and the largest in 2018, an increasing approximately linear trend is visible, which can be attributed to the rising sea surface temperatures caused by climate change. It is important to consider this non-stationarity when analyzing the data.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "In previous tutorials, we assumed that the probability density function (PDF) shape remains constant in time. In other words, the precipitation values are derived from the same distribution regardless of the timeframe. However, in today's world heavily influenced by climate change, we cannot assume that the PDF remains stable. For instance, global temperature is increasing, which causes a shift in the distribution's location. Additionally, local precipitation patterns are becoming more variable, leading to a widening of the distribution. Moreover, extreme events are becoming more severe, resulting in thicker tails of the distribution. We refer to this phenomenon as **non-stationarity**.\n", "\n", "To further investigate this, we can group our data into three 30-year periods known as \"**climate normals**\". We can create a first record for the 1931 to 1960 period, a second for 1961 to 1990, and a third for 1991 to 2020, respectively. By plotting the histogram of each dataset within the same frame, we can gain a more comprehensive understanding of the changes over time.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# 1931-1960\n", "data_period1 = data.iloc[0:30]\n", "# 1961-1990\n", "data_period2 = data.iloc[30:60]\n", "# 1990-2020\n", "data_period3 = data.iloc[60:90]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# plot the histograms for each climate normal identified above\n", "fig, ax = plt.subplots()\n", "colors =[\"C0\",\"C1\",\"C2\"]\n", "\n", "# loop over all periods of 30 years\n", "for ind, yr in enumerate(range(30,120,30)):\n", " sns.histplot(\n", " # select climate normal periods from dataset\n", " data.iloc[yr-30:yr].ssh,\n", " bins=np.arange(-200, 650, 50),\n", " color=colors[ind],\n", " element=\"step\",\n", " alpha=0.5,\n", " # estimate PDF via kernel density estimation\n", " kde=True,\n", " label=f\"{1901+yr}-{1930+yr}\",\n", " ax=ax,\n", ")\n", "\n", "# aesthetics\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Sea Surface Height (mm)\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "
\n", " Click here for a description of the plot \n", "Histogram of the annual maximum sea surface height in millimeters and a respective kernel density estimate of the PDF for three 30-year periods of the Washington DC dataset. The mean increases and the distribution becomes wider over time.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Let's also calculate the moments of each climate normal period:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 561, "status": "ok", "timestamp": 1680721216790, "user": { "displayName": "Yosmely Tamira Bermudez Gutierrez", "userId": "07776907551108334395" }, "user_tz": 180 }, "tags": [] }, "outputs": [], "source": [ "# setup pandas dataframe\n", "periods_stats = pd.DataFrame(index=[\"Mean\", \"Standard Deviation\", \"Skew\"])\n", "\n", "# add info for each climate normal period\n", "periods_stats[\"1931-1960\"] = [\n", " data_period1.ssh.mean(),\n", " data_period1.ssh.std(),\n", " data_period1.ssh.skew(),\n", "]\n", "periods_stats[\"1961-1990\"] = [\n", " data_period2.ssh.mean(),\n", " data_period2.ssh.std(),\n", " data_period2.ssh.skew(),\n", "]\n", "periods_stats[\"1991-2020\"] = [\n", " data_period3.ssh.mean(),\n", " data_period3.ssh.std(),\n", " data_period3.ssh.skew(),\n", "]\n", "\n", "periods_stats = periods_stats.T\n", "periods_stats" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "The mean increases as well as the standard deviation. Conversely, the skewness remains relatively stable over time, just decreasing slightly. This observation indicates that the dataset is non-stationary. To visualize the overall shape of the distribution changes, we can fit a Generalized Extreme Value (GEV) distribution to the data for each period and plot the corresponding PDF." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "tags": [] }, "outputs": [], "source": [ "# 1931-1960\n", "params_period1 = gev.fit(data_period1.ssh.values, 0)\n", "shape_period1, loc_period1, scale_period1 = params_period1\n", "\n", "# 1961-1990\n", "params_period2 = gev.fit(data_period2.ssh.values, 0)\n", "shape_period2, loc_period2, scale_period2 = params_period2\n", "\n", "# 1991-2020\n", "params_period3 = gev.fit(data_period3.ssh.values, 0)\n", "shape_period3, loc_period3, scale_period3 = params_period3" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# initialize list to store climate normals\n", "data_periods = []\n", "\n", "fig, ax = plt.subplots()\n", "# initialize array of possible annual maximum sea surface heights in mm\n", "x = np.linspace(-200, 600, 1000)\n", "\n", "# loop over three 30-year periods\n", "for ind, yr in enumerate(range(30,120,30)):\n", " # select climate normal periods from dataset\n", " data_periods.append(data.iloc[yr-30:yr])\n", " # fit the GEV distribution parameter\n", " shape_period, loc_period, scale_period = gev.fit(data_periods[ind].ssh.values, 0)\n", " # plot GEV distribution with these parameters\n", " ax.plot(\n", " x,\n", " gev.pdf(x, shape_period, loc=loc_period, scale=scale_period),\n", " color=colors[ind],\n", " linewidth=3,\n", " # add climate normal period label\n", " label=f\"{1901+yr}-{1930+yr}\"\n", ")\n", "# plot aesthetics\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Sea Surface Height (mm)\")\n", "ax.set_ylabel(\"Density\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "data_periods = []\n", "for ind, yr in enumerate(range(30,120,30)):\n", " data_periods.append(data.iloc[yr-30:yr])\n", "\n", "fig, ax = plt.subplots()\n", "x = np.linspace(-200, 600, 1000)\n", "\n", "for ind, yr in enumerate(range(30,120,30)):\n", " shape_period, loc_period, scale_period = gev.fit(data_periods[ind].ssh.values, 0)\n", " ax.plot(\n", " x,\n", " gev.pdf(x, shape_period, loc=loc_period, scale=scale_period),\n", " c=colors[ind],\n", " lw=3,\n", " label=f\"{1901+yr}-{1930+yr}\"\n", ")\n", "\n", "ax.legend()\n", "ax.set_xlabel(\"Annual Maximum Sea Surface Height (mm)\")\n", "ax.set_ylabel(\"Density\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now, let's examine the changes in the GEV parameters. This analysis will provide valuable insights into how we can incorporate non-stationarity into our model in one of the upcoming tutorials." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Question 1\n", "\n", "1. Look at the plot above. Just by visual inspection, describe how the distribution changes between the periods. Which parameters of the GEV distribution do you think are responsible? How and why? " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\"\"\"\n", "1. There is a relatively large increase in the location parameter over time, and the scale parameter from 1961-1990 and 1991-2020 also increases. This is likely due to sea level rise. Furthermore, the variability around the mean is increasing. The shape parameter also increases over this time but as commented in the lecture this parameter is not usually reliably estimated from such a short time span.\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": [ "## Coding Exercise 1\n", "\n", "1. Compare the location, scale and shape parameters of the fitted distribution for the three time periods. How do they change? Compare with your answers to the question above." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "# setup dataframe with titles for each parameter\n", "parameters = pd.DataFrame(index=[\"Location\", \"Scale\", \"Shape\"])\n", "\n", "# add in 1931-1960 parameters\n", "parameters[\"1931-1960\"] = ...\n", "\n", "# add in 1961-1990 parameters\n", "parameters[\"1961-1990\"] = ...\n", "\n", "# add in 1991-202 parameters\n", "parameters[\"1991-2020\"] = ...\n", "\n", "# transpose the dataset so the time periods are rows\n", "parameters = ...\n", "\n", "# round the values for viewing\n", "_ = ...\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "executionInfo": { "elapsed": 9, "status": "ok", "timestamp": 1680721280152, "user": { "displayName": "Yosmely Tamira Bermudez Gutierrez", "userId": "07776907551108334395" }, "user_tz": 180 }, "tags": [] }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# setup dataframe with titles for each parameter\n", "parameters = pd.DataFrame(index=[\"Location\", \"Scale\", \"Shape\"])\n", "\n", "# add in 1931-1960 parameters\n", "parameters[\"1931-1960\"] = [loc_period1, scale_period1, shape_period1]\n", "parameters[\"1931-1960\"] = gev.fit(data_period1.ssh.values, 0)\n", "\n", "# add in 1961-1990 parameters\n", "parameters[\"1961-1990\"] = [loc_period2, scale_period2, shape_period2]\n", "parameters[\"1961-1990\"] = gev.fit(data_period2.ssh.values, 0)\n", "\n", "# add in 1991-202 parameters\n", "parameters[\"1991-2020\"] = [loc_period3, scale_period3, shape_period3]\n", "parameters[\"1991-2020\"] = gev.fit(data_period3.ssh.values, 0)\n", "\n", "# transpose the dataset so the time periods are rows\n", "parameters = parameters.T\n", "\n", "# round the values for viewing\n", "_ = parameters.round(4) # .astype('%.2f')" ] }, { "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": [ "# Summary\n", "In this tutorial, you focused on the analysis of annual maximum sea surface heights in Washington DC, specifically considering the impact of **non-stationarity** due to climate change. You've learned how to analyze time series data across different 30-year \"**climate normal**\" periods, and how to evaluate changes in statistical moments and parameter values over time to identify this non-stationarity. By segmenting our data into climate normal periods, you were able to compare changes in sea level trends over time and understand the increasing severity of extreme events. " ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Resources\n", "\n", "Original data from this tutorial can be found [here](https://climexp.knmi.nl/getsealev.cgi?id=someone@somewhere&WMO=360&STATION=WASHINGTON_DC&extraargs=). Note the data used in this tutorial were preprocessed to extract the annual maximum values." ] } ], "metadata": { "colab": { "collapsed_sections": [], "include_colab_link": true, "name": "W2D3_Tutorial5", "provenance": [ { "file_id": "1q9gfBO_7bG3xrQfyFisZAKq8lK_RN2BX", "timestamp": 1680619878447 } ], "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 }