Regional and Local Sea Level Change

Regional and Local Sea Level Change#

../../../_images/84078410f877729fb9b6e6b2fbb4947a8714c9a6fc5c59b6b93d2c44cf1080e7.png

Since the start of continuous satellite measurements of global sea level in 1993, the absolute mean sea level in the vicinity of the Malakal tide gauge has risen by cm ( in). This compares to a relative value of cm ( in) measured by the tide gauge over this same time period. These values correspond to satellite-derived and station-derived rates of mm ( in) per year and mm ( in) per year, respectively. Differences in the rates of absolute (satellite-measured) and relative (tide gauge measured) sea level rise is largely attributable to vertical land motion. The sea level near Palau is rising faster than the global average, with satellite data showing a global mean sea level (GMSL) rise of 4.4 mm (0.17 in) per year (Willis et al., 2023).

In this notebook, we’ll be creating a table, a map, and a time series plot of absolute and regional sea level change at the Malakal, Palau tide gauge station from 1993-2022. Absolute sea level, typically measured by satellite altimetry, refers to the height of the sea surface relative to a reference ellipsoid. Here, we’ll use the global ocean gridded L4 Sea Level Anomalies (SLA) available from Copernicus, which is the sea surface height (SSH) minus the mean sea surface (MSS), where the MSS is the temporal mean of SSH over a given period. Relative sea level is measured by a tide gauge, and is the sea level relative to land at that location. Differences between the two measurements can arise from vertical land motion, regional oceanographic conditions like currents, and changes to the gravitational field (affecting the geoid).

Download Files: Map | Time Series Plot | Table

Setup#

First, we need to import all the necessary libraries.

# import necessary libraries
import numpy as np
import xarray as xr
import datetime as dt
from pathlib import Path
import pandas as pd
import os
import os.path as op
import sys

# data retrieval libraries
import requests
from urllib.request import urlretrieve #used for downloading files
import json
import copernicusmarine

# data processing libraries
from scipy import stats
import seaborn as sns
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

from myst_nb import glue #used for figure numbering when exporting to LaTeX

sys.path.append("../../../functions")
from data_downloaders import download_oni_index

Next, we’ll establish our directory for saving and reading data, and our output files.

data_base_dir = Path('../../../data')
path_figs = "../../../matrix_cc/figures"
data_dir = Path(data_base_dir,'sea_level')


#CHANGE THIS TO YOUR PATH!!
# output_dir = Path('/Users/laurac/Library/Mobile Documents/com~apple~CloudDocs/Projects/CC_indicators/CC_indicators/data/output')  # EDIT THIS TO YOUR PATH
output_dir = Path('/Users/laurac/Library/Mobile Documents/com~apple~CloudDocs/Projects/CC_indicators/CC_indicators/data/')  # EDIT THIS TO YOUR PATH

# Create the output directory if it does not exist
output_dir.mkdir(exist_ok=True)
data_dir.mkdir(exist_ok=True)


# Also let's just make our figure sizes the same:
plt.rcParams['figure.figsize'] = [10, 4]  # Set a default figure size for the notebook

Retrieve the Tide Station Data#

Assess Station Data Quality#

for the POR (1993-2022 …consult w/Ayesha and John)

Plot the hourly time series#

For a quick glance at the data

# Set the style of the plot
sns.set_style("whitegrid")

# Set the style of the plot
sns.set_style("whitegrid")

# Create a color palette
palette = sns.color_palette("Paired")

# Create the figure and the axes
fig, ax = plt.subplots()

# Plot the data

# plot tide gauge data
ax.scatter(rsl_daily['time'],100*rsl_daily, label='RSL', 
           color=palette[1], alpha=1, s= 5)

# Set the title, labels, and limits
ax.set(
    title=f'Tide Gauge ({station_name}) ({start_date_str} to {end_date_str})',
    xlabel='Year',
    ylabel='Water Level (cm)',
    ylim=[-50, 50],
    xlim=[start_date, end_date]
);

../../../_images/783db0f284ab0b59152332d207ec1d2cb7ba2a8665eea26c8c1fa4a12b7e7e5e.png

Calculate rate and magnitude of change#

Calculate values for both the Trend (rate of change) and Magnitude of Change

# calculate the rate of change for the tide gauge
trend_mag_rsl, trend_line_rsl, trend_rate_rsl, p_value_rsl, trend_err_rsl = process_trend_with_nan(rsl_daily)

