In [ ]:
import polars as pl  # Import Polars
import numpy as np
import requests
import io
import plotly.graph_objects as go
import plotly.express as px
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
In [29]:
# --- 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.")
--- Monthly Puerto Rico Analysis ---
Head of the monthly data:
shape: (5, 5)
┌───────────┬────────┬──────────────┬────────────────────┬─────────────┐
│ time      ┆ sl     ┆ missing_days ┆ flag_for_attention ┆ sl_anomaly  │
│ ---       ┆ ---    ┆ ---          ┆ ---                ┆ ---         │
│ f64       ┆ f64    ┆ i64          ┆ i64                ┆ f64         │
╞═══════════╪════════╪══════════════╪════════════════════╪═════════════╡
│ 1962.2917 ┆ 6997.0 ┆ 0            ┆ 0                  ┆ -93.321321  │
│ 1962.375  ┆ 6966.0 ┆ 0            ┆ 0                  ┆ -124.321321 │
│ 1962.4583 ┆ 6981.0 ┆ 0            ┆ 0                  ┆ -109.321321 │
│ 1962.5417 ┆ 7003.0 ┆ 0            ┆ 0                  ┆ -87.321321  │
│ 1962.625  ┆ 7030.0 ┆ 0            ┆ 0                  ┆ -60.321321  │
└───────────┴────────┴──────────────┴────────────────────┴─────────────┘

Summary of the monthly data:
shape: (9, 6)
┌────────────┬───────────┬─────────────┬──────────────┬────────────────────┬─────────────┐
│ statistic  ┆ time      ┆ sl          ┆ missing_days ┆ flag_for_attention ┆ sl_anomaly  │
│ ---        ┆ ---       ┆ ---         ┆ ---          ┆ ---                ┆ ---         │
│ str        ┆ f64       ┆ f64         ┆ f64          ┆ f64                ┆ f64         │
╞════════════╪═══════════╪═════════════╪══════════════╪════════════════════╪═════════════╡
│ count      ┆ 753.0     ┆ 666.0       ┆ 753.0        ┆ 753.0              ┆ 666.0       │
│ null_count ┆ 0.0       ┆ 87.0        ┆ 0.0          ┆ 0.0                ┆ 87.0        │
│ mean       ┆ 1993.625  ┆ 7090.321321 ┆ 3.549801     ┆ 0.0                ┆ 1.2837e-13  │
│ std        ┆ 18.126389 ┆ 67.152763   ┆ 18.419548    ┆ 0.0                ┆ 67.152763   │
│ min        ┆ 1962.2917 ┆ 6908.0      ┆ 0.0          ┆ 0.0                ┆ -182.321321 │
│ 25%        ┆ 1977.9583 ┆ 7042.0      ┆ 0.0          ┆ 0.0                ┆ -48.321321  │
│ 50%        ┆ 1993.625  ┆ 7089.0      ┆ 0.0          ┆ 0.0                ┆ -1.321321   │
│ 75%        ┆ 2009.2917 ┆ 7137.0      ┆ 0.0          ┆ 0.0                ┆ 46.678679   │
│ max        ┆ 2024.9583 ┆ 7286.0      ┆ 99.0         ┆ 0.0                ┆ 195.678679  │
└────────────┴───────────┴─────────────┴──────────────┴────────────────────┴─────────────┘
<Figure size 1000x600 with 0 Axes>
<Figure size 1000x600 with 0 Axes>
--- Yearly Puerto Rico Analysis ---
Head of the yearly data:
shape: (5, 5)
┌────────┬────────┬──────────────┬────────────────────┬────────────┐
│ time   ┆ sl     ┆ missing_days ┆ flag_for_attention ┆ sl_anomaly │
│ ---    ┆ ---    ┆ ---          ┆ ---                ┆ ---        │
│ f64    ┆ f64    ┆ str          ┆ i64                ┆ f64        │
╞════════╪════════╪══════════════╪════════════════════╪════════════╡
│ 1963.0 ┆ 7038.0 ┆ N            ┆ 0                  ┆ -52.019608 │
│ 1964.0 ┆ 7071.0 ┆ N            ┆ 0                  ┆ -19.019608 │
│ 1965.0 ┆ 7012.0 ┆ N            ┆ 0                  ┆ -78.019608 │
│ 1966.0 ┆ 7044.0 ┆ N            ┆ 0                  ┆ -46.019608 │
│ 1967.0 ┆ 7043.0 ┆ N            ┆ 0                  ┆ -47.019608 │
└────────┴────────┴──────────────┴────────────────────┴────────────┘

