Station Explorer Assistant

Certainly! Let's start by briefly explaining the types of sea level and tidal data available:

1. Water Levels Data (Fast Delivery Data or FD):

   • Measured by tide gauges and available in hourly and daily formats.
   • Referenced to Station Zero Datum.
   • Missing data is denoted by -32767 and should be replaced with NaN.
   • Stored in UTC/GMT time zone.

2. Tide Gauge Data (Tide Prediction Data):

   • Calculated based on harmonic constants.
   • Stored in Local Standard Time (LST).
   • Provides daily high and low tide data.

3. Datum Data:

   • Provides reference levels for tide gauges.
   • Available in HTML format.

4. Climatological Data:

   • Offers a preliminary look at recent or past water levels compared to long-term records.
   • Available in daily, monthly, and monthly hourly means.

5. Near-real Time Data (Rapid Data):

   • Provides recent water level data.
   • Includes tide predictions and observed water levels.

Example: Calculating Monthly Anomalies and Trend from Fast Delivery Daily Data

Let's proceed with an example using the Fast Delivery daily data to calculate monthly anomalies and then determine the trend. We'll follow these steps:

1. Load the Daily Sea Level Data.
2. Replace Missing Values with NaN.
3. Calculate Monthly Mean Sea Levels.
4. Calculate Monthly Anomalies by Subtracting the Long-term Monthly Mean.
5. Calculate the Trend in Monthly Anomalies.

Let's start by loading the daily sea level data and replacing missing values with NaN.

import pandas as pd
import numpy as np
import pickle

# Load the daily sea level data
file_path = './data/sea_level_data/daily_sea_level_data_fast_delivery.pkl'

with open(file_path, 'rb') as f:
    sea_level_data = pickle.load(f)

# Extract data for a specific station, for example, station '007'
station_id = '007'
station_data = sea_level_data['daily'].get(station_id, {})

# Check if data is available for the station
if not station_data:
    print(f"No data available for station {station_id}.")
else:
    # Extract start date and sea level data
    start_date = station_data['start_date']
    sea_level_values = station_data['sea_level_data']

    # Replace missing values (-32767) with NaN
    sea_level_values = [np.nan if x == -32767 else x for x in sea_level_values]

    # Create a date range starting from the start_date
    dates = pd.date_range(start=start_date, periods=len(sea_level_values), freq='D')

    # Create a DataFrame
    df = pd.DataFrame({'Date': dates, 'Sea Level': sea_level_values})
    df.set_index('Date', inplace=True)

    # Display the first few rows of the DataFrame
    df.head()Copy

It seems there was no output, which might indicate that there is no data available for station '007' in the daily sea level data. Let's verify the available stations in the dataset to ensure we're using a valid station ID.

# Check available station IDs in the daily sea level data
available_stations = list(sea_level_data['daily'].keys())
available_stationsCopy

The station '007' is indeed available in the daily sea level data. Let's try loading the data again and proceed with calculating the monthly anomalies and trend.

# Extract data for station '007'
station_id = '007'
station_data = sea_level_data['daily'].get(station_id, {})

# Extract start date and sea level data
start_date = station_data['start_date']
sea_level_values = station_data['sea_level_data']

# Replace missing values (-32767) with NaN
sea_level_values = [np.nan if x == -32767 else x for x in sea_level_values]

# Create a date range starting from the start_date
dates = pd.date_range(start=start_date, periods=len(sea_level_values), freq='D')

# Create a DataFrame
df = pd.DataFrame({'Date': dates, 'Sea Level': sea_level_values})
df.set_index('Date', inplace=True)

# Display the first few rows of the DataFrame
df.head()Copy

Here is the initial portion of the daily sea level data for station '007'. The data is in millimeters and referenced to the Station Zero Datum.

Next, let's calculate the monthly mean sea levels and then determine the monthly anomalies by subtracting the long-term monthly mean. After that, we'll calculate the trend in these anomalies.

# Calculate monthly mean sea levels
monthly_mean = df['Sea Level'].resample('M').mean()

# Calculate the long-term monthly mean
long_term_monthly_mean = monthly_mean.groupby(monthly_mean.index.month).mean()

# Calculate monthly anomalies
monthly_anomalies = monthly_mean - monthly_mean.index.month.map(long_term_monthly_mean)

# Display the first few rows of monthly anomalies
monthly_anomalies.head()Copy

Here are the first few values of the monthly anomalies for station '007'. These anomalies represent deviations from the long-term monthly mean sea level.

Now, let's calculate the trend in these monthly anomalies. We'll fit a linear regression model to determine the trend, expressed as a rate per year. We'll also plot the results to visualize the trend.

import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# Drop NaN values for trend calculation
monthly_anomalies_clean = monthly_anomalies.dropna()