print(f'The trend magnitude for the tide gauge is {100*trend_mag_rsl:.2f} cm.')
print(f'The trend rate for the tide gauge is {100*trend_rate_rsl:.2f} cm/year.')

The trend magnitude for the tide gauge is 12.26 cm.
The trend rate for the tide gauge is 0.41 cm/year.

Plot time series#

Now we’ll combine our daily relative sea level time series, the monthly mean, and a trend line into the same plot. We will label it with the Relative Sea Level Sea Level Change (magnitude in cm) and rate of change (cm/yr) at the Tide Station. The “zero” line is referenced to the [1993,2012] period, following the altimetry.

# make an rsl monthly mean for plotting
rsl_monthly = rsl_daily.resample(time='1ME').mean().squeeze()

# Set the style of the plot
sns.set_style("whitegrid")

# Create a color palette
palette = sns.color_palette("Set1")

# Plot the timeseries
# Create the figure and the axes
fig, ax = plt.subplots()
# plot altimetry data
ax.scatter(rsl_daily['time'], 100*rsl_daily, 
           label='Tide Gauge', color=palette[1], alpha=0.2, s= 5)
ax.plot(rsl_daily['time'], 100*trend_line_rsl, 
        label='Tide Gauge Trend', color=palette[1], linestyle='--')
# plot the monthly mean sea level
ax.plot(rsl_monthly['time'], 100*rsl_monthly, 
        label='Tide Gauge', color=palette[3])


# Set the title and labels
ax.set_title(f'Tide Gauge ({station_name}) ({start_date_str} to {end_date_str})')
ax.set_xlabel('Year')
ax.set_ylabel('Height (cm)')

# Set the y limits
ax.set_ylim([-50, 50])

# Set the x limits
ax.set_xlim([start_date, end_date])

# Add a legend
trendmag_str = f'Δ Sea Level: {100*trend_mag_rsl:.2f} cm, Trend: {100*trend_rate_rsl:.2f} cm/year'
# Add text in a white box to bottom right of plot
ax.text(0.95, 0.05, trendmag_str, transform=ax.transAxes, 
        fontsize=14, verticalalignment='bottom', horizontalalignment='right',
        bbox=dict(facecolor='white', alpha=0.5))

Text(0.95, 0.05, 'Δ Sea Level: 12.26 cm, Trend: 0.41 cm/year')

../../../_images/5df6f988d5065c995156d98d17ceb9308e7f6369e14577719a75d28c5b8389e1.png

Combining both sources#

Create a Table#

That compares the results of the absolute and relative time series.

# Constants
DATA_SOURCE_ALTIMETRY = 'CMEMS SSH L4 0.125 deg (SLA)'
DATA_SOURCE_TIDE_GAUGE = 'UHSLC RQDS'
TIME_PERIOD = f'{start_date_str} to {end_date_str}'

# Calculated values 
trend_mmyr_altimetry = 1000 * trend_rate_asl.values
trend_mmyr_tide_gauge = 1000 * trend_rate_rsl.values
delta_sea_level_altimetry = 100 * trend_mag_asl.values
delta_sea_level_tide_gauge = 100 * trend_mag_rsl.values


# Create DataFrame
SL_magnitude_results = pd.DataFrame({
    'Trend (mm/yr)': [trend_mmyr_altimetry, trend_mmyr_tide_gauge],
    'Trend (in/yr)': [trend_mmyr_altimetry * 0.0393701, trend_mmyr_tide_gauge * 0.0393701],
    'Δ Sea Level (cm)': [delta_sea_level_altimetry, delta_sea_level_tide_gauge],
    'Δ Sea Level (in)': [delta_sea_level_altimetry * 0.393701, delta_sea_level_tide_gauge * 0.393701],
    'Latitude': [sla_nearest_lat, rsl['lat'].values[0]],
    'Longitude': [sla_nearest_lon, rsl['lon'].values[0]],
    'Time_Period': [TIME_PERIOD, TIME_PERIOD],
    'Data_Source': [DATA_SOURCE_ALTIMETRY, DATA_SOURCE_TIDE_GAUGE],
}, index=['Altimetry', 'Tide Gauge'])

# Save to CSV
output_file_path = output_dir / 'SL_magnitude_results.csv'

# Use the path for operations, e.g., saving a DataFrame
SL_magnitude_results.to_csv(output_file_path)

