Velocity model–ADCP comparison

import subprocess
import os
import pandas as pd
import netCDF4
import numpy as np
import glob
import time
import matplotlib.pyplot as plt
import copy
import xarray as xr
from datetime import datetime, timedelta 
import dask
from scipy.interpolate import griddata
#from ocean_c_lab_tools import *
#from celluloid import Camera 
#import PyCO2SYS as csys
import seawater as sw
from roms_regrid import *

PROJECT_ROOT = "/anvil/projects/x-ees250129/x-uheede/C-Star-in-Hvalfjordur"
MODEL_ROOT = "/anvil/scratch/x-uheede/Iceland2_MARBL_2024_60m"
GRID_ROOT = "/home/x-uheede/S/Iceland2_MARBL_2024_60m/P_INPUT"

obs_adcp_path = f"{PROJECT_ROOT}/data/staged/seanoe_113246/2026-04-21/adcp_mooring_qc/Mooring_data_Final_.nc/HV??1*"
model_grid_path = f"{GRID_ROOT}/Iceland2_grid.nc"
model_data_path = f"{MODEL_ROOT}/Iceland2_MARBL_2024_avg.20240[7]????????.nc"

target_depth_levels = [10]  # Specify depth levels of interest
thinner = 6  # Sample one point per hour to save compute time

from roms_tools import Grid, ROMSOutput
grid = Grid.from_file(
    model_grid_path
)

#Only run this cell if grid is made in MATLAB
grid.update_vertical_coordinate(N=vert_levels, theta_s=theta_s_model, theta_b=theta_b_model, hc=hc_model, verbose=False)

import xarray as xr
import numpy as np
import glob

# =========================
# --- FIND FILES ---
# =========================
files = sorted(glob.glob(obs_adcp_path))

datasets = {}

for f in files:
    name = f.split("/")[-1].replace(".nc", "")
    ds = xr.open_dataset(f)

    datasets[name] = ds

def apply_qc(ds):

    good_velocity = ds["WECT_flag"] <= 1
    good_depth = ds["DEPH_QC_flag"] <= 1
    
    physical_depth = (ds["DEPH"] > 0) & (ds["DEPH"] < 150)

    # Combined mask
    good_total = good_velocity & good_depth & physical_depth

    # 2. Apply the mask to your variables
    ds["WECT_clean"] = ds["WECT"].where(good_total)
    ds["NSCT_clean"] = ds["NSCT"].where(good_total)
    ds["VCSP_clean"] = ds["VCSP"].where(good_total)

    return ds

datasets = {k: apply_qc(v) for k, v in datasets.items()}

def standardize(ds):

    ds = ds.rename({
        "TIME": "time",
        "DEPTH": "depth",
        "LATITUDE": "lat",
        "LONGITUDE": "lon",
        "WECT_clean": "u",
        "NSCT_clean": "v"
    })

    # flatten singleton coords
    ds = ds.assign_coords(
        lat=ds["lat"].values[0],
        lon=ds["lon"].values[0] + 360  # match model grid if needed
    )

    return ds

datasets = {k: standardize(v) for k, v in datasets.items()}

ds=datasets['HVIN1_ADCP_64m_20240404_20240814']
ds.DEPH.values

from roms_tools import Grid, ROMSOutput

import xarray as xr
import numpy as np

# Load ROMS output using your pattern
roms_output = ROMSOutput(
    grid=grid,
    path=[
        model_data_path,
    ],
    use_dask=True,
)

ds = roms_output.regrid(var_names=["u", "v"],depth_levels=target_depth_levels)

u_rg = ds['u'].thin({'time': thinner}).load()
v_rg = ds['v'].thin({'time': thinner}).load()

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from windrose import WindroseAxes

# Define the number of stations
station_labels = ['HVIN', 'HVNA', 'HVSA', 'HVNV', 'HVSV']
num_stations = 5

# Define the observation subset time range (2024)
obs_time_range = slice('2024-07-05', '2024-07-31')

# Define the model subset time range (2023)
model_time_range = slice('2024-07-05', '2024-07-31')

# Create figure for timeseries plots
fig, axs = plt.subplots(
    nrows=num_stations, ncols=2, figsize=(12, 3 * num_stations), sharex=False
)
fig.suptitle("Timeseries Comparison: Observations vs Model", fontsize=14)