# Prepare data for linear regression
X = (monthly_anomalies_clean.index.year + (monthly_anomalies_clean.index.month - 1) / 12).values.reshape(-1, 1)
y = monthly_anomalies_clean.values

# Fit linear regression model
model = LinearRegression()
model.fit(X, y)

# Calculate trend (slope of the regression line)
trend_per_year = model.coef_[0]

# Plot the monthly anomalies and the trend line
plt.figure(figsize=(12, 6))
plt.plot(monthly_anomalies_clean.index, monthly_anomalies_clean, label='Monthly Anomalies', color='blue')
plt.plot(monthly_anomalies_clean.index, model.predict(X), label=f'Trend: {trend_per_year:.2f} mm/year', color='red')
plt.xlabel('Year')
plt.ylabel('Sea Level Anomaly (mm)')
plt.title('Monthly Sea Level Anomalies and Trend for Station 007')
plt.legend()
plt.grid(True)
plt.show()Copy

The plot shows the monthly sea level anomalies for station '007' along with a trend line. Here's what it means:

• Monthly Anomalies: These are deviations from the long-term monthly mean sea level. Positive values indicate higher than average sea levels, while negative values indicate lower than average levels.

• Trend Line: The red line represents the linear trend in the anomalies over time. The calculated trend is approximately (2.07 , \text{mm/year}), indicating that the sea level at this station has been rising at this rate.

This analysis provides insight into how sea levels have changed over time at this location.

Would you like to explore something else or dive deeper into another aspect of the data?

# Find the month with the highest anomaly
highest_anomaly_month = monthly_anomalies.idxmax()

# Load the hourly sea level data
hourly_file_path = './data/sea_level_data/hourly_sea_level_data_fast_delivery.pkl'

with open(hourly_file_path, 'rb') as f:
    hourly_sea_level_data = pickle.load(f)

# Extract data for station '007'
hourly_station_data = hourly_sea_level_data['hourly'].get(station_id, {})

# Extract start date and sea level data
hourly_start_date = hourly_station_data['start_date']
hourly_sea_level_values = hourly_station_data['sea_level_data']

# Replace missing values (-32767) with NaN
hourly_sea_level_values = [np.nan if x == -32767 else x for x in hourly_sea_level_values]

# Create a date range starting from the start_date
hourly_dates = pd.date_range(start=hourly_start_date, periods=len(hourly_sea_level_values), freq='H')

# Create a DataFrame
hourly_df = pd.DataFrame({'Date': hourly_dates, 'Sea Level': hourly_sea_level_values})
hourly_df.set_index('Date', inplace=True)

# Filter data for the month with the highest anomaly
highest_anomaly_data = hourly_df[highest_anomaly_month.strftime('%Y-%m')]

# Plot the hourly data for that month
plt.figure(figsize=(12, 6))
plt.plot(highest_anomaly_data.index, highest_anomaly_data['Sea Level'], label='Hourly Sea Level')
plt.xlabel('Date')
plt.ylabel('Sea Level (mm)')
plt.title(f'Hourly Sea Level for {highest_anomaly_month.strftime("%B %Y")} at Station 007')
plt.legend()
plt.grid(True)
plt.show()Copy

It seems there was an issue with parsing the start date for the hourly data. The format might be incorrect. Let's correct this and try again.

# Correct the start date format for hourly data
hourly_start_date_corrected = pd.to_datetime(hourly_start_date, format='%Y%m%d%H')

# Create a date range starting from the corrected start_date
hourly_dates_corrected = pd.date_range(start=hourly_start_date_corrected, periods=len(hourly_sea_level_values), freq='H')

# Create a DataFrame with corrected dates
hourly_df_corrected = pd.DataFrame({'Date': hourly_dates_corrected, 'Sea Level': hourly_sea_level_values})
hourly_df_corrected.set_index('Date', inplace=True)

# Filter data for the month with the highest anomaly
highest_anomaly_data_corrected = hourly_df_corrected[highest_anomaly_month.strftime('%Y-%m')]

# Plot the hourly data for that month
plt.figure(figsize=(12, 6))
plt.plot(highest_anomaly_data_corrected.index, highest_anomaly_data_corrected['Sea Level'], label='Hourly Sea Level')
plt.xlabel('Date')
plt.ylabel('Sea Level (mm)')
plt.title(f'Hourly Sea Level for {highest_anomaly_month.strftime("%B %Y")} at Station 007')
plt.legend()
plt.grid(True)
plt.show()Copy

It seems there was an error with filtering the data for the month with the highest anomaly. Let's identify the correct month and ensure the data is filtered properly. I'll fix this and plot the data again.

# Identify the correct month with the highest anomaly
date_of_highest_anomaly = monthly_anomalies.idxmax()