Summary of the yearly data:
shape: (9, 6)
┌────────────┬───────────┬─────────────┬──────────────┬────────────────────┬─────────────┐
│ statistic  ┆ time      ┆ sl          ┆ missing_days ┆ flag_for_attention ┆ sl_anomaly  │
│ ---        ┆ ---       ┆ ---         ┆ ---          ┆ ---                ┆ ---         │
│ str        ┆ f64       ┆ f64         ┆ str          ┆ f64                ┆ f64         │
╞════════════╪═══════════╪═════════════╪══════════════╪════════════════════╪═════════════╡
│ count      ┆ 62.0      ┆ 51.0        ┆ 62           ┆ 62.0               ┆ 51.0        │
│ null_count ┆ 0.0       ┆ 11.0        ┆ 0            ┆ 0.0                ┆ 11.0        │
│ mean       ┆ 1993.5    ┆ 7090.019608 ┆ null         ┆ 0.0                ┆ 1.7833e-14  │
│ std        ┆ 18.041619 ┆ 46.252131   ┆ null         ┆ 0.0                ┆ 46.252131   │
│ min        ┆ 1963.0    ┆ 6977.0      ┆ N            ┆ 0.0                ┆ -113.019608 │
│ 25%        ┆ 1978.0    ┆ 7058.0      ┆ null         ┆ 0.0                ┆ -32.019608  │
│ 50%        ┆ 1994.0    ┆ 7085.0      ┆ null         ┆ 0.0                ┆ -5.019608   │
│ 75%        ┆ 2009.0    ┆ 7122.0      ┆ null         ┆ 0.0                ┆ 31.980392   │
│ max        ┆ 2024.0    ┆ 7193.0      ┆ Y            ┆ 0.0                ┆ 102.980392  │
└────────────┴───────────┴─────────────┴──────────────┴────────────────────┴─────────────┘
<Figure size 1000x600 with 0 Axes>
<Figure size 1000x600 with 0 Axes>
In [30]:
# --- Comparing Other Stations ---
print("\n--- Comparing Other Stations ---")

# Call the function for each city
cities = {
    1001: "San Juan, Puerto Rico",
    1545: "Tokyo, Japan",
    43: "Mumbai, India",
    12: "New York, New York, USA",
    723: "Dalian, China",
    331: "Takoradi, Ghana",
    245: "Los Angeles, California, USA",
    369: "Kolkata, India",
    832: "Buenos Aires, Argentina",
    196: "Sydney, Australia"
}

for dataset_id, city_name in cities.items():
    fig = analyze_city(dataset_id, city_name)
    if fig:
        fig.show()

print("\n--- Analysis Complete ---")
--- Comparing Other Stations ---

--- Analyzing San Juan, Puerto Rico (Dataset ID: 1001) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing Tokyo, Japan (Dataset ID: 1545) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing Mumbai, India (Dataset ID: 43) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing New York, New York, USA (Dataset ID: 12) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing Dalian, China (Dataset ID: 723) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing Takoradi, Ghana (Dataset ID: 331) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing Los Angeles, California, USA (Dataset ID: 245) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing Kolkata, India (Dataset ID: 369) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing Buenos Aires, Argentina (Dataset ID: 832) ---
<Figure size 1000x500 with 0 Axes>
--- Analyzing Sydney, Australia (Dataset ID: 196) ---
<Figure size 1000x500 with 0 Axes>
--- Analysis Complete ---
In [28]:
analyze_city(196, "Sydney, Australia")
--- Analyzing Sydney, Australia (Dataset ID: 196) ---
<Figure size 1000x500 with 0 Axes>