for i in range(num_stations):

    # Subset observation data for 2024
    obs = list(datasets.values())[i]
    obs_subset = obs.sel(time=obs_time_range)
    # --- use physical depths ---
    depths = obs_subset['DEPH'].values
    # find bin closest to 10 m
    depth_index = np.abs(depths - 10).argmin()
    u_obs = obs_subset['u'].isel(depth=depth_index)[::6]
    obs_lon = obs.lon
    obs_lat = obs.lat

    # Extract model data
    u_sel = u_rg.sel(lon=obs_lon.values, method='nearest').sel(lat=obs_lat.values, method='nearest')
    u_mod = u_sel.squeeze()

    # ----------------------------
    # Compute shared y-limits
    # ----------------------------
    ymin = min(u_obs.min().item(), u_mod.min().item())
    ymax = max(u_obs.max().item(), u_mod.max().item())

    # Add small padding
    yrange = ymax - ymin
    ymin -= 0.05 * yrange
    ymax += 0.05 * yrange

    # Plot observations (2024)
    axs[i, 0].plot(obs_subset.time[::6], u_obs, label="Observed U", color="blue")
    axs[i, 0].set_title(f"Station {station_labels[i]} - Observations")

    axs[i, 0].tick_params(axis='x', rotation=45)
    axs[i, 0].set_ylim([ymin, ymax])
    axs[i, 0].legend()

    # Plot model (2023)
    axs[i, 1].plot(u_rg['time'], u_mod, label="Model U", color="green")
    axs[i, 1].set_title(f"Station {station_labels[i]} - Model")
    axs[i, 1].tick_params(axis='x', rotation=45)
    axs[i, 1].set_ylim([ymin, ymax])
    axs[i, 1].legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from windrose import WindroseAxes

# -----------------------------
# Select station 4
# -----------------------------
station_labels = [10, 8, 6, 4, 2]
station_idx = station_labels.index(4)

obs_time_range = slice('2024-07-05', '2024-07-31')

obs = list(datasets.values())[station_idx]
obs_subset = obs.sel(time=obs_time_range)

# --- use physical depth (DEPH) ---
depths = obs_subset['DEPH'].values
depth_index = np.abs(depths - 10).argmin()

u_obs = obs_subset['u'].isel(depth=depth_index)[::6]
v_obs = obs_subset['v'].isel(depth=depth_index)[::6]

obs_lon = obs.lon
obs_lat = obs.lat

u_sel = u_rg.sel(lon=obs_lon.values, lat=obs_lat.values, method='nearest')
v_sel = v_rg.sel(lon=obs_lon.values, lat=obs_lat.values, method='nearest')

u_mod = u_sel.squeeze()
v_mod = v_sel.squeeze()

# -----------------------------
# Wind speed & direction
# -----------------------------
direction_obs = (np.arctan2(-u_obs.values, -v_obs.values) * 180 / np.pi + 360) % 360
speed_obs = np.sqrt(u_obs.values**2 + v_obs.values**2)

direction_model = (np.arctan2(-u_mod.values, -v_mod.values) * 180 / np.pi + 360) % 360
speed_model = np.sqrt(u_mod.values**2 + v_mod.values**2)

# -----------------------------
# Shared y-limits
# -----------------------------
ymin = min(u_obs.min().item(), u_mod.min().item())
ymax = max(u_obs.max().item(), u_mod.max().item())
yrange = ymax - ymin
ymin -= 0.05 * yrange
ymax += 0.05 * yrange

# -----------------------------
# Figure + GridSpec
# -----------------------------
fig = plt.figure(figsize=(16, 10))
fig.suptitle("Mooring HVNV: Observations vs Model", fontsize=16)

gs = GridSpec(
    nrows=2, ncols=2,
    height_ratios=[1, 2],   # ⬅️ top short, bottom tall
    hspace=0.35,
    wspace=0.25
)

# ---- Timeseries (wide & short) ----
ax_ts_obs = fig.add_subplot(gs[0, 0])
ax_ts_mod = fig.add_subplot(gs[0, 1])

ax_ts_obs.plot(obs_subset.time[::6], u_obs, color='blue', label='Observed U')
ax_ts_obs.set_title("HVNV – Obs (Timeseries)")
ax_ts_obs.set_ylim(ymin, ymax)
ax_ts_obs.set_ylabel('m/s')
ax_ts_obs.tick_params(axis='x', rotation=45)
ax_ts_obs.legend()

