# --- Helper Functions ---
def load_and_clean_sea_level_data(data_url, is_monthly=True):
"""
Loads sea level data from a PSMSL URL and performs initial cleaning using Polars.
Args:
data_url (str): The URL to the PSMSL data file.
is_monthly (bool): True if the data is monthly, False if annual.
Affects how 'missing_days' is interpreted and dtypes.
Returns:
polars.DataFrame: A DataFrame with cleaned sea level data and anomalies.
"""
try:
response = requests.get(data_url)
response.raise_for_status() # Raise an exception for HTTP errors
except requests.exceptions.RequestException as e:
print(f"Error fetching data from {data_url}: {e}")
return pl.DataFrame() # Return empty Polars DataFrame on error
data = io.StringIO(response.text)
col_names = ["time", "sl", "missing_days", "flag_for_attention"]
# Define dtypes explicitly for Polars read_csv for robustness
if is_monthly:
col_dtypes = {
"time": pl.Float64,
"sl": pl.Float64,
"missing_days": pl.Int64, # Numeric for monthly
"flag_for_attention": pl.Int64
}
else:
col_dtypes = {
"time": pl.Float64,
"sl": pl.Float64,
"missing_days": pl.Utf8, # String for annual
"flag_for_attention": pl.Int64
}
sea_level_data = pl.read_csv(
data,
separator=";",
has_header=False,
null_values="-99999", # Specifies how missing values in 'sl' are read
new_columns=col_names,
schema_overrides=col_dtypes
)
# Handle missing values based on 'missing_days' and 'flag_for_attention'
# Use Polars expressions with .with_columns() for conditional logic
if is_monthly:
sea_level_data = sea_level_data.with_columns(
pl.when((pl.col('missing_days') != 0) | (pl.col('flag_for_attention') != 0))
.then(None) # Set 'sl' to NaN if conditions met
.otherwise(pl.col('sl'))
.alias('sl')
)
else:
sea_level_data = sea_level_data.with_columns(
pl.when((pl.col('missing_days') != 'N') | (pl.col('flag_for_attention') != 0))
.then(None) # Set 'sl' to NaN if conditions met
.otherwise(pl.col('sl'))
.alias('sl')
)
# Calculate sea level anomaly
# pl.col('sl').mean() automatically skips null values
sea_level_data = sea_level_data.with_columns(
(pl.col('sl') - pl.col('sl').mean()).alias('sl_anomaly')
)
return sea_level_data
def sl_rmse(params, time, sea_levels):
"""
Calculates the Root Mean Squared Error (RMSE) for a second-order polynomial fit.
Args:
params (list/array): Coefficients [a, b, c, t_0] for the polynomial
a * (time - t_0)^2 + b * (time - t_0) + c.
time (array): Time values.
sea_levels (array): Observed sea level anomaly values.
Returns:
float: The calculated RMSE.
"""
a, b, c, t_0 = params
est_sl = a * (time - t_0)**2 + b * (time - t_0) + c
resids = sea_levels - est_sl
rmse = np.sqrt(np.mean(resids**2))
return rmse
def analyze_city(dataset_number, city_name, plot_acf_residuals=True):
"""
Performs sea level trend analysis for a given city using annual data.
Includes data loading, polynomial regression, and bootstrapping.
Args:
dataset_number (int): The PSMSL dataset ID for the city.
city_name (str): The name of the city for plot titles.
plot_acf_residuals (bool): Whether to plot autocorrelation for residuals.
Returns:
plotly.graph_objects.Figure: The Plotly figure with bootstrap realizations.
Returns None if no valid data is found.
"""
print(f"\n--- Analyzing {city_name} (Dataset ID: {dataset_number}) ---")
data_url = f"https://psmsl.org/data/obtaining/rlr.annual.data/{dataset_number}.rlrdata"
sea_level_data = load_and_clean_sea_level_data(data_url, is_monthly=False)
# Drop rows where 'sl_anomaly' is null
sea_level_data_filtered = sea_level_data.drop_nulls(subset=['sl_anomaly'])
if sea_level_data_filtered.is_empty():
print(f"No valid data points found for {city_name}. Skipping analysis.")
return None
# Extract data as NumPy arrays for optimization and plotting
t = sea_level_data_filtered.get_column('time').to_numpy()
obs_sl = sea_level_data_filtered.get_column('sl_anomaly').to_numpy()
start_params = [0, 0, 0, 1962] # Initial guess for a, b, c, t_0
# Optimize the sl_rmse function (Nelder-Mead is R's optim default for unconstrained problems)
optim_fit = minimize(sl_rmse, start_params, args=(t, obs_sl), method='Nelder-Mead')
best_a, best_b, best_c, best_t_0 = optim_fit.x
best_est_sl = best_a * (t - best_t_0)**2 + best_b * (t - best_t_0) + best_c
resids = obs_sl - best_est_sl
if plot_acf_residuals:
plt.figure(figsize=(10, 5))
plot_acf(resids, lags=10, title=f"Autocorrelation of Original Residuals - {city_name}")
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.grid(True)
plt.show()
# Bootstrapping
n_boot = 1000
proj_sl = np.full((n_boot, len(t)), np.nan)
for i in range(n_boot):
# Bootstrap residuals by resampling with replacement
sl_boot = best_est_sl + np.random.choice(resids, size=len(resids), replace=True)
# Optimize for the bootstrapped data
optim_boot = minimize(sl_rmse, start_params, args=(t, sl_boot), method='Nelder-Mead')
boot_a_i, boot_b_i, boot_c_i, boot_t_0_i = optim_boot.x
proj_sl[i, :] = boot_a_i * (t - boot_t_0_i)**2 + boot_b_i * (t - boot_t_0_i) + boot_c_i
# Calculate 95% confidence interval from bootstrap realizations
lower_bound = np.percentile(proj_sl, 2.5, axis=0)
upper_bound = np.percentile(proj_sl, 97.5, axis=0)
# Create the Plotly plot for bootstrap realizations
fig = go.Figure()
# Add 95% confidence interval
fig.add_trace(go.Scatter(x=np.concatenate([t, t[::-1]]),
y=np.concatenate([upper_bound, lower_bound[::-1]]),
fill='toself',
fillcolor='rgba(200,200,200,0.5)',
line=dict(color='rgba(255,255,255,0)'),
name='95% Bootstrap Trend CI',
showlegend=True))
# Add original data
fig.add_trace(go.Scatter(x=t, y=obs_sl, mode='markers', name='Original Data',
marker=dict(color='black'),
hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
# Add best-fit polynomial
fig.add_trace(go.Scatter(x=t, y=best_est_sl, mode='lines', name='Best-Fit Polynomial',
line=dict(color='red', width=2),
hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig.update_layout(title=f'Sea level over time -- {city_name}',
xaxis_title='Time (yr)',
yaxis_title='Sea level anomaly (mm)',
plot_bgcolor='white',
xaxis=dict(
showgrid=True,
gridcolor='lightgray',
zeroline=True,
zerolinecolor = 'lightgray'
),
yaxis=dict(
showgrid=True,
gridcolor='lightgray',
zeroline=True,
zerolinecolor = 'lightgray'
))
return fig
# --- Monthly Puerto Rico Analysis ---
print("--- Monthly Puerto Rico Analysis ---")
data_url_monthly = "https://psmsl.org/data/obtaining/rlr.monthly.data/1001.rlrdata"
sea_level_data_monthly = load_and_clean_sea_level_data(data_url_monthly, is_monthly=True)
print("Head of the monthly data:")
print(sea_level_data_monthly.head()) # Removed .to_string()
print("\nSummary of the monthly data:")
print(sea_level_data_monthly.describe()) # Removed .to_string()
# Interactive Histogram of Sea Level Anomalies
sea_level_data_monthly_plot = sea_level_data_monthly.drop_nulls(subset=['sl_anomaly'])
if not sea_level_data_monthly_plot.is_empty():
# Convert to pandas for plotly express, as it often expects pandas DataFrames
fig_hist_monthly = px.histogram(sea_level_data_monthly_plot, x="sl_anomaly",
nbins=int(sea_level_data_monthly_plot.get_column('sl_anomaly').max() - sea_level_data_monthly_plot.get_column('sl_anomaly').min()) // 10,
title="Interactive Histogram of Sea Level Anomalies (Monthly)",
labels={"sl_anomaly": "Sea Level Anomaly (mm)"})
fig_hist_monthly.update_layout(yaxis_title="Frequency")
fig_hist_monthly.show()
# Interactive Sea Level Anomaly Over Time
if not sea_level_data_monthly.drop_nulls(subset=['sl_anomaly']).is_empty():
fig_timeseries_monthly = px.line(sea_level_data_monthly, x="time", y="sl_anomaly",
title="Interactive Sea Level Anomaly Over Time (Monthly)",
labels={"time": "Year", "sl_anomaly": "Sea Level Anomaly (mm)"},
markers=True)
fig_timeseries_monthly.update_traces(hovertemplate='Year: %{x:.4f}<br>Anomaly: %{y:.2f} mm<extra></extra>')
fig_timeseries_monthly.show()
# Running Polynomial Regression (Monthly)
sea_level_data_monthly_filtered = sea_level_data_monthly.drop_nulls(subset=['sl_anomaly'])
if not sea_level_data_monthly_filtered.is_empty():
t_monthly = sea_level_data_monthly_filtered.get_column('time').to_numpy()
obs_sl_monthly = sea_level_data_monthly_filtered.get_column('sl_anomaly').to_numpy()
start_params = [0, 0, 0, 1962] # Initial guess for a, b, c, t_0
optim_fit_monthly = minimize(sl_rmse, start_params, args=(t_monthly, obs_sl_monthly), method='Nelder-Mead')
best_a_monthly, best_b_monthly, best_c_monthly, best_t_0_monthly = optim_fit_monthly.x
best_est_sl_monthly = best_a_monthly * (t_monthly - best_t_0_monthly)**2 + \
best_b_monthly * (t_monthly - best_t_0_monthly) + best_c_monthly
resids_monthly = obs_sl_monthly - best_est_sl_monthly
fig_poly_monthly = go.Figure()
fig_poly_monthly.add_trace(go.Scatter(x=t_monthly, y=obs_sl_monthly, mode='lines', name='Data',
line=dict(color='red'),
hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig_poly_monthly.add_trace(go.Scatter(x=t_monthly, y=best_est_sl_monthly, mode='lines', name='2nd-Order Polynomial',
line=dict(color='blue'),
hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig_poly_monthly.update_layout(title='Sea Level Anomaly and 2nd-Order Polynomial Fit (Monthly)',
xaxis_title='Time (yr)',
yaxis_title='Sea level anomaly (mm)')
fig_poly_monthly.show()
# Residual plot (Monthly)
plt.figure(figsize=(10, 6))
plt.plot(t_monthly, resids_monthly, label='Residuals', linewidth=1)
plt.xlabel('Time (yr)')
plt.ylabel('Residuals (mm)')
plt.title('Residuals of Monthly Sea Level Anomaly Fit')
plt.grid(True)
plt.show()
# Autocorrelation of original residuals (Monthly)
plt.figure(figsize=(10, 6))
plot_acf(resids_monthly, lags=10, title="Autocorrelation of Original Residuals (Monthly)")
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.grid(True)
plt.show()
# Bootstrapping (Monthly)
n_boot = 1000
proj_sl_monthly = np.full((n_boot, len(t_monthly)), np.nan)
for i in range(n_boot):
sl_boot = best_est_sl_monthly + np.random.choice(resids_monthly, size=len(resids_monthly), replace=True)
optim_boot = minimize(sl_rmse, start_params, args=(t_monthly, sl_boot), method='Nelder-Mead')
boot_a_i, boot_b_i, boot_c_i, boot_t_0_i = optim_boot.x
proj_sl_monthly[i, :] = boot_a_i * (t_monthly - boot_t_0_i)**2 + \
boot_b_i * (t_monthly - boot_t_0_i) + boot_c_i
# Calculate 95% confidence interval from bootstrap realizations for monthly data
lower_bound_monthly = np.percentile(proj_sl_monthly, 2.5, axis=0)
upper_bound_monthly = np.percentile(proj_sl_monthly, 97.5, axis=0)
fig_bootstrap_monthly = go.Figure()
# Add 95% confidence interval for monthly data
fig_bootstrap_monthly.add_trace(go.Scatter(x=np.concatenate([t_monthly, t_monthly[::-1]]),
y=np.concatenate([upper_bound_monthly, lower_bound_monthly[::-1]]),
fill='toself',
fillcolor='rgba(200,200,200,0.3)',
line=dict(color='rgba(255,255,255,0)'),
name='95% CI',
showlegend=True))
fig_bootstrap_monthly.add_trace(go.Scatter(x=t_monthly, y=obs_sl_monthly, mode='lines', name='Original Data',
line=dict(color='black', width=2),
hovertemplate='Year: %{x:.4f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig_bootstrap_monthly.add_trace(go.Scatter(x=t_monthly, y=best_est_sl_monthly, mode='lines', name='Best-Fit Polynomial',
line=dict(color='red', width=2),
hovertemplate='Year: %{x:.4f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig_bootstrap_monthly.update_layout(title='95% Confidence Interval of Second-Order Polynomials (Monthly)',
xaxis_title='Time (yr)',
yaxis_title='Sea level anomaly (mm)')
fig_bootstrap_monthly.show()
# Autocorrelation of bootstrap residuals (Monthly)
boot_resids_sample_monthly = np.random.choice(resids_monthly, size=len(resids_monthly), replace=True)
plt.figure(figsize=(10, 6))
plot_acf(boot_resids_sample_monthly, lags=10, title="Autocorrelation of Bootstrap Residuals Replicate (Monthly)")
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.grid(True)
plt.show()
else:
print("Monthly Puerto Rico analysis skipped due to no valid data.")
# --- Yearly Puerto Rico Analysis ---
print("\n--- Yearly Puerto Rico Analysis ---")
data_url_yearly = "https://psmsl.org/data/obtaining/rlr.annual.data/1001.rlrdata"
sea_level_data_yearly = load_and_clean_sea_level_data(data_url_yearly, is_monthly=False)
print("Head of the yearly data:")
print(sea_level_data_yearly.head()) # Removed .to_string()
print("\nSummary of the yearly data:")
print(sea_level_data_yearly.describe()) # Removed .to_string()
# Interactive Histogram of Sea Level Anomalies
sea_level_data_yearly_plot = sea_level_data_yearly.drop_nulls(subset=['sl_anomaly'])
if not sea_level_data_yearly_plot.is_empty():
fig_hist_yearly = px.histogram(sea_level_data_yearly_plot, x="sl_anomaly", # Convert to pandas for plotly express
nbins=int(sea_level_data_yearly_plot.get_column('sl_anomaly').max() - sea_level_data_yearly_plot.get_column('sl_anomaly').min()) // 10,
title="Interactive Histogram of Sea Level Anomalies (Yearly)",
labels={"sl_anomaly": "Sea Level Anomaly (mm)"})
fig_hist_yearly.update_layout(yaxis_title="Frequency")
fig_hist_yearly.show()
# Interactive Sea Level Anomaly Over Time
if not sea_level_data_yearly.drop_nulls(subset=['sl_anomaly']).is_empty():
fig_timeseries_yearly = px.line(sea_level_data_yearly, x="time", y="sl_anomaly", # Convert to pandas for plotly express
title="Interactive Sea Level Anomaly Over Time (Yearly)",
labels={"time": "Year", "sl_anomaly": "Sea Level Anomaly (mm)"},
markers=True)
fig_timeseries_yearly.update_traces(hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.0f} mm<extra></extra>')
fig_timeseries_yearly.show()
# Running Polynomial Regression (Yearly)
sea_level_data_yearly_filtered = sea_level_data_yearly.drop_nulls(subset=['sl_anomaly'])
if not sea_level_data_yearly_filtered.is_empty():
t_yearly = sea_level_data_yearly_filtered.get_column('time').to_numpy()
obs_sl_yearly = sea_level_data_yearly_filtered.get_column('sl_anomaly').to_numpy()
start_params = [0, 0, 0, 1962] # Initial guess for a, b, c, t_0
optim_fit_yearly = minimize(sl_rmse, start_params, args=(t_yearly, obs_sl_yearly), method='Nelder-Mead')
best_a_yearly, best_b_yearly, best_c_yearly, best_t_0_yearly = optim_fit_yearly.x
best_est_sl_yearly = best_a_yearly * (t_yearly - best_t_0_yearly)**2 + \
best_b_yearly * (t_yearly - best_t_0_yearly) + best_c_yearly
resids_yearly = obs_sl_yearly - best_est_sl_yearly
fig_poly_yearly = go.Figure()
fig_poly_yearly.add_trace(go.Scatter(x=t_yearly, y=obs_sl_yearly, mode='lines', name='Data',
line=dict(color='red'),
hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig_poly_yearly.add_trace(go.Scatter(x=t_yearly, y=best_est_sl_yearly, mode='lines', name='2nd-Order Polynomial',
line=dict(color='blue'),
hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig_poly_yearly.update_layout(title='Sea Level Anomaly and 2nd-Order Polynomial Fit (Yearly)',
xaxis_title='Time (yr)',
yaxis_title='Sea level anomaly (mm)')
fig_poly_yearly.show()
# Residual plot (Yearly)
plt.figure(figsize=(10, 6))
plt.plot(t_yearly, resids_yearly, label='Residuals', linewidth=1)
plt.xlabel('Time (yr)')
plt.ylabel('Residuals (mm)')
plt.title('Residuals of Yearly Sea Level Anomaly Fit')
plt.grid(True)
plt.show()
# Autocorrelation of original residuals (Yearly)
plt.figure(figsize=(10, 6))
plot_acf(resids_yearly, lags=10, title="Autocorrelation of Original Residuals (Yearly)")
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.grid(True)
plt.show()
# Bootstrapping (Yearly)
n_boot = 1000
proj_sl_yearly = np.full((n_boot, len(t_yearly)), np.nan)
for i in range(n_boot):
sl_boot = best_est_sl_yearly + np.random.choice(resids_yearly, size=len(resids_yearly), replace=True)
optim_boot = minimize(sl_rmse, start_params, args=(t_yearly, sl_boot), method='Nelder-Mead')
boot_a_i, boot_b_i, boot_c_i, boot_t_0_i = optim_boot.x
proj_sl_yearly[i, :] = boot_a_i * (t_yearly - boot_t_0_i)**2 + \
boot_b_i * (t_yearly - boot_t_0_i) + boot_c_i
# Calculate 95% confidence interval from bootstrap realizations for yearly data
lower_bound_yearly = np.percentile(proj_sl_yearly, 2.5, axis=0)
upper_bound_yearly = np.percentile(proj_sl_yearly, 97.5, axis=0)
fig_bootstrap_yearly = go.Figure()
# Add 95% confidence interval for yearly data
fig_bootstrap_yearly.add_trace(go.Scatter(x=np.concatenate([t_yearly, t_yearly[::-1]]),
y=np.concatenate([upper_bound_yearly, lower_bound_yearly[::-1]]),
fill='toself',
fillcolor='rgba(200,200,200,0.3)',
line=dict(color='rgba(255,255,255,0)'),
name='95% Bootstrap Trend CI',
showlegend=True))
fig_bootstrap_yearly.add_trace(go.Scatter(x=t_yearly, y=obs_sl_yearly, mode='lines', name='Original Data',
line=dict(color='black', width=2),
hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig_bootstrap_yearly.add_trace(go.Scatter(x=t_yearly, y=best_est_sl_yearly, mode='lines', name='Best-Fit Polynomial',
line=dict(color='red', width=2),
hovertemplate='Year: %{x:.0f}<br>Anomaly: %{y:.2f} mm<extra></extra>'))
fig_bootstrap_yearly.update_layout(title='95% Confidence Interval of Second-Order Polynomials (Yearly)',
xaxis_title='Time (yr)',
yaxis_title='Sea level anomaly (mm)')
fig_bootstrap_yearly.show()
# Autocorrelation of bootstrap residuals (Yearly)
boot_resids_sample_yearly = np.random.choice(resids_yearly, size=len(resids_yearly), replace=True)
plt.figure(figsize=(10, 6))
plot_acf(boot_resids_sample_yearly, lags=10, title="Autocorrelation of Bootstrap Residuals Replicate (Yearly)")
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.grid(True)
plt.show()
else:
print("Yearly Puerto Rico analysis skipped due to no valid data.")