Open In Colab   Open in Kaggle

Tutorial 5: Non-stationarity in Historical Records#

Week 2, Day 4, Extremes & Variability

Content creators: Matthias Aengenheyster, Joeri Reinders

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

Content editors: Paul Heubel, Jenna Pearson, Chi Zhang, Ohad Zivan

Production editors: Wesley Banfield, Paul Heubel, Jenna Pearson, Chi Zhang, Ohad Zivan

Our 2024 Sponsors: NFDI4Earth, CMIP

Tutorial Objectives#

Estimated timing of tutorial: 25 minutes

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.

By the end of this tutorial, you will gain the following abilities:

  • Analyze time series data during different climate normal periods.

  • Evaluate changes in statistical moments and parameter values over time to identify non-stationarity.

Setup#

# imports
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from scipy import stats
from scipy.stats import genextreme as gev
import os
import pooch
import tempfile

Figure Settings#

Hide code cell source
# @title Figure Settings
import ipywidgets as widgets  # interactive display

%config InlineBackend.figure_format = 'retina'
plt.style.use(
    "https://raw.githubusercontent.com/neuromatch/climate-course-content/main/cma.mplstyle"
)

Helper functions#

Hide code cell source
# @title Helper functions


def pooch_load(filelocation=None, filename=None, processor=None):
    shared_location = "/home/jovyan/shared/Data/tutorials/W2D4_ClimateResponse-Extremes&Variability"  # this is different for each day
    user_temp_cache = tempfile.gettempdir()

    if os.path.exists(os.path.join(shared_location, filename)):
        file = os.path.join(shared_location, filename)
    else:
        file = pooch.retrieve(
            filelocation,
            known_hash=None,
            fname=os.path.join(user_temp_cache, filename),
            processor=processor,
        )

    return file

Video 1: Non-Stationarity#

Section 1#

Let’s inspect the annual maximum sea surface height data and create a plot over time.

# download file: 'WashingtonDCSSH1930-2022.csv'

filename_WashingtonDCSSH1 = "WashingtonDCSSH1930-2022.csv"
url_WashingtonDCSSH1 = "https://osf.io/4zynp/download"

data = pd.read_csv(
    pooch_load(url_WashingtonDCSSH1, filename_WashingtonDCSSH1), index_col=0
).set_index("years")
Downloading data from 'https://osf.io/4zynp/download' to file '/tmp/WashingtonDCSSH1930-2022.csv'.
SHA256 hash of downloaded file: d1725e1e5c6ccd7561644771c9ad17ab07e553c05b4156a38cb6780aaca36edb
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
data
ssh
years
1931 -140
1932 -112
1933 -6
1934 -51
1935 -94
... ...
2018 517
2019 299
2020 257
2021 272
2022 306

92 rows × 1 columns

data.ssh.plot(
    linestyle="-", marker=".", xlabel="Time (years)", ylabel="Annual Maximum Sea Surface Height (mm)", grid=True
)
<Axes: xlabel='Time (years)', ylabel='Annual Maximum Sea Surface Height (mm)'>
../../../_images/9d1a81a326f3d8758adaf7fad7729490bfca338bd622bdd8900779dc500c1b86.png
Click here for a description of the plot 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.

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.

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.

# 1931-1960
data_period1 = data.iloc[0:30]
# 1961-1990
data_period2 = data.iloc[30:60]
# 1990-2020
data_period3 = data.iloc[60:90]
# plot the histograms for each climate normal identified above
fig, ax = plt.subplots()
colors =["C0","C1","C2"]

# loop over all periods of 30 years
for ind, yr in enumerate(range(30,120,30)):
    sns.histplot(
        # select climate normal periods from dataset
        data.iloc[yr-30:yr].ssh,
        bins=np.arange(-200, 650, 50),
        color=colors[ind],
        element="step",
        alpha=0.5,
        # estimate PDF via kernel density estimation
        kde=True,
        label=f"{1901+yr}-{1930+yr}",
        ax=ax,
)