# panel label
ax_ts_obs.text(0.01, 0.95, 'a)', transform=ax_ts_obs.transAxes,
               fontsize=14, fontweight='bold', va='top')

ax_ts_mod.plot(u_rg['time'], u_mod, color='green', label='Model U')
ax_ts_mod.set_title("HVNV – Model (Timeseries)")
ax_ts_mod.set_ylim(ymin, ymax)
ax_ts_mod.set_ylabel('m/s')
ax_ts_mod.tick_params(axis='x', rotation=45)
ax_ts_mod.legend()

# panel label
ax_ts_mod.text(0.01, 0.95, 'b)', transform=ax_ts_mod.transAxes,
               fontsize=14, fontweight='bold', va='top')

# -----------------------------
# Wind roses (larger area)
# -----------------------------
ax_wr_obs = WindroseAxes(fig, [0.08, 0.06, 0.38, 0.45])
fig.add_axes(ax_wr_obs)

ax_wr_obs.bar(direction_obs, speed_obs,
              normed=True, opening=0.8, edgecolor='white',
              bins=np.linspace(0, 0.6, 7))
ax_wr_obs.set_title("HVNV – Obs (Current Rose)")
ax_wr_obs.set_legend(
    title="Current\nSpeed (m s$^{-1}$)",
    loc="lower right",
    bbox_to_anchor=(1.2, 0.0)
)
ax_wr_obs.set_rticks([5, 10, 15, 20, 25, 30])
ax_wr_obs.set_yticklabels(['5%', '10%', '15%', '20%', '25%', '30%'])

# panel label
ax_wr_obs.text(0.02, 0.98, 'c)', transform=ax_wr_obs.transAxes,
               fontsize=14, fontweight='bold', va='top')

ax_wr_mod = WindroseAxes(fig, [0.54, 0.06, 0.38, 0.45])
fig.add_axes(ax_wr_mod)

ax_wr_mod.bar(direction_model, speed_model,
              normed=True, opening=0.8, edgecolor='white',
              bins=np.linspace(0, 0.6, 7))
ax_wr_mod.set_title("HVNV – Model (Current Rose)")


ax_wr_mod.set_legend(
    title="Current\nSpeed (m s$^{-1}$)",
    loc="lower right",
    bbox_to_anchor=(1.2, 0.0)
)
ax_wr_mod.set_rticks([5, 10, 15, 20, 25, 30])
ax_wr_mod.set_yticklabels(['5%', '10%', '15%', '20%', '25%', '30%'])

# panel label
ax_wr_mod.text(0.02, 0.98, 'd)', transform=ax_wr_mod.transAxes,
               fontsize=14, fontweight='bold', va='top')

plt.show()

import numpy as np
import matplotlib.pyplot as plt
from windrose import WindroseAxes
station_labels = ['HVIN', 'HVNA', 'HVSA', 'HVNV', 'HVSV']

# Define number of stations
num_stations = 4

# Create a figure with 2 columns (Obs & Model) and num_stations rows
fig = plt.figure(figsize=(14, 7 * num_stations),layout='tight')

fig.suptitle("Wind Roses: Observations (2024) vs Model (2024)", fontsize=16)