# glue trend_rates to the notebook
# glue("trend_rate_rsl", trend_rate_rsl, display=False)
# glue("trend_rate_asl", trend_rate_asl, display=False)




SL_magnitude_results

	Trend (mm/yr)	Trend (in/yr)	Δ Sea Level (cm)	Δ Sea Level (in)	Latitude	Longitude	Time_Period	Data_Source
Altimetry	5.779614	0.227544	17.336469	6.825385	7.3125	134.4375	1993-01-01 to 2022-12-31	CMEMS SSH L4 0.125 deg (SLA)
Tide Gauge	4.086637	0.160891	12.258233	4.826078	7.33	134.462997	1993-01-01 to 2022-12-31	UHSLC RQDS

The cell below will save all these variables for printing:

Create a Map#

Now we’ll combine a both the tide gauge and altimetry sources into a map that includes the absolute change with the addition of an icon depicting the magnitude of relative change at the tide station.

fig, ax, crs,cmap = plot_map(vmin,vmax,xlims,ylims)

# plot the trend*100 for centimeters
trend_mag_cmems_cm = trend_mag_cmems * 100

# plot a map of the magnitude of SL change in centimeters
trend_mag_cmems_cm.plot(ax=ax, transform=crs,
                         vmin=vmin, vmax=vmax, cmap=cmap, add_colorbar=True, 
                         cbar_kwargs={'label': 'Δ Sea Level (cm, 1993-2022)'},
)

# add the EEZ
ax.plot(palau_eez[:, 0], palau_eez[:, 1], transform=crs, color='black', linewidth=2)
# add the tide gauge location with black outlined dot, colored by the sea level value
ax.scatter(rsl['lon'], rsl['lat'], transform=crs, s=200, 
           c=100*trend_mag_rsl, vmin=vmin, vmax=vmax, cmap=cmap,
           linewidth=0.5, edgecolor='black')

plot_zebra_frame(ax, lw=5, segment_length=2, crs=crs)

glue("mag_fig", fig, display=False)

# save the figure
output_file_path = output_dir / 'SL_magnitude_map.png'
fig.savefig(output_file_path, dpi=300, bbox_inches='tight')

../../../_images/a518c3f34b2e6e73602389cf455e08916263ef0f8bdd0c4e1552c40839ac69e8.png — Fig. 3 Map of absolute and relative sea level change from the altimetry and tide gauge record near Malakal, Palau station from 1993-01-01 to 2022-12-31.#

Create a Time series plot#

Finally we will combine both tide gauge and altimetry sources into a time series plot that includes both Absolute and Relative Time Series.

# Set the style of the plot
sns.set_style("whitegrid")

# Create a color palette
palette = sns.color_palette("Set1")

# Create the figure and the axes
fig, ax = plt.subplots()

# plot altimetry data
labelSat = f'Altimetry Trend ({1000*trend_rate_asl:.2f} mm/year)'
ax.scatter(sla_nearest['time'], 100*sla_nearest, 
           label='Altimetry', color=palette[0], alpha=0.2, s=5)
ax.plot(sla_nearest['time'], 100*trend_line_asl, 
        label=labelSat, color=palette[0], linestyle='-')

# plot tide gauge data
labelTG = f'Tide Gauge Trend ({1000*trend_rate_rsl:.2f} mm/year)'
ax.scatter(rsl_daily['time'], 100*rsl_daily, 
           label='Tide Gauge', color=palette[1], alpha=0.2, s=10)
ax.plot(rsl_daily['time'], 100*trend_line_rsl, 
        label=labelTG, color=palette[1], linestyle='-')

# Set the title and labels
title = (
    f'Altimetry ({lat_str}, {lon_str}) vs \n'
    f'Tide Gauge ({station_name}) ({start_date_str} to {end_date_str})'
)
ax.set_title(title)
ax.set_xlabel('Year')
ax.set_ylabel('Height (cm, MHHW)')

# Set the y and x limits
ax.set_ylim([-40, 40])
ax.set_xlim([start_date, end_date])

# Add a legend below the plot
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=2)

glue("trend_fig", fig, display=False)

# save the figure
output_file_path = output_dir / 'SL_magnitude_timeseries.png'
fig.savefig(output_file_path, dpi=300, bbox_inches='tight')

plt.savefig(op.join(path_figs, 'F10_SeaLevel_trends.png'), dpi=300, bbox_inches='tight')