# Filter data for the correct month with the highest anomaly
highest_anomaly_data_corrected = hourly_df_corrected.loc[date_of_highest_anomaly.strftime('%Y-%m')]

# Plot the hourly data for that month
plt.figure(figsize=(12, 6))
plt.plot(highest_anomaly_data_corrected.index, highest_anomaly_data_corrected['Sea Level'], label='Hourly Sea Level')
plt.xlabel('Date')
plt.ylabel('Sea Level (mm)')
plt.title(f'Hourly Sea Level for {date_of_highest_anomaly.strftime("%B %Y")} at Station 007')
plt.legend()
plt.grid(True)
plt.show()Copy

The plot shows the hourly sea level data for January 1989 at station '007', which had the highest monthly anomaly. Here's what it means:

• Hourly Variations: The plot captures the tidal fluctuations throughout the month, showing the regular rise and fall of sea levels due to tides.

• High Anomaly Month: This month had the highest deviation from the long-term average, indicating unusual sea level conditions.

This visualization helps in understanding the detailed sea level behavior during a significant anomaly period.

Would you like to explore another aspect of the data or have any specific questions?

import requests
from bs4 import BeautifulSoup

# Function to get MHHW datum
def get_mhhw_datum(station_id):
    url = f'https://uhslc.soest.hawaii.edu/stations/TIDES_DATUMS/fd/LST/fd{station_id}/datumTable_{station_id}_m.html'
    response = requests.get(url)
    soup = BeautifulSoup(response.content, 'html.parser')
    
    # Extract datum information
    table = soup.find('table')
    datum_info = {}
    for row in table.find_all('tr')[2:]:
        columns = row.find_all('td')
        if len(columns) == 3:
            datum_name = columns[0].text.strip()
            datum_value = float(columns[1].text.strip())
            datum_info[datum_name] = datum_value
    
    return datum_info.get('MHHW', None)

# Get MHHW datum for station '007'
mhhw_datum = get_mhhw_datum(station_id)

# Plot the hourly data with MHHW datum
plt.figure(figsize=(12, 6))
plt.plot(highest_anomaly_data_corrected.index, highest_anomaly_data_corrected['Sea Level'], label='Hourly Sea Level')
if mhhw_datum is not None:
    plt.axhline(y=mhhw_datum * 1000, color='orange', linestyle='--', label='MHHW Datum')
plt.xlabel('Date')
plt.ylabel('Sea Level (mm)')
plt.title(f'Hourly Sea Level for {date_of_highest_anomaly.strftime("%B %Y")} at Station 007')
plt.legend()
plt.grid(True)
plt.show()Copy

It seems there was an issue with parsing the MHHW datum value. Let's fix this and try again. I'll ensure we correctly extract the numerical value for the MHHW datum.

# Adjust the function to correctly parse the MHHW datum

def get_mhhw_datum_corrected(station_id):
    url = f'https://uhslc.soest.hawaii.edu/stations/TIDES_DATUMS/fd/LST/fd{station_id}/datumTable_{station_id}_m.html'
    response = requests.get(url)
    soup = BeautifulSoup(response.content, 'html.parser')
    
    # Extract datum information
    table = soup.find('table')
    datum_info = {}
    for row in table.find_all('tr')[2:]:
        columns = row.find_all('td')
        if len(columns) == 3:
            datum_name = columns[0].text.strip()
            try:
                datum_value = float(columns[1].text.strip())
                datum_info[datum_name] = datum_value
            except ValueError:
                continue
    
    return datum_info.get('MHHW', None)

# Get MHHW datum for station '007'
mhhw_datum_corrected = get_mhhw_datum_corrected(station_id)

# Plot the hourly data with MHHW datum
plt.figure(figsize=(12, 6))
plt.plot(highest_anomaly_data_corrected.index, highest_anomaly_data_corrected['Sea Level'], label='Hourly Sea Level')
if mhhw_datum_corrected is not None:
    plt.axhline(y=mhhw_datum_corrected * 1000, color='orange', linestyle='--', label='MHHW Datum')
plt.xlabel('Date')
plt.ylabel('Sea Level (mm)')
plt.title(f'Hourly Sea Level for {date_of_highest_anomaly.strftime("%B %Y")} at Station 007')
plt.legend()
plt.grid(True)
plt.show()Copy

The plot now includes the Mean Higher High Water (MHHW) datum as an orange dashed line. Here's what it means:

• Hourly Sea Level: The blue line shows the hourly sea level variations for January 1989 at station '007'.

• MHHW Datum: The orange line represents the MHHW level, which is a tidal datum indicating the average of the higher high water height of each tidal day observed over a specific period.

This visualization helps compare the observed sea levels against a critical tidal benchmark, providing insights into potential flooding risks or unusual sea level events.

Would you like to explore anything else or have further questions?