From 91e59f8645acc8b798173d5c8014dadc72155289 Mon Sep 17 00:00:00 2001 From: Allie Date: Thu, 10 Sep 2026 13:11:35 -0600 Subject: [PATCH 1/9] Rename h5outputsReader.py to h5_reader.py --- .../{run_seastack/h5outputsReader.py => features/h5_reader.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename data/demos/{run_seastack/h5outputsReader.py => features/h5_reader.py} (100%) diff --git a/data/demos/run_seastack/h5outputsReader.py b/data/demos/features/h5_reader.py similarity index 100% rename from data/demos/run_seastack/h5outputsReader.py rename to data/demos/features/h5_reader.py From abfa897d43b7bcc9094d324f519bb7d550ce8ebd Mon Sep 17 00:00:00 2001 From: Allie Date: Thu, 10 Sep 2026 13:31:19 -0600 Subject: [PATCH 2/9] Create h5_reader_README.md --- data/demos/features/h5_reader_README.md | 254 ++++++++++++++++++++++++ 1 file changed, 254 insertions(+) create mode 100644 data/demos/features/h5_reader_README.md diff --git a/data/demos/features/h5_reader_README.md b/data/demos/features/h5_reader_README.md new file mode 100644 index 0000000..fbd0ce3 --- /dev/null +++ b/data/demos/features/h5_reader_README.md @@ -0,0 +1,254 @@ + +# H5 output reader + +A flexible Python tool for quick visualizing HDF5 SEA-Stack simulation outputs. + +It supports: +- Body motion (position, velocity, acceleration) +- Joint reactions (force & torque) +- Translational spring-damper actuators (TSDAs) outputs (energy, power, reaction forces) + +With capabilities such as +- Multiple devices plotting +- Various naming conventions +- Filtering via CLI +- File introspection (`--list`) + +--------- + +## Basic command structure: + +run h5 output reader python file in terminal +- located in data/demos/run_seastack/features + - either cd data/demos/run_seastack/features or direct to data/demos/run_seastack/features/h5outputsReader.py + +```bash +python h5outputsReader.py --field [options] --bodies [options] --dof [options] --joints [options] --tsda [options] --list --mooring [index] --mooring_path +``` + +## Features + +### `--field` +Select what data to plot. + +**Body fields:** +- position +- velocity +- acceleration + +**Joint fields:** +- joint_force +- joint_torque + +**TSDA fields:** +- energy +- power +- tsda_force + +--- + +### `--bodies` +Select body numbers to plot on same figure (default = all) + +**Used with:** body fields + +**Example:** +```bash +--bodies 1 2 4 +``` + +--- + +### `--dof` +Select degrees of freedom to plot on separate figures (default = heave) + +Options: +- surge +- sway +- heave +- roll +- pitch +- yaw + +**Used with:** body fields (not acceleration) + +**Example:** +```bash +--dof heave +``` + +or for multiple plots: +```bash +--dof surge sway heave +``` + +--- + +### `--joints` +Select joint names to plot + +**Used with:** joint fields + +**Tip:** use `--list` to see available joint names + +**Example:** +```bash +--joints float_plate_joint +``` + +can plot multiple joint parameters on same plot: +```bash +--joints joint_12 joint_34 +``` + +--- + +### `--tsdas` +Select TSDA names to plot + +**Used with:** TSDA fields + +**Example:** +```bash +--tsdas ram_23_stbd ram_45_stbd +``` + +**Special option:** +```bash +--tsdas all +``` + +- Plots the sum of all TSDAs (for energy and power fields specification) +- Will prints: + - total energy + - average power + - peak power + +--- + +### `--list` +List all available joint and TSDA names in the file for force and torque data + +**Example:** +```bash +--list +``` + +--- + +### `--mooring` +Load and visualize mooring output data from .out files +Supports: +- force Fx, Fy, Fz [on body] +- moment Mx, My, Mz [on body] +- tension [on line] + +With the same `--mooring` command and [integer] for respective field + +Mooring directory can be specified with `--mooring_path` + +**Example:** +```bash +--mooring tension 1 3 5 --mooring_path +``` +--mooring [index] --mooring_path + + +--- + +## Full example usage + +### Plot heave position (all bodies) of 5sa demo - bodies not specified + +```bash +python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field position --dof heave +``` +image + +--- + +### Plot sway velocity (of selected bodies) of 5sa model +(note: order bodies are listed can change plotting order!) + +```bash +python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field velocity --dof sway --bodies 3 5 +``` +image + +--- + +### Plot total TSDA power of RM3 + +```bash +python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field power +``` +image + +--- + +### List available joint force outputs for the 5sa +```bash +python h5outputsReader.py 5sa/regular_waves/outputs/results.regular.h5 --list +``` + +image + +--- + +### Plot joint torque for RM3 + +```bash +python h5outputsReader.py ../rm3/regular_waves/outputs/results.regular.h5 --field joint_torque --joints float_plate_joint +``` +image + +--- + +### Plot selected TSDA forces (between body 1 and 2) of 5sa (tsda names from list function output) + +```bash +python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field tsda_force --tsdas ram_12_top ram_12_bottom ram_12_stbd ram_12_port +``` + +image + +--- + +### Plot selected TSDA energy of 5sa (tsda names from list function output) + +```bash +python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field energy --tsdas ram_23_stbd ram_45_stbd +``` +image + +--- + +### Plot total TSDA power of 5sa + +```bash +python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field energy --tsdas all +``` + +image + +image + +--- + +### Plot moment on body 2 of 5sa + +```bash +python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --mooring moment 2 +``` +image + +--- + +### Plot tension on lines of RM3 specifying path + +```bash +python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --mooring tension --mooring_path C:\Users\ariley\source\repos\sea-stack\build\data\demos\run_seastack\rm3\assets\mooring +``` +image + +--- From 0e40891aba6f518f64a0a379bbc69c699f0fcd08 Mon Sep 17 00:00:00 2001 From: Allie Date: Thu, 10 Sep 2026 13:32:52 -0600 Subject: [PATCH 3/9] Update README to reference h5_reader.py --- data/demos/features/h5_reader_README.md | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/data/demos/features/h5_reader_README.md b/data/demos/features/h5_reader_README.md index fbd0ce3..513e571 100644 --- a/data/demos/features/h5_reader_README.md +++ b/data/demos/features/h5_reader_README.md @@ -20,10 +20,10 @@ With capabilities such as run h5 output reader python file in terminal - located in data/demos/run_seastack/features - - either cd data/demos/run_seastack/features or direct to data/demos/run_seastack/features/h5outputsReader.py + - either cd data/demos/run_seastack/features or direct to data/demos/run_seastack/features/h5_reader.py ```bash -python h5outputsReader.py --field [options] --bodies [options] --dof [options] --joints [options] --tsda [options] --list --mooring [index] --mooring_path +python h5_reader.py --field [options] --bodies [options] --dof [options] --joints [options] --tsda [options] --list --mooring [index] --mooring_path ``` ## Features @@ -161,7 +161,7 @@ Mooring directory can be specified with `--mooring_path` ### Plot heave position (all bodies) of 5sa demo - bodies not specified ```bash -python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field position --dof heave +python h5_reader.py ../5sa/regular_waves/outputs/results.regular.h5 --field position --dof heave ``` image @@ -171,7 +171,7 @@ python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --fiel (note: order bodies are listed can change plotting order!) ```bash -python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field velocity --dof sway --bodies 3 5 +python h5_reader.py ../5sa/regular_waves/outputs/results.regular.h5 --field velocity --dof sway --bodies 3 5 ``` image @@ -180,7 +180,7 @@ python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --fiel ### Plot total TSDA power of RM3 ```bash -python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field power +python h5_reader.py ../5sa/regular_waves/outputs/results.regular.h5 --field power ``` image @@ -188,7 +188,7 @@ python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --fie ### List available joint force outputs for the 5sa ```bash -python h5outputsReader.py 5sa/regular_waves/outputs/results.regular.h5 --list +python h5_reader.py 5sa/regular_waves/outputs/results.regular.h5 --list ``` image @@ -198,7 +198,7 @@ python h5outputsReader.py 5sa/regular_waves/outputs/results.regular.h5 --list ### Plot joint torque for RM3 ```bash -python h5outputsReader.py ../rm3/regular_waves/outputs/results.regular.h5 --field joint_torque --joints float_plate_joint +python h5_reader.py ../rm3/regular_waves/outputs/results.regular.h5 --field joint_torque --joints float_plate_joint ``` image @@ -207,7 +207,7 @@ python h5outputsReader.py ../rm3/regular_waves/outputs/results.regular.h5 --fie ### Plot selected TSDA forces (between body 1 and 2) of 5sa (tsda names from list function output) ```bash -python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field tsda_force --tsdas ram_12_top ram_12_bottom ram_12_stbd ram_12_port +python h5_reader.py ../5sa/regular_waves/outputs/results.regular.h5 --field tsda_force --tsdas ram_12_top ram_12_bottom ram_12_stbd ram_12_port ``` image @@ -217,7 +217,7 @@ python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --fiel ### Plot selected TSDA energy of 5sa (tsda names from list function output) ```bash -python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field energy --tsdas ram_23_stbd ram_45_stbd +python h5_reader.py ../5sa/regular_waves/outputs/results.regular.h5 --field energy --tsdas ram_23_stbd ram_45_stbd ``` image @@ -226,7 +226,7 @@ python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --fie ### Plot total TSDA power of 5sa ```bash -python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --field energy --tsdas all +python h5_reader.py ../5sa/regular_waves/outputs/results.regular.h5 --field energy --tsdas all ``` image @@ -238,7 +238,7 @@ python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --fie ### Plot moment on body 2 of 5sa ```bash -python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --mooring moment 2 +python h5_reader.py ../5sa/regular_waves/outputs/results.regular.h5 --mooring moment 2 ``` image @@ -247,7 +247,7 @@ python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --moo ### Plot tension on lines of RM3 specifying path ```bash -python h5outputsReader.py ../5sa/regular_waves/outputs/results.regular.h5 --mooring tension --mooring_path C:\Users\ariley\source\repos\sea-stack\build\data\demos\run_seastack\rm3\assets\mooring +python h5_reader.py ../5sa/regular_waves/outputs/results.regular.h5 --mooring tension --mooring_path C:\Users\ariley\source\repos\sea-stack\build\data\demos\run_seastack\rm3\assets\mooring ``` image From e9dbc91190bc7d1232ae3168b503277a66da4598 Mon Sep 17 00:00:00 2001 From: Allie Date: Thu, 10 Sep 2026 13:38:43 -0600 Subject: [PATCH 4/9] Add files via upload --- data/demos/features/MHKiTWaveYamlWriter.py | 3452 ++++++++++++++++++++ 1 file changed, 3452 insertions(+) create mode 100644 data/demos/features/MHKiTWaveYamlWriter.py diff --git a/data/demos/features/MHKiTWaveYamlWriter.py b/data/demos/features/MHKiTWaveYamlWriter.py new file mode 100644 index 0000000..5ada5da --- /dev/null +++ b/data/demos/features/MHKiTWaveYamlWriter.py @@ -0,0 +1,3452 @@ + +import os +import sys +import argparse +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import yaml +from scipy import integrate +from scipy.optimize import minimize, minimize_scalar +from mhkit import wave +from mhkit.wave.io import ndbc + +# Redirected stdout defaults to cp1252 on Windows, which cannot encode the gamma +# and superscript characters used in the printed diagnostics. +sys.stdout.reconfigure(encoding="utf-8", errors="replace") + + +#~~~~~~~~~~~~~~ Command-line Arguments ~~~~~~~~~~~~~~ + +parser = argparse.ArgumentParser( + description="Update YAML hydro file with NDBC buoy wave data" +) +parser.add_argument( + "path", nargs="?", default=None, + help="Path to the YAML hydro file to update (e.g., 5sa/regular_waves/5sa_regular.hydro.yaml)" +) +parser.add_argument( + "--buoy", + type=str, + default=None, + help="NDBC buoy number (e.g., 46050). Required unless --elevation_file points at an existing eta file." +) +parser.add_argument( + "--type", + type=str, + default="regular", + choices=["regular", "irregular"], + help="Wave type: 'regular' (default) or 'irregular'" +) +parser.add_argument( + "--year", + nargs="*", + default=None, + help="Year(s) to use for wave statistics. Options: no args (latest year), from:YYYY (years >= YYYY), YYYY (single year), YYYY YYYY... (multiple years). Default: 2025 or latest available." +) +parser.add_argument( + "--plot_spectrum", + nargs="*", + default=None, + help="Plot wave spectrum. Optional: specify a year for monthly spectra (e.g., --plot_spectrum 2020), or one or more months to overlay (e.g., --plot_spectrum jan feb mar), or omit for the --year flag" +) +parser.add_argument( + "--plot_heatmap", + nargs="*", + default=None, + help="Plot sea state frequency distribution (Hm0 vs Te) for the selection. Optional: a year (e.g., --plot_heatmap 2020) or one or more months (e.g., --plot_heatmap jan feb), or omit to use --year / --date." +) +parser.add_argument( + "--plot_wavestats", + nargs="*", + default=None, + help="Plot statistics (median + IQR) of Hm0, Te, Tp, energy flux J and steepness Sm. The resolution follows the selection: hourly for a day (--date), daily for up to ~2 months, monthly for longer. Optional: a year (e.g., --plot_wavestats 2020) or one or more months (e.g., --plot_wavestats jan feb), or omit to use --year / --date." +) +parser.add_argument( + "--plot_wavedirection", + nargs="?", + const="elevation", + default=None, + choices=["elevation", "energy", "spread"], + help="Plot directional wave spectrum (elevation/energy/spread). Default: elevation" +) +parser.add_argument( + "--spectrum", + type=str, + default=None, + choices=["pm", "pierson_moskowitz", "jonswap", "js", "custom"], + help="Spectrum type for irregular waves (pm, pierson_moskowitz, jonswap, js, custom). Only used with --type irregular. 'custom' generates a surface elevation time series from the measured buoy spectrum and writes an eta file next to the YAML." +) +parser.add_argument( + "--spectrum_file", + type=str, + default=None, + help="Read the custom spectrum from a text file instead of downloading it from a buoy, so --buoy is not needed. Accepts a bare name, a path relative to the YAML directory, or an absolute path. Two layouts are recognised: two rows (frequencies then spectral densities) or two columns (frequency, density) - blank lines and lines starting with # or %% are ignored. Units are Hz and m^2/Hz. Implies --type irregular --spectrum custom." +) +parser.add_argument( + "--elevation_duration", + type=float, + default=None, + help="Simulation duration in seconds, written to simulation.end_time for any spectrum type. With --spectrum custom it also sets the length of the generated elevation record. Default: 600 (10 minutes) for the eta record; the simulation YAML is left alone unless this flag is given. Ignored when --elevation_file points at an existing record." +) +parser.add_argument( + "--elevation_dt", + type=float, + default=None, + help="Solver time step in seconds, written to simulation.time_step for any spectrum type. With --spectrum custom it also sets the sample rate of the generated elevation record. Default: the existing simulation.time_step (falls back to 0.05 if no simulation YAML is found), so the record is sampled at the solver step." +) +parser.add_argument( + "--elevation_nfreq", + type=int, + default=512, + help="[deprecated] Ignored. The spectrum grid is fixed by the record: df = 1/(duration+dt), up to the Nyquist frequency, because that is the grid surface_elevation's ifft assumes." +) +parser.add_argument( + "--elevation_file", + type=str, + default=None, + help="[--spectrum custom] Eta file: a bare name, a path relative to the YAML directory, or an absolute path. If the file ALREADY EXISTS it is used as-is - no buoy data is downloaded, no record is generated and the file is not rewritten. Otherwise it is the name of the eta file written next to the YAML. Default: _eta.txt" +) +parser.add_argument( + "--eta_method", + type=str, + default="dft", + choices=["dft", "irf"], + help="[--spectrum custom] How sea-stack consumes the eta file. 'dft' (default) extracts discrete wave components, giving a spatial free surface so the GUI wireframe renders. 'irf' uses EtaTableWaveField, a point time series with no spatial variation - the body moves but no wave surface is drawn." +) +parser.add_argument( + "--seed", + type=int, + default=None, + help="Random seed for wave generation (optional)" +) +parser.add_argument( + "--ramp_time", + type=float, + default=None, + help="Ramp time in seconds for wave generation (optional)" +) +parser.add_argument( + "--height", + type=float, + default=None, + help="Override calculated wave height in meters (optional)" +) +parser.add_argument( + "--period", + type=float, + default=None, + help="Override calculated wave period in seconds (optional)" +) +parser.add_argument( + "--depth", + type=float, + default=None, + help="Override calculated water depth in meters (optional)" +) +parser.add_argument( + "--nofill", + action="store_true", + help="Remove fill from directional spectrum plot (show contours only)" +) +parser.add_argument( + "--plot_spectrum_pm", "--plot_spectrum_PM", + dest="plot_spectrum_pm", + nargs="*", + default=None, + help="Plot wave spectrum against Pierson-Moskowitz estimate. Optional: a year (e.g., --plot_spectrum_pm 2020) or one or more months to overlay (e.g., --plot_spectrum_pm jan feb), or omit to use --year flag" +) +parser.add_argument( + "--plot_spectrum_js", "--plot_spectrum_JS", + "--plot_spectrum_jonswap", "--plot_spectrum_JONSWAP", + dest="plot_spectrum_js", + nargs="*", + default=None, + help="Plot wave spectrum against fitted JONSWAP spectrum. Optional: a year (e.g., --plot_spectrum_js 2020) or one or more months to overlay (e.g., --plot_spectrum_js jan feb), or omit to use --year flag" +) +parser.add_argument( + "--gamma", + nargs="*", + type=float, + action='append', + help="Manual gamma (peak enhancement factor) values to plot alongside fitted JONSWAP. Can be used as: --gamma 1.0 1.5 or --gamma 1.0 --gamma 1.5 (e.g., --gamma 0.5 1.0 2.5)" +) +parser.add_argument( + "--plot_wind", + nargs="*", + default=None, + help="Plot wind speed, gusts and direction. With --date: daily 10-minute interval time series. Otherwise median + IQR at a resolution that follows the selection: daily for up to ~2 months, monthly for longer. Optional: a year (e.g., --plot_wind 2020) or one or more months (e.g., --plot_wind jan feb). Does not update the YAML file." +) +parser.add_argument( + "--date", + nargs="*", + default=None, + help="Extract data for one or more dates. Format: YYYY-MM-DD (e.g., 2024-06-15), MM-DD-YYYY (e.g., 06-15-2024), or MM-DD with --year (e.g., --date 06-15 --year 2024). Multiple dates are overlaid on one plot (e.g., --date 01-01 06-15). Without --time, each date is a daily average. With --time, each date/time combination is a single measurement." +) +parser.add_argument( + "--time", + nargs="*", + default=None, + help="[Optional] One or more times for single measurement extraction. Format: HH:MM in 24-hour format (e.g., 18:30), or a range HH:MM-HH:MM (e.g., 05:20-16:40) which averages every measurement in that window. 24:00 means end of day. Multiple times or ranges are overlaid (e.g., --time 00:00-12:00 12:00-24:00). Each date is combined with each time. NDBC wave data is sampled at :10 and :40 past each hour. If omitted with --date, daily averages are used." +) + +args = parser.parse_args() + +user_set_spectrum = args.spectrum is not None +user_set_type = args.type != "regular" + +if args.spectrum_file: + if args.spectrum is None: + args.spectrum = "custom" + elif args.spectrum != "custom": + parser.error("--spectrum_file requires --spectrum custom") + if args.type != "irregular": + args.type = "irregular" + if args.date or args.time: + parser.error("--date / --time select a buoy measurement and cannot be used with --spectrum_file") + if args.plot_wind is not None: + parser.error("wind plots need buoy data and cannot be used with --spectrum_file") + +if args.buoy is None and not (args.spectrum == "custom" and (args.elevation_file or args.spectrum_file)): + parser.error("--buoy is required, unless --spectrum custom is reading an existing " + "--elevation_file or a --spectrum_file") + +# Flatten gamma list if it was provided (handles both --gamma 1.0 1.5 and --gamma 1.0 --gamma 1.5) +if args.gamma: + args.gamma = [g for sublist in args.gamma for g in sublist] if any(isinstance(i, list) for i in args.gamma) else args.gamma +else: + args.gamma = [] + +# --date alone gives a daily average per date; --date + --time gives a single raw +# measurement per date/time pair. Multiple values form a cross product, overlaid. +single_measurement_date = None +daily_start = None +daily_end = None +extract_single = False +extract_daily = False +all_wind_data = None +extract_targets = [] # list of dicts: {'kind': 'single'|'daily', 'timestamp': ..., 'label': ...} + +if args.date is not None and len(args.date) > 0: + def _parse_one_date(date_str): + """Parse YYYY-MM-DD, MM-DD-YYYY or MM-DD (with --year) into (year, month, day).""" + parts = date_str.split('-') + y = m = d = None + + if len(parts) == 3: + try: + a, b, c = (int(p) for p in parts) + except ValueError: + parser.error(f"Invalid date format: {date_str}. Use YYYY-MM-DD, MM-DD-YYYY or MM-DD") + # The 4-digit field tells the two orderings apart + if len(parts[0]) == 4: + y, m, d = a, b, c + elif len(parts[2]) == 4: + y, m, d = c, a, b + else: + parser.error(f"Invalid date format: {date_str}. The year must be 4 digits, " + f"as YYYY-MM-DD or MM-DD-YYYY") + elif len(parts) == 2: + if args.year is None or len(args.year) == 0: + parser.error("--date in MM-DD format requires --year to be specified") + try: + m, d = map(int, parts) + y = int(args.year[0]) if isinstance(args.year[0], str) else args.year[0] + except (ValueError, IndexError): + parser.error(f"Invalid date or year format: {date_str} with year {args.year}") + else: + parser.error(f"Invalid date format: {date_str}. Use YYYY-MM-DD, MM-DD-YYYY or MM-DD") + + try: + pd.Timestamp(year=y, month=m, day=d) + except ValueError as e: + parser.error(f"Invalid date: {date_str} ({e})") + return y, m, d + + def _parse_one_time(time_str): + """Parse HH:MM into (hour, minute); 24:00 is accepted as end of day.""" + try: + h, mi = map(int, time_str.split(':')) + except ValueError: + parser.error(f"Invalid time format: {time_str}. Use HH:MM (24-hour format)") + if (h, mi) == (24, 0): + return 23, 59 + if not (0 <= h < 24 and 0 <= mi < 60): + parser.error(f"Invalid time format: {time_str}. Use HH:MM (24-hour format), " + f"or 24:00 for end of day") + return h, mi + + def _parse_time_token(token): + """Parse HH:MM or HH:MM-HH:MM into ('point', (h, mi)) or ('range', start, end).""" + if '-' in token: + lo, hi = token.split('-', 1) + start, end = _parse_one_time(lo), _parse_one_time(hi) + if end <= start: + parser.error(f"Invalid time range: {token}. The end time must be after the start") + return ('range', start, end) + return ('point', _parse_one_time(token)) + + parsed_dates = [_parse_one_date(d) for d in args.date] + + if args.time is not None and len(args.time) > 0: + # Cross product of dates x times; a range collapses to one averaged spectrum + parsed_times = [_parse_time_token(t) for t in args.time] + for (y, m, d) in parsed_dates: + for spec in parsed_times: + if spec[0] == 'range': + (h1, mi1), (h2, mi2) = spec[1], spec[2] + start = pd.Timestamp(year=y, month=m, day=d, hour=h1, minute=mi1) + end = pd.Timestamp(year=y, month=m, day=d, hour=h2, minute=mi2) + extract_targets.append({ + 'kind': 'range', + 'timestamp': start, + 'end': end, + 'label': f"{start:%Y-%m-%d %H:%M}-{end:%H:%M}", + }) + else: + h, mi = spec[1] + ts = pd.Timestamp(year=y, month=m, day=d, hour=h, minute=mi) + extract_targets.append({ + 'kind': 'single', + 'timestamp': ts, + 'label': ts.strftime('%Y-%m-%d %H:%M'), + }) + extract_single = True + single_measurement_date = extract_targets[0]['timestamp'] + mode_word = "Time range" if all(t['kind'] == 'range' for t in extract_targets) else "Single measurement" + if len(extract_targets) == 1: + print(f"{mode_word} mode: {extract_targets[0]['label']}") + else: + print(f"{mode_word} mode: {len(extract_targets)} selections") + for t in extract_targets: + print(f" - {t['label']}") + print(f" (All plotting and writing operations will use these measurements)\n") + else: + # Daily average mode: one target per date + for (y, m, d) in parsed_dates: + ts = pd.Timestamp(year=y, month=m, day=d, hour=0, minute=0) + extract_targets.append({ + 'kind': 'daily', + 'timestamp': ts, + 'label': ts.strftime('%Y-%m-%d'), + }) + extract_daily = True + daily_start = extract_targets[0]['timestamp'] + daily_end = daily_start.replace(hour=23, minute=59) + if len(extract_targets) == 1: + print(f"Daily average mode: {extract_targets[0]['label']} (00:00 to 23:59)") + else: + print(f"Daily average mode: {len(extract_targets)} days") + for t in extract_targets: + print(f" - {t['label']}") + print(f" (All plotting and writing operations will use the daily average(s))\n") + + + +# nargs="*" hands back a list; downstream code expects None/True/str, so keep both +def normalize_plot_arg(value): + """Returns (normalized_arg, values_list). normalized_arg is None/True/str.""" + if value is None: + return None, [] + if isinstance(value, list): + if len(value) == 0: + return True, [] + if len(value) == 1: + return value[0], [value[0]] + return value[0], list(value) + return value, [value] if isinstance(value, str) else [] + +args.plot_spectrum, plot_spectrum_values = normalize_plot_arg(args.plot_spectrum) +args.plot_spectrum_pm, plot_spectrum_pm_values = normalize_plot_arg(args.plot_spectrum_pm) +args.plot_spectrum_js, plot_spectrum_js_values = normalize_plot_arg(args.plot_spectrum_js) +args.plot_wavestats, plot_wavestats_values = normalize_plot_arg(args.plot_wavestats) +args.plot_heatmap, plot_heatmap_values = normalize_plot_arg(args.plot_heatmap) +args.plot_wind, plot_wind_values = normalize_plot_arg(args.plot_wind) + + +def parse_month_list(values, min_count=2): + """Parse a list of month tokens. Returns list of month numbers, or None if not all months. + + min_count guards how many values are required before the list branch is used, so + flags that already have dedicated single-month handling keep it. + """ + if not values or len(values) < min_count: + return None + months = [] + for v in values: + num, is_month = parse_month(v) + if not is_month: + return None + months.append(num) + return months + + +# Plot-only flags do not by themselves imply intent to rewrite the YAML +plot_only_options = ( + args.plot_spectrum is not None or + args.plot_heatmap is not None or + args.plot_wavestats is not None or + args.plot_wavedirection is not None or + args.plot_spectrum_pm is not None or + args.plot_spectrum_js is not None or + args.plot_wind is not None +) + +# Check if user explicitly requested a YAML update (not just plot) +# These are parameters that indicate intent to update YAML +explicit_update_intent = ( + user_set_spectrum or + args.depth is not None or + user_set_type +) + +# If plot-only and no explicit update intent, skip YAML update (height/period used only for plotting) +# If no plotting options, then height/period alone means update YAML +skip_yaml_update = plot_only_options and not explicit_update_intent +should_update_yaml = not skip_yaml_update + +# Validate that path was provided, but only if updating YAML (not for plot-only) +if args.path is None and not skip_yaml_update: + parser.error("Please provide a path to the YAML hydro file") + +yaml_file_path = args.path if args.path is not None else "" + +# Fail fast on a bad YAML path rather than after downloading buoy data and +# generating a full elevation record +if yaml_file_path and not skip_yaml_update: + if not os.path.exists(yaml_file_path): + target_dir = os.path.dirname(os.path.abspath(yaml_file_path)) + print(f"ERROR: YAML file not found: {yaml_file_path}") + if not os.path.isdir(target_dir): + print(f" The directory does not exist either: {target_dir}") + print(f" Create the directory and place a .hydro.yaml file in it first.") + else: + existing = sorted(f for f in os.listdir(target_dir) if f.endswith('.yaml')) + if existing: + print(f" YAML files present in {target_dir}:") + for name in existing: + print(f" - {name}") + else: + print(f" No .yaml files found in {target_dir}") + print(" This script updates an existing hydro YAML; it does not create one.") + exit(1) + +#~~~~~~~~~~~~~~ Existing eta file -> YAML (no buoy data needed) ~~~~~~~~~~~~~~~~ + +# An --elevation_file that already exists is used exactly as it is: nothing is +# downloaded, nothing is generated and the file is never rewritten. +existing_eta_path = None +if args.spectrum == "custom" and args.elevation_file and yaml_file_path: + _yaml_dir = os.path.dirname(os.path.abspath(yaml_file_path)) + _candidate = args.elevation_file if os.path.isabs(args.elevation_file) else \ + os.path.normpath(os.path.join(_yaml_dir, args.elevation_file)) + if not os.path.isfile(_candidate) and not os.path.splitext(_candidate)[1]: + _candidate += ".txt" + if os.path.isfile(_candidate): + existing_eta_path = _candidate + +if existing_eta_path: + print("=" * 60) + print("Existing elevation file -> YAML (buoy data not used)") + print("=" * 60) + print(f"eta file: {existing_eta_path}") + + if args.elevation_duration is not None or args.elevation_dt is not None: + print(" NOTE: --elevation_duration / --elevation_dt are ignored; the record is used as-is") + + times, etas = [], [] + with open(existing_eta_path, 'r') as f: + for line in f: + line = line.strip() + if not line or line[0] in '#%': + continue + parts = line.replace(':', ' ').replace(',', ' ').split() + try: + times.append(float(parts[0])) + etas.append(float(parts[1])) + except (IndexError, ValueError): + if not times: + continue # header row + print(f"ERROR: cannot parse '{line}'") + print(" Expected 'time:elevation' (':', ',' or whitespace separated)") + exit(1) + + if len(times) < 2: + print(f"ERROR: {existing_eta_path} has fewer than 2 samples") + exit(1) + + eta_values = np.asarray(etas) + n_samples = len(times) + duration = float(times[-1] - times[0]) + dt = duration / (n_samples - 1) + print(f" {n_samples} samples, {duration:.2f} s at dt={dt:g} s") + print(f" Elevation range {eta_values.min():.3f} to {eta_values.max():.3f} m, " + f"Hm0 (4*sigma) = {4.0 * np.std(eta_values):.4f} m") + + # nfrequencies must span the default 0.001 - 1.0 Hz eta band: the C++ default of + # 1000 strides the Fourier grid, silently dropping bins and attenuating the record. + df_natural = 1.0 / (n_samples * dt) + k_min = max(1, int(np.ceil(0.001 / df_natural))) + k_max = int(np.floor(min(1.0, 1.0 / (2.0 * dt)) / df_natural)) + eta_nfrequencies = max(1, k_max - k_min + 1) + print(f" Fourier grid: df = {df_natural:.3e} Hz -> {eta_nfrequencies} bins " + f"across the default 0.001 - 1.0 Hz eta band") + if args.eta_method == "dft" and eta_nfrequencies > 4000: + print(f" WARNING: {eta_nfrequencies} wave components is a lot. The DFT is O(nf*N) and") + print(f" the GUI evaluates every component per surface vertex per frame.") + + eta_ref = os.path.relpath(existing_eta_path, _yaml_dir).replace('\\', '/') + + waves_block = {'type': 'irregular', 'eta_file': eta_ref} + if args.eta_method == "irf": + waves_block['method'] = "irf_convolution" + else: + waves_block['nfrequencies'] = int(eta_nfrequencies) + + ramp = float(args.ramp_time) if args.ramp_time is not None else 60.0 + if ramp > 0.0: + waves_block['ramp_duration'] = ramp + waves_block['ramp_type'] = "cosine" + if args.depth is not None: + waves_block['depth'] = float(args.depth) + + with open(yaml_file_path, 'r') as f: + yaml_data = yaml.safe_load(f) or {} + if 'waves' not in (yaml_data.get('hydrodynamics') or {}): + print(f"ERROR: 'hydrodynamics.waves' section not found in {yaml_file_path}") + exit(1) + yaml_data['hydrodynamics']['waves'] = waves_block + + class IndentDumper(yaml.SafeDumper): + def increase_indent(self, flow=False, indentless=False): + return super().increase_indent(flow, False) + + with open(yaml_file_path, 'w') as f: + yaml.dump(yaml_data, f, Dumper=IndentDumper, default_flow_style=False, sort_keys=False) + + print(f"\nUpdated hydro YAML: {yaml_file_path}") + for key, value in waves_block.items(): + print(f" {key}: {value}") + + # The wave field returns zero elevation outside the table range, so the + # companion simulation end_time must not exceed the record. + base_name = os.path.basename(yaml_file_path) + sim_stem = base_name[:-len('.hydro.yaml')] if base_name.endswith('.hydro.yaml') \ + else os.path.splitext(base_name)[0] + sim_yaml_path = os.path.join(_yaml_dir, f"{sim_stem}.simulation.yaml") + if not os.path.exists(sim_yaml_path): + import glob + matches = sorted(glob.glob(os.path.join(_yaml_dir, "*.simulation.yaml"))) + sim_yaml_path = matches[0] if len(matches) == 1 else None + + if sim_yaml_path is None: + print(f"\nNOTE: no companion *.simulation.yaml found; end_time not updated.") + else: + with open(sim_yaml_path, 'r') as f: + sim_lines = f.read().splitlines() + + sim_step = None + in_simulation = False + updated = False + old_end_time = None + for i, line in enumerate(sim_lines): + stripped = line.strip() + if not stripped or stripped.startswith('#'): + continue + indent = len(line) - len(line.lstrip()) + if indent == 0: + in_simulation = stripped.rstrip(':') == 'simulation' + continue + if not in_simulation: + continue + key = stripped.split(':')[0].strip() + if key == 'time_step': + try: + sim_step = float(stripped.split(':', 1)[1].split('#')[0].strip()) + except ValueError: + pass + elif key == 'end_time': + old_end_time = stripped.split(':', 1)[1].strip() + sim_lines[i] = f"{' ' * indent}end_time: {duration:.1f}" + updated = True + + if updated: + with open(sim_yaml_path, 'w', newline='\n') as f: + f.write("\n".join(sim_lines) + "\n") + print(f"\nUpdated simulation YAML: {sim_yaml_path}") + print(f" end_time: {old_end_time} -> {duration:.1f} s (matches eta record length)") + else: + print(f"\nWARNING: no 'end_time' key found under 'simulation:' in {sim_yaml_path}") + + if sim_step and dt > sim_step: + print(f" WARNING: eta dt ({dt} s) is coarser than simulation.time_step ({sim_step} s)") + + exit(0) + +#~~~~~~~~~~~from buoy number~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +# # use if SSLcertificateError / VPN blocks from downloading data from NDBC website +os.environ["REQUESTS_CA_BUNDLE"] = r"C:\Users\ariley\OneDrive - NREL\Documents\nrel_root_ca.cer" +os.environ["SSL_CERT_FILE"] = r"C:\Users\ariley\OneDrive - NREL\Documents\nrel_root_ca.cer" + +#~~~~~~~~~~~~~~~INPUTS~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +buoy_number = args.buoy +yaml_file_path = args.path +wave_type = args.type + +#~~~~~~~~~~~~~~ Custom spectrum from a file (no buoy data) ~~~~~~~~~~~~~~~~~~~~~ + +def load_spectrum_file(path): + """Read a two-row or two-column frequency/spectral-density file.""" + rows = [] + with open(path, 'r') as f: + for line in f: + line = line.split('#')[0].split('%')[0].strip() + if not line: + continue + try: + rows.append([float(tok) for tok in line.replace(',', ' ').split()]) + except ValueError: + if not rows: + continue # header row + raise ValueError(f"cannot parse numbers from line: {line!r}") + + if len(rows) == 2 and len(rows[0]) == len(rows[1]) and len(rows[0]) > 2: + freq, dens = np.asarray(rows[0]), np.asarray(rows[1]) + elif len(rows) > 2 and all(len(r) == 2 for r in rows): + arr = np.asarray(rows) + freq, dens = arr[:, 0], arr[:, 1] + else: + raise ValueError( + "unrecognised layout. Expected either two rows (frequencies, then " + "spectral densities) or two columns (frequency, density)") + + if np.any(np.diff(freq) <= 0): + order = np.argsort(freq) + freq, dens = freq[order], dens[order] + if np.any(freq < 0) or np.any(dens < 0): + raise ValueError("frequencies and spectral densities must be non-negative") + return freq, dens + + +spectrum_from_file = bool(args.spectrum_file) +custom_freq = custom_S = None +spectrum_label = None +# Titles say "NDBC Buoy " for downloaded data and the file name otherwise +source_title = f"NDBC Buoy {buoy_number}" + +if spectrum_from_file: + _yaml_dir_for_spec = os.path.dirname(os.path.abspath(yaml_file_path)) if yaml_file_path else os.getcwd() + spectrum_file_path = args.spectrum_file if os.path.isabs(args.spectrum_file) else \ + os.path.normpath(os.path.join(_yaml_dir_for_spec, args.spectrum_file)) + if not os.path.isfile(spectrum_file_path): + # Also accept a path relative to the working directory + alt = os.path.abspath(args.spectrum_file) + if os.path.isfile(alt): + spectrum_file_path = alt + if not os.path.isfile(spectrum_file_path): + print(f"ERROR: spectrum file not found: {args.spectrum_file}") + print(f" Looked in {_yaml_dir_for_spec} and {os.getcwd()}") + exit(1) + + try: + custom_freq, custom_S = load_spectrum_file(spectrum_file_path) + except (ValueError, OSError) as e: + print(f"ERROR: could not read spectrum file {spectrum_file_path}: {e}") + exit(1) + + spectrum_label = os.path.basename(spectrum_file_path) + source_title = spectrum_label + _m0 = float(integrate.trapezoid(custom_S, custom_freq)) + print("=" * 60) + print("Custom spectrum from file (buoy data not used)") + print("=" * 60) + print(f"File: {spectrum_file_path}") + print(f" {len(custom_freq)} bins, {custom_freq.min():.4f} - {custom_freq.max():.4f} Hz") + print(f" Hm0 = {4.0 * np.sqrt(_m0):.4f} m") + +# ========== Define spectrum helper functions EARLY (before YAML update) ========== + +def plot_heading(label, model=""): + """Plot title head, without repeating the source when it already is the label.""" + head = label if label == source_title else f"{source_title} - {label}" + return f"{head} - {model}" if model else head + + +# JONSWAP is only defined for gamma >= 1 by the standard definition; 0.9 is the +# practical lower limit allowed here (empirical best fits sit just below 1.0). +GAMMA_MIN = 0.9 +GAMMA_MAX = 10.0 + +def parse_month(month_input): + """Parse month name or number. Returns (month_number, is_month) tuple.""" + month_map = { + 'jan': 1, 'january': 1, + 'feb': 2, 'february': 2, + 'mar': 3, 'march': 3, + 'apr': 4, 'april': 4, + 'may': 5, + 'jun': 6, 'june': 6, + 'jul': 7, 'july': 7, + 'aug': 8, 'august': 8, + 'sep': 9, 'september': 9, + 'oct': 10, 'october': 10, + 'nov': 11, 'november': 11, + 'dec': 12, 'december': 12 + } + + if isinstance(month_input, str): + lower_input = month_input.lower() + if lower_input in month_map: + return month_map[lower_input], True + try: + month_num = int(month_input) + if 1 <= month_num <= 12: + return month_num, True + except ValueError: + pass + return None, False + + +def pm_spectrum(freq, hm0, tp): + """ + Unified Pierson-Moskowitz spectrum function. + + Parameters: + - freq: frequency array [Hz] + - hm0: significant wave height [m] + - tp: peak period [s] + + Returns: spectrum S(f) normalized so that 4*sqrt(integral(S)) = hm0 + """ + freq = np.asarray(freq) + fp = 1.0 / tp # peak frequency + + # Avoid division by zero + safe_freq = np.where(freq > 1e-10, freq, 1e-10) + + # Unnormalized PM shape: f^(-5) * exp(-5/4 * (fp/f)^4) + exponent = -5.0/4.0 * (fp / safe_freq)**4 + S_unnormalized = (safe_freq**(-5)) * np.exp(exponent) + + # Integrate to get unnormalized zeroth moment + m0_unnormalized = integrate.trapezoid(S_unnormalized, freq) + + # Scale so that 4*sqrt(m0) = hm0 (significant wave height definition) + target_m0 = (hm0 / 4.0) ** 2.0 + if m0_unnormalized > 0: + S = S_unnormalized * (target_m0 / m0_unnormalized) + else: + S = S_unnormalized + + return S + + +def jonswap_spectrum(freq, hm0, tp, gamma=3.3): + """ + Unified JONSWAP spectrum function. + + JONSWAP = C(gamma) * S_PM(f) * gamma^[exp(-(f-fp)²/(2σ²fp²))] + + where C(gamma) = integral(S_PM) / integral(S_PM * gamma^peak_factor) + + Parameters: + - freq: frequency array [Hz] + - hm0: significant wave height [m] + - tp: peak period [s] + - gamma: peak enhancement factor (1-5, typical 3.3) + + Returns: spectrum S(f) normalized so that 4*sqrt(integral(S)) = hm0 + """ + freq = np.asarray(freq) + fp = 1.0 / tp # peak frequency + + # Get base PM spectrum + S_pm = pm_spectrum(freq, hm0, tp) + + # Peak enhancement factor: gamma^[-(f - fp)^2 / (2*sigma^2 * fp^2)] + sigma = np.where(freq <= fp, 0.07, 0.09) + exponent = -((freq - fp)**2.0) / (2.0 * (sigma**2.0) * (fp**2.0)) + exponent = np.clip(exponent, -100, 100) + + # Compute gamma^exponent + if gamma > 0: + log_gamma = np.log(gamma) + peak_factor = np.exp(np.clip(exponent * log_gamma, -100, 100)) + else: + peak_factor = np.ones_like(freq) + + # Compute C(gamma) normalization constant to maintain m0 when peak is enhanced + safe_freq = np.where(freq > 1e-10, freq, 1e-10) + S_pm_unnormalized = (safe_freq**(-5.0)) * np.exp(-5.0/4.0 * (fp / safe_freq)**4.0) + + integral_pm_unnormalized = integrate.trapezoid(S_pm_unnormalized, freq) + integral_pm_enhanced = integrate.trapezoid(S_pm_unnormalized * peak_factor, freq) + + if integral_pm_enhanced > 0 and np.isfinite(integral_pm_enhanced): + c_gamma = integral_pm_unnormalized / integral_pm_enhanced + else: + c_gamma = 1.0 + + # Apply peak enhancement to normalized PM spectrum + # Then renormalize to target Hm0 + S_jonswap_before_norm = c_gamma * S_pm * peak_factor + + # Renormalize to target Hm0 + m0_jonswap = integrate.trapezoid(S_jonswap_before_norm, freq) + target_m0 = (hm0 / 4.0) ** 2.0 + + if m0_jonswap > 0 and np.isfinite(m0_jonswap): + S = S_jonswap_before_norm * (target_m0 / m0_jonswap) + else: + S = S_jonswap_before_norm + + # Ensure no NaNs or infs + S = np.nan_to_num(S, nan=0.0, posinf=0.0, neginf=0.0) + + return S + + +def jonswap_spectrum_custom(freq, Tp, Hs, gamma=3.3, alpha=0.0081): + """ + Legacy wrapper - calls unified jonswap_spectrum for backwards compatibility. + """ + return jonswap_spectrum(freq, Hs, Tp, gamma=gamma) + + +def pierson_moskowitz_spectrum(freq, hm0, tp): + """ + Legacy wrapper - calls unified pm_spectrum for backwards compatibility. + """ + return pm_spectrum(freq, hm0, tp) + + +def fit_jonswap_gamma(measured_spectrum, freq, Tp, Hs, gamma_initial=3.3, fit_tp=True): + """ + Fit a JONSWAP spectrum to a measured spectrum by least squares. + + Fits gamma (peak enhancement), and by default also refines Tp (peak period). + Refining Tp matters because the measured Tp comes from the single highest + spectral bin, which is noisy and can land on the wrong lobe of a broad or + double-peaked sea. Hs is held fixed so the fitted spectrum keeps the same + total energy (m0) as the measurement. + + gamma is constrained to [GAMMA_MIN, GAMMA_MAX]; JONSWAP is undefined below 1 + by the standard definition. + + Returns: (gamma_fit, error, tp_fit, hit_bound) + """ + freq_array = np.asarray(freq) + Tp = float(Tp) + Hs = float(Hs) + + if hasattr(measured_spectrum, 'values'): + measured_values = measured_spectrum.values + else: + measured_values = np.asarray(measured_spectrum) + + def sse(gamma, tp): + """Sum of squared error between measured spectrum and JONSWAP.""" + if not np.isfinite(gamma) or not np.isfinite(tp) or tp <= 0: + return 1e12 + js_spectrum = jonswap_spectrum(freq_array, Hs, tp, gamma=float(gamma)) + return float(np.sum((measured_values - js_spectrum) ** 2)) + + if fit_tp: + # Allow the peak period to move +/-30% around the measured value + tp_lo, tp_hi = 0.7 * Tp, 1.3 * Tp + + def objective(params): + gamma, tp = float(params[0]), float(params[1]) + if not (GAMMA_MIN <= gamma <= GAMMA_MAX) or not (tp_lo <= tp <= tp_hi): + return 1e12 + return sse(gamma, tp) + + best = None + # Multi-start so the optimizer does not settle in a local minimum + for g0 in (1.0, 2.0, 3.3, 5.0): + for t0 in (Tp, 0.9 * Tp, 1.1 * Tp): + res = minimize(objective, x0=[g0, t0], method='Nelder-Mead', + options={'maxiter': 800, 'xatol': 1e-6, 'fatol': 1e-9}) + if best is None or res.fun < best.fun: + best = res + + gamma_fit = float(np.clip(best.x[0], GAMMA_MIN, GAMMA_MAX)) + tp_fit = float(np.clip(best.x[1], tp_lo, tp_hi)) + final_error = sse(gamma_fit, tp_fit) + else: + # Gamma-only fit on a bounded 1-D problem + res = minimize_scalar(lambda g: sse(g, Tp), bounds=(GAMMA_MIN, GAMMA_MAX), + method='bounded', options={'xatol': 1e-6}) + gamma_fit = float(np.clip(res.x, GAMMA_MIN, GAMMA_MAX)) + tp_fit = Tp + final_error = sse(gamma_fit, tp_fit) + + # Flag when the fit is pinned to a gamma bound (means JONSWAP cannot match the shape) + hit_bound = (abs(gamma_fit - GAMMA_MIN) < 1e-3) or (abs(gamma_fit - GAMMA_MAX) < 1e-3) + + return gamma_fit, final_error, tp_fit, hit_bound + +# ========== End spectrum helper functions ========== + +def find_simulation_yaml(hydro_path): + """Locate the companion *.simulation.yaml next to a hydro YAML.""" + if not hydro_path: + return None + directory = os.path.dirname(os.path.abspath(hydro_path)) + base = os.path.basename(hydro_path) + stem = base[:-len('.hydro.yaml')] if base.endswith('.hydro.yaml') \ + else os.path.splitext(base)[0] + + candidate = os.path.join(directory, f"{stem}.simulation.yaml") + if os.path.exists(candidate): + return candidate + + import glob + matches = sorted(glob.glob(os.path.join(directory, "*.simulation.yaml"))) + if len(matches) == 1: + return matches[0] + if len(matches) > 1: + print(f"\nWARNING: multiple *.simulation.yaml files found in {directory}; " + f"expected {stem}.simulation.yaml") + return None + + +def update_simulation_yaml(sim_path, end_time=None, time_step=None, reason=""): + """Rewrite end_time / time_step in place so comments and 1e-4 notation survive.""" + wanted = {} + if end_time is not None: + wanted['end_time'] = f"{float(end_time):.1f}" + if time_step is not None: + wanted['time_step'] = f"{float(time_step):g}" + if not wanted: + return + + with open(sim_path, 'r') as f: + lines = f.read().splitlines() + + in_simulation = False + changed = {} + for i, line in enumerate(lines): + stripped = line.strip() + if not stripped or stripped.startswith('#'): + continue + indent = len(line) - len(line.lstrip()) + if indent == 0: + in_simulation = stripped.rstrip(':') == 'simulation' + continue + if not in_simulation: + continue + key = stripped.split(':')[0].strip() + if key in wanted: + old = stripped.split(':', 1)[1].split('#')[0].strip() + lines[i] = f"{' ' * indent}{key}: {wanted[key]}" + changed[key] = (old, wanted[key]) + + if changed: + with open(sim_path, 'w', newline='\n') as f: + f.write("\n".join(lines) + "\n") + print(f"\nUpdated simulation YAML: {sim_path}") + for key, (old, new) in changed.items(): + print(f" {key}: {old} -> {new}{reason if key == 'end_time' else ''}") + for key in wanted: + if key not in changed: + print(f"\nWARNING: no '{key}' key found under 'simulation:' in {sim_path}") + + +water_depth = None +if spectrum_from_file: + water_depth = float(args.depth) if args.depth is not None else None +else: + try: + buoy_metadata = ndbc.get_buoy_metadata(buoy_number) + except Exception as e: + print(f"ERROR: Could not fetch buoy metadata for {buoy_number}: {e}") + exit(1) + + # Extract water depth from metadata + if 'Water depth' in buoy_metadata: + depth_str = buoy_metadata['Water depth'] + # Parse depth from string like "149 m" + try: + water_depth = float(depth_str.split()[0]) + except (ValueError, IndexError): + print(f"Warning: Could not parse water depth from '{depth_str}'") + +if water_depth: + print(f"\nExtracted Water Depth: {water_depth} m") + +#~~~~~~~~~~~~~~ Parse year specification ~~~~~~~~~~~~~~ + +def parse_year_spec(year_args, available_years): + """ + Parse --year arguments to determine which years to use + + Args: + year_args: List of arguments from --year (or None if not provided) + available_years: List of available years + + Returns: + Tuple of (years_list, is_default) + """ + if year_args is None or len(year_args) == 0: + # Default: use latest year (2025 if available, otherwise most recent) + if 2025 in available_years: + return [2025], True + else: + return [available_years[-1]], True + + years_selected = [] + + for arg in year_args: + if arg.startswith('from:'): + # Range: from:2020 means 2020 and onwards + threshold_year = int(arg[5:]) + years_selected.extend([y for y in available_years if y >= threshold_year]) + else: + # Single year + year = int(arg) + if year in available_years: + years_selected.append(year) + + return sorted(list(set(years_selected))), False # Remove duplicates and sort + + +parameter = "swden" +ndbc_data = {} + +if spectrum_from_file: + available_years = [] + selected_years = [] + is_default = False +else: + ndbc_available_data = ndbc.available_data(parameter, buoy_number) + available_years = sorted([int(y) for y in ndbc_available_data["year"].values]) + print(f"Available years: {available_years}") + + if not available_years: + # Buoys operated by partners (e.g. CDIP/Scripps) can have a live realtime feed on + # ndbc.noaa.gov while having nothing in https://www.ndbc.noaa.gov/data/historical/swden/ + print(f"ERROR: NDBC has no historical spectral wave density (swden) archive for buoy {buoy_number}.") + print(f" The buoy may still show data on the website - that is the realtime feed, which") + print(f" this script cannot use for wave statistics. Check for swden files at:") + print(f" https://www.ndbc.noaa.gov/station_history.php?station={buoy_number}") + print(f" Pick a buoy with a swden archive (e.g. 46050, 46042, 46086).") + exit(1) + + # Parse year specification + selected_years, is_default = parse_year_spec(args.year, available_years) + + # A YYYY-MM-DD passed to --date already pins the year, so don't fall back to "latest". + if extract_targets and (args.year is None or len(args.year) == 0): + date_years = sorted({t['timestamp'].year for t in extract_targets}) + unavailable = [y for y in date_years if y not in available_years] + if unavailable: + print(f"WARNING: year(s) {', '.join(map(str, unavailable))} from --date are not in the " + f"swden archive for buoy {buoy_number}.") + selected_years = [y for y in date_years if y in available_years] + is_default = False + + # Only warn about a defaulted year when no plot flag already pins one + plot_spec_has_year = args.plot_spectrum is not None and isinstance(args.plot_spectrum, str) + plot_spec_pm_has_year = args.plot_spectrum_pm is not None and isinstance(args.plot_spectrum_pm, str) + plot_spec_js_has_year = args.plot_spectrum_js is not None and isinstance(args.plot_spectrum_js, str) + + if is_default and not plot_spec_has_year and not plot_spec_pm_has_year and not plot_spec_js_has_year: + print(f"WARNING: No year specified, using latest available: {selected_years}") + else: + if not is_default: + print(f"Using years: {selected_years}") + + filenames = ndbc_available_data["filename"] + try: + raw_data = ndbc.request_data(parameter, filenames) + except Exception as e: + print(f"ERROR: Could not fetch NDBC data: {e}") + exit(1) + for year in raw_data: + year_data = raw_data[year] + swden_data = ndbc.to_datetime_index(parameter, year_data) + swden_data = swden_data.T + swden_data = swden_data.sort_index() + ndbc_data[year] = swden_data + +wind_data = {} +if args.plot_wind is not None: + print("\nFetching wind data (stdmet)...") + try: + wind_parameter = "stdmet" + wind_ndbc_available = ndbc.available_data(wind_parameter, buoy_number) + wind_filenames = wind_ndbc_available["filename"] + wind_raw = ndbc.request_data(wind_parameter, wind_filenames) + + for year in wind_raw: + year_wind_data = wind_raw[year] + stdmet_data = ndbc.to_datetime_index(wind_parameter, year_wind_data) + stdmet_data = stdmet_data.T + stdmet_data = stdmet_data.sort_index() + wind_data[year] = stdmet_data + print(f" Wind data loaded for years: {list(wind_data.keys())}") + except Exception as e: + print(f" WARNING: Could not fetch wind data: {e}") + +# Handle date/time extraction, building one labeled spectrum per requested target +extracted_spectra = [] # list of (label, pd.Series) in the order requested + +if extract_single: + for target in extract_targets: + ts = target['timestamp'] + year_key = str(ts.year) + if year_key not in ndbc_data: + print(f"ERROR: Year {year_key} not found in available data") + exit(1) + + year_data = ndbc_data[year_key] + + if target['kind'] == 'range': + end_ts = target['end'] + window = [c for c in year_data.columns if ts <= c <= end_ts] + window_spectra = [year_data[c].dropna() for c in window] + window_spectra = [s for s in window_spectra if len(s) > 0] + if not window_spectra: + print(f"ERROR: No measurements between {ts:%Y-%m-%d %H:%M} and {end_ts:%H:%M}") + exit(1) + print(f" Found {len(window_spectra)} measurements in {target['label']}. Using average.") + extracted_spectra.append( + (target['label'], pd.concat(window_spectra, axis=1).mean(axis=1))) + continue + + try: + spectrum_clean = year_data[ts].dropna() + if len(spectrum_clean) == 0: + raise KeyError(ts) + print(f" Found measurement at {ts.strftime('%Y-%m-%d %H:%M')} with {len(spectrum_clean)} frequencies") + extracted_spectra.append((target['label'], spectrum_clean)) + except KeyError: + # Fall back to the closest available time within 1 hour + available_times = year_data.columns + closest_time = min(available_times, key=lambda x: abs((x - ts).total_seconds())) + time_diff = abs((closest_time - ts).total_seconds()) + + if time_diff < 3600: + print(f" Exact time not found for {ts.strftime('%Y-%m-%d %H:%M')}. " + f"Using closest available: {closest_time.strftime('%Y-%m-%d %H:%M')} ({time_diff/60:.0f} min difference)") + spectrum_clean = year_data[closest_time].dropna() + extracted_spectra.append((closest_time.strftime('%Y-%m-%d %H:%M'), spectrum_clean)) + else: + print(f"ERROR: No data near {ts.strftime('%Y-%m-%d %H:%M')}. " + f"Closest available: {closest_time.strftime('%Y-%m-%d %H:%M')} ({time_diff/3600:.1f} hours away)") + exit(1) + +elif extract_daily: + for target in extract_targets: + day_ts = target['timestamp'] + year_key = str(day_ts.year) + if year_key not in ndbc_data: + print(f"ERROR: Year {year_key} not found in available data") + exit(1) + + year_data = ndbc_data[year_key] + + # Collect every measurement on this calendar date + daily_spectra = [] + for timestamp, spectrum in year_data.items(): + if day_ts.date() == timestamp.date(): + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + daily_spectra.append(spectrum_clean) + + if len(daily_spectra) == 0: + print(f"ERROR: No data found for {day_ts.strftime('%Y-%m-%d')}") + exit(1) + + daily_avg = pd.concat(daily_spectra, axis=1).mean(axis=1) + print(f" Found {len(daily_spectra)} measurements for {day_ts.strftime('%Y-%m-%d')}. Using daily average.") + extracted_spectra.append((target['label'], daily_avg)) + +if extracted_spectra: + # all_swden_data drives the wave statistics and YAML update; keep every extracted + # spectrum so stats reflect the full selection + all_swden_data = pd.concat( + [s.rename(label) for label, s in extracted_spectra], axis=1 + ) + +# Extract wind data for the same measurement/day if wind flag is set +all_wind_data = None +if args.plot_wind is not None and len(wind_data) > 0: + if extract_single: + year_key = str(single_measurement_date.year) + if year_key in wind_data: + year_wind = wind_data[year_key] + try: + wind_measurement = year_wind[single_measurement_date] + all_wind_data = pd.DataFrame({single_measurement_date: wind_measurement}) + except KeyError: + # Try closest wind measurement + available_times = year_wind.columns + closest_time = min(available_times, key=lambda x: abs((x - single_measurement_date).total_seconds())) + time_diff = abs((closest_time - single_measurement_date).total_seconds()) + if time_diff < 3600: + wind_measurement = year_wind[closest_time] + all_wind_data = pd.DataFrame({closest_time: wind_measurement}) + + elif extract_daily: + year_key = str(daily_start.year) + if year_key in wind_data: + year_wind = wind_data[year_key] + # Filter for this day only + daily_winds = [] + for timestamp, wind_row in year_wind.items(): + if daily_start.date() == timestamp.date(): + wind_clean = wind_row.dropna() + if len(wind_clean) > 0: + daily_winds.append(wind_clean) + + if len(daily_winds) > 0: + daily_wind_avg = pd.concat(daily_winds, axis=1).mean(axis=1) + all_wind_data = pd.DataFrame({daily_start: daily_wind_avg}) + +# A file-supplied spectrum stands in for a single extracted measurement, so the +# statistics, plotting and eta-generation paths below need no special casing. +if spectrum_from_file: + _series = pd.Series(custom_S, index=pd.Index(custom_freq, name="Frequency")) + all_swden_data = pd.DataFrame({spectrum_label: _series}) + extracted_spectra = [(spectrum_label, _series)] + extract_single = True + +# Extract and print min/max frequencies only for irregular waves +if wave_type == "irregular" and (extract_single or extract_daily): + # Get first spectrum from extracted data + spectrum = all_swden_data.iloc[:, 0] + spectrum = spectrum.dropna() + + # Calculate min and max frequencies + freq_index = spectrum.index.values + min_freq = freq_index.min() + max_freq = freq_index.max() + print(f"\nSpectrum Frequency Range:") + print(f" Min frequency: {min_freq:.6f} Hz") + print(f" Max frequency: {max_freq:.6f} Hz") + +# Concatenate data from selected years only (if not in extract_single/daily mode) +if not extract_single and not extract_daily: + selected_year_keys = [str(y) for y in selected_years] + selected_swden_data = [ndbc_data[year_key] for year_key in selected_year_keys if year_key in ndbc_data] + + if not selected_swden_data: + print(f"ERROR: No data found for selected years {selected_years}") + exit(1) + + all_swden_data = pd.concat(selected_swden_data, axis=1, sort=False) + all_swden_data = all_swden_data.sort_index(axis=0) + all_swden_data = all_swden_data.sort_index(axis=1) +else: + # all_swden_data already set in single measurement extraction above + pass + +Hm0 = wave.resource.significant_wave_height( + all_swden_data) +Te = wave.resource.energy_period( + all_swden_data) +Tp = wave.resource.peak_period( + all_swden_data) + +print("\nWave Statistics") +print("----------------") +print( + f"Mean Significant Wave Height Hm0: " + f"{Hm0.mean():.2f} m") +print( + f"Mean Energy Period Te: " + f"{Te.mean():.2f} s") +print( + f"Mean Peak Period Tp: " + f"{Tp.mean():.2f} s") + +Hs_rep = Hm0.mean() +Te_rep = Te.mean() + +# Display wind statistics if requested +if args.plot_wind is not None and len(wind_data) > 0: + print("\nWind Statistics") + print("----------------") + + # Use extracted wind data if available (single/daily mode) or fetch from selected years + if all_wind_data is not None: + wind_df = all_wind_data + else: + # Concatenate wind data from selected years + selected_wind_keys = [str(y) for y in selected_years] + selected_wind_data = [wind_data[year_key] for year_key in selected_wind_keys if year_key in wind_data] + + if selected_wind_data: + wind_df = pd.concat(selected_wind_data, axis=1, sort=False) + else: + wind_df = None + + if wind_df is not None: + # NDBC stdmet columns: WSPD (wind speed m/s), WDIR (direction deg), GST (gusts m/s) + # Missing values use sentinels: 99.0 for speed/gust, 999 for direction + wind_speed = wind_df.loc['WSPD'] if 'WSPD' in wind_df.index else None + wind_direction = wind_df.loc['WDIR'] if 'WDIR' in wind_df.index else None + wind_gusts = wind_df.loc['GST'] if 'GST' in wind_df.index else None + + if wind_speed is not None: + valid_speed = wind_speed[wind_speed < 90] + if len(valid_speed) > 0: + print(f"Mean Wind Speed: {valid_speed.mean():.2f} m/s") + + if wind_direction is not None: + valid_dir = wind_direction[wind_direction <= 360] + if len(valid_dir) > 0: + print(f"Mean Wind Direction: {valid_dir.mean():.1f} deg") + + if wind_gusts is not None: + valid_gusts = wind_gusts[wind_gusts < 90] + if len(valid_gusts) > 0: + print(f"Mean Wind Gusts: {valid_gusts.mean():.2f} m/s") + +#~~~~~~~~~~~~~~ Custom spectrum -> surface elevation time series ~~~~~~~~~~~~~~ + +# EtaTableWaveField / ComponentSampler::BuildFromEtaFile require one "time:elevation" +# pair per line, colon delimited, no header and no blank lines. +eta_file_name = None + +if wave_type == "irregular" and args.spectrum == "custom" and not skip_yaml_update: + from mhkit.wave.resource import surface_elevation + + print("\n" + "=" * 60) + print("Custom spectrum -> surface elevation") + print("=" * 60) + + # Pick the source spectrum: an extracted measurement/day if one was requested, + # otherwise the average over the selected years + if extracted_spectra: + source_label, source_spectrum = extracted_spectra[0] + if len(extracted_spectra) > 1: + print(f"NOTE: {len(extracted_spectra)} spectra selected; using the first ({source_label})") + else: + source_spectrum = all_swden_data.mean(axis=1).dropna() + source_label = f"Average of {', '.join(map(str, selected_years))}" + + freq_measured = source_spectrum.index.astype(float).values + S_measured = source_spectrum.values.astype(float) + + print(f"Source spectrum: {source_label}") + print(f" Measured bins: {len(freq_measured)} " + f"({freq_measured.min():.4f} - {freq_measured.max():.4f} Hz, non-uniform spacing)") + + # The companion simulation YAML sets the eta sample rate via time_step, and its + # end_time must match the record: the wave field is zero outside the table range. + yaml_dir = os.path.dirname(os.path.abspath(yaml_file_path)) if yaml_file_path else os.getcwd() + sim_yaml_path = find_simulation_yaml(yaml_file_path) + + sim_time_step = None + if sim_yaml_path: + with open(sim_yaml_path, 'r') as f: + sim_data = yaml.safe_load(f) or {} + raw_step = (sim_data.get('simulation') or {}).get('time_step') + if raw_step is not None: + sim_time_step = float(raw_step) + + # Build the elevation time series + duration = float(args.elevation_duration) if args.elevation_duration is not None else 600.0 + if args.elevation_dt is not None: + dt = float(args.elevation_dt) + elif sim_time_step and sim_time_step > 0.0: + dt = sim_time_step + print(f"\nEta dt = {dt} s (from simulation.time_step in " + f"{os.path.basename(sim_yaml_path)})") + else: + dt = 0.05 + + # Only a concern when the solver step is left alone; --elevation_dt rewrites it below + if args.elevation_dt is None and sim_time_step and dt > sim_time_step: + print(f"\nWARNING: eta dt ({dt} s) is coarser than simulation.time_step " + f"({sim_time_step} s). The solver steps faster than the record is sampled.") + + n_samples = int(round(duration / dt)) + 1 + time_vec = np.arange(n_samples) * dt + + # The frequency grid is dictated by the time vector, NOT chosen freely. + # surface_elevation's ifft calls np.fft.irfftn(..., s=[n_samples]), which reads bin + # i as i/(n_samples*dt) and ignores the index on S. Any other spacing silently + # rescales every frequency (a 512-bin 0-0.485 Hz grid over 120 s stretched the + # spectrum 8.8x, turning 17 s swell into 2 s chop). + df_uniform = 1.0 / (n_samples * dt) + n_freq = n_samples // 2 + 1 + freq_uniform = np.arange(n_freq) * df_uniform + f_max = float(freq_measured.max()) + + nyquist = 1.0 / (2.0 * dt) + if f_max > nyquist: + print(f" WARNING: spectrum extends to {f_max:.3f} Hz but the Nyquist frequency for " + f"dt={dt} s is {nyquist:.3f} Hz. Energy above Nyquist will alias; " + f"reduce --elevation_dt to resolve it.") + + # Interpolate; energy outside the measured band is zero (no extrapolation) + S_uniform = np.interp(freq_uniform, freq_measured, S_measured, left=0.0, right=0.0) + S_uniform[freq_uniform < freq_measured.min()] = 0.0 + + # Preserve total variance (m0) so Hm0 of the realization matches the measurement + m0_measured = float(integrate.trapezoid(S_measured, freq_measured)) + m0_uniform = float(integrate.trapezoid(S_uniform, freq_uniform)) + if m0_uniform > 0: + S_uniform *= (m0_measured / m0_uniform) + + hm0_measured = 4.0 * np.sqrt(m0_measured) + hm0_uniform = 4.0 * np.sqrt(float(integrate.trapezoid(S_uniform, freq_uniform))) + n_in_band = int(np.count_nonzero(S_uniform > 0.0)) + print(f" Record grid: {n_freq} bins (0 - {freq_uniform[-1]:.4f} Hz, df = {df_uniform:.6f} Hz)") + print(f" Bins carrying measured energy: {n_in_band}") + print(f" Hm0 measured = {hm0_measured:.4f} m | interpolated = {hm0_uniform:.4f} m") + if n_in_band < 30: + print(f" WARNING: only {n_in_band} bins fall inside the measured band. The sea state will") + print(f" look repetitive; lengthen --elevation_duration for a finer df.") + + # Plot measured vs interpolated spectrum + plt.figure(figsize=(12, 6)) + plt.plot(freq_measured, S_measured, 'o-', label=f"Measured ({len(freq_measured)} bins)", + linewidth=2.0, markersize=4, color='steelblue', alpha=0.9) + plt.plot(freq_uniform, S_uniform, '-', label=f"Interpolated (df = {df_uniform:.5f} Hz)", + linewidth=1.5, color='coral', alpha=0.85) + plt.xlim(0.0, f_max) + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"{source_title} - Measured vs Interpolated Spectrum\n{source_label}", fontsize=13) + plt.legend(fontsize=10, loc='upper right') + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + + S_df = pd.DataFrame({0: S_uniform}, index=pd.Index(freq_uniform, name="Frequency")) + seed = args.seed if args.seed is not None else 42 + + print(f"\nGenerating elevation: {duration:.0f} s at dt={dt} s ({n_samples} samples), seed={seed}") + eta = surface_elevation(S_df, time_vec, seed=seed) + eta_values = np.asarray(eta).squeeze() + + hm0_realized = 4.0 * np.std(eta_values) + print(f" Elevation range: {eta_values.min():.3f} to {eta_values.max():.3f} m") + print(f" Realized Hm0 (4*sigma) = {hm0_realized:.4f} m (target {hm0_measured:.4f} m)") + + # Write the eta file next to the YAML + if args.elevation_file: + eta_file_name = args.elevation_file.replace('\\', '/') + else: + yaml_stem = os.path.splitext(os.path.basename(yaml_file_path))[0] if yaml_file_path else "waves" + yaml_stem = yaml_stem.replace('.hydro', '') + eta_file_name = f"{yaml_stem}_eta.txt" + + eta_path = os.path.join(yaml_dir, eta_file_name) + # Do not create directories; the target must already exist + eta_parent = os.path.dirname(eta_path) + if eta_parent and not os.path.isdir(eta_parent): + print(f"\nERROR: directory does not exist: {eta_parent}") + print(f" Create it first, or choose a different --elevation_file.") + exit(1) + + with open(eta_path, 'w', newline='\n') as f: + # No header, no comments, no trailing blank line: the C++ reader parses + # every line with (ss >> t >> ':' >> eta) and throws on anything else + f.write("\n".join(f"{t:.6f}:{e:.6f}" for t, e in zip(time_vec, eta_values))) + f.write("\n") + + print(f"\nWrote eta file: {eta_path}") + print(f" {n_samples} samples, format 'time:elevation'") + + # BuildFromEtaFile only evaluates the orthogonal Fourier frequencies k/(N*dt) and + # strides, dropping energy, when nfrequencies is below the in-band bin count. + t_window = n_samples * dt + df_natural = 1.0 / t_window + k_min = max(1, int(np.ceil(0.001 / df_natural))) + k_max = int(np.floor(min(1.0, 1.0 / (2.0 * dt)) / df_natural)) + eta_nfrequencies = max(1, k_max - k_min + 1) + + if args.eta_method == "dft": + print(f" nfrequencies: {eta_nfrequencies} (every Fourier bin below 1 Hz; " + f"df = {df_natural:.3e} Hz)") + if eta_nfrequencies > 4000: + print(f" WARNING: {eta_nfrequencies} wave components is a lot. The DFT is O(nf*N)") + print(f" and the GUI evaluates every component per surface vertex per frame.") + print(f" Shorten --elevation_duration or increase --elevation_dt to reduce it.") + else: + print(f" method: irf_convolution (point elevation table, no spatial free surface)") + + # Keep the companion simulation YAML in sync: the wave field returns zero + # elevation outside the table range, so end_time must not exceed the record. + if sim_yaml_path is None: + print(f"\nNOTE: no companion *.simulation.yaml found; timing not updated.") + else: + update_simulation_yaml(sim_yaml_path, end_time=duration, + time_step=args.elevation_dt, + reason=" s (matches eta record length)") + + +#~~~~~~~~~~~~~~ Update YAML file with calculated wave parameters ~~~~~~~~~~~~~~ + +if not skip_yaml_update and os.path.exists(yaml_file_path): + print(f"\n\nUpdating YAML file: {yaml_file_path}") + + with open(yaml_file_path, 'r') as f: + yaml_data = yaml.safe_load(f) + + if 'hydrodynamics' in yaml_data and 'waves' in yaml_data['hydrodynamics']: + # With an eta file the sea state comes entirely from the record, so the waves + # block is rebuilt rather than layered on top of what a previous run left. + custom_eta = (wave_type == "irregular" and args.spectrum == "custom") + + yaml_data['hydrodynamics']['waves']['type'] = wave_type + + if not custom_eta: + # Use override values if provided, otherwise use calculated values + height_value = args.height if args.height is not None else Hs_rep + period_value = args.period if args.period is not None else Te_rep + + yaml_data['hydrodynamics']['waves']['height'] = float(round(height_value, 2)) + yaml_data['hydrodynamics']['waves']['period'] = float(round(period_value, 2)) + + # Handle spectrum: only for irregular waves + if wave_type == "irregular": + if args.spectrum == "custom": + # dft -> BuildFromEtaFile gives discrete components and a spatial free + # surface, so the GUI wireframe renders. irf -> EtaTableWaveField is a + # point time series: kinematics are zero and no wireframe is drawn. + waves_block = {'type': wave_type} + + if args.eta_method == "irf": + waves_block['method'] = "irf_convolution" + waves_block['eta_file'] = eta_file_name + else: + # Legacy DFT path: waves.method is omitted entirely. + waves_block['eta_file'] = eta_file_name + # nfrequencies must span the default 0.001 - 1.0 Hz eta band; the + # C++ default of 1000 strides the grid and attenuates the record. + waves_block['nfrequencies'] = int(eta_nfrequencies) + + ramp = float(args.ramp_time) if args.ramp_time is not None else 60.0 + if ramp > 0.0: + waves_block['ramp_duration'] = ramp + waves_block['ramp_type'] = "cosine" + + if args.depth is not None: + # Only when asked for. With no depth key the DFT path assumes deep + # water and EtaTableWaveField falls back to the H5 water depth. + waves_block['depth'] = float(args.depth) + + yaml_data['hydrodynamics']['waves'] = waves_block + elif args.spectrum: + # Normalize spectrum name to lowercase single char or full name + spectrum_type = args.spectrum.lower() + if spectrum_type in ["pm", "pierson_moskowitz"]: + yaml_data['hydrodynamics']['waves']['spectrum'] = "pierson_moskowitz" + # For PM, use peak period (Tp) if not overridden + if args.period is None: + yaml_data['hydrodynamics']['waves']['period'] = float(round(Tp.mean(), 2)) + elif spectrum_type in ["jonswap", "js"]: + yaml_data['hydrodynamics']['waves']['spectrum'] = "jonswap" + # For JONSWAP, use peak period (Tp) if not overridden + if args.period is None: + yaml_data['hydrodynamics']['waves']['period'] = float(round(Tp.mean(), 2)) + + # Check if gamma was provided as override + if args.gamma: + # Use first gamma value as override + yaml_data['hydrodynamics']['waves']['gamma'] = float(round(args.gamma[0], 3)) + else: + # Calculate fitted gamma for JONSWAP + # IMPORTANT: Average spectrum first, then calculate Hm0/Tp (same as plotting) + avg_spectrum = all_swden_data.mean(axis=1) + avg_spectrum_clean = avg_spectrum.dropna() + + if len(avg_spectrum_clean) > 0: + # Use period override if provided, otherwise calculate from averaged spectrum + if args.period is not None: + Tp_mean = float(args.period) + else: + Tp_mean = float(wave.resource.peak_period(avg_spectrum_clean)) + + # Use height override if provided, otherwise calculate from averaged spectrum + if args.height is not None: + Hs_mean = float(args.height) + else: + Hs_mean = float(wave.resource.significant_wave_height(avg_spectrum_clean)) + + freq_array = avg_spectrum_clean.index.astype(float).values + gamma_fit, error, tp_fit, hit_bound = fit_jonswap_gamma(avg_spectrum_clean, freq_array, Tp_mean, Hs_mean) + yaml_data['hydrodynamics']['waves']['gamma'] = float(round(gamma_fit, 3)) + else: + # Fallback to default gamma if no spectrum data + yaml_data['hydrodynamics']['waves']['gamma'] = 3.3 + else: + # Default spectrum for irregular waves is PM + yaml_data['hydrodynamics']['waves']['spectrum'] = "pierson_moskowitz" + yaml_data['hydrodynamics']['waves']['period'] = float(round(Tp.mean(), 2)) + else: + # Regular waves: deterministic, so spectrum shape keys do not apply + for stale_key in ('spectrum', 'gamma'): + if stale_key in yaml_data['hydrodynamics']['waves']: + del yaml_data['hydrodynamics']['waves'][stale_key] + + # Drop eta-import leftovers whenever a parametric spectrum is in use, so a + # previous --spectrum custom run does not keep overriding the sea state. + # ramp keys are re-added below when --ramp_time asks for them. + if args.spectrum != "custom": + for stale_key in ('eta_file', 'eta_file_path', 'method', 'nfrequencies', + 'ramp_duration', 'ramp_type'): + if stale_key in yaml_data['hydrodynamics']['waves']: + del yaml_data['hydrodynamics']['waves'][stale_key] + + + if not custom_eta: + # Add water depth if available + if args.depth is not None: + yaml_data['hydrodynamics']['waves']['depth'] = float(args.depth) + elif water_depth: + yaml_data['hydrodynamics']['waves']['depth'] = float(water_depth) + + # Add optional parameters if provided + if args.seed is not None: + yaml_data['hydrodynamics']['waves']['seed'] = args.seed + + if args.ramp_time is not None: + # sea-stack reads waves.ramp_duration / waves.ramp_type; there is no + # waves.ramp_time key. + yaml_data['hydrodynamics']['waves']['ramp_duration'] = float(args.ramp_time) + yaml_data['hydrodynamics']['waves']['ramp_type'] = "cosine" + + print(f"Updated wave parameters:") + print(f" Type: {yaml_data['hydrodynamics']['waves']['type']}") + if custom_eta: + pass + elif args.height is not None: + print(f" Height: {yaml_data['hydrodynamics']['waves']['height']} m (OVERRIDE)") + else: + print(f" Height: {yaml_data['hydrodynamics']['waves']['height']} m (from {Hs_rep:.4f} m)") + + # Period handling - varies by spectrum type + if custom_eta: + wb = yaml_data['hydrodynamics']['waves'] + print(f" eta_file: {wb['eta_file']}") + if 'method' in wb: + print(f" method: {wb['method']}") + else: + print(f" nfrequencies: {wb['nfrequencies']}") + if 'ramp_duration' in wb: + print(f" ramp: {wb['ramp_duration']} s ({wb['ramp_type']})") + if 'depth' in wb: + print(f" depth: {wb['depth']} m") + elif wave_type == "irregular" and args.spectrum: + spectrum_type = args.spectrum.lower() + if spectrum_type in ["pm", "pierson_moskowitz"]: + if args.period is not None: + print(f" Period: {yaml_data['hydrodynamics']['waves']['period']} s (OVERRIDE)") + else: + print(f" Period: {yaml_data['hydrodynamics']['waves']['period']} s (Tp, from {Tp.mean():.4f} s)") + elif spectrum_type in ["jonswap", "js"]: + if args.period is not None: + print(f" Period: {yaml_data['hydrodynamics']['waves']['period']} s (OVERRIDE)") + else: + print(f" Period: {yaml_data['hydrodynamics']['waves']['period']} s (Tp, from {Tp.mean():.4f} s)") + if 'gamma' in yaml_data['hydrodynamics']['waves']: + print(f" Spectrum: JONSWAP") + gamma_source = "override" if args.gamma else "fitted" + print(f" Gamma: {yaml_data['hydrodynamics']['waves']['gamma']} ({gamma_source} peak enhancement factor)") + elif wave_type == "irregular": + # Default PM for irregular + print(f" Period: {yaml_data['hydrodynamics']['waves']['period']} s (Tp, from {Tp.mean():.4f} s, default PM)") + else: + # Regular waves + if args.period is not None: + print(f" Period: {yaml_data['hydrodynamics']['waves']['period']} s (OVERRIDE)") + else: + print(f" Period: {yaml_data['hydrodynamics']['waves']['period']} s (from {Te_rep:.4f} s)") + + if 'spectrum' in yaml_data['hydrodynamics']['waves']: + print(f" Spectrum: {yaml_data['hydrodynamics']['waves']['spectrum']}") + if not custom_eta: + if 'depth' in yaml_data['hydrodynamics']['waves']: + if args.depth is not None: + print(f" Depth: {yaml_data['hydrodynamics']['waves']['depth']} m (OVERRIDE)") + else: + print(f" Depth: {yaml_data['hydrodynamics']['waves']['depth']} m") + if args.seed is not None: + print(f" Seed: {yaml_data['hydrodynamics']['waves']['seed']}") + if args.ramp_time is not None: + print(f" Ramp: {yaml_data['hydrodynamics']['waves']['ramp_duration']} s (cosine)") + + # with open(yaml_file_path, 'w') as f: + # yaml.dump(yaml_data, f, default_flow_style=False, sort_keys=False) + + class IndentDumper(yaml.SafeDumper): + def increase_indent(self, flow=False, indentless=False): + return super().increase_indent(flow, False) + + with open(yaml_file_path, "w") as f: + yaml.dump( + yaml_data, + f, + Dumper=IndentDumper, + default_flow_style=False, + sort_keys=False, + ) + + print(f"\nYAML file updated successfully!") + else: + print("ERROR: 'hydrodynamics.waves' section not found in YAML file") +elif not skip_yaml_update: + # Trying to update but file not found + print(f"ERROR: YAML file not found at {yaml_file_path}") + exit(1) + +# Notify user if YAML was not updated +if skip_yaml_update: + print(f"\nNo YAML update requested (plot only). Skipping YAML file update.") + +# --elevation_duration / --elevation_dt set the simulation timing for every spectrum +# type. On the custom path they also size the eta record and were applied with it. +if (not skip_yaml_update and args.spectrum != "custom" + and (args.elevation_duration is not None or args.elevation_dt is not None)): + _sim_path = find_simulation_yaml(yaml_file_path) + if _sim_path is None: + print("\nNOTE: no companion *.simulation.yaml found; timing not updated.") + else: + update_simulation_yaml(_sim_path, end_time=args.elevation_duration, + time_step=args.elevation_dt) + + +#~~~~~~~~~~~~~~ Plot spectrum if requested ~~~~~~~~~~~~~~ + +# Handle spectrum plotting if requested +if args.plot_spectrum is not None: + # If in extract_single or extract_daily mode, plot every extracted spectrum + if extract_single or extract_daily: + plt.figure(figsize=(12, 6)) + colors = plt.cm.tab10(np.linspace(0, 1, max(len(extracted_spectra), 1))) + + for idx, (label, spectrum) in enumerate(extracted_spectra): + plt.plot( + spectrum.index.astype(float), + spectrum.values, + label=label, + linewidth=2.5, + color=colors[idx], + alpha=0.9 + ) + + mode_name = "Measurement" if extract_single else "Daily Average" + if len(extracted_spectra) == 1: + title = plot_heading(extracted_spectra[0][0]) + else: + title = f"{source_title} - {len(extracted_spectra)} {mode_name}s" + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(title, fontsize=14) + if len(extracted_spectra) > 1: + plt.legend(fontsize=10, loc='upper right') + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + print(f"Spectrum plot displayed") + + # One or more months requested: overlay each month's average spectrum + elif parse_month_list(plot_spectrum_values, min_count=1): + month_nums = parse_month_list(plot_spectrum_values, min_count=1) + month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + if not selected_years: + print("No valid years found to plot") + else: + plot_year = selected_years[0] + print(f"\nPlotting {len(month_nums)} monthly averaged spectra for {plot_year}...") + + if str(plot_year) not in ndbc_data: + print(f"Year {plot_year} not found in available data") + else: + year_data = ndbc_data[str(plot_year)] + plt.figure(figsize=(12, 6)) + colors = plt.cm.tab10(np.linspace(0, 1, max(len(month_nums), 1))) + plotted_any = False + + for idx, month_num in enumerate(month_nums): + month_spectra = [] + for timestamp, spectrum in year_data.items(): + if timestamp.month == month_num: + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + month_spectra.append(spectrum_clean) + + if not month_spectra: + print(f" No data for {month_names[month_num-1]} {plot_year}") + continue + + month_avg = pd.concat(month_spectra, axis=1).mean(axis=1) + month_avg_series = pd.Series(month_avg.values, index=month_avg.index) + + hm0 = float(wave.resource.significant_wave_height(month_avg_series)) + te = float(wave.resource.energy_period(month_avg_series)) + tp = float(wave.resource.peak_period(month_avg_series)) + + label = month_names[month_num-1] + plt.plot(month_avg.index.astype(float), month_avg.values, + label=label, linewidth=2.5, color=colors[idx], alpha=0.9) + + print(f"{label} {plot_year}: Hm0={hm0:6.2f} m | Te={te:6.2f} s | Tp={tp:6.2f} s") + plotted_any = True + + if plotted_any: + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Monthly Averaged Spectra ({plot_year})", fontsize=14) + plt.legend(fontsize=10, loc='upper right') + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + print(f"Spectrum plot displayed") + else: + plt.close() + print("No month data available to plot") + + # Normal plotting logic (when not in extract mode) + elif isinstance(args.plot_spectrum, str): + # Plot monthly spectra for specified year + year_for_monthly = args.plot_spectrum + print(f"\nPlotting monthly averaged spectra for year {year_for_monthly}...") + + if year_for_monthly in ndbc_data: + year_data = ndbc_data[year_for_monthly] + + # Group data by month + monthly_spectra = {} + for timestamp, spectrum in year_data.items(): + month = timestamp.month + if month not in monthly_spectra: + monthly_spectra[month] = [] + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + monthly_spectra[month].append(spectrum_clean) + + # Plot monthly spectra + plt.figure(figsize=(12, 6)) + colors = plt.cm.viridis(np.linspace(0, 1, len(monthly_spectra))) + month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + + print(f"\nMonthly Wave Metrics for {year_for_monthly}:") + print("-" * 50) + + for idx, (month, spectra_list) in enumerate(sorted(monthly_spectra.items())): + if len(spectra_list) > 0: + # Average spectra for this month + month_avg = pd.concat(spectra_list, axis=1).mean(axis=1) + + # Calculate wave metrics for this month + month_avg_series = pd.Series(month_avg.values, index=month_avg.index) + hm0 = wave.resource.significant_wave_height(month_avg_series) + te = wave.resource.energy_period(month_avg_series) + tp = wave.resource.peak_period(month_avg_series) + + print(f"{month_names[month-1]:>3}: Hm0={hm0:6.2f} m | Te={te:6.2f} s | Tp={tp:6.2f} s") + + plt.plot( + month_avg.index.astype(float), + month_avg.values, + label=f"{month_names[month-1]} (n={len(spectra_list)})", + linewidth=2, + color=colors[idx], + alpha=0.8 + ) + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Monthly Averaged Wave Spectra ({year_for_monthly})", fontsize=14) + plt.legend(fontsize=10, loc='upper right') + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + print(f"Monthly spectrum plot displayed") + else: + print(f"Year {year_for_monthly} not found in available data") + + elif args.plot_spectrum is True: + if selected_years: + plt.figure(figsize=(12, 6)) + colors = plt.cm.viridis(np.linspace(0, 1, len(selected_years))) + + all_year_avg_spectra = [] + + for idx, year in enumerate(selected_years): + year_key = str(year) + if year_key in ndbc_data: + year_data = ndbc_data[year_key] + year_avg_spectrum = year_data.mean(axis=1) + year_avg_spectrum = year_avg_spectrum.dropna() + all_year_avg_spectra.append(year_avg_spectrum) + plt.plot( + year_avg_spectrum.index.astype(float), + year_avg_spectrum.values, + label=f"Year {year}", + linewidth=2, + color=colors[idx], + alpha=0.8 + ) + + if len(all_year_avg_spectra) > 1: + overall_avg_spectrum = pd.concat(all_year_avg_spectra, axis=1).mean(axis=1) + plt.plot( + overall_avg_spectrum.index.astype(float), + overall_avg_spectrum.values, + label="Overall Average", + linewidth=2.5, + linestyle="--", + color="red", + alpha=0.9 + ) + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Year-Averaged Wave Spectra", fontsize=14) + plt.legend(fontsize=10) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + else: + print("No valid years found to plot") + + +#~~~~~~~~~~~~~~ Plot spectrum vs Pierson-Moskowitz if requested ~~~~~~~~~~~~~~ + +if args.plot_spectrum_pm is not None: + # If in extract_single or extract_daily mode, plot each extracted spectrum with its PM estimate + if extract_single or extract_daily: + plt.figure(figsize=(12, 6)) + colors = plt.cm.tab10(np.linspace(0, 1, max(len(extracted_spectra), 1))) + multi = len(extracted_spectra) > 1 + + for idx, (label, spectrum) in enumerate(extracted_spectra): + spectrum_series = pd.Series(spectrum.values, index=spectrum.index) + + hm0 = float(wave.resource.significant_wave_height(spectrum_series)) + te = float(wave.resource.energy_period(spectrum_series)) + tp = float(wave.resource.peak_period(spectrum_series)) + + # Use override values if provided + hm0_plot = float(args.height) if args.height is not None else hm0 + tp_plot = float(args.period) if args.period is not None else tp + + freq = spectrum.index.astype(float).values + + measured_label = f"{label} (measured)" if multi else "Measured" + plt.plot( + freq, + spectrum.values, + label=measured_label, + linewidth=2.5, + color=colors[idx], + alpha=0.9 + ) + + pm_spec = pierson_moskowitz_spectrum(freq, hm0_plot, tp_plot) + pm_label = f"{label} PM" if multi else f"Pierson-Moskowitz (Hm0={hm0_plot:.2f} m, Tp={tp_plot:.2f} s)" + plt.plot( + freq, + pm_spec, + label=pm_label, + linewidth=2.0, + color=colors[idx], + alpha=0.75, + linestyle='--' + ) + + error = float(np.sum((spectrum.values - pm_spec) ** 2)) + print(f"\n{label}") + print(f"Hm0={hm0:.2f} m | Te={te:.2f} s | Tp={tp:.2f} s") + print(f"Pierson-Moskowitz: Error={error:.6e}") + + mode_name = "Measurement" if extract_single else "Daily Average" + if multi: + title = f"{source_title} - {len(extracted_spectra)} {mode_name}s vs Pierson-Moskowitz" + else: + title = plot_heading(extracted_spectra[0][0], "Pierson-Moskowitz") + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(title, fontsize=14) + plt.legend(fontsize=9, loc='upper right', ncol=2 if multi else 1) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + print(f"\nPierson-Moskowitz spectrum plot displayed") + + # Multiple months requested: overlay each month's average with its PM estimate + elif parse_month_list(plot_spectrum_pm_values): + month_nums = parse_month_list(plot_spectrum_pm_values) + month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + if not selected_years: + print("No valid years found to plot") + else: + plot_year = selected_years[0] + print(f"\nPlotting spectrum vs Pierson-Moskowitz for {len(month_nums)} months of {plot_year}...") + + if str(plot_year) not in ndbc_data: + print(f"Year {plot_year} not found in available data") + else: + year_data = ndbc_data[str(plot_year)] + plt.figure(figsize=(12, 6)) + colors = plt.cm.tab10(np.linspace(0, 1, max(len(month_nums), 1))) + plotted_any = False + + for idx, month_num in enumerate(month_nums): + month_spectra = [] + for timestamp, spectrum in year_data.items(): + if timestamp.month == month_num: + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + month_spectra.append(spectrum_clean) + + if not month_spectra: + print(f" No data for {month_names[month_num-1]} {plot_year}") + continue + + month_avg = pd.concat(month_spectra, axis=1).mean(axis=1) + month_avg_series = pd.Series(month_avg.values, index=month_avg.index) + + hm0 = float(wave.resource.significant_wave_height(month_avg_series)) + te = float(wave.resource.energy_period(month_avg_series)) + tp = float(wave.resource.peak_period(month_avg_series)) + + hm0_plot = float(args.height) if args.height is not None else hm0 + tp_plot = float(args.period) if args.period is not None else tp + + freq = month_avg.index.astype(float).values + label = month_names[month_num-1] + + plt.plot(freq, month_avg.values, label=f"{label} (measured)", + linewidth=2.5, color=colors[idx], alpha=0.9) + + pm_spec = pierson_moskowitz_spectrum(freq, hm0_plot, tp_plot) + plt.plot(freq, pm_spec, label=f"{label} PM", + linewidth=2.0, color=colors[idx], alpha=0.75, linestyle='--') + + error = float(np.sum((month_avg.values - pm_spec) ** 2)) + print(f"{label} {plot_year}: Hm0={hm0:6.2f} m | Te={te:6.2f} s | Tp={tp:6.2f} s | Error={error:.6e}") + plotted_any = True + + if plotted_any: + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Monthly Measured vs Pierson-Moskowitz ({plot_year})", + fontsize=14) + plt.legend(fontsize=9, loc='upper right', ncol=2) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + print(f"\nPierson-Moskowitz spectrum plot displayed") + else: + plt.close() + print("No month data available to plot") + + # Check if year or month was provided as argument + elif isinstance(args.plot_spectrum_pm, str): + # Determine if input is a month or year + month_num, is_month = parse_month(args.plot_spectrum_pm) + + if is_month: + # Plot specific month for the selected year + month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + if selected_years: + plot_year = selected_years[0] + print(f"\nPlotting spectrum vs Pierson-Moskowitz for {month_names[month_num-1]} {plot_year}...") + + if str(plot_year) in ndbc_data: + year_data = ndbc_data[str(plot_year)] + + # Filter for just this month + month_spectra = [] + for timestamp, spectrum in year_data.items(): + if timestamp.month == month_num: + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + month_spectra.append(spectrum_clean) + + if month_spectra: + # Average spectra for this month + month_avg = pd.concat(month_spectra, axis=1).mean(axis=1) + month_avg_series = pd.Series(month_avg.values, index=month_avg.index) + + # Calculate wave metrics + hm0 = wave.resource.significant_wave_height(month_avg_series) + te = wave.resource.energy_period(month_avg_series) + tp = wave.resource.peak_period(month_avg_series) + + # Use override values if provided + hm0_plot = args.height if args.height is not None else hm0 + tp_plot = args.period if args.period is not None else tp + + print(f"{month_names[month_num-1]} {plot_year}: Hm0={hm0:.2f} m | Te={te:.2f} s | Tp={tp:.2f} s") + if args.height is not None or args.period is not None: + print(f" (Using overrides: Hm0={hm0_plot} m, Tp={tp_plot} s)") + + # Plot measured vs PM + plt.figure(figsize=(12, 6)) + + plt.plot( + month_avg.index.astype(float), + month_avg.values, + label=f"{month_names[month_num-1]} (measured)", + linewidth=2.5, + color='steelblue', + alpha=0.9 + ) + + # Calculate and plot PM spectrum + freq = month_avg.index.astype(float) + pm_curve = pierson_moskowitz_spectrum(freq, hm0_plot, tp_plot) + + plt.plot( + freq, + pm_curve, + label=f"{month_names[month_num-1]} (PM)", + linewidth=1.5, + color='coral', + alpha=0.7, + linestyle='--' + ) + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - {month_names[month_num-1]} {plot_year} Measured vs Pierson-Moskowitz Spectrum", fontsize=14) + plt.legend(fontsize=10, loc='upper right') + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + else: + print(f"No data found for {month_names[month_num-1]} {plot_year}") + else: + print(f"Year {plot_year} not found in available data") + else: + print("No years found to plot") + else: + # Plot monthly spectra vs PM for specified year + year_for_pm = args.plot_spectrum_pm + print(f"\nPlotting spectrum vs Pierson-Moskowitz for year {year_for_pm}...") + + if year_for_pm in ndbc_data: + year_data = ndbc_data[year_for_pm] + + # Group data by month + monthly_spectra = {} + for timestamp, spectrum in year_data.items(): + month = timestamp.month + if month not in monthly_spectra: + monthly_spectra[month] = [] + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + monthly_spectra[month].append(spectrum_clean) + + # Plot monthly spectra vs PM + plt.figure(figsize=(14, 8)) + colors = plt.cm.viridis(np.linspace(0, 1, len(monthly_spectra))) + month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + + print(f"\nMonthly Spectrum Comparison for {year_for_pm}:") + print("-" * 70) + + for idx, (month, spectra_list) in enumerate(sorted(monthly_spectra.items())): + if len(spectra_list) > 0: + # Average spectra for this month + month_avg = pd.concat(spectra_list, axis=1).mean(axis=1) + month_avg_series = pd.Series(month_avg.values, index=month_avg.index) + + # Calculate wave metrics + hm0 = wave.resource.significant_wave_height(month_avg_series) + te = wave.resource.energy_period(month_avg_series) + tp = wave.resource.peak_period(month_avg_series) + + # Use override values if provided + hm0_plot = args.height if args.height is not None else hm0 + tp_plot = args.period if args.period is not None else tp + + print(f"{month_names[month-1]:>3}: Hm0={hm0:6.2f} m | Te={te:6.2f} s | Tp={tp:6.2f} s") + if args.height is not None or args.period is not None: + print(f" (Using overrides: Hm0={hm0_plot} m, Tp={tp_plot} s)") + + # Plot measured spectrum + plt.plot( + month_avg.index.astype(float), + month_avg.values, + label=f"{month_names[month-1]} (measured)", + linewidth=2.5, + color=colors[idx], + alpha=0.9 + ) + + # Calculate and plot PM spectrum (using Tp, not Te) + freq = month_avg.index.astype(float) + pm_curve = pierson_moskowitz_spectrum(freq, hm0_plot, tp_plot) + + # Diagnostic: check spectrum properties + # Integrate to verify Hm0 matches + df = freq[1] - freq[0] if len(freq) > 1 else 1.0 + m0_measured = integrate.trapezoid(month_avg.values, freq) + m0_pm = integrate.trapezoid(pm_curve, freq) + hm0_from_m0_measured = 4 * np.sqrt(m0_measured) + hm0_from_m0_pm = 4 * np.sqrt(m0_pm) + + print(f" Diagnostics: m0_meas={m0_measured:.4f}, Hm0_calc={hm0_from_m0_measured:.2f} vs input={hm0:.2f}") + print(f" PM diagnostics: m0_pm={m0_pm:.4f}, Hm0_pm={hm0_from_m0_pm:.2f}, peak_meas={month_avg.max():.4f}, peak_pm={pm_curve.max():.4f}") + + plt.plot( + freq, + pm_curve, + label=f"{month_names[month-1]} (PM)", + linewidth=1.5, + color=colors[idx], + alpha=0.5, + linestyle='--' + ) + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Monthly Measured vs Pierson-Moskowitz Spectra ({year_for_pm})", fontsize=14) + plt.legend(fontsize=9, loc='upper right', ncol=2) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + else: + print(f"Year {year_for_pm} not found in available data") + + elif args.plot_spectrum_pm is True: + # Plot year-averaged spectra vs PM + if selected_years: + plt.figure(figsize=(14, 8)) + colors = plt.cm.viridis(np.linspace(0, 1, len(selected_years))) + + all_year_avg_spectra = [] + + print(f"\nYear-Averaged Spectrum Comparison:") + print("-" * 70) + + for idx, year in enumerate(selected_years): + year_key = str(year) + if year_key in ndbc_data: + year_data = ndbc_data[year_key] + year_avg_spectrum = year_data.mean(axis=1) + year_avg_spectrum = year_avg_spectrum.dropna() + year_avg_series = pd.Series(year_avg_spectrum.values, index=year_avg_spectrum.index) + + # Calculate wave metrics + hm0 = wave.resource.significant_wave_height(year_avg_series) + te = wave.resource.energy_period(year_avg_series) + tp = wave.resource.peak_period(year_avg_series) + + # Use override values if provided + hm0_plot = args.height if args.height is not None else hm0 + tp_plot = args.period if args.period is not None else tp + + print(f"Year {year}: Hm0={hm0:6.2f} m | Te={te:6.2f} s | Tp={tp:6.2f} s") + if args.height is not None or args.period is not None: + print(f" (Using overrides: Hm0={hm0_plot} m, Tp={tp_plot} s)") + + all_year_avg_spectra.append(year_avg_spectrum) + + # Plot measured spectrum + plt.plot( + year_avg_spectrum.index.astype(float), + year_avg_spectrum.values, + label=f"Year {year} (measured)", + linewidth=2.5, + color=colors[idx], + alpha=0.9 + ) + + # Calculate and plot PM spectrum (using Tp, not Te) + freq = year_avg_spectrum.index.astype(float) + pm_curve = pierson_moskowitz_spectrum(freq, hm0_plot, tp_plot) + plt.plot( + freq, + pm_curve, + label=f"Year {year} (PM)", + linewidth=1.5, + color=colors[idx], + alpha=0.5, + linestyle='--' + ) + + if len(all_year_avg_spectra) > 1: + overall_avg_spectrum = pd.concat(all_year_avg_spectra, axis=1).mean(axis=1) + overall_avg_series = pd.Series(overall_avg_spectrum.values, index=overall_avg_spectrum.index) + + hm0_overall = wave.resource.significant_wave_height(overall_avg_series) + te_overall = wave.resource.energy_period(overall_avg_series) + tp_overall = wave.resource.peak_period(overall_avg_series) + + # Use override values if provided + hm0_overall_plot = args.height if args.height is not None else hm0_overall + tp_overall_plot = args.period if args.period is not None else tp_overall + + print(f"Overall: Hm0={hm0_overall:6.2f} m | Te={te_overall:6.2f} s | Tp={tp_overall:6.2f} s") + if args.height is not None or args.period is not None: + print(f" (Using overrides: Hm0={hm0_overall_plot} m, Tp={tp_overall_plot} s)") + + plt.plot( + overall_avg_spectrum.index.astype(float), + overall_avg_spectrum.values, + label="Overall Average (measured)", + linewidth=3, + linestyle="-", + color="black", + alpha=0.9 + ) + + # PM for overall average (using Tp, not Te) + freq = overall_avg_spectrum.index.astype(float) + pm_overall = pierson_moskowitz_spectrum(freq, hm0_overall_plot, tp_overall_plot) + plt.plot( + freq, + pm_overall, + label="Overall Average (PM)", + linewidth=2, + linestyle="--", + color="black", + alpha=0.5 + ) + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Year-Averaged Measured vs Pierson-Moskowitz Spectra", fontsize=14) + plt.legend(fontsize=10, loc='upper right', ncol=2) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + else: + print("No valid years found to plot") + + +#~~~~~~~~~~~~~~ Plot JONSWAP spectrum if requested ~~~~~~~~~~~~~~ + +if args.plot_spectrum_js is not None: + # If in extract_single or extract_daily mode, plot each extracted spectrum with its own JONSWAP fit + if extract_single or extract_daily: + plt.figure(figsize=(12, 6)) + colors = plt.cm.tab10(np.linspace(0, 1, max(len(extracted_spectra), 1))) + multi = len(extracted_spectra) > 1 + + for idx, (label, spectrum) in enumerate(extracted_spectra): + spectrum_series = pd.Series(spectrum.values, index=spectrum.index) + + # Calculate wave metrics + hm0 = float(wave.resource.significant_wave_height(spectrum_series)) + te = float(wave.resource.energy_period(spectrum_series)) + tp = float(wave.resource.peak_period(spectrum_series)) + + # Use override values if provided + hm0_plot = float(args.height) if args.height is not None else hm0 + tp_plot = float(args.period) if args.period is not None else tp + + # Fit JONSWAP gamma (and refine Tp, since the measured peak bin is noisy) + freq = spectrum.index.astype(float).values + gamma_fit, error, tp_fit, hit_bound = fit_jonswap_gamma( + spectrum_series, freq, tp_plot, hm0_plot) + + measured_label = f"{label} (measured)" if multi else "Measured" + plt.plot( + freq, + spectrum.values, + label=measured_label, + linewidth=2.5, + color=colors[idx], + alpha=0.9 + ) + + # Generate fitted JONSWAP spectrum + js_spectrum = jonswap_spectrum_custom(freq, tp_fit, hm0_plot, gamma=gamma_fit) + + js_label = (f"{label} JONSWAP (γ={gamma_fit:.3f}, Tp={tp_fit:.2f} s)" if multi + else f"JONSWAP fit (γ={gamma_fit:.3f}, Tp={tp_fit:.2f} s)") + plt.plot( + freq, + js_spectrum, + label=js_label, + linewidth=2.0, + color=colors[idx], + alpha=0.75, + linestyle='--' + ) + + print(f"\n{label}") + print(f"Hm0={hm0:.2f} m | Te={te:.2f} s | Tp={tp:.2f} s") + print(f"Fitted JONSWAP: gamma={gamma_fit:.3f}, Tp={tp_fit:.2f} s (measured {tp_plot:.2f} s), Error={error:.6e}") + if hit_bound: + print(f" NOTE: gamma is pinned at the {GAMMA_MIN}-{GAMMA_MAX} bound. The measured") + print(f" spectrum is likely broad or multi-peaked (swell + wind sea), which a") + print(f" single-peak JONSWAP cannot represent. Check Te vs Tp: Te > Tp usually") + print(f" indicates a second low-frequency peak.") + + # Manual gamma overlays, drawn at the measured Hm0/Tp (or the overrides) + if args.gamma: + gamma_colors = ['lime', 'orange', 'purple', 'brown', 'pink', 'gray', 'olive', 'cyan'] + gamma_linestyles = [':', '-.', '--', '-'] + for g_idx, gamma_val in enumerate(args.gamma): + gamma_val = float(gamma_val) + if gamma_val < GAMMA_MIN: + print(f" ERROR: gamma={gamma_val:.3f} < {GAMMA_MIN} - outside valid range (skipping)") + continue + + js_manual = jonswap_spectrum_custom(freq, tp_plot, hm0_plot, gamma=gamma_val) + if np.any(np.isnan(js_manual)) or np.any(np.isinf(js_manual)): + print(f" WARNING: gamma={gamma_val:.3f} produced NaN/inf values") + continue + + error_vs_fitted = float(np.sum((js_spectrum - js_manual) ** 2)) + error_vs_measured = float(np.sum((spectrum.values - js_manual) ** 2)) + manual_label = (f"{label} manual γ={gamma_val:.3f}" if multi + else f"Manual γ={gamma_val:.3f}") + plt.plot( + freq, + js_manual, + label=f"{manual_label} (err={error_vs_measured:.2e})", + linewidth=1.6, + color=gamma_colors[g_idx % len(gamma_colors)], + alpha=0.9, + linestyle=gamma_linestyles[g_idx % len(gamma_linestyles)] + ) + print(f" Manual gamma={gamma_val:.3f}: Error vs measured = {error_vs_measured:.6e}, " + f"vs fitted = {error_vs_fitted:.6e}") + + mode_name = "Measurement" if extract_single else "Daily Average" + if multi: + title = f"{source_title} - {len(extracted_spectra)} {mode_name}s vs Fitted JONSWAP" + else: + title = plot_heading(extracted_spectra[0][0], "Fitted JONSWAP") + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(title, fontsize=14) + plt.legend(fontsize=9, loc='upper right', ncol=2 if multi else 1) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + print(f"\nJONSWAP spectrum plot displayed") + + # Multiple months requested: overlay each month's average with its own JONSWAP fit + elif parse_month_list(plot_spectrum_js_values): + month_nums = parse_month_list(plot_spectrum_js_values) + month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + if not selected_years: + print("No valid years found to plot") + else: + plot_year = selected_years[0] + print(f"\nPlotting spectrum vs fitted JONSWAP for {len(month_nums)} months of {plot_year}...") + + if str(plot_year) not in ndbc_data: + print(f"Year {plot_year} not found in available data") + else: + year_data = ndbc_data[str(plot_year)] + plt.figure(figsize=(12, 6)) + colors = plt.cm.tab10(np.linspace(0, 1, max(len(month_nums), 1))) + plotted_any = False + + for idx, month_num in enumerate(month_nums): + month_spectra = [] + for timestamp, spectrum in year_data.items(): + if timestamp.month == month_num: + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + month_spectra.append(spectrum_clean) + + if not month_spectra: + print(f" No data for {month_names[month_num-1]} {plot_year}") + continue + + month_avg = pd.concat(month_spectra, axis=1).mean(axis=1) + month_avg_series = pd.Series(month_avg.values, index=month_avg.index) + + hm0 = float(wave.resource.significant_wave_height(month_avg_series)) + te = float(wave.resource.energy_period(month_avg_series)) + tp = float(wave.resource.peak_period(month_avg_series)) + + hm0_plot = float(args.height) if args.height is not None else hm0 + tp_plot = float(args.period) if args.period is not None else tp + + freq = month_avg.index.astype(float).values + gamma_fit, error, tp_fit, hit_bound = fit_jonswap_gamma( + month_avg_series, freq, tp_plot, hm0_plot) + + label = month_names[month_num-1] + plt.plot(freq, month_avg.values, label=f"{label} (measured)", + linewidth=2.5, color=colors[idx], alpha=0.9) + + js_spectrum = jonswap_spectrum_custom(freq, tp_fit, hm0_plot, gamma=gamma_fit) + plt.plot(freq, js_spectrum, + label=f"{label} JONSWAP (γ={gamma_fit:.3f})", + linewidth=2.0, color=colors[idx], alpha=0.75, linestyle='--') + + print(f"{label} {plot_year}: Hm0={hm0:6.2f} m | Te={te:6.2f} s | Tp={tp:6.2f} s | " + f"gamma={gamma_fit:.3f} | Tp_fit={tp_fit:6.2f} s | Error={error:.6e}") + plotted_any = True + + if plotted_any: + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Monthly Measured vs Fitted JONSWAP ({plot_year})", + fontsize=14) + plt.legend(fontsize=9, loc='upper right', ncol=2) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + print(f"\nJONSWAP spectrum plot displayed") + else: + plt.close() + print("No month data available to plot") + + # Normal plotting logic (when not in extract mode) + elif isinstance(args.plot_spectrum_js, str): + # Determine if input is a month or year + month_num, is_month = parse_month(args.plot_spectrum_js) + + if is_month: + # Plot specific month for the selected year + month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + if selected_years: + plot_year = selected_years[0] + print(f"\nPlotting spectrum vs fitted JONSWAP for {month_names[month_num-1]} {plot_year}...") + + if str(plot_year) in ndbc_data: + year_data = ndbc_data[str(plot_year)] + + # Filter for just this month + month_spectra = [] + for timestamp, spectrum in year_data.items(): + if timestamp.month == month_num: + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + month_spectra.append(spectrum_clean) + + if month_spectra: + # Average spectra for this month + month_avg = pd.concat(month_spectra, axis=1).mean(axis=1) + month_avg_series = pd.Series(month_avg.values, index=month_avg.index) + + # Calculate wave metrics (ensure they're Python floats, not numpy types) + hm0 = float(wave.resource.significant_wave_height(month_avg_series)) + te = float(wave.resource.energy_period(month_avg_series)) + tp = float(wave.resource.peak_period(month_avg_series)) + + print(f"{month_names[month_num-1]} {plot_year}: Hm0={hm0:.2f} m | Te={te:.2f} s | Tp={tp:.2f} s") + + # Fit JONSWAP gamma to MEASURED spectrum (NOT using overrides) + freq = month_avg.index.astype(float).values + gamma_fit, error, tp_fit, hit_bound = fit_jonswap_gamma(month_avg_series, freq, tp, hm0) + print(f"Fitted JONSWAP: gamma={gamma_fit:.3f} | Tp={tp_fit:.2f} s | Error={error:.6e}") + + # Plot measured vs JONSWAP + plt.figure(figsize=(12, 6)) + + plt.plot( + freq, + month_avg.values, + label=f"{month_names[month_num-1]} (measured)", + linewidth=2.5, + color='steelblue', + alpha=0.9 + ) + + # Generate fitted JONSWAP spectrum using MEASURED frequency grid + js_spectrum = jonswap_spectrum_custom(freq, tp, hm0, gamma=gamma_fit) + + plt.plot( + freq, + js_spectrum, + label=f"{month_names[month_num-1]} (JONSWAP fit, γ={gamma_fit:.3f})", + linewidth=2.8, + color='coral', + alpha=0.9, + linestyle='--' + ) + + # Plot manual gammas if provided (support multiple) + # These use override values if provided, otherwise measured values + if args.gamma: + hm0_manual = float(args.height) if args.height is not None else hm0 + tp_manual = float(args.period) if args.period is not None else tp + print(f" Manual gamma: Using Hm0={hm0_manual:.3f} m, Tp={tp_manual:.3f} s (from {'override' if args.height is not None else 'measured data'})") + gamma_colors = ['lime', 'orange', 'purple', 'brown', 'pink', 'gray', 'olive', 'cyan'] + gamma_linestyles = [':', '-.', '--', '-', ':', '-.', '--', '-'] + print(f"\nManual gamma vs fitted comparison ({month_names[month_num-1]}):") + print(f"DEBUG: args.gamma = {args.gamma}, type = {type(args.gamma)}") + + for idx, gamma_val in enumerate(args.gamma): + print(f"DEBUG: Processing gamma_val = {gamma_val}, type = {type(gamma_val)}") + if isinstance(gamma_val, list): + continue + gamma_val = float(gamma_val) + + # JONSWAP constrained to gamma >= GAMMA_MIN for physical validity + if gamma_val < GAMMA_MIN: + print(f" ERROR: gamma={gamma_val:.3f} < {GAMMA_MIN} - outside valid range (skipping)") + continue + + try: + js_manual = jonswap_spectrum_custom(freq, tp_manual, hm0_manual, gamma=gamma_val) + # Check for NaN or inf values + if np.any(np.isnan(js_manual)) or np.any(np.isinf(js_manual)): + print(f" WARNING: γ={gamma_val:.3f} produced NaN/inf values") + continue + error_vs_fitted = np.sum((js_spectrum - js_manual) ** 2) + error_vs_measured = np.sum((month_avg.values - js_manual) ** 2) + + color = gamma_colors[idx % len(gamma_colors)] + linestyle = gamma_linestyles[idx % len(gamma_linestyles)] + + plt.plot( + freq, + js_manual, + label=f"{month_names[month_num-1]} (manual γ={gamma_val:.3f}, err={error_vs_fitted:.2e})", + linewidth=1.8, + color=color, + alpha=0.95, + linestyle=linestyle + ) + print(f" Manual γ={gamma_val:.3f}: Error vs fitted = {error_vs_fitted:.6e}, Error vs measured = {error_vs_measured:.6e}") + except Exception as e: + print(f" ERROR with γ={gamma_val:.3f}: {e}") + import traceback + traceback.print_exc() + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - {month_names[month_num-1]} {plot_year} Measured vs Fitted JONSWAP Spectrum", fontsize=14) + plt.legend(fontsize=10, loc='upper right', ncol=2) + plt.grid(True, alpha=0.3) + plt.tight_layout() + + # Save to file before showing + plot_filename = f"spectrum_jonswap_{month_names[month_num-1].lower()}_{plot_year}.png" + plt.savefig(plot_filename, dpi=150) + print(f"✓ Plot saved: {plot_filename}") + + plt.show() + plt.close() + else: + print(f"No data found for {month_names[month_num-1]} {plot_year}") + else: + print(f"Year {plot_year} not found in available data") + else: + print("No years found to plot") + else: + # Plot monthly spectra vs fitted JONSWAP for specified year + year_for_js = args.plot_spectrum_js + print(f"\nPlotting spectrum vs fitted JONSWAP for year {year_for_js}...") + + if year_for_js in ndbc_data: + year_data = ndbc_data[year_for_js] + + # Group data by month + monthly_spectra = {} + for timestamp, spectrum in year_data.items(): + month = timestamp.month + if month not in monthly_spectra: + monthly_spectra[month] = [] + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + monthly_spectra[month].append(spectrum_clean) + + # Plot monthly spectra vs JONSWAP + plt.figure(figsize=(14, 8)) + colors = plt.cm.viridis(np.linspace(0, 1, len(monthly_spectra))) + month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + + print(f"\nMonthly Spectrum Comparison for {year_for_js}:") + print("-" * 80) + + for idx, (month, spectra_list) in enumerate(sorted(monthly_spectra.items())): + if len(spectra_list) > 0: + # Average spectra for this month + month_avg = pd.concat(spectra_list, axis=1).mean(axis=1) + month_avg_series = pd.Series(month_avg.values, index=month_avg.index) + + # Calculate wave metrics (ensure they're Python floats, not numpy types) + hm0 = float(wave.resource.significant_wave_height(month_avg_series)) + te = float(wave.resource.energy_period(month_avg_series)) + tp = float(wave.resource.peak_period(month_avg_series)) + + # Use override values if provided + hm0_plot = float(args.height) if args.height is not None else hm0 + tp_plot = float(args.period) if args.period is not None else tp + + # Fit JONSWAP gamma (peak enhancement) + freq = month_avg.index.astype(float).values + gamma_fit, error, tp_fit, hit_bound = fit_jonswap_gamma(month_avg_series, freq, tp_plot, hm0_plot) + + print(f"{month_names[month-1]:>3}: Hm0={hm0:6.2f} m | Te={te:6.2f} s | Tp={tp:6.2f} s | gamma={gamma_fit:.3f} | Tp_fit={tp_fit:6.2f} s | Error={error:.6e}") + if args.height is not None or args.period is not None: + print(f" (Using overrides: Hm0={hm0_plot} m, Tp={tp_plot} s)") + + # Plot measured spectrum + plt.plot( + freq, + month_avg.values, + label=f"{month_names[month-1]} (measured)", + linewidth=2.5, + color=colors[idx], + alpha=0.9 + ) + + # Generate and plot fitted JONSWAP spectrum + js_spectrum = jonswap_spectrum_custom(freq, tp_plot, hm0_plot, gamma=gamma_fit) + + plt.plot( + freq, + js_spectrum, + label=f"{month_names[month-1]} (γ={gamma_fit:.3f})", + linewidth=1.5, + color=colors[idx], + alpha=0.5, + linestyle='--' + ) + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Monthly Measured vs Fitted JONSWAP Spectra ({year_for_js})", fontsize=14) + plt.legend(fontsize=9, loc='upper right', ncol=2) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + else: + print(f"Year {year_for_js} not found in available data") + + elif args.plot_spectrum_js is True: + # Plot year-averaged spectra vs fitted JONSWAP + if selected_years: + plt.figure(figsize=(14, 8)) + colors = plt.cm.viridis(np.linspace(0, 1, len(selected_years))) + + all_year_avg_spectra_js = [] + + print(f"\nYear-Averaged Spectrum Comparison:") + print("-" * 80) + + for idx, year in enumerate(selected_years): + year_key = str(year) + if year_key in ndbc_data: + year_data = ndbc_data[year_key] + + # Collect all spectra for this year + year_spectra = [] + for timestamp, spectrum in year_data.items(): + spectrum_clean = spectrum.dropna() + if len(spectrum_clean) > 0: + year_spectra.append(spectrum_clean) + + if year_spectra: + # Year average + year_avg = pd.concat(year_spectra, axis=1).mean(axis=1) + year_avg_series = pd.Series(year_avg.values, index=year_avg.index) + all_year_avg_spectra_js.append(year_avg_series) + + # Calculate metrics (ensure they're Python floats, not numpy types) + hm0 = float(wave.resource.significant_wave_height(year_avg_series)) + te = float(wave.resource.energy_period(year_avg_series)) + tp = float(wave.resource.peak_period(year_avg_series)) + + # Fit JONSWAP gamma to MEASURED spectrum (NOT using overrides) + freq = year_avg.index.astype(float).values + gamma_fit, error, tp_fit, hit_bound = fit_jonswap_gamma(year_avg_series, freq, tp, hm0) + + print(f"Year {year}: Hm0={hm0:6.2f} m | Te={te:6.2f} s | Tp={tp:6.2f} s | gamma={gamma_fit:.3f} | Tp_fit={tp_fit:6.2f} s | Error={error:.6e}") + + # Plot measured spectrum + plt.plot( + freq, + year_avg.values, + label=f"Year {year} (measured)", + linewidth=2.5, + color=colors[idx], + alpha=0.9 + ) + + # Generate and plot fitted JONSWAP spectrum using MEASURED Hm0/Tp + js_spectrum = jonswap_spectrum_custom(freq, tp, hm0, gamma=gamma_fit) + plt.plot( + freq, + js_spectrum, + label=f"Year {year} (γ={gamma_fit:.3f}, fitted)", + linewidth=2.0, + color=colors[idx], + alpha=0.6, + linestyle='--' + ) + + # Plot manual gammas for each year if provided (support multiple) + # These use override values if provided, otherwise measured values + if args.gamma: + hm0_manual = float(args.height) if args.height is not None else hm0 + tp_manual = float(args.period) if args.period is not None else tp + print(f" Year {year} manual gamma: Using Hm0={hm0_manual:.3f} m, Tp={tp_manual:.3f} s (from {'override' if args.height is not None else 'measured data'})") + gamma_linestyles = ['-', '-.'] + for g_idx, gamma_val in enumerate(args.gamma): + gamma_val = float(gamma_val) + + # JONSWAP constrained to gamma >= GAMMA_MIN for physical validity + if gamma_val < GAMMA_MIN: + print(f" WARNING: gamma={gamma_val:.3f} < {GAMMA_MIN} - outside valid range (skipping)") + continue + + js_manual = jonswap_spectrum_custom(freq, tp_manual, hm0_manual, gamma=gamma_val) + # Calculate error vs fitted spectrum + error_vs_fitted = np.sum((js_spectrum - js_manual) ** 2) + linestyle = gamma_linestyles[g_idx % len(gamma_linestyles)] + plt.plot( + freq, + js_manual, + label=f"Year {year} (manual γ={gamma_val:.3f}, err={error_vs_fitted:.2e})", + linewidth=1.2, + color=colors[idx], + alpha=0.7, + linestyle=linestyle + ) + print(f" Manual γ={gamma_val:.3f}: Error vs fitted = {error_vs_fitted:.6e}") + + if len(all_year_avg_spectra_js) > 1: + overall_avg_spectrum_js = pd.concat(all_year_avg_spectra_js, axis=1).mean(axis=1) + overall_avg_series_js = pd.Series(overall_avg_spectrum_js.values, index=overall_avg_spectrum_js.index) + + hm0_overall = float(wave.resource.significant_wave_height(overall_avg_series_js)) + te_overall = float(wave.resource.energy_period(overall_avg_series_js)) + tp_overall = float(wave.resource.peak_period(overall_avg_series_js)) + + # Fit to MEASURED overall spectrum (NOT using overrides) + freq = overall_avg_spectrum_js.index.astype(float).values + gamma_fit, error, tp_fit, hit_bound = fit_jonswap_gamma(overall_avg_series_js, freq, tp_overall, hm0_overall) + + print(f"Overall: Hm0={hm0_overall:6.2f} m | Te={te_overall:6.2f} s | Tp={tp_overall:6.2f} s | gamma={gamma_fit:.3f} | Tp_fit={tp_fit:6.2f} s | Error={error:.6e}") + + plt.plot( + freq, + overall_avg_spectrum_js.values, + label="Overall Average (measured)", + linewidth=3, + linestyle="-", + color="black", + alpha=0.9 + ) + + # JONSWAP for overall average using MEASURED Hm0/Tp + js_overall = jonswap_spectrum_custom(freq, tp_overall, hm0_overall, gamma=gamma_fit) + + # Check peak frequency of fitted spectrum + fp_overall = 1.0 / tp_overall + peak_idx_fit = np.argmax(js_overall) + peak_freq_fit = freq[peak_idx_fit] + peak_height_fit = js_overall[peak_idx_fit] + + plt.plot( + freq, + js_overall, + label=f"Overall Average (γ={gamma_fit:.3f}, fitted)", + linewidth=2, + linestyle="--", + color="black", + alpha=0.5 + ) + print(f" Fitted γ={gamma_fit:.3f}: Peak at f={peak_freq_fit:.6f} Hz (expect {fp_overall:.6f}), height={peak_height_fit:.4f} m²/Hz") + + # Plot manual gammas if provided (support multiple) + # These use override values if provided, otherwise measured values + if args.gamma: + hm0_manual = float(args.height) if args.height is not None else hm0_overall + tp_manual = float(args.period) if args.period is not None else tp_overall + gamma_colors = ['green', 'blue', 'red', 'purple', 'orange', 'brown', 'pink', 'gray'] + gamma_linestyles = [':', '-.', '--', '-', ':', '-.', '--', '-'] + fp_expected = 1.0 / tp_manual + print(f"\nManual gamma vs fitted comparison (overall):") + print(f" Using Hm0={hm0_manual:.3f} m, Tp={tp_manual:.3f} s (from {'override' if args.height is not None else 'measured data'})") + print(f" Expected peak frequency: fp = 1/{tp_manual:.3f} = {fp_expected:.6f} Hz") + print(f"DEBUG: args.gamma = {args.gamma}, type = {type(args.gamma)}") + for idx, gamma_val in enumerate(args.gamma): + print(f"DEBUG: Processing gamma_val = {gamma_val}, type = {type(gamma_val)}") + if isinstance(gamma_val, list): + continue + gamma_val = float(gamma_val) + + # JONSWAP constrained to gamma >= GAMMA_MIN for physical validity + if gamma_val < GAMMA_MIN: + print(f" ERROR: gamma={gamma_val:.3f} < {GAMMA_MIN} - outside valid range (skipping)") + continue + + try: + js_manual = jonswap_spectrum_custom(freq, tp_manual, hm0_manual, gamma=gamma_val) + # Check for NaN or inf values + if np.any(np.isnan(js_manual)) or np.any(np.isinf(js_manual)): + print(f" WARNING: γ={gamma_val:.3f} produced NaN/inf values") + continue + error_vs_fitted = np.sum((js_overall - js_manual) ** 2) + + # Find the peak frequency of this spectrum + peak_idx = np.argmax(js_manual) + peak_freq = freq[peak_idx] + peak_height = js_manual[peak_idx] + + color = gamma_colors[idx % len(gamma_colors)] + linestyle = gamma_linestyles[idx % len(gamma_linestyles)] + plt.plot( + freq, + js_manual, + label=f"Overall Average (manual γ={gamma_val:.3f}, err={error_vs_fitted:.2e})", + linewidth=2, + linestyle=linestyle, + color=color, + alpha=0.5 + ) + print(f" γ={gamma_val:.3f}: Peak at f={peak_freq:.6f} Hz (expect {fp_expected:.6f}), height={peak_height:.4f} m²/Hz, error={error_vs_fitted:.6e}") + print(f" Manual γ={gamma_val:.3f}: Error vs fitted = {error_vs_fitted:.6e}") + except Exception as e: + print(f" ERROR with γ={gamma_val:.3f}: {e}") + import traceback + traceback.print_exc() + + plt.xlabel("Frequency [Hz]", fontsize=12) + plt.ylabel("Spectral Density [m²/Hz]", fontsize=12) + plt.title(f"NDBC Buoy {buoy_number} - Year-Averaged Measured vs Fitted JONSWAP Spectra", fontsize=14) + plt.legend(fontsize=10, loc='upper right', ncol=2) + plt.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + plt.close() + else: + print("No valid years found to plot") + + +#~~~~~~~~~~~~~~ Plot monthly wave power statistics if requested ~~~~~~~~~~~~~~ + +def select_records(flag_values): + """Raw swden records matching the current --year / --date / --time selection. + + A month list or year passed to the flag itself overrides the global year choice. + Returns (DataFrame, label) on success or (None, reason). + """ + months = parse_month_list(flag_values, min_count=1) + year_arg = None + if months is None and len(flag_values) == 1: + try: + candidate = int(flag_values[0]) + if candidate > 12: + year_arg = candidate + except ValueError: + pass + + years = [year_arg] if year_arg else list(selected_years) + frames = [ndbc_data[str(y)] for y in years if str(y) in ndbc_data] + if not frames: + return None, f"no data found for {years}" + + raw = pd.concat(frames, axis=1, sort=False) + + def keep(ts): + if months and ts.month not in months: + return False + if not extract_targets: + return True + for t in extract_targets: + if t['kind'] == 'range': + if t['timestamp'] <= ts <= t['end']: + return True + elif t['kind'] == 'daily': + if ts.date() == t['timestamp'].date(): + return True + elif abs((ts - t['timestamp']).total_seconds()) <= 1800: + return True + return False + + cols = sorted([c for c in raw.columns if keep(c)]) + if len(cols) < 2: + return None, f"only {len(cols)} record(s) in the selection" + + names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', + 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + if extract_targets: + labels = [t['label'] for t in extract_targets] + label = labels[0] if len(labels) == 1 else f"{labels[0]} ... {labels[-1]}" + elif months: + label = f"{', '.join(names[m - 1] for m in months)} {years[0]}" + else: + label = ', '.join(map(str, years)) + return raw[cols], label + + +if args.plot_wavestats is not None: + if spectrum_from_file: + print("\n--plot_wavestats needs buoy records and is not available with --spectrum_file") + else: + window, sel_label = select_records(plot_wavestats_values) + if window is None: + print(f"\nCannot plot wave statistics: {sel_label}") + else: + idx = pd.DatetimeIndex(window.columns) + + hm0_w = np.asarray(wave.resource.significant_wave_height(window)).ravel() + tz_w = np.asarray(wave.resource.average_zero_crossing_period(window)).ravel() + metrics = { + 'Hm0': hm0_w, + 'Te': np.asarray(wave.resource.energy_period(window)).ravel(), + 'Tp': np.asarray(wave.resource.peak_period(window)).ravel(), + } + if water_depth: + metrics['J'] = np.asarray( + wave.resource.energy_flux(window, water_depth)).ravel() / 1000.0 + # Wave steepness Sm = Hm0 / ((g/2pi) * Tz^2) + metrics['Sm'] = hm0_w / ((9.81 / (2 * np.pi)) * (tz_w ** 2)) + + data_clean = pd.DataFrame(metrics, index=idx).dropna() + + if len(data_clean) < 2: + print("\nNot enough valid data for wave statistics") + else: + clean_idx = data_clean.index + span = clean_idx[-1] - clean_idx[0] + + # Resolution follows the span of the selection + month_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', + 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + if span <= pd.Timedelta(days=2): + group_key, xlabel, res = clean_idx.hour, "Hour of Day (UTC)", "Hourly" + elif span <= pd.Timedelta(days=62): + group_key, xlabel, res = clean_idx.normalize(), "Date", "Daily" + else: + group_key, xlabel, res = clean_idx.month, "Month", "Monthly" + + print(f"\nGenerating {res.lower()} wave statistics for {sel_label} " + f"({len(data_clean)} records)...") + + data_group = data_clean.groupby(group_key) + medians = data_group.median() + quartiles = data_group.describe() + + QoIs = list(data_clean.columns) + fig, axs = plt.subplots(len(QoIs), 1, figsize=(10, 10), sharex=True) + if len(QoIs) == 1: + axs = [axs] + + units = {'Hm0': 'm', 'Te': 's', 'Tp': 's', 'J': 'kW/m', 'Sm': '-'} + for i, QoI in enumerate(QoIs): + axs[i].plot(medians.index, medians[QoI], marker="o", linewidth=2, + markersize=6, color="steelblue") + axs[i].fill_between( + medians.index, + quartiles[QoI, "25%"], + quartiles[QoI, "75%"], + alpha=0.3, + color="steelblue", + label="25%-75% IQR" + ) + axs[i].grid(True, alpha=0.3) + axs[i].set_ylabel(f"{QoI} [{units.get(QoI, '-')}]", fontsize=11) + + mx, mn = medians[QoI].max(), medians[QoI].min() + mx_at, mn_at = medians[QoI].idxmax(), medians[QoI].idxmin() + if res == "Monthly": + mx_at, mn_at = month_labels[mx_at - 1], month_labels[mn_at - 1] + elif res == "Daily": + mx_at, mn_at = f"{mx_at:%Y-%m-%d}", f"{mn_at:%Y-%m-%d}" + else: + mx_at, mn_at = f"{mx_at:02d}:00", f"{mn_at:02d}:00" + print(f"\n{QoI}:") + print(f" Max: {mx:.4f} at {mx_at}") + print(f" Min: {mn:.4f} at {mn_at}") + + if i == 0: + axs[i].set_title( + f"{source_title} - {res} Wave Statistics ({sel_label})", + fontsize=13) + if i == len(QoIs) - 1: + axs[i].set_xlabel(xlabel, fontsize=11) + + if res == "Monthly": + axs[-1].set_xticks(range(1, 13)) + axs[-1].set_xticklabels(month_labels) + elif res == "Hourly": + axs[-1].set_xticks(range(0, 24, 2)) + else: + fig.autofmt_xdate() + + plt.tight_layout() + plt.show() + plt.close() + + print(f"\n{res} wave statistics plot displayed") + + +#~~~~~~~~~~~~~~ Plot wind data if requested ~~~~~~~~~~~~~~ + +if args.plot_wind is not None: + if len(wind_data) == 0: + print("\nNo wind data available to plot") + else: + # A month or year given to the flag overrides the global year selection + wind_months = parse_month_list(plot_wind_values, min_count=1) + wind_year = None + if wind_months is None and len(plot_wind_values) == 1: + try: + candidate = int(plot_wind_values[0]) + if candidate > 12: + wind_year = candidate + except ValueError: + pass + + wind_years = [wind_year] if wind_year else list(selected_years) + frames = [wind_data[str(y)].T for y in wind_years if str(y) in wind_data] + + if not frames: + print(f"\nNo wind data found for {wind_years}") + else: + wind_ts = pd.concat(frames, axis=0).sort_index() + + # Mask NDBC missing-value sentinels + for col, limit in (('WSPD', 90), ('GST', 90), ('WDIR', 400)): + if col in wind_ts.columns: + wind_ts[col] = wind_ts[col].where(wind_ts[col] < limit) + + if wind_months: + wind_ts = wind_ts[wind_ts.index.month.isin(wind_months)] + + _names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', + 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + wind_label = (f"{', '.join(_names[m - 1] for m in wind_months)} {wind_years[0]}" + if wind_months else ', '.join(map(str, wind_years))) + + series_specs = [ + ('WSPD', 'Wind Speed [m/s]', 'steelblue'), + ('GST', 'Wind Gusts [m/s]', 'darkorange'), + ('WDIR', 'Wind Direction [deg]', 'seagreen'), + ] + series_specs = [s for s in series_specs if s[0] in wind_ts.columns] + + if extract_single or extract_daily: + # Daily view: raw 10-minute interval time series for the requested day + target_day = (single_measurement_date if extract_single else daily_start).date() + day_ts = wind_ts[wind_ts.index.date == target_day] + + if len(day_ts) == 0: + print(f"\nNo wind data found for {target_day}") + else: + print(f"\nGenerating daily wind plot for {target_day} ({len(day_ts)} 10-minute records)...") + fig, axs = plt.subplots(len(series_specs), 1, figsize=(12, 8), sharex=True) + if len(series_specs) == 1: + axs = [axs] + + for i, (col, label, color) in enumerate(series_specs): + axs[i].plot(day_ts.index, day_ts[col], marker='o', markersize=3, + linewidth=1.2, color=color) + axs[i].set_ylabel(label, fontsize=11) + axs[i].grid(True, which='major', alpha=0.3) + axs[i].grid(True, which='minor', axis='x', alpha=0.15, linestyle=':') + if col == 'WDIR': + axs[i].set_ylim(0, 360) + axs[i].set_yticks([0, 90, 180, 270, 360]) + if extract_single: + axs[i].axvline(single_measurement_date, color='red', + linestyle='--', linewidth=1, alpha=0.7) + + axs[0].set_title( + f"NDBC Buoy {buoy_number} - Wind at 10-Minute Intervals ({target_day})", + fontsize=13) + axs[-1].set_xlabel("Hour of Day (UTC)", fontsize=11) + + # Tick every 2 hours (minor ticks hourly) across the full 24-hour day + day_start = pd.Timestamp(target_day) + axs[-1].set_xlim(day_start, day_start + pd.Timedelta(days=1)) + axs[-1].xaxis.set_major_locator(mdates.HourLocator(interval=2)) + axs[-1].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M')) + axs[-1].xaxis.set_minor_locator(mdates.HourLocator(interval=1)) + plt.setp(axs[-1].get_xticklabels(), rotation=0, ha='center') + + plt.tight_layout() + plt.show() + plt.close() + print("Daily wind plot displayed") + else: + # Median + IQR at a resolution that follows the selection + idx = wind_ts.index + span = idx[-1] - idx[0] + month_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', + 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + if span <= pd.Timedelta(days=62): + group_key, xlabel, res = idx.normalize(), "Date", "Daily" + else: + group_key, xlabel, res = idx.month, "Month", "Monthly" + + print(f"\nGenerating {res.lower()} wind statistics for {wind_label} " + f"({len(wind_ts)} records)...") + + grouped = wind_ts.groupby(group_key) + medians = grouped.median() + quartiles = grouped.describe() + + fig, axs = plt.subplots(len(series_specs), 1, figsize=(12, 10), sharex=True) + if len(series_specs) == 1: + axs = [axs] + + for i, (col, label, color) in enumerate(series_specs): + axs[i].plot(medians.index, medians[col], marker="o", linewidth=2, + markersize=5, color=color, label="Median") + axs[i].fill_between( + medians.index, + quartiles[col, "25%"], + quartiles[col, "75%"], + alpha=0.3, + color=color, + label="25%-75% IQR" + ) + axs[i].set_ylabel(label, fontsize=11) + axs[i].grid(True, alpha=0.3) + if col == 'WDIR': + axs[i].set_ylim(0, 360) + axs[i].set_yticks([0, 90, 180, 270, 360]) + if i == 0: + axs[i].legend(fontsize=9, loc='upper right') + + mx, mn = medians[col].max(), medians[col].min() + mx_at, mn_at = medians[col].idxmax(), medians[col].idxmin() + if res == "Monthly": + mx_at, mn_at = month_labels[mx_at - 1], month_labels[mn_at - 1] + else: + mx_at, mn_at = f"{mx_at:%Y-%m-%d}", f"{mn_at:%Y-%m-%d}" + print(f"\n{col}:") + print(f" Max median: {mx:.2f} at {mx_at}") + print(f" Min median: {mn:.2f} at {mn_at}") + print(f" Overall mean/std: {wind_ts[col].mean():.2f} / {wind_ts[col].std():.2f}") + + axs[0].set_title( + f"{source_title} - {res} Wind Distribution ({wind_label})", fontsize=13) + axs[-1].set_xlabel(xlabel, fontsize=11) + if res == "Monthly": + axs[-1].set_xticks(range(1, 13)) + axs[-1].set_xticklabels(month_labels) + axs[-1].set_xlim(0.4, 12.6) + else: + fig.autofmt_xdate() + plt.tight_layout() + plt.show() + plt.close() + print(f"\n{res} wind statistics plot displayed") + + +#~~~~~~~~~~~~~~ Plot directional wave spectrum if requested ~~~~~~~~~~~~~~ + +if args.plot_wavedirection is not None: + if selected_years: + # Use the first selected year for directional spectrum + dir_year = selected_years[0] + plot_type = args.plot_wavedirection # elevation, energy, or spread + print(f"\nFetching directional wave data for year {dir_year}...") + + try: + # Check if directional data is available + dir_available = ndbc.available_data("swdir", buoy_number) + if len(dir_available) == 0: + print(f"No directional wave data available for buoy {buoy_number}") + else: + # Request directional data for the selected year + dir_data_all = ndbc.request_directional_data(buoy_number, dir_year) + + if len(dir_data_all) == 0: + print(f"No directional data found for year {dir_year}") + else: + # Drop NaN values before selecting or averaging + dir_data_clean = dir_data_all.dropna(dim='date', how='all') + + # --date / --time narrow the average to those records; the swdir + # archive is yearly, so the whole year is fetched either way. + period_label = f"Year-Averaged {dir_year}" + if extract_targets and len(dir_data_clean.date) > 0: + dir_dates = pd.DatetimeIndex(dir_data_clean.date.values) + labels = [t['label'] for t in extract_targets] + if extract_single: + mask = np.zeros(len(dir_dates), dtype=bool) + for target in extract_targets: + if target['kind'] == 'range': + mask |= ((dir_dates >= target['timestamp']) & + (dir_dates <= target['end'])) + continue + ts = target['timestamp'] + nearest = dir_dates[np.argmin(np.abs(dir_dates - ts))] + diff = abs((nearest - ts).total_seconds()) + if diff < 3600: + mask |= (dir_dates == nearest) + print(f" Directional record at {nearest:%Y-%m-%d %H:%M} " + f"({diff/60:.0f} min from requested {ts:%Y-%m-%d %H:%M})") + else: + print(f" WARNING: no directional record within 1 h of " + f"{ts:%Y-%m-%d %H:%M}") + else: + target_days = {t['timestamp'].normalize() for t in extract_targets} + mask = dir_dates.normalize().isin(target_days) + + if mask.any(): + dir_data_clean = dir_data_clean.isel(date=np.flatnonzero(mask)) + period_label = labels[0] if len(labels) == 1 else \ + f"{len(labels)} selections ({labels[0]} ... {labels[-1]})" + print(f"Averaging {int(mask.sum())} directional records for {period_label}...") + else: + print(f"WARNING: no directional data for the requested date(s); " + f"falling back to the {dir_year} average.") + else: + print(f"Averaging directional spectrum across year {dir_year}...") + + if len(dir_data_clean.date) == 0: + print(f"No valid directional data for year {dir_year}") + else: + # Average the spectral components across all times + dir_data_mean = dir_data_clean.mean(dim='date') + + # Create directional spectrum from averaged data (2 degree resolution) + directions = np.arange(0, 360, 2.0) + raw_spectrum = ndbc.create_directional_spectrum(dir_data_mean, directions) + + if plot_type == "elevation": + # Standard elevation spectrum (2D wave spectrum S(f,θ)) [m²/Hz/deg] + directional_spectrum = raw_spectrum + plot_title = f"{period_label} Elevation Spectrum" + colorbar_label = "m²/Hz/deg" + elif plot_type == "energy": + # Energy spectrum - multiply by frequency to get energy distribution [J/Hz/deg] + directional_spectrum = raw_spectrum.copy() + freq_values = raw_spectrum.frequency.values + # Convert to energy by multiplying by frequency + directional_spectrum.values = raw_spectrum.values * freq_values[:, np.newaxis] + plot_title = f"{period_label} Energy Spectrum" + colorbar_label = "J/Hz/deg" + else: # spread + # Spreading function - normalize each frequency to sum to 1 [1/Hz/deg] + directional_spectrum = raw_spectrum.copy() + spectrum_values = raw_spectrum.values + for freq_idx in range(spectrum_values.shape[0]): + freq_sum = spectrum_values[freq_idx, :].sum() + if freq_sum > 0: + directional_spectrum.values[freq_idx, :] = spectrum_values[freq_idx, :] / freq_sum + plot_title = f"{period_label} Spreading Function" + colorbar_label = "1/Hz/deg" + + # Plot directional spectrum + print(f"Plotting {plot_type} spectrum for {period_label}...") + wave.graphics.plot_directional_spectrum(directional_spectrum) + + # Set nautical convention: 0° north (up), 90° east (right) + ax = plt.gca() + ax.set_theta_zero_location('N') + ax.set_theta_direction(-1) # Clockwise + + # Add colorbar label with units + cbar = ax.collections[0].colorbar if hasattr(ax.collections[0], 'colorbar') and ax.collections[0].colorbar else None + if cbar is None: + # Try to find colorbar in figure + for cbar_obj in ax.figure.axes: + if hasattr(cbar_obj, 'yaxis') and cbar_obj.yaxis.get_label_text() == '': + cbar_obj.set_ylabel(colorbar_label) + else: + cbar.set_label(colorbar_label) + + # Apply nofill if requested + if args.nofill: + # Remove filled contours, keep only contour lines with color mapping + ax.clear() + # Re-plot with contour lines only (no fill), but with colors + theta = np.deg2rad(directions) + freq = np.arange(directional_spectrum.values.shape[0]) * 0.01 # Approximate frequency bins + theta_mesh, freq_mesh = np.meshgrid(theta, freq) + # Create contour lines with colormap - no fill + contours = ax.contour(theta_mesh, freq_mesh, directional_spectrum.values, levels=12, cmap='viridis', linewidths=1.5) + ax.clabel(contours, inline=False, fontsize=0) # Hide labels but keep contours + ax.set_theta_zero_location('N') + ax.set_theta_direction(-1) + ax.set_title(plot_title) + else: + plt.title(plot_title) + + plt.show() + plt.close() + print(f"Directional spectrum plot displayed") + except Exception as e: + print(f"Error plotting directional spectrum: {e}") + else: + print("No years found for directional spectrum") + + +#~~~~~~~~~~~~~~ Plot environmental contour heatmap if requested ~~~~~~~~~~~~~~ + +if args.plot_heatmap is not None: + if spectrum_from_file: + print("\n--plot_heatmap needs buoy records and is not available with --spectrum_file") + else: + window, sel_label = select_records(plot_heatmap_values) + if window is None: + print(f"\nCannot plot heatmap: {sel_label}") + else: + df_metrics = pd.DataFrame({ + 'Hm0': np.asarray(wave.resource.significant_wave_height(window)).ravel(), + 'Te': np.asarray(wave.resource.energy_period(window)).ravel(), + }).dropna() + + if len(df_metrics) < 2: + print("\nNot enough valid data for the heatmap") + else: + print(f"\nGenerating sea state heatmap for {sel_label} " + f"({len(df_metrics)} records)...") + + # Finer bins when the selection is small, so a day is not one blob + if len(df_metrics) < 200: + Hm0_bin_size, Te_bin_size = 0.25, 0.5 + else: + Hm0_bin_size, Te_bin_size = 0.5, 1.0 + + # Bin over the observed range so empty rows/columns are not drawn + def _edges(series, step): + lo = max(0.0, np.floor(float(series.min()) / step) * step) + hi = np.ceil(float(series.max()) / step) * step + return np.arange(lo, hi + step, step) + + Hm0_bins = _edges(df_metrics['Hm0'], Hm0_bin_size) + Te_bins = _edges(df_metrics['Te'], Te_bin_size) + hist, Hm0_edges, Te_edges = np.histogram2d( + df_metrics['Hm0'].values, + df_metrics['Te'].values, + bins=[Hm0_bins, Te_bins] + ) + + fig, ax = plt.subplots(figsize=(12, 8)) + im = ax.pcolormesh( + Te_edges, + Hm0_edges, + hist, + shading='auto', + cmap='YlOrRd' + ) + ax.set_xlabel("Energy Period, Te [s]", fontsize=12) + ax.set_ylabel("Significant Wave Height, Hm0 [m]", fontsize=12) + ax.set_title( + f"{source_title} - Sea State Frequency Distribution ({sel_label})", + fontsize=14 + ) + ax.set_xlim(Te_edges[0], Te_edges[-1]) + ax.set_ylim(Hm0_edges[0], Hm0_edges[-1]) + cbar = plt.colorbar(im, ax=ax) + cbar.set_label("Count (Occurrences)", fontsize=11) + plt.tight_layout() + plt.show() + plt.close() + print("Sea state heatmap displayed") From fe21d5cc40e43a4eefdd09332ed311a538a536e2 Mon Sep 17 00:00:00 2001 From: Allie Date: Thu, 10 Sep 2026 13:39:10 -0600 Subject: [PATCH 5/9] Rename MHKiTWaveYamlWriter.py to MHKiT_writer.py --- data/demos/features/{MHKiTWaveYamlWriter.py => MHKiT_writer.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename data/demos/features/{MHKiTWaveYamlWriter.py => MHKiT_writer.py} (100%) diff --git a/data/demos/features/MHKiTWaveYamlWriter.py b/data/demos/features/MHKiT_writer.py similarity index 100% rename from data/demos/features/MHKiTWaveYamlWriter.py rename to data/demos/features/MHKiT_writer.py From 747cad41651f59c08fce46b9453f8f22a7d2edbe Mon Sep 17 00:00:00 2001 From: Allie Date: Thu, 10 Sep 2026 13:42:30 -0600 Subject: [PATCH 6/9] Create MHKiT_writer_README.md with usage instructions Add MHKiT Wave YAML Writer documentation detailing features, command structure, and usage examples. --- data/demos/features/MHKiT_writer_README.md | 451 +++++++++++++++++++++ 1 file changed, 451 insertions(+) create mode 100644 data/demos/features/MHKiT_writer_README.md diff --git a/data/demos/features/MHKiT_writer_README.md b/data/demos/features/MHKiT_writer_README.md new file mode 100644 index 0000000..26685d1 --- /dev/null +++ b/data/demos/features/MHKiT_writer_README.md @@ -0,0 +1,451 @@ +# MHKiT Wave YAML Writer + +A Python tool for analyzing real ocean wave data and writing it into SEA-Stack YAML input files. + +It pulls spectral wave density (`swden`) records from NDBC buoys via MHKiT, computes sea-state statistics, and writes the result into device `*.hydro.yaml` file + +It supports: +- Regular waves (height + period from measured data) +- Irregular waves with a parametric spectrum (Pierson-Moskowitz or JONSWAP with a fitting gamma feature) +- Irregular waves from a measured spectrum, written as a surface elevation (eta) file +- Custom spectra read from a local text file + +With capabilities such as +- Direct data access to the various NDBC buoys +- Selecting data by year, month, date, or an individual 30-minute measurement +- Comparing measured spectra against PM and fitted-JONSWAP estimates +- Override capabilities of wave height, period, peak enhancement, water depth, and simulation parameters +- Plotting-only for wave characteristic analysis such as: + - environmental contours + - joint distributions + - monthly wave statistics + - directional spectra + - wind data + +--------- + +## Basic command structure: + +run the python file in terminal +- located in data/demos/run_seastack/features + - either cd data/demos/run_seastack/features or direct to data/demos/run_seastack/features/MHKiT_writer.py + +```bash +python MHKiT_writer.py --buoy +``` + +The **positional path** is the hydro YAML to update. It must already exist; this script edits an existing `hydrodynamics.waves` block, it does not create one. + +--- + +## How it works + +``` + NDBC buoy ──┐ + ├──► select years / months / date / measurement ──► data analysis plotting + spectrum file┘ │ + ▼ + Hm0 , Te , Tp , water depth + │ + ┌──────────────────────┼──────────────────────┐ + ▼ ▼ ▼ + regular waves parametric spectrum custom spectrum + height + period pm / jonswap eta file + nfrequencies + │ │ │ + └──────────────────────┼──────────────────────┘ + ▼ + .hydro.yaml (waves block) + .simulation.yaml (end_time) +``` + +When plot flags are called (without yaml definitions like `--type`, `--spectrum`, or `--depth`), then the yaml files will not be updated and the code is *plot-only* for data analysis. + +--- + +### `--buoy` +NDBC buoy number to download spectral wave data from, can be found from website: https://www.ndbc.noaa.gov/ + +Required unless `--spectrum_file` or an existing `--elevation_file` is supplied. + +**Note:** not every buoy has a historical `swden` archive. Partner-operated buoys (CDIP/Scripps, e.g. 46258) show live data on the NDBC website but have no spectral archive to download. Some buoys does not have continuous year-long data. + +--- + +### `--year` +Select which year(s) of buoy data to use (default = latest available) for longer time-scale data analysis + +- `--year 2024` — a single year +- `--year 2022 2023 2024` — several years averaged +- `--year from:2020` — that year and everything after +- omit for the latest available year + +--- + +### `--date` / `--time` +Narrow the selection to specific days or individual measurements. + +- `--date 01-01-2024` → daily average +- `--date 01-01 01-02 01-03` — list of dates with `--year 2024` +- `--date` and `--time` → single raw measurement at date/time pair +- `--time 05:10 05:40` — list of times +- `--time 14:00-17:00` — time range +- multiple list values are overlaid on one plot + +NDBC samples at `:10` and `:40` past the hour; the closest measurement within an hour is used and reported. + +--- + +## Data analysis features + + +### `--plot_spectrum` +Plot the measured wave spectrum + +- `--plot_spectrum` → with `--date` and/or `--time` plots corresponding spectrum +- `--plot_spectrum year` → every month of that year overlaid +- `--plot_spectrum month` → individual or list months with `--year` + +**Year breakdown example:** +```bash +python MHKiT_writer.py -buoy 46050 --plot_spectrum 2020 +``` + +image + + +**Multiple years example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_spectrum --year from:2010 +``` + +image + + +**List of months example:** +```bash +python MHKiT_writer.py -buoy 46050 --plot_spectrum jan feb march apr may --year 2020 +``` + +image + + +**Exact date/times example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_spectrum --date 01-01-2020 --time 05:00 06:00 07:00 +``` + +image + + +--- + +### `--plot_spectrum_pm` +Plot the measured spectrum against a Pierson-Moskowitz estimate at the same Hm0 and Tp, and print the sum of squared residuals fit error. + +**Aliases:** `--plot_spectrum_PM` + +Accepts the same bare / year / months / `--date` forms as `--plot_spectrum`. + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_spectrum_pm 2020 +``` + +image + +--- + +### `--plot_spectrum_js` +Plot the measured spectrum against a **fitted** JONSWAP spectrum. Gamma and Tp are least-squares fitted to the measurement (Hs held fixed so total energy matches). + +**Aliases:** `--plot_spectrum_JS`, `--plot_spectrum_jonswap`, `--plot_spectrum_JONSWAP` + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_spectrum_js --date 01-01 05-01 --time 00:00-12:00 12:00-24:00 --year 2024 +``` + +image + +Technically JONSWAP spectrum is only valid for gamma > 0.1, so if the printed gamma pins below, then the sea is broad or double-peaked (swell + wind sea) and a single-peak JONSWAP cannot represent it. + +--- + +### `--gamma` +Overlay manual gamma values alongside the fitted JONSWAP, with the error against both the fit and the measurement. + +**Used with:** `--plot_spectrum_js` or to overwrite yaml + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_spectrum_js --date 01-03-2022 --gamma 1.0 1.1 1.2 +``` + +image + + + +--- + +### `--plot_wavestats` +Median + interquartile range of Hm0, Te, Tp, energy flux J and steepness Sm, with the max/min month printed for each. + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wavestats --year 2020 +``` + +image + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wavestats jan --year 2020 +``` + +image + + +--- + +### `--plot_heatmap` +Sea-state frequency distribution — a 2D histogram of Hm0 against Te over the selected time. + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_heatmap --year 2020 +``` + +image + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_heatmap jan --year 2024 +``` + +image + + +--- + +### `--plot_wavedirection` +Directional wave spectrum on a polar plot, nautical convention (0° north, clockwise). + +Options: +- elevation (default) +- energy +- spread + +Add `--nofill` for contour lines only. + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wavedirection --year 2020 +``` + +image + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wavedirection spread --date 01-01-2020 01-02-2020 01-02-2020 --nofill +``` + +image + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wavedirection spread --year 2024 --nofill +``` + +image + + +--- + +### `--plot_wind` +Wind speed, direction and gusts from the `stdmet` record. + +- with `--date` → 10-minute interval series for that day +- without → monthly box-and-whisker with mean ± 1 std + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wind --date 01-01-2020 +``` + +image + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wind jan --year 2024 +``` + +image + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wind jan 2023 +``` + +image + + +--- + +## Writing to Yaml input files + +### Path to yaml file name +This path tells the script where the yaml file of the device is. It will only edit the waves information in the hydro.yaml file, and sometimes the simulation.yaml if necessary. + +```bash +python MHKiT_writer.py +``` + +### Override flags +To force a value to write to the yaml, instead of using the measured one. + +- `--height` — significant wave height [m] +- `--period` — wave period [s] +- `--depth` — water depth [m] (otherwise taken from buoy metadata) +- `--seed` — random seed for eta generation (default 42) +- `--ramp_time` — ramp duration [s], cosine ramp (default 60) +- `--gamma` — peak enhancement factor, overrides for plotting and yaml writing + + +**Override height, period, and depth:** +```bash +--height 3.0 --period 9.0 --depth 100 +``` + + +### `--type` +Wave type to define YAML hydro (default = `regular`) + +- regular +- irregular + +**Write year averaged regular waves:** +```bash +python MHKiT_writer.py ../5sa/regular_waves/5sa_regular.hydro.yaml --buoy 46050 --year 2024 +``` + +--- + +### `--spectrum` +Spectrum type for irregular waves + +Options: +- `pm` / `pierson_moskowitz` — parametric PM, uses Tp +- `js` / `jonswap` — parametric JONSWAP, gamma fitted to the measured spectrum associated with data +- `custom` — build a surface elevation record from the measured spectrum, or spectrum file, and write an eta file + +**Used with:** `--type irregular` + +**Jonswap fit and write:** +```bash +python MHKiT_writer.py <../path_to_hydro_yaml/case.hydro.yaml> --type irregular --spectrum jonswap --date 01-01-2020 --gamma 1.2 +``` +This pulls spectrum data for the day 01-01-2020, fits a jonswap function with height, period, and gamma, overrides the gamma with 1.2, then writes this all to the input hydro.yaml file to generate this type of irregular wave. + +--- + +### `--spectrum_file` +Read the spectrum from a local text file instead of downloading it. **`--buoy` is not needed.** + +Implies `--type irregular --spectrum custom`. + +Two layouts are accepted — blank lines and `#` / `%` comments are ignored. Units are Hz and m²/Hz. + +**Two rows** (frequencies, then spectral densities): +``` +.0200 .0325 .0375 .0425 .0475 ... +0.00 0.00 0.00 0.00 0.00 ... +``` + +**Two columns** (frequency, density): +``` +# f S +0.0200 0.00 +0.0325 0.00 +0.0375 0.00 +``` + +Path may be a bare name (resolved next to the YAML), a relative path, or absolute. + +**Example to write an elevation time-series from custom spectrum file:** +```bash +python MHKiT_writer.py <../path_to_hydro_yaml/case.hydro.yaml> --spectrum_file custom_spectrum.txt +``` + + +**Fit a JONSWAP curve to a custom spectrum without writing to yaml:** +```bash +python MHKiT_writer.py ../5sa/custom_waves/5sa_custom.hydro.yaml --spectrum_file test_spectrum.txt --plot_spectrum_js +``` + +image + + +--- + +### `--elevation_duration` / `--elevation_dt` +Length and time step of the generated eta record (default = 600 s; dt from the companion `simulation.time_step`) + +**Works with:** any `--type` and `--spectrum`, but defines time length and intervals for `--spectrum custom` generated eta files + +**Important:** duration *is* frequency resolution. The eta record's frequency grid is fixed by the time vector as + +``` +df = 1 / (n_samples * dt) +``` + +so a longer record resolves the spectrum more finely but produces more wave components and a bigger file. At `dt = 0.01 s`, 600 s gives 600 components and a 1.2 MB eta file; 10800 s gives 10790 components and 22 MB. + +**Example:** +```bash +python MHKiT_writer.py --type irregular --spectrum jonswap --elevation_duration 1200 --elevation_dt 0.01 +``` + +--- + +### `--elevation_file` +Name of the eta file to write, **or** an existing eta file to consume as-is. + +If the file already exists, nothing is downloaded or generated — the record is used exactly as it is and only the YAMLs are updated. Default name is `_eta.txt`. + +To use existing elevation time-series file: +**Example:** +```bash +--elevation_file downloaded_eta.txt +``` + +To use write an elevation time-series file from buoy data: +**Example:** +```bash +--buoy 46050 --date 2024-10-12 --type irregular --spectrum custom --elevation_duration 100 --elevation_file newer_eta.txt +``` +image + +Interpolates the buoy spectral data and generates a time-series elevation txt file + +image + +And changes the yaml file to read from the new .txt eta file generated + + +--- + +### `--eta_method` +How SEA-Stack consumes the eta file (default = `dft`) + +- `dft` — extracts discrete wave components, giving a spatial free surface, so the GUI wireframe renders +- `irf` — `EtaTableWaveField`, a point time series; the body moves but no wave surface is drawn + + +--- + +## Small notes + +- **`nfrequencies` is computed** SEA-Stack's DFT eta import only evaluates the orthogonal Fourier frequencies `k/(N·dt)`, and strides — silently dropping energy — if `nfrequencies` is below the in-band bin count. The script sizes it to cover the whole default 0.001–1.0 Hz eta band. +- **Hm0 in the file header vs the statistics block** can differ by a percent or two because the header uses a plain trapezoid integration whereas the statistics use MHKiT's `significant_wave_height`. The eta record targets the former. +- **SSL / VPN.** If NDBC downloads fail with a certificate error, the `REQUESTS_CA_BUNDLE` and `SSL_CERT_FILE` lines near the top of the script point at a local root CA and may need updating for your machine. +- **Plot flags do not write YAML.** Add `--spectrum`, `--type irregular` if you want a plot *and* an update in the same run. From 6587646bd52013c6e5c13b9761e27888d8dc0664 Mon Sep 17 00:00:00 2001 From: Allie Date: Mon, 14 Sep 2026 09:52:48 -0600 Subject: [PATCH 7/9] Update MHKiT_writer_README.md --- data/demos/features/MHKiT_writer_README.md | 115 +++++++++++++++++++++ 1 file changed, 115 insertions(+) diff --git a/data/demos/features/MHKiT_writer_README.md b/data/demos/features/MHKiT_writer_README.md index 26685d1..a6c8553 100644 --- a/data/demos/features/MHKiT_writer_README.md +++ b/data/demos/features/MHKiT_writer_README.md @@ -21,6 +21,7 @@ With capabilities such as - monthly wave statistics - directional spectra - wind data + --------- @@ -441,6 +442,120 @@ How SEA-Stack consumes the eta file (default = `dft`) - `irf` — `EtaTableWaveField`, a point time series; the body moves but no wave surface is drawn +--- + +## Directional wave partitions + +NDBC transmits the directional distribution as four moments per frequency (`swdir` α₁, `swdir2` α₂, `swr1` r₁, `swr2` r₂) alongside `swden`. These are reconstructed into a full 2D spectrum S(f,θ), split into distinct wave systems, and each system is fitted to the JONSWAP + cos-2s form that SEA-Stack accepts as `waves.partitions`. + +The reconstruction uses the **maximum entropy method** (Lygre & Krogstad) rather than the truncated Fourier series. With only two directional harmonics, the Fourier form rings negative opposite the peak whenever r₁ + r₂ > 0.5 — inventing energy from directions that carry none. MEM is non-negative by construction and resolves twin directional peaks. + +Partitions are then found by a **steepest-ascent watershed** over the (frequency, direction) plane, wrapping in θ. Every cell climbs to a local maximum; cells sharing a maximum form one wave system. This finds systems in 2D, so a swell and a wind sea arriving from different directions at overlapping frequencies are separated — which a frequency-only split cannot do. + +Per partition: + +| YAML key | Derived from | +|---|---| +| `Hs` | 4√(energy integrated over the partition's cells) | +| `Tp` | 1 / frequency of the partition's peak cell | +| `gamma` | least-squares JONSWAP fit to the partition's frequency spectrum | +| `direction` | local maximum of the directional marginal, parabolic sub-bin refined | +| `s` | least-squares cos-2s fit to the measured D(θ), weighted by energy | + +**Direction convention:** NDBC reports the direction waves come *from*, clockwise from true North. SEA-Stack `direction` is the direction waves travel *toward*, counter-clockwise from +X. The conversion is `(270 − α) mod 360`, and both values are printed. + + + +### `--plot_wavedirection_cos` +Compare the measured directional spectrum against the cos-2s reconstruction SEA-Stack would generate. **Never writes the YAML**, even alongside `--partitions` or the custom-spectrum flags. + +- `--plot_wavedirection_cos` → uses the recommended number of partitions +- `--plot_wavedirection_cos 4` → forces 4 partitions + +Produces three polar panels — measured (MEM), reconstructed (N × JONSWAP × cos-2s), and their difference — plus the per-partition directional distributions with the `n_theta` sampling marked. Prints the RMS error and where the largest residual sits. + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wavedirection_cos --date 01-03-2022 +``` + + + +--- + +### `--partitions` +Split the measured directional spectrum into wave systems and write them as `waves.partitions`. Implies `--type irregular`. Plots the same comparison as `--plot_wavedirection_cos`. + +- `--partitions` → uses the recommended count +- `--partitions 3` → forces 3 partitions + +The recommendation counts systems holding at least 5% of the total energy, capped at 6. The full energy ladder is printed so you can see whether your chosen number agrees with the data: + +``` +Watershed found 14 raw system(s); 2 hold >=5% of the energy -> recommending 2 partition(s) + energy share of the strongest: 63.3%, 31.9%, 1.7%, 1.0%, 0.7%, 0.4% +``` + +**Example:** +```bash +python MHKiT_writer.py ../5sa/bimodal/5sa_bimodal.hydro.yaml --buoy 46050 --partitions --date 01-03-2022 +``` + + + +Produces: + +```yaml + waves: + type: irregular + discretization: + n_omega: 64 + n_theta: 21 + seed: 42 + partitions: + - spectrum: jonswap + Hs: 5.647 + Tp: 12.121 + gamma: 1.048 + direction: 82.4 + spreading: + type: cos2s + s: 23.91 + - spectrum: jonswap + Hs: 3.863 + Tp: 11.429 + gamma: 1.024 + direction: 359.4 + spreading: + type: cos2s + s: 11.04 +``` + +--- + +### `--n_omega` / `--n_theta` +Discretization written to `waves.discretization` (defaults 64 and 21). + +**Used with:** `--partitions` + +`s` controls how concentrated each lobe is — higher `s` is a narrower beam — while `n_theta` controls how finely SEA-Stack samples it. Total wave components are `n_omega × n_theta × partitions`, all evaluated per surface vertex per frame, so this is the main cost knob. + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --partitions 2 --n_omega 64 --n_theta 31 +``` + +--- + +### Notes on directional partitioning + +- **Use `--date` or `--time`, not a year.** Averaging thousands of records smooths every local maximum away — a full year of buoy 46050 collapses to a single system, while one day resolves 9–14. The script warns when more partitions are requested than the spectrum contains, since merging can combine systems but never split one. +- **More partitions is not automatically better.** Beyond the systems that hold real energy, extra partitions fragment a single physical wave train and each one still costs `n_omega × n_theta` components. +- **cos-2s cannot match MEM exactly.** Both are normalised so ∫D dθ = 1, but the measured distribution is zero outside its partition while cos-2s has tails around the whole circle. A peak match of 90–97% is the structural limit, not a fitting failure. +- **Swell is inner, wind sea is outer.** Radius is frequency, so low-frequency swell sits near the centre with a tight lobe (high `s`), and locally generated wind sea sits toward the rim with a broad lobe (low `s`) aligned with the local wind. Cross-check against `--plot_wind` for the same date. +- Needs buoy directional data (`swdir` / `swr1` / `swr2`); not available with `--spectrum_file`. + + --- ## Small notes From a3f10021afcf978d2b404d4303bcb069896653d4 Mon Sep 17 00:00:00 2001 From: Allie Date: Mon, 14 Sep 2026 12:00:57 -0600 Subject: [PATCH 8/9] Update MHKiT_writer_README.md --- data/demos/features/MHKiT_writer_README.md | 119 +++++++-------------- 1 file changed, 41 insertions(+), 78 deletions(-) diff --git a/data/demos/features/MHKiT_writer_README.md b/data/demos/features/MHKiT_writer_README.md index a6c8553..4178241 100644 --- a/data/demos/features/MHKiT_writer_README.md +++ b/data/demos/features/MHKiT_writer_README.md @@ -1,8 +1,8 @@ # MHKiT Wave YAML Writer -A Python tool for analyzing real ocean wave data and writing it into SEA-Stack YAML input files. +A Python tool for analyzing real ocean wave data and writing characteristics into SEA-Stack YAML input files. -It pulls spectral wave density (`swden`) records from NDBC buoys via MHKiT, computes sea-state statistics, and writes the result into device `*.hydro.yaml` file +It pulls spectral wave density (`swden`) records from NDBC buoys via MHKiT, computes sea-state statistics, and writes the result into device `*.hydro.yaml` file. For custom wave directional distributions, it pulls four moments per frequency (`swdir` α₁, `swdir2` α₂, `swr1` r₁, `swr2` r₂) from MKHiT alongside `swden`. These are reconstructed into a full 2D spectrum S(f,θ), and fitted to the JONSWAP + cos-2s form that SEA-Stack accepts. The reconstruction uses the **maximum entropy method** (MEM) (Lygre & Krogstad, 1986). It supports: - Regular waves (height + period from measured data) @@ -260,6 +260,25 @@ python MHKiT_writer.py --buoy 46050 --plot_wavedirection spread --year 2024 --n image +--- + +### `--plot_wavedirection_cos` +Compare the measured directional spectrum against the cos-2s reconstruction for SEA-Stack. Use with date/time/year **Does not rewrites the YAML** + +- `--plot_wavedirection_cos` → uses the recommended number of partitions +- `--plot_wavedirection_cos jan --year 2020` → averages over the month for 2020 +- `--plot_wavedirection_cos --partitions 3` → uses three partitions + +Produces three polar panels — measured (MEM), reconstructed (N × JONSWAP × cos-2s), and their RMS error difference, as well as the per-partition directional distributions. + +**Example:** +```bash +python MHKiT_writer.py --buoy 46050 --plot_wavedirection_cos --date 01-03-2022 +``` + +image + + --- ### `--plot_wind` @@ -444,118 +463,60 @@ How SEA-Stack consumes the eta file (default = `dft`) --- -## Directional wave partitions - -NDBC transmits the directional distribution as four moments per frequency (`swdir` α₁, `swdir2` α₂, `swr1` r₁, `swr2` r₂) alongside `swden`. These are reconstructed into a full 2D spectrum S(f,θ), split into distinct wave systems, and each system is fitted to the JONSWAP + cos-2s form that SEA-Stack accepts as `waves.partitions`. - -The reconstruction uses the **maximum entropy method** (Lygre & Krogstad) rather than the truncated Fourier series. With only two directional harmonics, the Fourier form rings negative opposite the peak whenever r₁ + r₂ > 0.5 — inventing energy from directions that carry none. MEM is non-negative by construction and resolves twin directional peaks. +### `--partitions` +Recommended number of partitions are identified using a local-maximum steepest-ascent watershed over the (frequency, direction) plane. This finds systems in 2D, so a swell and a wind sea arriving from different directions at overlapping frequencies are separated. Forcing a count higher than the number of local maxima switches to residual peeling (fit the strongest JONSWAP × cos2s system, subtract it, fit the next from what is left). This is written to the hydro.yaml under waves.partitions and plots the measured vs reconstructed wave directions same as --plot_wavedirection_cos. -Partitions are then found by a **steepest-ascent watershed** over the (frequency, direction) plane, wrapping in θ. Every cell climbs to a local maximum; cells sharing a maximum form one wave system. This finds systems in 2D, so a swell and a wind sea arriving from different directions at overlapping frequencies are separated — which a frequency-only split cannot do. +- `--partitions` → uses the recommended count and initiates writing to the yaml +- `--partitions 3` → forces 3 partitions Per partition: - | YAML key | Derived from | |---|---| | `Hs` | 4√(energy integrated over the partition's cells) | -| `Tp` | 1 / frequency of the partition's peak cell | -| `gamma` | least-squares JONSWAP fit to the partition's frequency spectrum | +| `Tp` | least-squares JONSWAP fit, seeded from the peak cell | +| `gamma` | same fit, over the full frequency axis | | `direction` | local maximum of the directional marginal, parabolic sub-bin refined | | `s` | least-squares cos-2s fit to the measured D(θ), weighted by energy | -**Direction convention:** NDBC reports the direction waves come *from*, clockwise from true North. SEA-Stack `direction` is the direction waves travel *toward*, counter-clockwise from +X. The conversion is `(270 − α) mod 360`, and both values are printed. - - - -### `--plot_wavedirection_cos` -Compare the measured directional spectrum against the cos-2s reconstruction SEA-Stack would generate. **Never writes the YAML**, even alongside `--partitions` or the custom-spectrum flags. - -- `--plot_wavedirection_cos` → uses the recommended number of partitions -- `--plot_wavedirection_cos 4` → forces 4 partitions - -Produces three polar panels — measured (MEM), reconstructed (N × JONSWAP × cos-2s), and their difference — plus the per-partition directional distributions with the `n_theta` sampling marked. Prints the RMS error and where the largest residual sits. **Example:** ```bash -python MHKiT_writer.py --buoy 46050 --plot_wavedirection_cos --date 01-03-2022 +python MHKiT_writer.py ../5sa/bimodal/5sa_bimodal.hydro.yaml --buoy 46050 --partitions --date 12-16 --year 2025 --time 00:00-12:00 ``` - - +image --- -### `--partitions` -Split the measured directional spectrum into wave systems and write them as `waves.partitions`. Implies `--type irregular`. Plots the same comparison as `--plot_wavedirection_cos`. +### `--spread` -- `--partitions` → uses the recommended count -- `--partitions 3` → forces 3 partitions +Overrides the fitted cos2s exponent s for each partition. One value per partition, ordered by decreasing energy. Higher s is a narrower. -The recommendation counts systems holding at least 5% of the total energy, capped at 6. The full energy ladder is printed so you can see whether your chosen number agrees with the data: +Used with: --partitions or --plot_wavedirection_cos -``` -Watershed found 14 raw system(s); 2 hold >=5% of the energy -> recommending 2 partition(s) - energy share of the strongest: 63.3%, 31.9%, 1.7%, 1.0%, 0.7%, 0.4% -``` +The value count must match the partition count or it errors out. Applies to the polar reconstruction, the D(θ) overlay curves and the YAML, so you can preview a value with --plot_wavedirection_cos before committing it with --partitions. **Example:** ```bash -python MHKiT_writer.py ../5sa/bimodal/5sa_bimodal.hydro.yaml --buoy 46050 --partitions --date 01-03-2022 -``` - - - -Produces: - -```yaml - waves: - type: irregular - discretization: - n_omega: 64 - n_theta: 21 - seed: 42 - partitions: - - spectrum: jonswap - Hs: 5.647 - Tp: 12.121 - gamma: 1.048 - direction: 82.4 - spreading: - type: cos2s - s: 23.91 - - spectrum: jonswap - Hs: 3.863 - Tp: 11.429 - gamma: 1.024 - direction: 359.4 - spreading: - type: cos2s - s: 11.04 +python MHKiT_writer.py --buoy 46050 --date 12-16 --time 00:00-12:00 --year 2025 --partitions 2 --spread 50 12 --plot_wavedirection_cos ``` +image + + --- ### `--n_omega` / `--n_theta` -Discretization written to `waves.discretization` (defaults 64 and 21). +Discretization for directional waves written to `waves.discretization` (defaults 64 and 21). **Used with:** `--partitions` -`s` controls how concentrated each lobe is — higher `s` is a narrower beam — while `n_theta` controls how finely SEA-Stack samples it. Total wave components are `n_omega × n_theta × partitions`, all evaluated per surface vertex per frame, so this is the main cost knob. +While `s` controls how concentrated the spread of each lobe is, n_theta bins span the full circle centered on each partition's direction, controlling how finely SEA-Stack samples it. A narrow lobe (high s) wastes most of its bins on near-zero directions, whereas n_omega bins are spaced linearly in Hz across 0.25/Tp_max to 4.0/Tp_min, taken from the partition peak periods. Total wave components are n_omega × n_theta × partitions, all evaluated per surface vertex per frame, so this is the main cost knob. Mixing a short-period wind sea with a long-period swell widens that range and coarsens resolution at the swell peak, and linear spacing makes the sea repeat every 1/df seconds — it is recommended to pick n_omega so that exceeds your simulation.end_time. **Example:** ```bash python MHKiT_writer.py --buoy 46050 --partitions 2 --n_omega 64 --n_theta 31 ``` ---- - -### Notes on directional partitioning - -- **Use `--date` or `--time`, not a year.** Averaging thousands of records smooths every local maximum away — a full year of buoy 46050 collapses to a single system, while one day resolves 9–14. The script warns when more partitions are requested than the spectrum contains, since merging can combine systems but never split one. -- **More partitions is not automatically better.** Beyond the systems that hold real energy, extra partitions fragment a single physical wave train and each one still costs `n_omega × n_theta` components. -- **cos-2s cannot match MEM exactly.** Both are normalised so ∫D dθ = 1, but the measured distribution is zero outside its partition while cos-2s has tails around the whole circle. A peak match of 90–97% is the structural limit, not a fitting failure. -- **Swell is inner, wind sea is outer.** Radius is frequency, so low-frequency swell sits near the centre with a tight lobe (high `s`), and locally generated wind sea sits toward the rim with a broad lobe (low `s`) aligned with the local wind. Cross-check against `--plot_wind` for the same date. -- Needs buoy directional data (`swdir` / `swr1` / `swr2`); not available with `--spectrum_file`. - - --- ## Small notes @@ -564,3 +525,5 @@ python MHKiT_writer.py --buoy 46050 --parti - **Hm0 in the file header vs the statistics block** can differ by a percent or two because the header uses a plain trapezoid integration whereas the statistics use MHKiT's `significant_wave_height`. The eta record targets the former. - **SSL / VPN.** If NDBC downloads fail with a certificate error, the `REQUESTS_CA_BUNDLE` and `SSL_CERT_FILE` lines near the top of the script point at a local root CA and may need updating for your machine. - **Plot flags do not write YAML.** Add `--spectrum`, `--type irregular` if you want a plot *and* an update in the same run. +- **cos-2s cannot match MEM exactly.** Both are normalised so ∫D dθ = 1, but the measured distribution is zero outside its partition while cos-2s has tails around the whole circle. A peak match of 90–97% is the structural limit, not a fitting failure. +- **Swell is inner, wind sea is outer.** Radius is frequency, so low-frequency swell sits near the center with a tight lobe, and locally generated wind sea sits toward the rim with a broad lobe aligned with the local wind. Cross-check wave_direction against `--plot_wind` for the same date. From 419d6f8867d48d6fd6b2656c2fedd3ec5001be0e Mon Sep 17 00:00:00 2001 From: Allie Date: Mon, 14 Sep 2026 13:56:17 -0600 Subject: [PATCH 9/9] Update MHKiT_writer.py --- data/demos/features/MHKiT_writer.py | 657 +++++++++++++++++++++++++++- 1 file changed, 654 insertions(+), 3 deletions(-) diff --git a/data/demos/features/MHKiT_writer.py b/data/demos/features/MHKiT_writer.py index 5ada5da..4ddf1d5 100644 --- a/data/demos/features/MHKiT_writer.py +++ b/data/demos/features/MHKiT_writer.py @@ -69,7 +69,13 @@ const="elevation", default=None, choices=["elevation", "energy", "spread"], - help="Plot directional wave spectrum (elevation/energy/spread). Default: elevation" + help="Plot the measured directional wave spectrum (elevation/energy/spread). Default: elevation" +) +parser.add_argument( + "--plot_wavedirection_cos", + nargs="*", + default=None, + help="Plot the measured directional distribution against the cos2s fit sea-stack would use. Optional: a 4-digit year and/or months (e.g., --plot_wavedirection_cos jan feb 2020), or omit to use --year / --date. Use --partitions to set the number of wave systems; combining the two plots that count without writing the YAML. Shows measured, reconstructed and difference as polar contours plus the per-partition directional distributions. Never writes the YAML." ) parser.add_argument( "--spectrum", @@ -84,6 +90,33 @@ default=None, help="Read the custom spectrum from a text file instead of downloading it from a buoy, so --buoy is not needed. Accepts a bare name, a path relative to the YAML directory, or an absolute path. Two layouts are recognised: two rows (frequencies then spectral densities) or two columns (frequency, density) - blank lines and lines starting with # or %% are ignored. Units are Hz and m^2/Hz. Implies --type irregular --spectrum custom." ) +parser.add_argument( + "--partitions", + nargs="?", + type=int, + const=0, + default=None, + help="Split the measured directional spectrum into wave systems and write them as sea-stack waves.partitions, each with its own Hs, Tp, gamma, direction and cos2s spreading s. Give a number to force it, or omit the number to use the count the watershed recommends (systems holding at least 5%% of the energy). Needs buoy directional data (swdir/swr); implies --type irregular." +) +parser.add_argument( + "--n_omega", + type=int, + default=64, + help="[--partitions] Number of frequency bins written to waves.discretization.n_omega. Default: 64" +) +parser.add_argument( + "--n_theta", + type=int, + default=21, + help="[--partitions] Number of directional bins written to waves.discretization.n_theta. Default: 21" +) +parser.add_argument( + "--spread", + nargs="+", + type=float, + default=None, + help="[--partitions / --plot_wavedirection_cos] Override the fitted cos2s exponent s. Give one value per partition, ordered by decreasing energy (e.g., --partitions 2 --spread 0.4 1.2). Higher s is a narrower, taller directional lobe. Applies to both the plot and the YAML." +) parser.add_argument( "--elevation_duration", type=float, @@ -212,6 +245,34 @@ parser.error("--buoy is required, unless --spectrum custom is reading an existing " "--elevation_file or a --spectrum_file") +if args.partitions is not None: + if args.spectrum_file: + parser.error("--partitions needs buoy directional data and cannot be used with --spectrum_file") + if args.buoy is None: + parser.error("--partitions requires --buoy") + if args.partitions < 0: + parser.error("--partitions must be 0 (auto) or more") + args.type = "irregular" + +if args.plot_wavedirection is not None: + if args.spectrum_file: + parser.error("--plot_wavedirection needs buoy directional data and cannot be used with --spectrum_file") + if args.buoy is None: + parser.error("--plot_wavedirection requires --buoy") + +if args.plot_wavedirection_cos is not None: + if args.spectrum_file: + parser.error("--plot_wavedirection_cos needs buoy directional data and cannot be used with --spectrum_file") + if args.buoy is None: + parser.error("--plot_wavedirection_cos requires --buoy") + +if args.spread is not None: + if args.partitions is None and args.plot_wavedirection_cos is None: + parser.error("--spread requires --partitions or --plot_wavedirection_cos") + bad_s = [v for v in args.spread if v <= 0] + if bad_s: + parser.error(f"--spread values must be greater than 0; got {bad_s}") + # Flatten gamma list if it was provided (handles both --gamma 1.0 1.5 and --gamma 1.0 --gamma 1.5) if args.gamma: args.gamma = [g for sublist in args.gamma for g in sublist] if any(isinstance(i, list) for i in args.gamma) else args.gamma @@ -388,6 +449,7 @@ def parse_month_list(values, min_count=2): args.plot_heatmap is not None or args.plot_wavestats is not None or args.plot_wavedirection is not None or + args.plot_wavedirection_cos is not None or args.plot_spectrum_pm is not None or args.plot_spectrum_js is not None or args.plot_wind is not None @@ -398,12 +460,19 @@ def parse_month_list(values, min_count=2): explicit_update_intent = ( user_set_spectrum or args.depth is not None or + args.partitions is not None or user_set_type ) # If plot-only and no explicit update intent, skip YAML update (height/period used only for plotting) # If no plotting options, then height/period alone means update YAML skip_yaml_update = plot_only_options and not explicit_update_intent + +# --plot_wavedirection_cos is inspection only: it never writes, even next to +# --partitions or the custom-spectrum flags. Use --partitions alone to write. +if args.plot_wavedirection_cos is not None: + skip_yaml_update = True + should_update_yaml = not skip_yaml_update # Validate that path was provided, but only if updating YAML (not for plot-only) @@ -710,6 +779,30 @@ def parse_month(month_input): return None, False +# --plot_wavedirection_cos mixes a partition count with an optional year and months. +# --plot_wavedirection_cos selects the period only; the partition count comes from +# --partitions. Runs after parse_month because it calls it. +dir_months = None +dir_year_override = None +if args.plot_wavedirection_cos is not None: + _months, _bad = [], [] + for tok in args.plot_wavedirection_cos: + try: + n = int(tok) + except (TypeError, ValueError): + n = None + if n is not None and n >= 1900: + dir_year_override = n + continue + num, is_month = parse_month(tok) + _months.append(num) if is_month else _bad.append(tok) + if _bad: + parser.error(f"--plot_wavedirection_cos: unrecognised value(s) {_bad}. Expected a " + f"4-digit year and/or months (jan, feb, ... or 1-12)") + dir_months = sorted(set(_months)) or None + args.plot_wavedirection_cos = True + + def pm_spectrum(freq, hm0, tp): """ Unified Pierson-Moskowitz spectrum function. @@ -889,6 +982,35 @@ def objective(params): return gamma_fit, final_error, tp_fit, hit_bound + +SPREAD_S_MIN = 0.5 +SPREAD_S_MAX = 100.0 + + +def fit_cos2s_spread(measured_D, thetas, theta_mean, d_theta): + """Fit the cos2s exponent s to a measured directional distribution. + + The residual is weighted by the measured density. Unweighted least squares is + dominated by the tails, which hold most of the sample points but little energy, + and that drags s down so the lobe fits too broad and too low at the peak. + + Returns: (s_fit, error, hit_bound) + """ + weight = measured_D / max(float(measured_D.sum()), 1e-30) + + def shape(s_val): + d = np.abs(np.cos(0.5 * (thetas - theta_mean))) ** (2.0 * s_val) + return d / max(d.sum() * d_theta, 1e-30) + + def sse(s_val): + return float(np.sum(weight * (shape(float(s_val)) - measured_D) ** 2)) + + res = minimize_scalar(sse, bounds=(SPREAD_S_MIN, SPREAD_S_MAX), + method='bounded', options={'xatol': 1e-3}) + s_fit = float(np.clip(res.x, SPREAD_S_MIN, SPREAD_S_MAX)) + hit_bound = (abs(s_fit - SPREAD_S_MIN) < 1e-2) or (abs(s_fit - SPREAD_S_MAX) < 1e-2) + return s_fit, sse(s_fit), hit_bound + # ========== End spectrum helper functions ========== def find_simulation_yaml(hydro_path): @@ -1308,6 +1430,487 @@ def parse_year_spec(year_args, available_years): if len(valid_gusts) > 0: print(f"Mean Wind Gusts: {valid_gusts.mean():.2f} m/s") +#~~~~~~~~~~~~~~ Directional partitions from the measured spectrum ~~~~~~~~~~~~~~ + +# NDBC transmits the directional distribution as a truncated Fourier series, +# D(f,th) = 1/pi * [1/2 + r1*cos(th-a1) + r2*cos(2*(th-a2))] +# and sea-stack accepts cos2s spreading per partition. For cos2s the first moment +# is r1 = s/(s+1), so s = r1/(1-r1) falls straight out of swr1 with no fitting. +wave_partitions = [] +partition_specs = [] + +if args.partitions is not None or args.plot_wavedirection_cos is not None: + # --partitions carries the count; 0 or absent means let the watershed recommend + n_bands = args.partitions if args.partitions else None + print("\n" + "=" * 60) + print(f"Directional partitions ({'auto' if n_bands is None else n_bands})") + print("=" * 60) + + dir_year = (dir_year_override if dir_year_override is not None + else (selected_years[0] if selected_years else available_years[-1])) + try: + dir_all = ndbc.request_directional_data(buoy_number, dir_year) + except Exception as e: + print(f"ERROR: could not fetch directional data for {dir_year}: {e}") + exit(1) + + dir_clean = dir_all.dropna(dim='date', how='all') + if len(dir_clean.date) == 0: + print(f"ERROR: no valid directional data for {dir_year}") + exit(1) + + dir_label = str(dir_year) + if dir_months: + month_names = ', '.join(pd.Timestamp(2000, m, 1).strftime('%b') for m in dir_months) + keep = np.isin(pd.DatetimeIndex(dir_clean.date.values).month, dir_months) + if not keep.any(): + print(f"ERROR: no directional records for {month_names} {dir_year}") + exit(1) + dir_clean = dir_clean.isel(date=np.flatnonzero(keep)) + dir_label = f"{month_names} {dir_year}" + + if extract_targets: + dd = pd.DatetimeIndex(dir_clean.date.values) + mask = np.zeros(len(dd), dtype=bool) + for t in extract_targets: + if t['kind'] == 'range': + mask |= (dd >= t['timestamp']) & (dd <= t['end']) + elif t['kind'] == 'daily': + mask |= (dd.normalize() == t['timestamp'].normalize()) + else: + near = dd[np.argmin(np.abs(dd - t['timestamp']))] + if abs((near - t['timestamp']).total_seconds()) < 3600: + mask |= (dd == near) + if mask.any(): + dir_clean = dir_clean.isel(date=np.flatnonzero(mask)) + labels = [t['label'] for t in extract_targets] + dir_label = labels[0] if len(labels) == 1 else f"{labels[0]} ... {labels[-1]}" + else: + print(f" WARNING: no directional records in the selection; using all of {dir_year}") + print(f" Source: {len(dir_clean.date)} directional records ({dir_label})") + + # Average the moments as vectors; averaging the angles directly is wrong + _rad = np.deg2rad + a1m = (dir_clean['swr1'] * np.cos(_rad(dir_clean['swdir']))).mean(dim='date').values + b1m = (dir_clean['swr1'] * np.sin(_rad(dir_clean['swdir']))).mean(dim='date').values + a2m = (dir_clean['swr2'] * np.cos(2 * _rad(dir_clean['swdir2']))).mean(dim='date').values + b2m = (dir_clean['swr2'] * np.sin(2 * _rad(dir_clean['swdir2']))).mean(dim='date').values + S_dir = dir_clean['swden'].mean(dim='date').values + freq_dir = dir_clean['frequency'].values + + ok = np.isfinite(S_dir) & np.isfinite(a1m) & np.isfinite(b1m) + freq_dir, S_dir = freq_dir[ok], S_dir[ok] + a1m, b1m, a2m, b2m = a1m[ok], b1m[ok], a2m[ok], b2m[ok] + + thetas_deg = np.arange(0.0, 360.0, 5.0) + thetas_rad = np.deg2rad(thetas_deg) + d_theta = np.deg2rad(5.0) + + def mem_spread(c1, c2, theta): + """Lygre & Krogstad maximum-entropy D(theta) from the first two moments. + + D = sigma / (2*pi*|1 - phi1*e^-i*th - phi2*e^-2i*th|^2) is non-negative by + construction and reproduces c1 and c2 exactly, so nothing is clipped here. + The 2-term Fourier reconstruction is the one that rings negative when + r1 + r2 > 0.5; that is a truncation artifact, not energy. + """ + denom = 1.0 - np.abs(c1) ** 2 + denom = np.where(np.abs(denom) < 1e-12, 1e-12, denom) + phi1 = (c1 - c2 * np.conj(c1)) / denom + phi2 = c2 - c1 * phi1 + sigma = np.real(1.0 - phi1 * np.conj(c1) - phi2 * np.conj(c2)) + + # sigma <= 0 means the moment set is not realisable by any distribution, + # which vector-averaging across time can in principle produce. + n_bad = int(np.sum(sigma <= 0.0)) + if n_bad: + print(f" WARNING: {n_bad} frequency bin(s) have an inconsistent moment set " + f"(MEM sigma <= 0); their spread is unreliable.") + sigma = np.clip(sigma, 1e-12, None) + + e1 = np.exp(-1j * theta)[None, :] + e2 = np.exp(-2j * theta)[None, :] + num = sigma[:, None] + den = np.abs(1.0 - phi1[:, None] * e1 - phi2[:, None] * e2) ** 2 + D = num / (2.0 * np.pi * np.clip(den, 1e-12, None)) + return D / np.clip(D.sum(axis=1, keepdims=True) * d_theta, 1e-30, None) + + c1_f = a1m + 1j * b1m + c2_f = a2m + 1j * b2m + D_mem = mem_spread(c1_f, c2_f, thetas_rad) + S2 = S_dir[:, None] * D_mem # m^2/Hz/rad on (frequency, direction) + + df_dir = np.gradient(freq_dir) + cell_E = S2 * df_dir[:, None] * d_theta + total_E = float(cell_E.sum()) + if total_E <= 0: + print("ERROR: directional spectrum carries no energy") + exit(1) + + def watershed_label(field): + """Steepest-ascent watershed over (frequency, direction), wrapping in theta. + + Each cell climbs to a local maximum; cells sharing a maximum form one wave + system, so systems are found in 2D rather than imposed as frequency bands. + """ + nf, nth = field.shape + lab = -np.ones((nf, nth), dtype=int) + peaks = [] + for flat in np.argsort(field, axis=None)[::-1]: + i, j = divmod(int(flat), nth) + if field[i, j] <= 0: + continue + best, best_val = None, field[i, j] + for di in (-1, 0, 1): + ii = i + di + if ii < 0 or ii >= nf: + continue + for dj in (-1, 0, 1): + if di == 0 and dj == 0: + continue + jj = (j + dj) % nth + if field[ii, jj] > best_val: + best, best_val = (ii, jj), field[ii, jj] + if best is None or lab[best] < 0: + lab[i, j] = len(peaks) + peaks.append((i, j)) + else: + lab[i, j] = lab[best] + return lab + + def merge_to(lab, n_keep): + """Merge the weakest system into its strongest neighbour until n_keep remain.""" + nf, nth = lab.shape + + def label_energy(l): + return float(cell_E[lab == l].sum()) + + alive = sorted(set(int(v) for v in np.unique(lab) if v >= 0)) + while len(alive) > n_keep: + weakest = min(alive, key=label_energy) + nbrs = set() + idx_i, idx_j = np.where(lab == weakest) + for i, j in zip(idx_i, idx_j): + for di in (-1, 0, 1): + ii = i + di + if ii < 0 or ii >= nf: + continue + for dj in (-1, 0, 1): + n = lab[ii, (j + dj) % nth] + if n >= 0 and n != weakest: + nbrs.add(int(n)) + target = (max(nbrs, key=label_energy) if nbrs + else max((a for a in alive if a != weakest), key=label_energy)) + lab[lab == weakest] = target + alive.remove(weakest) + return lab, alive + + def unit_model(tp, gamma, th0_rad, s_val): + """JONSWAP x cos2s on the (frequency, direction) grid, scaled to m0 = 1.""" + S_f = jonswap_spectrum(freq_dir, 4.0, tp, gamma=gamma) # hm0 = 4 -> m0 = 1 + D_t = np.abs(np.cos(0.5 * (thetas_rad - th0_rad))) ** (2.0 * s_val) + D_t /= max(D_t.sum() * d_theta, 1e-30) + return S_f[:, None] * D_t[None, :] + + def peel_systems(field, n_want): + """Locate n_want systems by repeatedly fitting a JONSWAP x cos2s and subtracting it. + + The watershed can only merge, never split, so it can never return more systems + than the field has local maxima - averaging a month of records smooths them into + one blob. Peeling finds the next system in whatever the previous fits leave + behind, which is what lets a forced --partitions count exceed that limit. + + Returns a list of fitted models, strongest first. + """ + w = df_dir[:, None] * d_theta + R = field.copy() + models = [] + for _ in range(n_want): + if R.max() <= 0: + break + pi0, pj0 = np.unravel_index(int(np.argmax(R)), R.shape) + tp0 = 1.0 / float(freq_dir[pi0]) + th0 = float(thetas_rad[pj0]) + + # m0 enters the model linearly, so only (Tp, gamma, s, theta) need searching + def fit_error(p): + tp, g, s_val, th = p + if not (0.5 * tp0 <= tp <= 2.0 * tp0) or not (GAMMA_MIN <= g <= GAMMA_MAX) \ + or not (SPREAD_S_MIN <= s_val <= SPREAD_S_MAX): + return 1e30 + M = unit_model(tp, g, th, s_val) + denom = float((w * M * M).sum()) + if denom <= 0: + return 1e30 + amp = float((w * R * M).sum()) / denom + if amp <= 0: + return 1e30 + return float((w * (R - amp * M) ** 2).sum()) + + best = None + for g0 in (1.5, 3.3, 6.0): + for s0 in (2.0, 8.0, 25.0): + res = minimize(fit_error, x0=[tp0, g0, s0, th0], method='Nelder-Mead', + options={'maxiter': 600, 'xatol': 1e-4, 'fatol': 1e-12}) + if best is None or res.fun < best.fun: + best = res + + M = unit_model(best.x[0], best.x[1], best.x[3], best.x[2]) + amp = float((w * R * M).sum()) / max(float((w * M * M).sum()), 1e-30) + if amp <= 0: + break + models.append(amp * M) + R = np.clip(R - amp * M, 0.0, None) + return models + + labels0 = watershed_label(S2) + raw_E = sorted((float(cell_E[labels0 == l].sum()) + for l in np.unique(labels0) if l >= 0), reverse=True) + significant = [e for e in raw_E if e >= 0.05 * total_E] + recommended = int(np.clip(len(significant), 1, 6)) + print(f" Watershed found {len(raw_E)} raw system(s); {len(significant)} hold " + f">=5% of the energy -> recommending {recommended} partition(s)") + print(" energy share of the strongest: " + + ", ".join(f"{100 * e / total_E:.1f}%" for e in raw_E[:6])) + if n_bands is None: + n_bands = recommended + + if n_bands > len(raw_E): + print(f" NOTE: {n_bands} partition(s) requested but the spectrum only contains " + f"{len(raw_E)} local maximum/maxima.") + print(f" Merging cannot split one, so the systems are peeled out instead: " + f"fit the strongest, subtract it, fit the next from what is left.") + models = peel_systems(S2, n_bands) + print(f" Peeled {len(models)} system(s) from the residual.") + else: + models = [] + + if len(models) >= 2: + # Each cell goes to whichever peeled model explains it best, so the partitions + # stay disjoint and every bit of measured energy is still accounted for + labels = np.argmax(np.stack(models), axis=0) + alive = sorted(set(int(v) for v in np.unique(labels))) + else: + if n_bands > len(raw_E): + print(f" Peeling did not separate the blob; falling back to " + f"{len(raw_E)} merged system(s).") + labels, alive = merge_to(labels0, n_bands) + alive = sorted(alive, key=lambda l: -float(cell_E[labels == l].sum())) + print(f" Using {len(alive)} partition(s)") + + if args.spread and len(args.spread) != len(alive): + print(f"ERROR: --spread has {len(args.spread)} value(s) but {len(alive)} partition(s) " + f"were found; give one value per partition") + exit(1) + + overlays = [] + for lbl in alive: + sel = (labels == lbl) + E_p = float(cell_E[sel].sum()) + if E_p <= 0: + continue + + m0 = E_p + Hs = 4.0 * np.sqrt(m0) + + # Partition peak cell sets Tp; its direction anchors the partition + masked = np.where(sel, S2, -np.inf) + pi_, pj_ = np.unravel_index(int(np.argmax(masked)), S2.shape) + Tp_bin = 1.0 / float(freq_dir[pi_]) + Tp = Tp_bin + + # Direction comes from the local maximum of the partition's directional + # marginal, not its circular mean: the mean is dragged off the peak whenever + # the lobe is skewed or carries a weak secondary shoulder. + w_cells = np.where(sel, cell_E, 0.0) + D_meas = w_cells.sum(axis=0) + D_meas /= max(D_meas.sum() * d_theta, 1e-30) + + n_th = len(thetas_deg) + jm = int(np.argmax(D_meas)) + y0, y1, y2 = D_meas[(jm - 1) % n_th], D_meas[jm], D_meas[(jm + 1) % n_th] + curv = y0 - 2.0 * y1 + y2 + shift = float(np.clip(0.5 * (y0 - y2) / curv, -0.5, 0.5)) if abs(curv) > 1e-30 else 0.0 + theta_from = float((thetas_deg[jm] + shift * (360.0 / n_th)) % 360.0) + direction_deg = float((270.0 - theta_from) % 360.0) + + # Circular mean retained for reference only + R = complex((w_cells * np.exp(1j * thetas_rad)[None, :]).sum() / E_p) + r1_bar = float(abs(R)) + theta_mean = float(np.degrees(np.angle(R)) % 360.0) + + s_fitted, s_err, s_hit_bound = fit_cos2s_spread( + D_meas, thetas_rad, np.deg2rad(theta_from), d_theta) + s_moment = float(np.clip(r1_bar / max(1.0 - r1_bar, 1e-6), + SPREAD_S_MIN, SPREAD_S_MAX)) + s_fit = float(args.spread[len(partition_specs)]) if args.spread else s_fitted + + # Frequency spectrum of this partition alone, for the JONSWAP fit + S_p = (np.where(sel, S2, 0.0) * d_theta).sum(axis=1) + band = np.flatnonzero(S_p > 0) + f_lo, f_hi = float(freq_dir[band[0]]), float(freq_dir[band[-1]]) + # Fit across the whole frequency axis, not just the band. jonswap_spectrum + # normalises m0 over whatever axis it is handed, so fitting on the band and + # reconstructing on freq_dir would scale the two curves differently; the + # zeros outside the band also make the fit pay for energy put in the tail. + if len(band) >= 4: + gamma_fit, _e, Tp, _hb = fit_jonswap_gamma(S_p, freq_dir, Tp_bin, Hs) + else: + gamma_fit = 3.3 + + partition_specs.append({ + 'spectrum': 'jonswap', + 'Hs': round(float(Hs), 3), + 'Tp': round(float(Tp), 3), + 'gamma': round(float(gamma_fit), 3), + 'direction': round(float(direction_deg), 1), + 'spreading': {'type': 'cos2s', 's': round(float(s_fit), 2)}, + }) + + # Measured directional distribution of this partition's energy + D_fit = np.abs(np.cos(0.5 * (thetas_rad - np.deg2rad(theta_from)))) ** (2.0 * s_fit) + D_fit /= max(D_fit.sum() * d_theta, 1e-30) + overlays.append((f_lo, f_hi, theta_from, s_fit, D_meas, D_fit)) + + print(f"\n Partition {len(partition_specs)}: {f_lo:.4f} - {f_hi:.4f} Hz, " + f"{100 * E_p / total_E:.1f}% of energy") + print(f" Hs = {Hs:.3f} m Tp = {Tp:.2f} s gamma = {gamma_fit:.3f}" + f" (peak bin Tp was {Tp_bin:.2f} s)") + print(f" peak direction: {theta_from:.1f} deg FROM (NDBC) " + f"-> {direction_deg:.1f} deg sea-stack (toward, CCW from +X)") + print(f" circular mean was {theta_mean:.1f} deg FROM " + f"({abs((theta_mean - theta_from + 180) % 360 - 180):.1f} deg off the peak)") + if args.spread: + print(f" cos2s: s = {s_fit:.2f} (--spread override; the fit gave " + f"{s_fitted:.2f}, moment-matched would be {s_moment:.2f})") + else: + print(f" Fitted cos2s: s = {s_fit:.2f}, Error = {s_err:.6e} " + f"(moment-matched s would be {s_moment:.2f})") + print(f" peak D: measured {D_meas.max():.3f}, cos2s {D_fit.max():.3f} " + f"({100 * D_fit.max() / max(D_meas.max(), 1e-30):.0f}% of measured)") + if s_hit_bound and not args.spread: + print(f" NOTE: s pinned at the {SPREAD_S_MIN}-{SPREAD_S_MAX} bound; the measured") + print(f" lobe is outside what a single cos2s can represent") + + if not partition_specs: + print("ERROR: could not extract any partition with energy") + exit(1) + + total_hs = np.sqrt(sum(p['Hs'] ** 2 for p in partition_specs)) + print(f"\n Combined Hs across partitions = {total_hs:.3f} m") + + # Only --partitions writes; --plot_wavedirection_cos is plot-only + wave_partitions = partition_specs if args.partitions is not None else [] + + if args.plot_wavedirection_cos is not None or args.partitions is not None: + n_theta_ss = int(args.n_theta) + _fine = np.linspace(-np.pi, np.pi, 4001) + + # Rebuild the full S(f,theta) sea-stack will generate: for each partition a + # JONSWAP in frequency times a cos2s in direction. This is what the solver + # actually sees, so gamma and Hs enter here, not in the D(theta) curves. + S_recon = np.zeros_like(S2) + for spec, (f_lo, f_hi, th0, s_val, D_m, D_c) in zip(partition_specs, overlays): + S_p = jonswap_spectrum(freq_dir, spec['Hs'], spec['Tp'], gamma=spec['gamma']) + D_p = np.abs(np.cos(0.5 * (thetas_rad - np.deg2rad(th0)))) ** (2.0 * s_val) + D_p /= max(D_p.sum() * d_theta, 1e-30) + S_recon += S_p[:, None] * D_p[None, :] + + diff = S_recon - S2 + vmax = float(np.nanmax(np.abs(diff))) or 1e-12 + smax = float(max(np.nanmax(S2), np.nanmax(S_recon))) or 1e-12 + + # Trim the frequency axis to where the energy actually is + e_cum = np.cumsum(S_dir * df_dir) + f_top = freq_dir[min(len(freq_dir) - 1, + int(np.searchsorted(e_cum, 0.995 * e_cum[-1])) + 3)] + + rel_rms = 100.0 * float(np.sqrt(np.mean(diff ** 2))) / smax + print(f"\n Reconstructed vs measured S(f,theta): RMS {rel_rms:.1f}% of peak") + i_w, j_w = np.unravel_index(int(np.argmax(np.abs(diff))), diff.shape) + print(f" largest error at {freq_dir[i_w]:.3f} Hz, {thetas_deg[j_w]:.0f} deg " + f"({'over' if diff[i_w, j_w] > 0 else 'under'}-predicted)") + + fig = plt.figure(figsize=(16, 9)) + TH, FR = np.meshgrid(thetas_rad, freq_dir) + colors_p = plt.cm.tab10(np.linspace(0, 1, max(len(overlays), 1))) + # contourf takes its levels from the data range, so vmin/vmax alone do + # nothing; explicit level arrays are what make the extremes saturate. + s2max = float(np.nanmax(S2)) or 1e-12 + srmax = float(np.nanmax(S_recon)) or 1e-12 + # One scale for both spectra so the panels are directly comparable + spec_levels = np.linspace(0.0, max(s2max, srmax), 25) + # The residual peak is often a single cell, which contourf renders as an + # invisible sliver. Saturate at the 99th percentile and draw it per cell. + vshow = float(np.nanpercentile(np.abs(diff), 99.0)) + if not np.isfinite(vshow) or vshow <= 0: + vshow = vmax + print(f" difference colour saturates at +/-{vshow:.1f} " + f"(true extreme {vmax:.1f})") + + panels = [ + ("Measured (MEM)", S2, 'viridis', spec_levels, True), + (f"Reconstructed ({len(partition_specs)} x JONSWAP x cos2s)", S_recon, + 'viridis', spec_levels, True), + ("Difference (reconstructed - measured)", diff, 'RdBu_r', None, False), + ] + for k, (title, field, cmap, levels, show_marks) in enumerate(panels): + ax = fig.add_subplot(2, 3, k + 1, projection='polar') + if levels is None: + cf = ax.pcolormesh(TH, FR, field, cmap=cmap, shading='auto', + vmin=-vshow, vmax=vshow) + else: + cf = ax.contourf(TH, FR, field, levels=levels, cmap=cmap, extend='both') + ax.set_theta_zero_location('N') + ax.set_theta_direction(-1) + ax.set_ylim(0, f_top) + ax.set_title(title, fontsize=10, pad=14) + ax.tick_params(labelsize=8) + # Radial line = partition mean direction. Left off the difference panel + # so nothing hides the residual. + if show_marks: + for pidx, ov in enumerate(overlays): + th0 = ov[2] + ax.plot([np.deg2rad(th0), np.deg2rad(th0)], [0.0, f_top], + linestyle='--', linewidth=1.4, color=colors_p[pidx], alpha=0.9) + plt.colorbar(cf, ax=ax, pad=0.10, shrink=0.8, extend='both', + label="S [m^2/Hz/rad]") + + # Directional distributions only: normalised, so no Hs or gamma dependence + axd = fig.add_subplot(2, 1, 2) + colors = colors_p + for i, (f_lo, f_hi, th0, s_val, D_m, D_c) in enumerate(overlays): + axd.plot(thetas_deg, D_m, linewidth=2.2, color=colors[i], + label=f"P{i + 1} measured ({f_lo:.3f}-{f_hi:.3f} Hz)") + axd.plot(thetas_deg, D_c, linewidth=1.8, color=colors[i], linestyle='--', + label=f"P{i + 1} cos2s s={s_val:.1f} @ {th0:.0f} deg") + d_th_ss = 2.0 * np.pi / n_theta_ss + C_s = 1.0 / np.trapezoid(np.abs(np.cos(0.5 * _fine)) ** (2.0 * s_val), _fine) + th_ss = np.deg2rad(th0) - np.pi + (np.arange(n_theta_ss) + 0.5) * d_th_ss + D_ss = C_s * np.abs(np.cos(0.5 * (th_ss - np.deg2rad(th0)))) ** (2.0 * s_val) + axd.plot(np.degrees(th_ss) % 360.0, D_ss, 'o', markersize=3.5, + color=colors[i], alpha=0.8) + axd.set_xlim(0, 360) + axd.set_xticks(range(0, 361, 30)) + axd.set_xlabel("Direction FROM [deg] (circles = sea-stack n_theta sampling)", + fontsize=10) + axd.set_ylabel("D [1/rad]", fontsize=10) + axd.set_title("Directional distribution per partition (normalised: no Hs or gamma)", + fontsize=11) + axd.grid(True, alpha=0.3) + axd.legend(fontsize=8, ncol=max(1, len(overlays)), loc='upper right') + + fig.suptitle(f"{source_title} - Directional Reconstruction ({dir_label}, " + f"{len(partition_specs)} partition(s), n_theta={n_theta_ss})\n" + f"polar: angle = direction FROM, radius = frequency [Hz]; " + f"dashed = partition mean direction", + fontsize=12) + plt.tight_layout() + plt.show() + plt.close() + print("\nDirectional reconstruction comparison displayed") + + #~~~~~~~~~~~~~~ Custom spectrum -> surface elevation time series ~~~~~~~~~~~~~~ # EtaTableWaveField / ComponentSampler::BuildFromEtaFile require one "time:elevation" @@ -1488,7 +2091,55 @@ def parse_year_spec(year_args, available_years): #~~~~~~~~~~~~~~ Update YAML file with calculated wave parameters ~~~~~~~~~~~~~~ -if not skip_yaml_update and os.path.exists(yaml_file_path): +if wave_partitions and not skip_yaml_update and os.path.exists(yaml_file_path): + print(f"\n\nUpdating YAML file: {yaml_file_path}") + + with open(yaml_file_path, 'r') as f: + yaml_data = yaml.safe_load(f) + + if 'waves' not in ((yaml_data or {}).get('hydrodynamics') or {}): + print("ERROR: 'hydrodynamics.waves' section not found in YAML file") + exit(1) + + # Rebuilt from scratch: height / period / spectrum / eta_file are all ignored + # by sea-stack once partitions are present. + waves_block = { + 'type': 'irregular', + 'discretization': {'n_omega': int(args.n_omega), 'n_theta': int(args.n_theta)}, + 'seed': int(args.seed) if args.seed is not None else 42, + 'partitions': wave_partitions, + } + ramp = float(args.ramp_time) if args.ramp_time is not None else 60.0 + if ramp > 0.0: + waves_block['ramp_duration'] = ramp + waves_block['ramp_type'] = "cosine" + if args.depth is not None: + waves_block['depth'] = float(args.depth) + elif water_depth: + waves_block['depth'] = float(water_depth) + + yaml_data['hydrodynamics']['waves'] = waves_block + + class IndentDumper(yaml.SafeDumper): + def increase_indent(self, flow=False, indentless=False): + return super().increase_indent(flow, False) + + # Serialise first: opening with "w" truncates, so a dump failure would destroy the file + rendered = yaml.dump(yaml_data, Dumper=IndentDumper, + default_flow_style=False, sort_keys=False) + with open(yaml_file_path, "w") as f: + f.write(rendered) + + print("Updated wave parameters:") + print(f" Type: irregular with {len(wave_partitions)} partition(s)") + print(f" discretization: n_omega={args.n_omega}, n_theta={args.n_theta} " + f"-> {args.n_omega * args.n_theta * len(wave_partitions)} wave components") + for i, p in enumerate(wave_partitions, 1): + print(f" partition {i}: Hs={p['Hs']} m, Tp={p['Tp']} s, gamma={p['gamma']}, " + f"direction={p['direction']} deg, s={p['spreading']['s']}") + print(f"\nYAML file updated successfully!") + +if not wave_partitions and not skip_yaml_update and os.path.exists(yaml_file_path): print(f"\n\nUpdating YAML file: {yaml_file_path}") with open(yaml_file_path, 'r') as f: @@ -1697,7 +2348,7 @@ def increase_indent(self, flow=False, indentless=False): print(f"\nYAML file updated successfully!") else: print("ERROR: 'hydrodynamics.waves' section not found in YAML file") -elif not skip_yaml_update: +elif not wave_partitions and not skip_yaml_update: # Trying to update but file not found print(f"ERROR: YAML file not found at {yaml_file_path}") exit(1)