for i in range(num_stations):

    # Compute observed wind speed and direction
    obs=list(datasets.values())[i]
    obs_subset = obs.sel(time=obs_time_range)
    depths = obs_subset['DEPH'].values
    depth_index = np.abs(depths - 10).argmin()
    u_obs = obs_subset['u'].sel(depth=depth_index)[::6]
    v_obs = obs_subset['v'].sel(depth=depth_index)[::6]
    obs_lon=obs.lon
    obs_lat=obs.lat
    direction_obs = (np.arctan2(-u_obs.squeeze().values, -v_obs.squeeze().values) * (180 / np.pi) + 360) % 360
    speed_obs = np.sqrt(u_obs.squeeze().values**2 + v_obs.squeeze().values**2)

    # Compute model wind speed and direction
    u_sel = u_rg.sel(lon=obs_lon.values, method='nearest').sel(lat=obs_lat.values, method='nearest')
    v_sel = v_rg.sel(lon=obs_lon.values, method='nearest').sel(lat=obs_lat.values, method='nearest')
    
    direction_model = (np.arctan2(-u_sel.squeeze().values, -v_sel.squeeze().values) * (180 / np.pi) + 360) % 360
    speed_model = np.sqrt(u_sel.squeeze().values**2 + v_sel.squeeze().values**2)

    # Define subplot positions using bbox (left, bottom, width, height)
    ax_obs = WindroseAxes(fig, [0.05, 0.75 - (i * 0.2), 0.4, 0.18])  # Left column
    fig.add_axes(ax_obs)
    ax_obs.bar(direction_obs, speed_obs, normed=True, opening=0.8, edgecolor='white', bins=np.linspace(0, 0.6, 7))
    ax_obs.set_title(f"Station {station_labels[i]} - Obs (2024)", fontsize=12)
    ax_obs.set_legend()
    ax_obs.set_rticks([5, 10, 15, 20, 25, 30])
    ax_obs.set_yticklabels(['5%', '10%', '15%', '20%', '25%', '30%'])
    
    ax_model = WindroseAxes(fig, [0.55, 0.75 - (i * 0.2), 0.4, 0.18])  # Right column
    fig.add_axes(ax_model)
    ax_model.bar(direction_model, speed_model, normed=True, opening=0.8, edgecolor='white', bins=np.linspace(0, 0.6, 7))
    ax_model.set_title(f"Station {station_labels[i]} - Model (2024)", fontsize=12)
    ax_model.set_legend()
    ax_model.set_rticks([5, 10, 15, 20, 25, 30])
    ax_model.set_yticklabels(['5%', '10%', '15%', '20%', '25%', '30%'])

plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks

import numpy as np
from scipy.signal import find_peaks


def compute_spectrum(time, signal):
    """Compute normalized power spectral density and identify spectral peaks."""

    # Remove NaNs
    mask = np.isfinite(signal) & np.isfinite(time)
    time = time[mask]
    signal = signal[mask]

    # Remove mean (important for spectral analysis)
    signal = signal - np.mean(signal)

    # FFT
    n = len(signal)
    dt = np.mean(np.diff(time))  # time step in hours
    freq = np.fft.rfftfreq(n, d=dt)
    spectrum = np.fft.rfft(signal)

    # Normalized PSD
    power = (np.abs(spectrum) ** 2) / (n * dt)

    # Peak detection
    peaks, _ = find_peaks(power, height=np.max(power) * 0.1)

    return freq, power, peaks
    
# Define reference time as a numpy datetime64 object
reference_time = np.datetime64('2000-01-01T00:00:00')    
# Example time series (Replace these with your observed and modeled data)
obs_time = (obs_subset.time[::6].load() - reference_time).astype('timedelta64[h]').astype(int).load().values
#obs_time = (obs_subset.time.load() - reference_time).astype('timedelta64[h]').astype(int).load().values
# Define reference time (2000-01-01T00:00:00)



model_time = (u_rg['time'].load() - reference_time).astype('timedelta64[h]').astype(int).load().values# Example time axis (modify as needed)
obs_time=(obs_time-obs_time[0])*10**(-0)*1/3600
model_time=(model_time-model_time[0])*10**(-0)*1/3600

observed_signal = u_obs.squeeze().values  # Example signal
modeled_signal = u_sel.squeeze().values # Example signal

# Compute spectra
freq_obs, power_obs, peaks_obs = compute_spectrum(obs_time, observed_signal)
freq_mod, power_mod, peaks_mod = compute_spectrum(model_time, modeled_signal)

# Plot results
plt.figure(figsize=(10, 5))

# Observed Spectrum
plt.plot(freq_obs, power_obs, label="Observed Spectrum", color='b')
plt.scatter(freq_obs[peaks_obs], power_obs[peaks_obs], color='r', label="Observed Peaks", zorder=3)

# Modeled Spectrum
plt.plot(freq_mod, power_mod, label="Modeled Spectrum", color='g', linestyle='dashed')
plt.scatter(freq_mod[peaks_mod], power_mod[peaks_mod], color='orange', label="Modeled Peaks", zorder=3)

#plt.text(0.01, 0.98, 'e)', 
#         transform=plt.gca().transAxes, 
#         fontsize=14, fontweight='bold', 
#         va='top', ha='left')

plt.xlabel("Frequency (hr$^{-1}$)")
plt.ylabel("Power spectral density (m$^2$ s$^{-2}$ hr)")

plt.xlim(0,0.2)
plt.legend()
plt.title("Spectral Peaks in Observed and Modeled Data")
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks

# Convert frequency to period (avoid division by zero)
period_obs = 1 / freq_obs[1:]  # Ignore first value to avoid division by zero
period_mod = 1 / freq_mod[1:]

# Plot results
plt.figure(figsize=(10, 5))

# Observed Spectrum
plt.plot(period_obs, power_obs[1:], label="Observed Spectrum", color='b')
#plt.scatter(period_obs[peaks_obs - 1], power_obs[peaks_obs][1:], color='r', label="Observed Peaks", zorder=3)

# Modeled Spectrum
plt.plot(period_mod, power_mod[1:], label="Modeled Spectrum", color='g', linestyle='dashed')
#plt.scatter(period_mod[peaks_mod - 1], power_mod[peaks_mod][1:], color='orange', label="Modeled Peaks", zorder=3)

#plt.xscale("log")  # Log scale to better show long periods
plt.xlabel("Period (hours)")
plt.ylabel("Power")
plt.legend()
plt.title("Spectral Peaks in Observed and Modeled Data")
plt.grid(which="both", linestyle="--", linewidth=0.5)
plt.gca().invert_xaxis()  # Invert x-axis so long periods are on the right
plt.xlim(1,50)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def bandpass_filter(time, signal, min_period=10, max_period=15):
    """Filter the signal to retain only frequencies corresponding to periods between min_period and max_period hours."""
    
    # Compute FFT
    n = len(signal)
    dt = np.mean(np.diff(time))  # Time step in hours
    freq = np.fft.rfftfreq(n, d=dt)  # Frequency axis (in Hz or cycles per hour)
    spectrum = np.fft.rfft(signal)   # Compute FFT

    # Convert periods to frequencies
    min_freq = 1 / max_period  # Lower cutoff frequency
    max_freq = 1 / min_period  # Upper cutoff frequency

    # Zero out frequencies outside the band
    spectrum[(freq < min_freq) | (freq > max_freq)] = 0

    # Inverse FFT to reconstruct filtered signal
    filtered_signal = np.fft.irfft(spectrum, n=n)

    return filtered_signal

# Convert time to hours since reference time
reference_time = np.datetime64('2000-01-01T00:00:00', 's')
obs_time = (obs_subset.time[::6].load() - reference_time).astype('timedelta64[h]').astype(int).values
model_time = (u_rg['time'] - reference_time).astype('timedelta64[h]').astype(int).load().values
obs_time=(obs_time-obs_time[0])*10**(-0)*1/3600
model_time=(model_time-model_time[0])*10**(-0)*1/3600

# Apply band-pass filter
u_obs_filtered = bandpass_filter(obs_time, u_obs.squeeze().values, min_period=10, max_period=15)
u_model_filtered = bandpass_filter(model_time, u_sel.squeeze().values, min_period=10, max_period=15)

# Create figure and subplots
fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(12, 8), sharex=True)

# Plot filtered observed signal
axs[0].plot(obs_time, u_obs_filtered, label="Filtered Observed (10-15h)", color='b')
axs[0].set_ylabel("Velocity (m/s)")
axs[0].set_ylim(-.41,.41)
axs[0].set_title("Filtered Observed u_obs (10-15h Period)")
axs[0].legend()
axs[0].grid()

# Plot filtered modeled signal
axs[1].plot(model_time, u_model_filtered, label="Filtered Modeled (10-15h)", color='g')
axs[1].set_xlabel("Time (hours)")
axs[1].set_ylabel("Velocity (m/s)")
axs[1].set_ylim(-.41,.41)
axs[1].set_title("Filtered Modeled u_model (10-15h Period)")
axs[1].legend()
axs[1].grid()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def highpass_filter(time, signal, cutoff_period=10):
    """Remove frequencies with periods longer than `cutoff_period` hours (high-pass filter)."""

    # Compute FFT
    n = len(signal)
    dt = np.mean(np.diff(time))  # Time step in hours
    freq = np.fft.rfftfreq(n, d=dt)  # Frequency axis (in Hz or cycles per hour)
    spectrum = np.fft.rfft(signal)   # Compute FFT

    # Define cutoff frequency
    cutoff_freq = 1 / cutoff_period  # Frequency corresponding to the longest period to keep

    # Zero out frequencies below the cutoff (i.e., remove long periods)
    spectrum[freq < cutoff_freq] = 0

    # Inverse FFT to reconstruct the high-pass filtered signal
    filtered_signal = np.fft.irfft(spectrum, n=n)

    return filtered_signal