# aesthetics
ax.legend()
ax.set_xlabel("Annual Maximum Sea Surface Height (mm)")
/opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/seaborn/_oldcore.py:1119: FutureWarning: use_inf_as_na option is deprecated and will be removed in a future version. Convert inf values to NaN before operating instead.
  with pd.option_context('mode.use_inf_as_na', True):
/opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/seaborn/_oldcore.py:1119: FutureWarning: use_inf_as_na option is deprecated and will be removed in a future version. Convert inf values to NaN before operating instead.
  with pd.option_context('mode.use_inf_as_na', True):
/opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/seaborn/_oldcore.py:1119: FutureWarning: use_inf_as_na option is deprecated and will be removed in a future version. Convert inf values to NaN before operating instead.
  with pd.option_context('mode.use_inf_as_na', True):
Text(0.5, 0, 'Annual Maximum Sea Surface Height (mm)')
../../../_images/fcd62f7f04179e978235c8b551aaa7e6b3384b25d61d22acaf2a9a55c51753fc.png
Click here for a description of the plot 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.

Let’s also calculate the moments of each climate normal period:

# setup pandas dataframe
periods_stats = pd.DataFrame(index=["Mean", "Standard Deviation", "Skew"])

# add info for each climate normal period
periods_stats["1931-1960"] = [
    data_period1.ssh.mean(),
    data_period1.ssh.std(),
    data_period1.ssh.skew(),
]
periods_stats["1961-1990"] = [
    data_period2.ssh.mean(),
    data_period2.ssh.std(),
    data_period2.ssh.skew(),
]
periods_stats["1991-2020"] = [
    data_period3.ssh.mean(),
    data_period3.ssh.std(),
    data_period3.ssh.skew(),
]

periods_stats = periods_stats.T
periods_stats
Mean Standard Deviation Skew
1931-1960 -9.966667 87.095066 0.922327
1961-1990 85.200000 87.953906 0.854553
1991-2020 216.633333 105.739264 0.701258

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.

# 1931-1960
params_period1 = gev.fit(data_period1.ssh.values, 0)
shape_period1, loc_period1, scale_period1 = params_period1

# 1961-1990
params_period2 = gev.fit(data_period2.ssh.values, 0)
shape_period2, loc_period2, scale_period2 = params_period2

# 1991-2020
params_period3 = gev.fit(data_period3.ssh.values, 0)
shape_period3, loc_period3, scale_period3 = params_period3
# initialize list to store climate normals
data_periods = []

fig, ax = plt.subplots()
# initialize array of possible annual maximum sea surface heights in mm
x = np.linspace(-200, 600, 1000)

# loop over three 30-year periods
for ind, yr in enumerate(range(30,120,30)):
    # select climate normal periods from dataset
    data_periods.append(data.iloc[yr-30:yr])
    # fit the GEV distribution parameter
    shape_period, loc_period, scale_period = gev.fit(data_periods[ind].ssh.values, 0)
    # plot GEV distribution with these parameters
    ax.plot(
        x,
        gev.pdf(x, shape_period, loc=loc_period, scale=scale_period),
        color=colors[ind],
        linewidth=3,
        # add climate normal period label
        label=f"{1901+yr}-{1930+yr}"
)
# plot aesthetics
ax.legend()
ax.set_xlabel("Annual Maximum Sea Surface Height (mm)")
ax.set_ylabel("Density")
Text(0, 0.5, 'Density')
../../../_images/3f525347b63d5f68678f9c67600582644816de0f8218219b6a07c9f373f25e10.png

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.

Question 1#

  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?

Click for solution

Coding Exercise 1#

  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.

# setup dataframe with titles for each parameter
parameters = pd.DataFrame(index=["Location", "Scale", "Shape"])

# add in 1931-1960 parameters
parameters["1931-1960"] = ...

# add in 1961-1990 parameters
parameters["1961-1990"] = ...

# add in 1991-202 parameters
parameters["1991-2020"] = ...

# transpose the dataset so the time periods are rows
parameters = ...

# round the values for viewing
_ = ...

Click for solution

Summary#

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.

Resources#

Original data from this tutorial can be found here. Note the data used in this tutorial were preprocessed to extract the annual maximum values.