ENSO Effects on Sea Level#

Next we’ll look at the variation of sea level with El Niño and La Niña.

# add nino
p_data = 'https://psl.noaa.gov/data/correlation/oni.data'
oni = download_oni_index(p_data)

/Users/laurac/Library/Mobile Documents/com~apple~CloudDocs/Projects/CC_indicators/CC_indicators/notebooks_historical/ocean/1_sea_level/../../../functions/data_downloaders.py:242: FutureWarning: The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead
  oni = pd.read_csv(

# Remove the trend from the tide gauge timeseries
detrended_rsl = rsl_daily - trend_line_rsl

# estimate the seasonal cycle (detrended_rsl is xarray)
seasonal_cycle = detrended_rsl.groupby('time.dayofyear').mean('time')

# Remove the seasonal cycle from the detrended tide gauge timeseries
deseasoned_rsl = detrended_rsl.groupby('time.dayofyear') - seasonal_cycle

# Resample the ONI index to daily frequency and align with the tide gauge timeseries
oni_daily = oni.resample('D').interpolate(method='linear')
oni_daily = oni_daily.reindex(deseasoned_rsl['time'], method='nearest')

# scatter plot of oni daily and detrended sea level
# Scatter plot of ONI daily and detrended sea level
fig, ax = plt.subplots(figsize=(8,6))
ax.scatter(oni_daily, deseasoned_rsl, color='gray', alpha=0.2, s=7)
ax.set_xlabel('ONI (°C)')
ax.set_ylabel('De-seasoned Sea Level (m)')

# Emphasize the 1997-1998 and 2015-2016 El Niño events using pandas indexing
nina_1998 = deseasoned_rsl.sel(time=slice('1998-08-01', '2001-01-30'))
nino_1997 = deseasoned_rsl.sel(time=slice('1997-06-01', '1998-04-30'))

# Plot as dashed lines
ax.scatter(oni_daily.loc['1998-08-01':'2001-01-30'], nina_1998, color='blue', label='1998-2001 La Niña',alpha=0.5)
ax.scatter(oni_daily.loc['1997-06-01':'1998-04-30'], nino_1997, color='red', label='1997-1998 El Niño',alpha=0.5)

# draw a correlation line
# Keep only valid (non-NaN) indices
deseasoned_rsl_series = deseasoned_rsl.to_series()
valid_mask = ~deseasoned_rsl_series.isna()
oni_valid = oni_daily.squeeze().loc[valid_mask]
deseasoned_valid = deseasoned_rsl_series.loc[valid_mask]

slope, intercept, r_value, p_value, std_err = stats.linregress(oni_valid, deseasoned_valid, )
x = np.linspace(-2,3,100)
y = slope * x + intercept
ax.plot(x, y, color='black', linestyle='-', label=f'{slope:.2f} m/ °C (r$^2$ = {r_value**2:.2f})')

ax.set_xlim(-2,3)
ax.legend(frameon=False)
# Add shading for different ONI ranges


ax.axvspan(-2, -1.5, color='blue', alpha=0.15, linewidth=0, zorder=0)
ax.axvspan(-1.5, -1, color='blue', alpha=0.10, linewidth=0, zorder=0)
ax.axvspan(-1, -0.5, color='blue', alpha=0.05, linewidth=0, zorder=0)
ax.axvspan(-0.5, 0.5, color='gray', alpha=0.05, linewidth=0, zorder=0)
ax.axvspan(0.5, 1, color='red', alpha=0.05, linewidth=0, zorder=0)
ax.axvspan(1, 1.5, color='red', alpha=0.10, linewidth=0, zorder=0)
ax.axvspan(1.5, 2, color='red', alpha=0.15, linewidth=0, zorder=0)
ax.axvspan(2, 3, color='red', alpha=0.20, linewidth=0, zorder=0)


# spines should be black not gray

for spine in ax.spines.values():
    spine.set_edgecolor('black')
    spine.set_linewidth(1)

# set all text to avenir
for text in ax.texts:
    text.set_fontname('Avenir')

# set all labels to avenir
for label in ax.get_xticklabels() + ax.get_yticklabels() + ax.get_legend().get_texts() :
    label.set_fontname('Avenir')

ax.xaxis.label.set_fontname('Avenir')
ax.yaxis.label.set_fontname('Avenir')

plt.show()

# save the figure
output_file_path = output_dir / 'SL_ONI_scatter.png'
fig.savefig(output_file_path, dpi=300, bbox_inches='tight')

../../../_images/26fcb0fed8c8010f5daca4f49c1768458c502d030991bb513c056ef844e47dc7.png

# Create a DataFrame to summarize the results
summary_data = {
    'Description': ['Slope of Correlation Line', 
                    'El Niño 1997-1998 Min Sea Level', 
                    'El Niño 1997-1998 Median Sea Level',
                    'La Niña 1998-2000 Max Sea Level', 
                    'La Niña 1998-2000 Median Sea Level',
                    'Correlation Coefficient (r)', 
                    'p-value'],
    'Value': [slope, 
              100*nino_1997.min().values,
              100*nino_1997.median().values, 
              100*nina_1998.max().values, 
              100*nina_1998.median().values,
              r_value, 
              p_value],
    'units': ['m/°C', 'cm', 'cm','cm', 'cm', '', '']
}

summary_df = pd.DataFrame(summary_data)

summary_df.attrs['description'] = "Correlation of Sea Level at Malakal with ONI 1983-2024."
summary_df.attrs['source'] = "NOAA / University of Hawaii Sea Level Center"
summary_df.attrs['notes'] = "Sea level data has been deseasoned and detrended."

# Display the summary table with attributes

# Function to format values with units
def format_with_units(value, unit):
    if pd.isna(value):  # Handle NaNs
        return "N/A"
    if unit:  # Add unit only if it's not an empty string
        return f"{value:.2f} {unit}"
    return f"{value:.2f}"  # No unit, just the value

# Apply formatting
summary_df['Formatted'] = summary_df.apply(lambda row: format_with_units(row['Value'], row['units']), axis=1)

# Keep only needed columns (remove 'Value' and 'units' separately)
formatted_df = summary_df[['Description', 'Formatted']].rename(columns={'Formatted': 'Value'})

# Remove index for clean display
formatted_df = formatted_df.reset_index(drop=True)



footer_text = f"Source: {summary_df.attrs['source']}"
notes_text = f"Notes: {summary_df.attrs['notes']}"
footer_df = pd.DataFrame({'Description': [footer_text, notes_text], 
                          'Value': ['','']})

# Combine summary table with footer
summary_table = pd.concat([formatted_df, footer_df], ignore_index=True)  # Removes index completely

# Apply styling
styled_df = summary_table.style.set_caption(
    summary_df.attrs['description']
).set_table_styles([
    {
        'selector': 'caption',
        'props': [('color', 'black'), ('font-size', '16px'), ('font-weight', 'bold')]
    },
    {
        'selector': 'tbody tr:last-child td, tbody tr:nth-last-child(2) td',  # Targets last two rows (footer)
        'props': [('font-style', 'italic'), ('font-align', 'center')]
    }
])



# Save the summary table to a CSV file
summary_file_path = output_dir / 'ENSO_SL_influence_summary.csv'
summary_table.to_csv(summary_file_path, index=False)


styled_df

Correlation of Sea Level at Malakal with ONI 1983-2024.
	Description	Value
0	Slope of Correlation Line	-0.10 m/°C
1	El Niño 1997-1998 Min Sea Level	-28.57 cm
2	El Niño 1997-1998 Median Sea Level	-17.74 cm
3	La Niña 1998-2000 Max Sea Level	35.22 cm
4	La Niña 1998-2000 Median Sea Level	13.34 cm
5	Correlation Coefficient (r)	-0.71
6	p-value	0.00
7	Source: NOAA / University of Hawaii Sea Level Center
8	Notes: Sea level data has been deseasoned and detrended.

Regional and Local Sea Level Change

Contents

Regional and Local Sea Level Change#

Setup#

Retrieve Data Sources#

Retrieve the location of the Malakal, Palau tide gauge#

Retrieve the EEZ boundary for Palau#

Retrieve altimetry data#

Retrieve Tide Gauge Data#

Process the data#

Clip#

Plot the timeseries#

Calculate change#

Plot a map#

Plot time series#

Retrieve the Tide Station Data#

Assess Station Data Quality#

Plot the hourly time series#

Calculate rate and magnitude of change#

Plot time series#

Combining both sources#

Create a Table#

Create a Map#

Create a Time series plot#

ENSO Effects on Sea Level#