# Apply high-pass filter (removing periods ≥ 10 hours)
u_obs_highpass = highpass_filter(obs_time, u_obs.squeeze().values, cutoff_period=10)
u_model_highpass = highpass_filter(model_time, u_sel.squeeze().values, cutoff_period=10)

# Create figure and subplots
fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(12, 8), sharex=True)

# Plot high-pass filtered observed signal
axs[0].plot(obs_time, u_obs_highpass, label="Observed (Periods < 10h)", color='b')
axs[0].set_ylabel("Velocity (m/s)")
axs[0].set_title("High-Pass Filtered Observed u_obs (Periods < 10h)")
axs[0].set_ylim(-.31,.31)
axs[0].legend()
axs[0].grid()
axs[0].set_xlabel("Time (hours)")

# Plot high-pass filtered modeled signal
axs[1].plot(model_time, u_model_highpass, label="Modeled (Periods < 10h)", color='g')
axs[1].set_xlabel("Time (hours)")
axs[1].set_ylabel("Velocity (m/s)")
axs[1].set_ylim(-.31,.31)
axs[1].set_title("High-Pass Filtered Modeled u_model (Periods < 10h)")
axs[1].legend()
axs[1].grid()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def lowpass_filter(time, signal, cutoff_period=15):
    """Remove frequencies with periods longer than `cutoff_period` hours (high-pass filter)."""

    # Compute FFT
    n = len(signal)
    dt = np.mean(np.diff(time))  # Time step in hours
    freq = np.fft.rfftfreq(n, d=dt)  # Frequency axis (in Hz or cycles per hour)
    spectrum = np.fft.rfft(signal)   # Compute FFT

    # Define cutoff frequency
    cutoff_freq = 1 / cutoff_period  # Frequency corresponding to the longest period to keep

    # Zero out frequencies below the cutoff (i.e., remove long periods)
    spectrum[freq > cutoff_freq] = 0

    # Inverse FFT to reconstruct the high-pass filtered signal
    filtered_signal = np.fft.irfft(spectrum, n=n)

    return filtered_signal

# Apply high-pass filter (removing periods ≥ 10 hours)
u_obs_lowpass = lowpass_filter(obs_time, u_obs.squeeze().values, cutoff_period=15)
u_model_lowpass = lowpass_filter(model_time, u_sel.squeeze().values, cutoff_period=15)

# Create figure and subplots
fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(12, 8), sharex=True)

# Plot high-pass filtered observed signal
axs[0].plot(obs_time, u_obs_lowpass, label="Observed (Periods > 15h)", color='b')
axs[0].set_ylabel("Velocity (m/s)")
axs[0].set_title("Low-Pass Filtered Observed u_obs (Periods > 15h)")
axs[0].legend()
axs[0].set_ylim(-.31,.31)
axs[0].grid()
axs[0].set_xlabel("Time (hours)")

# Plot high-pass filtered modeled signal
axs[1].plot(model_time, u_model_lowpass, label="Modeled (Periods > 15h)", color='g')
axs[1].set_xlabel("Time (hours)")
axs[1].set_ylabel("Velocity (m/s)")
axs[1].set_ylim(-.31,.31)
axs[1].set_title("Low-Pass Filtered Modeled u_model (Periods > 15h)")
axs[1].legend()
axs[1].grid()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def bandpass_filter(time, signal, low_cutoff=10, high_cutoff=15):
    """Apply a bandpass filter to retain frequencies between `low_cutoff` and `high_cutoff` hours."""
    n = len(signal)
    dt = np.mean(np.diff(time))
    freq = np.fft.rfftfreq(n, d=dt)
    spectrum = np.fft.rfft(signal)

    # Define cutoff frequencies
    low_freq = 1 / high_cutoff
    high_freq = 1 / low_cutoff

    # Zero out frequencies outside the band
    spectrum[(freq < low_freq) | (freq > high_freq)] = 0

    return np.fft.irfft(spectrum, n=n)

def highpass_filter(time, signal, cutoff_period=10):
    """Remove frequencies with periods longer than `cutoff_period` hours (retain high frequencies)."""
    n = len(signal)
    dt = np.mean(np.diff(time))
    freq = np.fft.rfftfreq(n, d=dt)
    spectrum = np.fft.rfft(signal)

    # Define cutoff frequency
    cutoff_freq = 1 / cutoff_period  

    # Zero out frequencies below the cutoff (remove long periods)
    spectrum[freq < cutoff_freq] = 0

    return np.fft.irfft(spectrum, n=n)

def lowpass_filter(time, signal, cutoff_period=15):
    """Remove frequencies with periods shorter than `cutoff_period` hours (retain low frequencies)."""
    n = len(signal)
    dt = np.mean(np.diff(time))
    freq = np.fft.rfftfreq(n, d=dt)
    spectrum = np.fft.rfft(signal)

    # Define cutoff frequency
    cutoff_freq = 1 / cutoff_period  

    # Zero out frequencies above the cutoff (remove short periods)
    spectrum[freq > cutoff_freq] = 0

    return np.fft.irfft(spectrum, n=n)

# Apply filters to observed and modeled data
u_obs_bandpass = bandpass_filter(obs_time, u_obs.squeeze().values, low_cutoff=10, high_cutoff=15)
u_model_bandpass = bandpass_filter(model_time, u_sel.squeeze().values, low_cutoff=10, high_cutoff=15)

u_obs_highpass = highpass_filter(obs_time, u_obs.squeeze().values, cutoff_period=10)
u_model_highpass = highpass_filter(model_time, u_sel.squeeze().values, cutoff_period=10)

u_obs_lowpass = lowpass_filter(obs_time, u_obs.squeeze().values, cutoff_period=15)
u_model_lowpass = lowpass_filter(model_time, u_sel.squeeze().values, cutoff_period=15)

# Prepare data for boxplot
data = [
    u_obs_bandpass, u_model_bandpass,
    u_obs_highpass, u_model_highpass,
    u_obs_lowpass, u_model_lowpass
]

labels = [
    "Obs Bandpass (10-15h)", "Model Bandpass (10-15h)",
    "Obs High-pass (<10h)", "Model High-pass (<10h)",
    "Obs Low-pass (>15h)", "Model Low-pass (>15h)"
]

# Define colors (Obs = Light Green, Model = Light Blue)
colors = ['lightgreen', 'lightblue', 'lightgreen', 'lightblue', 'lightgreen', 'lightblue']

# Create boxplot
fig, ax = plt.subplots(figsize=(10, 6))
boxplot = ax.boxplot(data, labels=labels, patch_artist=True, medianprops=dict(color="red"))

# Apply colors to boxes
for patch, color in zip(boxplot['boxes'], colors):
    patch.set_facecolor(color)
    patch.set_alpha(0.6)  # Transparency for better readability

ax.set_ylabel("Velocity (m/s)")
ax.set_title("Boxplot of Filtered Velocity Components")

plt.xticks(rotation=20)  # Rotate labels for better readability
plt.grid(axis='y', linestyle='--', alpha=0.6)

plt.show()

import xarray as xr
import numpy as np
target_depth_levels=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,25,30,35,40,45,50]
# Load ROMS output using your pattern
roms_output = ROMSOutput(
    grid=grid,
    path=[
        model_data_path,
    ],
    use_dask=True,
)

ds = roms_output.regrid(var_names=["u", "v"],depth_levels=target_depth_levels)

u_mean=ds['u'].thin({'time': 48}).mean('time').load()
v_mean=ds['v'].thin({'time': 48}).mean('time').load()

import matplotlib.pyplot as plt

# Create figure and subplots
fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(20, 5), sharey=True)

# Loop through each dataset and create a subplot
for i, (key, obs) in enumerate(datasets.items()):
    ax = axes[i]  # Select subplot
    obs=list(datasets.values())[i]
    obs_lon=obs.lon
    obs_lat=obs.lat
    # Plot observed profile
    ax.plot(obs['u'].mean('time').squeeze(), obs['DEPH'], label="Observed", color='b')

    # Extract nearest model data
    u_sel = u_mean.sel(lon=obs_lon.values, method='nearest').sel(lat=obs_lat.values, method='nearest')

    # Plot modeled profile
    ax.plot(u_sel.squeeze(), u_mean.depth, label="Modeled", color='r')

    # Reverse y-axis
    ax.invert_yaxis()

    # Titles and labels
    ax.set_title(f"Station {i+1}")
    if i == 0:
        ax.set_ylabel("Depth (m)")
    ax.set_xlabel("Velocity (m/s)")
    ax.legend()

plt.tight_layout()
plt.show()