diff --git a/docs/files/Spectrum_plot.png b/docs/files/Spectrum_plot.png deleted file mode 100644 index edf4bda9..00000000 Binary files a/docs/files/Spectrum_plot.png and /dev/null differ diff --git a/docs/files/wave_spectra_2d_direction.png b/docs/files/wave_spectra_2d_direction.png new file mode 100644 index 00000000..37a7b931 Binary files /dev/null and b/docs/files/wave_spectra_2d_direction.png differ diff --git a/docs/files/wave_spectra_2d_monthly_mean.png b/docs/files/wave_spectra_2d_monthly_mean.png new file mode 100644 index 00000000..9d1f7309 Binary files /dev/null and b/docs/files/wave_spectra_2d_monthly_mean.png differ diff --git a/docs/files/wave_spectrum_1d_month.png b/docs/files/wave_spectrum_1d_month.png new file mode 100644 index 00000000..da061d3e Binary files /dev/null and b/docs/files/wave_spectrum_1d_month.png differ diff --git a/docs/files/wave_spectrum_1d_months.png b/docs/files/wave_spectrum_1d_months.png new file mode 100644 index 00000000..53f37576 Binary files /dev/null and b/docs/files/wave_spectrum_1d_months.png differ diff --git a/docs/files/wave_spectrum_2d.png b/docs/files/wave_spectrum_2d.png new file mode 100644 index 00000000..b091a0b7 Binary files /dev/null and b/docs/files/wave_spectrum_2d.png differ diff --git a/docs/files/wave_spectrum_diana.png b/docs/files/wave_spectrum_diana.png new file mode 100644 index 00000000..11493dee Binary files /dev/null and b/docs/files/wave_spectrum_diana.png differ diff --git a/docs/index.rst b/docs/index.rst index 196bde2a..f09a678c 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1576,22 +1576,108 @@ CCA profiles Table :header-rows: 1 :file: files/table_cca_profiles.csv -Directional Wave Energy Spectrum Plot + +Wave Spectrum Plots ----------------------------------------- +Monthly Mean 1D Wave Spectrum + +.. code-block:: python + + plots.plot_spectra_1d( + data, + var = 'SPEC', + period = None, + month = None, + method = 'mean', + output_file = 'wave_spectrum_1d_months.png' + ) + +.. image:: files/wave_spectrum_1d_months.png + :width: 500 + +1D Spectrum Mean for a Specific Month Across Several Years + +.. code-block:: python + + plots.plot_spectra_1d( + data, + var = 'SPEC', + period = None, + month = 1, + method = 'mean', + output_file = 'wave_spectrum_1d_month.png' + ) + +.. image:: files/wave_spectrum_1d_month.png + :width: 500 + +2D Wave Spectrum at Time of Maximum Hm0 in Selected Period + +.. code-block:: python + + plots.plot_spectrum_2d( + data, + var = 'SPEC', + period = ('2021-01-01T00', '2024-12-31T23'), + month = None, + method = 'hm0_max', + plot_type = 'pcolormesh', + output_file = 'wave_spectrum_2d.png' + ) + +.. image:: files/wave_spectrum_2d.png + :width: 500 + +Diana Wave Spectrum with Swell and Windsea Partitions, Averaged Over Times with Hm0 ≥ 99th Percentile + +.. code-block:: python + + plots.plot_diana_spectrum( + data, + var = 'SPEC', + period = ('2022-01-01T00', '2022-07-31T23'), + month = None, + method = 'top_1_percent_mean', + partition=True, + plot_type = 'pcolormesh' + freq_mask=True, + output_file = 'wave_spectrum_diana.png' + ) + +.. image:: files/wave_spectrum_diana.png + :width: 500 + +2D Monthly Mean Wave Spectra + +.. code-block:: python + + plots.plot_spectra_2d( + data, + var = 'SPEC', + method = 'mean', + plot_type = 'contour', + cbar = 'multiple', + output_file = 'wave_spectra_2d_monthly_mean.png' + ) + +.. image:: files/wave_spectra_2d_monthly_mean.png + :width: 500 + +Directional Wave Spectra Averaged by Peak Wave Direction in 30° Sectors + .. code-block:: python - plots.plot_spectrum( - spectrum, - frequencies, - directions, - spec_unit='m**2 s', - radius='frequency', - log_radius=False, - output_file='Spectrum_plot.png') + plots.plot_spectra_2d( + data, + var = 'SPEC', + method = 'direction' + plot_type = 'contour', + cbar = 'single', + output_file = 'wave_spectra_2d_direction.png' ) -.. image:: files/Spectrum_plot.png +.. image:: files/wave_spectra_2d_direction.png :width: 500 Tidal Levels Table diff --git a/environment.yml b/environment.yml index 33edd997..63a5b1e0 100644 --- a/environment.yml +++ b/environment.yml @@ -18,4 +18,6 @@ dependencies: - netcdf4 - pydap - setuptools - - virocon \ No newline at end of file + - metpy + - virocon + diff --git a/metocean_stats/__pycache__/__init__.cpython-312.pyc b/metocean_stats/__pycache__/__init__.cpython-312.pyc deleted file mode 100644 index a2d0a2c9..00000000 Binary files a/metocean_stats/__pycache__/__init__.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/plots/__init__.py b/metocean_stats/plots/__init__.py index 812e44f4..9c6a32c2 100644 --- a/metocean_stats/plots/__init__.py +++ b/metocean_stats/plots/__init__.py @@ -10,3 +10,4 @@ from .extreme import * from .general import * from .verification import * +from .spectra import * diff --git a/metocean_stats/plots/__pycache__/__init__.cpython-312.pyc b/metocean_stats/plots/__pycache__/__init__.cpython-312.pyc deleted file mode 100644 index 05fb4de9..00000000 Binary files a/metocean_stats/plots/__pycache__/__init__.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/plots/__pycache__/dir.cpython-312.pyc b/metocean_stats/plots/__pycache__/dir.cpython-312.pyc deleted file mode 100644 index 900a00c0..00000000 Binary files a/metocean_stats/plots/__pycache__/dir.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/plots/__pycache__/extreme.cpython-312.pyc b/metocean_stats/plots/__pycache__/extreme.cpython-312.pyc deleted file mode 100644 index 97d32018..00000000 Binary files a/metocean_stats/plots/__pycache__/extreme.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/plots/__pycache__/general.cpython-312.pyc b/metocean_stats/plots/__pycache__/general.cpython-312.pyc deleted file mode 100644 index 1cf193ed..00000000 Binary files a/metocean_stats/plots/__pycache__/general.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/plots/__pycache__/verification.cpython-312.pyc b/metocean_stats/plots/__pycache__/verification.cpython-312.pyc deleted file mode 100644 index fab75075..00000000 Binary files a/metocean_stats/plots/__pycache__/verification.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/plots/dir.py b/metocean_stats/plots/dir.py index e7a835b1..9a700962 100644 --- a/metocean_stats/plots/dir.py +++ b/metocean_stats/plots/dir.py @@ -138,82 +138,3 @@ def monthly_var_rose(data, return fig -def plot_spectrum_2d(ds, var='SPEC', radius='frequency', log_radius=False, plot_type='contourf', dir_letters=False, output_file='Spectrum_plot.png'): - ''' - This function plots a directional wave spectrum - - Parameters - ---------- - ds : xarray dataset - Should contain variable var with two (frequencies, directions). - var : string - Name of the spectrum variable - radius : string - Should be 'period'/'frequency' to plot the periods/frequencies in the radial direction - log_radius : Boolean - True to get the radial axis on a logarithmic scale - plot_type : string - Can be 'pcolormesh' or 'contourf' - dir_letters : Boolean - Displays the directions as letters instead of degrees if True' - output_file : string - Name of the figure file (xxx.pdf or xxx.png) - - Returns - ------- - Figure matplotlib - - Authors - ------- - Function written by clio-met - ''' - spectrum=ds[var] - dims=list(spectrum.dims) - # Get the coordinates - frequencies=spectrum.coords[dims[0]].to_numpy() - directions=spectrum.coords[dims[1]].to_numpy() - # The last direction should be the same as the first - if directions[0]!=directions[-1]: - if directions[0]<360: - dirc=np.concatenate([directions,directions[0:1]+360]) - else: - dirc=np.concatenate([directions,directions[0:1]]) - spec=np.concatenate([spectrum[:,:],spectrum[:,0:1]],axis=1) - - # Color map with 10 colors - cmap = mpl.cm.hot_r(np.linspace(0,1,10)) - cmap = mpl.colors.ListedColormap(cmap) - - fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}) - ax.set_theta_zero_location('N') - ax.set_theta_direction(-1) - if radius=='period': - rad_data=1/frequencies - elif radius=='frequency': - rad_data=frequencies - - if plot_type=='pcolormesh': - c=ax.pcolormesh(np.radians(dirc), rad_data, spec, cmap='hot_r', shading='auto') - elif plot_type=='contourf': - c=ax.contourf(np.radians(dirc), rad_data, spec, cmap=cmap) - else: - print('plot_type should pcolormesh or contourf') - sys.exit() - - ax.grid(True) - if log_radius==True: - ax.set_rscale('log') - #ax.set_rlim(0,10) - if dir_letters==True: - ticks_loc = ax.get_xticks().tolist() - ax.xaxis.set_major_locator(mticker.FixedLocator(ticks_loc)) - ax.set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW']) - plt.colorbar(c, label='['+spectrum.units+']', pad=0.1) - plt.tight_layout() - if output_file is not None: - plt.savefig(output_file, dpi=200, bbox_inches='tight') - plt.close() - - return fig - - diff --git a/metocean_stats/plots/spectra.py b/metocean_stats/plots/spectra.py new file mode 100644 index 00000000..0a8d7951 --- /dev/null +++ b/metocean_stats/plots/spectra.py @@ -0,0 +1,1009 @@ +import pandas as pd +import xarray as xr +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.colors as mcol +import matplotlib.cm as cm +import matplotlib.ticker as mticker +import matplotlib.patches as patches +import re +from ..stats import spec_funcs +from ..tables import spectra + +def plot_spectra_1d(data, var='SPEC', period=None, month=None, method='mean', output_file='monthly_spectra_1d.png'): + ''' + Plots 1D spectra aggregated by month, for a specific month across multiple years, or for a single timestamp. + + - If no month is specified: + Plots 12 monthly spectra (one subplot per month) plus an overall average spectrum + computed from the monthly aggregates. + - If a month is specified: + Plots one spectrum per year for that month, with an additional average spectrum + across all selected years. + - If `period` contains a single timestamp: + Plots only the spectrum for that specific time. + + Parameters + - data : xarray.Dataset + Dataset containing 2D wave spectra over time. + - var : str, optional, default='SPEC' + Name of the spectral variable to analyze. + - period : tuple of two str, optional, default=None + A tuple specifying the desired time range. + - (start_time, end_time): Filters between start_time and end_time, e.g., ('2021-01-01T00', '2021-12-31T23'). + - (start_time): Filters to a single timestamp. + - None: Uses the full time range available in data. + Both start_time and end_time may be strings or datetime-like objects. + - month : int or None, optional, default=None + If set (1 = January, ..., 12 = December), plots only that month. + If None, plots all monthly spectra. + - method : str, optional, default='mean' + Aggregation method: + - 'mean' : Average over time. + - 'top_1_percent_mean' : Average over times where Hm0 ≥ 99th percentile. + - 'hm0_max' : Use time step with maximum Hm0. + - 'hm0_top3_mean' : Average of Hm0 over the three time steps with the highest values. + output_file : str, optional, default='monthly_spectra_1d.png' + Filename for the saved plot. + + Returns + ------- + - fig : matplotlib.figure.Figure + The figure object containing the plot. + ''' + + filtered_data, period_label = spec_funcs.filter_period(data=data[var], period=period) # Filters the dataset to the specified time period. + + df = spectra.table_monthly_freq_1dspectrum(filtered_data, var=var, month=month, method=method, average_over='month', output_file=None) # DataFrame containing one aggregated 1D spectrum per month and a final row for the overall mean spectrum based on the selected method. + + freq = [col for col in df.columns if col not in ['Year', 'Month', 'Month_no', 'Hm0']] # Select all columns except 'Month', 'Month_no', and 'Hm0' + spec_values = df[freq].to_numpy() + hm0 = df['Hm0'].tolist() + + fig, ax = plt.subplots() + + color_list = [ + "#051BE2", + "#079FF7", + "#77C4F1", + "#117733", + "#A3C26B", + "#00EBAC", + "#EED039", + "#E69F00", + "#D55E00", + "#CC79A7", + "#999999", + "#990303", + ] + + dashed_lines = {"#E69F00", "#079FF7", "#A3C26B"} + dotted_lines = {"#00EBAC", "#77C4F1", "#990303", "#051BE2"} + + # Build (color, linestyle) pairs + line_styles_12 = [] + for color in color_list: + linestyle = '--' if color in dashed_lines else (':' if color in dotted_lines else '-') + line_styles_12.append((color, linestyle)) + + # Title + parts = period_label.split(" to ") + year_label = parts[0][:4] if len(parts) == 1 or parts[0][:4] == parts[1][:4] else f"{parts[0][:4]}-{parts[1][:4]}" + + + ######## Plot of monthly aggregation across all years ######## + if month is None and (period is None or isinstance(period, tuple)): + months = df['Month'].tolist() + + for i in range(len(months)): + color, ls = ('black', '-') if i == len(months) -1 else line_styles_12[i] + if period is None or isinstance(period, tuple): + ax.plot(freq, spec_values[i, :], label=rf"{months[i]} ($H_{{m_0}}$={round(hm0[i], 1)} m)", color=color, linestyle=ls, linewidth=1.5, markevery=5) + ax.set_title(f'Period: {year_label}, method = {method}', fontsize=14) + + ######## Single timestamp plot ######## + elif isinstance(period, str): + ax.plot(freq, spec_values[0, :], label=rf"{period_label} ($H_{{m_0}}$={round(hm0[0], 1)} m)", color='blue', linewidth=1.5) + + ######## Plot of yearly aggregation for a specific month ######## + else: + print(df) + for i, label_val in enumerate(df['Year']): + spec = spec_values[i] + label = f"{label_val} ($H_{{m_0}}$={hm0[i]:.1f} m)" + color, ls = ('black', '-') if i == len(df['Year']) - 1 else line_styles_12[i] + ax.plot(freq, spec, linewidth=1.5, label=label, color=color, linestyle=ls) + ax.set_title(f'Period: {year_label}, month = {df.index.name}, method = {method}', fontsize=14) + + + ax.legend(fontsize=8) + ax.set_xlabel('Frequency [Hz]', fontsize=12) + ax.set_ylabel(r'Variance Density [m$^2$/Hz]', fontsize=12) + ax.set_xlim(0, None) + ax.set_ylim(0, None) + ax.tick_params(axis='both', labelsize=12) + plt.tight_layout() + + if output_file is not None: + plt.savefig(output_file, dpi=200, bbox_inches='tight') + plt.close() + return fig + +def plot_spectrum_2d(data,var='SPEC', period=None, month = None, method='mean', plot_type='pcolormesh', + freq_mask = False, radius='frequency', dir_letters=False, output_file='spectrum_2d.png'): + ''' + Plots a 2D directional spectrum for a specific month or time period. + + Parameters: + - data : xarray.Dataset + Dataset containing 2D wave spectra over time. + - var : str, optional, default='SPEC' + Name of the spectral variable to plot. + - period : tuple of str, optional, default=None + A tuple specifying the desired time range. + - (start_time, end_time): Filters between start_time and end_time, e.g., ('2021-01-01T00', '2021-12-31T23'). + - (start_time): Filters to a single timestamp. + - None: Uses the full time range available in data. + Both start_time and end_time may be strings or datetime-like objects. + - month : int or None, optional, default=None + Filter data by calendar month (1 = January, ..., 12 = December). + If None, no filtering by month is applied. + - method : str, optional, default='mean' + Aggregation method: + - 'mean' : Average over time. + - 'top_1_percent_mean' : Average over times when Hm0 ≥ 99th percentile. + - 'hm0_max' : Use time step with maximum Hm0. + - 'hm0_top3_mean' : Average of Hm0 over the three time steps with the highest values. + - plot_type : str, optional, default='pcolormesh' + 2D spectrum plotting method. Choose 'pcolormesh' or 'contour'. + - freq_mask : bool, optional, default=False + If True, mask low-energy frequencies. + - radius : str, optional, default='frequency' + Units for radial axis, period or frequency. + - dir_letters : bool, optional, default=False + Use compass labels instead of degrees. + - output_file : str, optional, default='spectrum_2d.png' + Filename for the saved plot. + + Returns + - fig: matplotlib.figure.Figure + The figure object containing the plot. + ''' + data = spec_funcs.standardize_wave_dataset(data) + + # Normalize and wrap directional spectral data to 0-360°. + spec, dirc = spec_funcs.wrap_directional_spectrum(data, var=var) + + # Filters the dataset to the specified time period. + filtered_data, period_label = spec_funcs.filter_period(data=spec, period=period) + + # Add hm0 to dataset if not present + if 'hm0' not in data: + hm0 = spec_funcs.integrated_parameters(filtered_data, var=var, params=['hm0']) + filtered_data['hm0'] = hm0 + + # Aggregate the data using the specified method + data_aggregated = spec_funcs.aggregate_spectrum(data=filtered_data, hm0=filtered_data['hm0'], method=method, month=month) + + fig, ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection="polar")) + cmap = plt.cm.hot_r + + frequencies=spec.coords['freq'] + + max_data_aggregated = np.nanmax(data_aggregated) + min_data_aggregated = np.nanmin(data_aggregated) + # If masking enabled, keep frequencies where data exceeds threshold (focusing on the most significant parts of the spectrum) + if freq_mask: + threshold = max(0.0001, min(0.03, 0.01 / max_data_aggregated))* max_data_aggregated # Set threshold dynamically: smaller (0.0001 * max) for large max values, larger (up to 0.03 * max) for small max values + threshold_val = min_data_aggregated + (threshold * (max_data_aggregated - min_data_aggregated)) # Convert to data value + frequencies = data_aggregated.freq[(data_aggregated >= threshold_val).any(dim='direction')] # Keep frequencies where data exceeds threshold in any direction + data_aggregated = data_aggregated.sel(freq=frequencies) + + # Set radial coordinate data based on chosen radius type (period or frequency) + if radius=='period': + rad_data=1/frequencies.values + elif radius=='frequency': + rad_data=frequencies.values + + # Plots 2D spectrum based on plot_type + if plot_type == 'contour': + step = np.round(max_data_aggregated / 10, 1) + levels = np.round(np.arange(0, max_data_aggregated + step, step), 1) + cp = ax.contourf(np.radians(dirc.values), rad_data, data_aggregated, cmap=cmap, levels=levels) + + elif plot_type == 'pcolormesh': + cp = ax.pcolormesh(np.radians(dirc.values), rad_data, data_aggregated, cmap=cmap, shading='auto') + + # Replace x-axis tick labels with compass directions (N, NE, E, ..., NW) if enabled + if dir_letters: + ticks_loc = ax.get_xticks().tolist() + ax.xaxis.set_major_locator(mticker.FixedLocator(ticks_loc)) + ax.set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW']) + + cbar = plt.colorbar(cp, pad=0.1, shrink=0.7) + cbar.set_label(r'Variance Density [$\mathrm{' + spec.units.replace('**', '^') + '}$]', fontsize=14) + cbar.ax.tick_params(labelsize=14) + + label = rf'$\mathbf{{Method:}}\ {method.replace("_", r"\_")}$' + '\n' if isinstance(period, tuple) else None + + # Labels + if month is None: + label2 = (rf"$\mathbf{{Period:}}$ {period_label}" + " " * 15 ) + else: + months=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'] + label2 = ( + rf"$\mathbf{{Year:}}$ {filtered_data.time.dt.year.min().item()} - {filtered_data.time.dt.year.max().item()}"+ "\n" + + rf"$\mathbf{{Month:}}$ {months[month-1]}") + + # Text based on selected method + if method == 'hm0_max': + period_label = pd.to_datetime(data_aggregated.time.values).to_pydatetime().strftime('%Y-%m-%dT%H') + 'Z' + label2 = (rf"$\mathbf{{Timestamp:}}$ {period_label}" + rf" (Hm0 = {np.round(data_aggregated['hm0'].values.item(), 1)})" + "\n") + + times = pd.to_datetime(filtered_data.time.values) + label2 += rf"$\mathbf{{Range:}}$ {(times.min().strftime('%Y-%m-%dT%H') + 'Z' if times.min() == times.max() else f'{times.min().strftime('%Y-%m-%dT%H') + 'Z'} to {times.max().strftime('%Y-%m-%dT%H') + 'Z'}')}" + + elif method == 'hm0_top3_mean': + hm0_top3_timesteps = [pd.to_datetime(t).strftime('%Y-%m-%dT%H') + 'Z' for t in data_aggregated.attrs['hm0_top3_timesteps']] + period_label = ', '.join(hm0_top3_timesteps[:3]) + + label2 = (rf"$\mathbf{{Timestamps:}}$ {period_label}" + " " * 15 + "\n") + times = pd.to_datetime(filtered_data.time.values) + label2 += rf"$\mathbf{{Range:}}$ {(times.min().strftime('%Y-%m-%dT%H') + 'Z' if times.min() == times.max() else f'{times.min().strftime('%Y-%m-%dT%H') + 'Z'} to {times.max().strftime('%Y-%m-%dT%H') + 'Z'}')}" + + label = label + label2 if label != None else label2 + + ax.text(-0.07,-0.15,label, + transform=ax.transAxes,fontsize=16 if method == 'hm0_top3_mean' else 18, ha='left') + + lat_str = ", ".join( + f"{abs(lat):.1f}°{'N' if lat >= 0 else 'S'}" # Format latitudes with N/S + for lat in np.unique(data['latitude'].round(1).values)) + + lon_str = ", ".join( + f"{abs(lon):.1f}°{'E' if lon >= 0 else 'W'}" # Format longitudes with E/W + for lon in np.unique(data['longitude'].round(1).values)) + + label_position = ( + rf"$\mathbf{{Position:}}$ Lon: {lon_str}, Lat: {lat_str}" # Position label + + "\u2007" * max(0, 27 - len(f"Lon: {lon_str}, Lat: {lat_str}")) + "\n" + ) + + ax.text(-0.07,-0.23,label_position, + transform=ax.transAxes,fontsize=16 if method == 'hm0_top3_mean' else 18, ha='left') + + # Set up polar plot style and tick appearance + ax.grid(True) + ax.set_rlabel_position(315) + ax.set_theta_zero_location('N') + ax.set_theta_direction(-1) + ax.tick_params(axis='y', labelsize=12) + ax.tick_params(axis='x', labelsize=12, pad=5) + plt.tight_layout(pad=0) + + if output_file is not None: + plt.savefig(output_file, dpi=200, bbox_inches='tight') + plt.close() + return fig + + +def plot_diana_spectrum(data, var='SPEC', period=None, month = None, method='mean', partition=False, plot_type='pcolormesh', freq_mask=False, + radius='frequency', dir_letters=False, bar='hm0', output_file='diana_spectrum.png'): + ''' + Generates a Diana plot showing the 2D wave spectrum with mean wave direction, and mean wind direction if available. + Includes 1D spectrum inset, position, period, method and Hm0. Mean swell and windsea directions are added if + partitioning is enabled (requires wind data). + + Parameters + - data : xarray.Dataset + Dataset with 2D directional-frequency wave spectra. + - var : str, optional, default='SPEC' + Spectral variable name. + - period : tuple of two str, optional, default=None + A tuple specifying the desired time range. + - (start_time, end_time): Filters between start_time and end_time, e.g., ('2021-01-01T00', '2021-12-31T23'). + - (start_time): Filters to a single timestamp. + - None: Uses the full time range available in data. + Both start_time and end_time may be strings or datetime-like objects. + - month : int, optional, default=None + Filter by month (1 = January, ..., 12 = December). If set, selects data from that month across all years. + If None, no filtering is applied. + - method : str, optional, default='mean' + Aggregation method: + - 'mean' : Average over time. + - 'top_1_percent_mean' : Average over times where Hm0 ≥ 99th percentile. + - 'hm0_max' : Use time step with maximum Hm0. + - 'hm0_top3_mean' : Average of Hm0 over the three time steps with the highest values. + - partition : bool, optional, default=False + If True, partitions the wave spectrum into wind sea and swell components, and computes the mean direction for each. + This is only applicable when wind direction data is available. + - plot_type : str, optional, default='pcolormesh' + 2D spectrum plotting method. Choose 'pcolormesh' or 'contour'. + - freq_mask : bool, optional, default=False + If True, mask low-energy frequencies. + - radius : str, optional, default='frequency' + Units for radial axis, period or frequency. + - dir_letters : bool, optional, default=False + Use compass labels instead of degrees. + - bar : str, optional, default='hm0' + Show colorbar for 'density' or significant wave height 'hm0'. + - output_file : str, optional, default='diana_spectrum.png' + Output filename for saving the figure. + + Returns + fig_main : matplotlib.figure.Figure + The generated plot figure. + ''' + + if partition and 'wind_direction' not in data: + raise ValueError("Wind data ('wind_direction') is required for partitioning. Change partition to False or add wind data.") + + valid_methods = ['mean', 'top_1_percent_mean', 'hm0_max', 'hm0_top3_mean'] + if method not in valid_methods: + raise ValueError(f"Invalid method '{method}'. Available methods are: {', '.join(valid_methods)}.") + + + ############ Data handling ############ + # Standardize the 2D wave spectrum dataset to match the WINDSURFER/NORA3 format. + data = spec_funcs.standardize_wave_dataset(data) + + # Filters the dataset to the specified time period. + filtered_data, period_label = spec_funcs.filter_period(data=data, period=period) + + if month is not None and month not in pd.to_datetime(filtered_data['time'].values).month: + raise ValueError(f"Month {month} is not present in the filtered data period ({period_label}).") + + # Add hm0 to dataset if not present + if 'hm0' not in data: + hm0 = spec_funcs.integrated_parameters(filtered_data, var=var, params=['hm0']) + filtered_data['hm0'] = hm0 + + # Normalize and wrap directional spectral data to 0-360°. + spec, dirc = spec_funcs.wrap_directional_spectrum(filtered_data, var=var) + + # Aggregate spectrum using the specified method + data_aggregated_spec = spec_funcs.aggregate_spectrum(data=spec, hm0=filtered_data['hm0'], var=var, method=method, month=month) + + # Aggregate the whole data set using the specified method + data_aggregated = spec_funcs.aggregate_spectrum(data=filtered_data, hm0=filtered_data['hm0'], var=None, method=method, month=month) + hm0_aggregated = np.round(data_aggregated['hm0'].values,1) + # Uncomment the following line to print calculated Hm0 + # print('Calculated Hm0: ', data_aggregated['hm0'].values) + + + ############ 2D Spectrum ############ + fig_main = plt.figure(figsize=(7.5, 6)) + ax_2D_spectra = fig_main.add_axes([0.25, 0.25, 0.6, 0.6], projection='polar') # [left, bottom, width, height] + + # Style preferences — change these to update all colors in the figure consistently + cmap = cm.ocean_r + color="#0463d7ff" + alpha = 0.22 + textcolor = 'black' + + frequencies=spec.coords['freq'] + + max_data_aggregated = np.nanmax(data_aggregated_spec) + min_data_aggregated = np.nanmin(data_aggregated_spec) + # If masking enabled, keep frequencies where data exceeds threshold (focusing on the most significant parts of the spectrum) + if freq_mask: + threshold = max(0.0001, min(0.03, 0.01 / max_data_aggregated))* max_data_aggregated # Set threshold dynamically: smaller (0.0001 * max) for large max values, larger (up to 0.03 * max) for small max values + threshold_val = min_data_aggregated + (threshold * (max_data_aggregated - min_data_aggregated)) # Convert to data value + frequencies = data_aggregated_spec.freq[(data_aggregated_spec >= threshold_val).any(dim='direction')] # Keep frequencies where data exceeds threshold in any direction + data_aggregated_spec = data_aggregated_spec.sel(freq=frequencies) # Filter data to those frequencies + + # Set radial coordinate data based on chosen radius type (period or frequency) + if radius=='period': + rad_data=1/frequencies.values + elif radius=='frequency': + rad_data=frequencies.values + + # Plots 2D spectrum based on plot_type + if plot_type == 'contour': + step = np.round(max_data_aggregated / 10, 1) + levels = np.round(np.arange(0, max_data_aggregated + step, step), 1) + cp = ax_2D_spectra.contourf(np.radians(dirc.values), rad_data, data_aggregated_spec, cmap=cmap, levels=levels) + + elif plot_type == 'pcolormesh': + cp = ax_2D_spectra.pcolormesh(np.radians(dirc.values), rad_data, data_aggregated_spec, cmap=cmap, shading='auto') + + # Show labels only on every second radial tick when radius is frequency + if radius == 'frequency': + ticks = ax_2D_spectra.get_yticks() + reduced_ticks = ticks[::2] # keep every 2nd tick + ax_2D_spectra.set_yticks(reduced_ticks) # Set fewer yticks + labels = [str(round(t, 2)) for t in reduced_ticks] # Set labels only on those ticks (rounded nicely) + ax_2D_spectra.set_yticklabels(labels) + + # Replace x-axis tick labels with compass directions (N, NE, E, ..., NW) if enabled + if dir_letters: + ticks_loc = ax_2D_spectra.get_xticks().tolist() + ax_2D_spectra.xaxis.set_major_locator(mticker.FixedLocator(ticks_loc)) + ax_2D_spectra.set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW']) + + # If variance selected, add variance density colorbar with label and ticks (dynamically adjusted based on plot type and data range) + if bar == 'density': + cbar_ax = fig_main.add_axes([0.87, 0.3, 0.03, 0.5]) + cbar = plt.colorbar(cp, cax=cbar_ax, pad=0.1, shrink=0.6) + cbar_ax.text(0.5, 1.05, r'Density [m$^2$/Hz]', ha='center', va='bottom', fontsize=12, transform=cbar_ax.transAxes) + cbar.ax.tick_params(labelsize=12) + + if plot_type == 'pcolormesh': + step = max(np.round(max_data_aggregated/5,1),0.1) # Define dynamic step size based on data range + ticks = np.arange(np.round(min_data_aggregated,1), np.round(max_data_aggregated+0.05,1), step) # Generate ticks from min to max with this step + cbar.set_ticks(ticks) + else: + cbar.set_ticks(levels[::2]) + + # Set up polar plot style and tick appearance + ax_2D_spectra.grid(True) + ax_2D_spectra.set_rlabel_position(315) + ax_2D_spectra.set_theta_zero_location('N') + ax_2D_spectra.set_theta_direction(-1) + ax_2D_spectra.tick_params(axis='y', labelsize=12) + ax_2D_spectra.tick_params(axis='x', labelsize=12, pad=5) + + + ############ 1D wave spectrum ############ + ax_1D_spectra = fig_main.add_axes([0.09, 0.15, 0.3, 0.3]) # [left, bottom, width, height] + + # Plot 1D wave spectrum and get the spectra dataframe + df = plot_1d_spectrum_on_ax(filtered_data, ax=ax_1D_spectra, var=var, month=month, color=color, alpha=alpha,method=method) + + # If 'hm0' bar plot is requested, create a vertical bar representing Hm0 beside the 2D spectrum + if bar == 'hm0': + cbar_ax = fig_main.add_axes([0.87, 0.3, 0.03, 0.5]) + cbar_ax.bar(0, hm0_aggregated, color=color, alpha=min(alpha+0.5,1)) + cbar_ax.set_ylim(0, 20) + if hm0_aggregated >20: + cbar_ax.set_ylim(0, 25) + cbar_ax.set_xlim(-0.5, 0.5) + cbar_ax.set_xticks([]) + cbar_ax.yaxis.tick_right() + cbar_ax.set_yticks(np.linspace(0, 20, 5)) + cbar_ax.tick_params(axis='y', labelsize=12) + cbar_ax.text(0, -1, r'H$_{m_0}$ [m]', ha='center', va='top', fontsize=12) + + ############ Wave direction arrow ############ + mean_dir_rad = spec_funcs.compute_mean_wave_direction(data=filtered_data, var=var, method=method, month=month) + + # Uncomment the following lines to print the mean wave direction + # mean_dir_deg = (450 - np.rad2deg(mean_dir_rad)) % 360 + # print('Mean wave direction (calculated):', mean_dir_deg.values) + + arrow_length = 0.35 # Total length of arrow + + dx = np.cos(mean_dir_rad) * arrow_length # Direction vector (unit length) + dy = np.sin(mean_dir_rad) * arrow_length + + x0, y0 = 0.5,0.4 + end = (x0 + dx, y0 + dy) + + ax_arrow = fig_main.add_axes([0.02, 0.55, 0.2, 0.3], frameon=False) # [left, bottom, width, height] + ax_arrow.set_axis_off() + ax_arrow.set_xlim(0, 1) + ax_arrow.set_ylim(0, 1) + ax_arrow.set_zorder(20) + + ax_arrow.annotate( + '', + xy=end, + xytext=(x0,y0), + xycoords='axes fraction', + textcoords='axes fraction', + arrowprops=dict(facecolor=textcolor, edgecolor='none', width=2, headwidth=8) + ) + + # If 'wind_direction' in data plot mean wind, swell wave and wind sea wave direction + if 'wind_direction' in data: + # Plot wind + plot_direction_arrow(ax_arrow, + direction=(filtered_data['wind_direction'].sel(height=10)), + speed_data=filtered_data['wind_speed'].sel(height=10), + hm0=filtered_data['hm0'], + method=method, + month=month, + x0=x0, y0=y0, + color=color, + label='wind', + draw_ticks=True, + style='line') + + if partition == True: + sp = spec_funcs.Spectral_Partition_wind(filtered_data, beta=1.3, method = method, month=month) + # Plot swell + plot_direction_arrow(ax_arrow, + direction=sp['mean_dir_swell_rad'], + hm0=sp['hm0'], + method=method, + month=month, + x0=x0, y0=y0, + color='green', + label='swell', + style='arrow') + + # Plot windsea + plot_direction_arrow(ax_arrow, + direction=sp['mean_dir_windsea_rad'], + hm0=sp['hm0'], + method=method, + month=month, + x0=x0, y0=y0, + color='#E69F00', + label='windsea', + style='arrow') + + + # Draw cross lines at the arrow base to show orientation relative to coordinate axes + line_len = 0.15 # Length of cross arms (in axes fraction) + angles = [np.radians(0), np.radians(90), np.radians(45), np.radians(135)] # Axes to display + + for angle in angles: + dx = np.cos(angle) * line_len + dy = np.sin(angle) * line_len + ax_arrow.plot([x0 - dx, x0 + dx], [y0 - dy, y0 + dy], color='gray', linestyle='-', linewidth=1) + + # Title of the arrow + if partition == False: + ax_arrow.text(0.05 if'wind_direction' in data else 0.2, 0.85, "Mean wave" + ('/' if'wind_direction' in data else ''), color='black', transform=ax_arrow.transAxes, fontsize=12, va='bottom') + if 'wind_direction' in data: + ax_arrow.text(0.72, 0.85, "wind", color=color, transform=ax_arrow.transAxes, fontsize=12, va='bottom') + ax_arrow.text(0.25, 0.75, "direction", color='black', transform=ax_arrow.transAxes, fontsize=12, va='bottom') + + elif partition == True: + x = 0.05 + y = 0.99 + + ax_arrow.text(x, y, "Mean wave/", color='black', transform=ax_arrow.transAxes, fontsize=12, va='bottom') + ax_arrow.text(x + 0.67, y, "swell", color='green', transform=ax_arrow.transAxes, fontsize=12, va='bottom') + ax_arrow.text(x + 0.95, y, "/", color='black', transform=ax_arrow.transAxes, fontsize=12, va='bottom') + ax_arrow.text(x + 0.1, y - 0.1, "windsea", color='#E69F00', transform=ax_arrow.transAxes, fontsize=12, va='bottom') + ax_arrow.text(x + 0.55, y - 0.1, "/", color='black', transform=ax_arrow.transAxes, fontsize=12, va='bottom') + ax_arrow.text(x + 0.6, y - 0.1, "wind", color=color, transform=ax_arrow.transAxes, fontsize=12, va='bottom') + ax_arrow.text(x + 0.2, y - 0.2, " direction", color='black', transform=ax_arrow.transAxes, fontsize=12, va='bottom') + + # Color of arrow box + rect = patches.Rectangle( + (0, 0), 1.1 if partition==True else 1, 1.1 if partition==True else 1, + transform=ax_arrow.transAxes, + facecolor=color, + alpha=alpha, + edgecolor='none' + ) + ax_arrow.add_patch(rect) + rect.set_clip_on(False) + + + ############ Title of the figure ############ + fig_main.text(0.2, 0.98, rf'$\mathbf{{Variance}}$ $\mathbf{{Density}}$ $\mathbf{{Spectrum}}$', + transform=fig_main.transFigure, + fontsize=24, + verticalalignment='top', + horizontalalignment='left', + bbox=dict(facecolor='none', edgecolor='none', boxstyle='square,pad=0.5')) + + + ############ Create label with position, period, method, and Hm0 for the plot ############ + lat_str = ", ".join( + f"{abs(lat):.1f}°{'N' if lat >= 0 else 'S'}" # Format latitudes with N/S + for lat in np.unique(data['latitude'].round(1).values) + ) + + lon_str = ", ".join( + f"{abs(lon):.1f}°{'E' if lon >= 0 else 'W'}" # Format longitudes with E/W + for lon in np.unique(data['longitude'].round(1).values) + ) + + label = ( + rf"$\mathbf{{Position:}}$ Lon: {lon_str}, Lat: {lat_str}" # Position label + + "\u2007" * max(0, 27 - len(f"Lon: {lon_str}, Lat: {lat_str}")) + "\n" + ) + + label += ( + rf"$\mathbf{{Method:}}$ {method}" # Method label + + "\u2007" * max(0, 31 - len(f"Method: {method}")) + ) + + fig_main.text(0.02, 0.025 if method == 'hm0_top3_mean' else 0.02, + label, + transform=fig_main.transFigure, + fontsize=10 if method == 'hm0_top3_mean' else 12, + color=textcolor, + horizontalalignment='left', + bbox=dict(facecolor='none', edgecolor='none', boxstyle='square,pad=0.5') + ) + + if method == 'hm0_max': + period_label = pd.to_datetime(data_aggregated_spec.time.values).to_pydatetime().strftime('%Y-%m-%dT%H') + 'Z' + elif method == 'hm0_top3_mean': + times = [pd.to_datetime(t).strftime('%Y-%m-%dT%H') + 'Z' for t in data_aggregated_spec.attrs['hm0_top3_timesteps']] + period_label = ', '.join(times[:3]) + + if month is None: + label2 = rf"$\mathbf{{Period:}}$ {period_label}" + "\n" # Show full period label if no month is selected + else: + min_year = filtered_data.time.dt.year.min().item() + max_year = filtered_data.time.dt.year.max().item() + year_str = f"{min_year}" if min_year == max_year else f"{min_year} - {max_year}" + month_name = df.loc[df['Month_no'] == month, 'Month'].values[0] + label2 = ( + rf"$\mathbf{{Year:}}$ {year_str}" # Show year range and selected month + + " " * 15 + + rf"$\mathbf{{Month:}}$ {month_name}" + "\n" + ) + label2 += (rf"$\mathbf{{H_{{m_0}}:}}$ {hm0_aggregated:.1f} m") # Hm0 label + + fig_main.text(0.4 if method == 'hm0_top3_mean' else 0.47, 0.025 if method == 'hm0_top3_mean' else 0.02, + label2, + transform=fig_main.transFigure, + fontsize=10 if method == 'hm0_top3_mean' else 12, + color=textcolor, + horizontalalignment='left', + bbox=dict(facecolor='none', edgecolor='none', boxstyle='square,pad=0.5') + ) + + + ############ Layout ############ + rect = patches.Rectangle( + (0, 0), 1, 0.1, + transform=fig_main.transFigure, + facecolor=color, + alpha=alpha, + edgecolor='none' + ) + fig_main.patches.append(rect) + + ax_2D_spectra.set_zorder(2) # Bring 2D spectrum plot to front + ax_1D_spectra.set_zorder(1) # Push 1D inset to back + + if output_file is not None: + plt.savefig(output_file, dpi=200, bbox_inches=None, pad_inches=0) + plt.close() + return fig_main + +def plot_spectra_2d(data,var='SPEC', period=None, method='monthly_mean', plot_type='pcolormesh', cbar='single', radius='frequency',dir_letters=False,output_file='spectra_2d.png'): + ''' + Creates a 12-panel plot of aggregated 2D spectra based on method. + + Parameters + - data : xarray.Dataset + Dataset with 2D directional-frequency wave spectra. + - var : str, optional, default='SPEC' + Spectral variable name. + - period : tuple of two str, optional, default=None + Time range to filter the data as (start_time, end_time), e.g., ('2021-01-01T00', '2021-12-31T23'). + If None, the full time range in the dataset is used. + Both start_time and end_time may be strings or datetime-like objects. + -method : str, optional, default='mean' + Aggregation method: + - 'monthly mean' : Average over month + - 'direction' : average of data with peak direction within the specified sector. + - plot_type : str, optional, default='pcolormesh' + 2D spectrum plotting method. Choose 'pcolormesh' or 'contour'. + - cbar : str, optional, default='single' + Number of colorbars + - 'single' : one colorbar for all 12 subplots + - 'multiple' : one colorbar per subplot + - radius : str, optional, default='frequency' + Units for radial axis, period or frequency. + - dir_letters : bool, optional, default=False + Use compass labels instead of degrees. + - output_file : str, optional, default='spectra_2d.png' + Output filename for saving the figure. + + Returns + fig_main : matplotlib.figure.Figure + The generated plot figure. + + Notes: + - If method = 'direction' each panel represents a 30° directional sector and displays the average 2D wave spectrum computed + from all time steps where the peak wave wave direction falls within that sector. The percentage + shown in each panel indicates the fraction of total spectra contributing to that sector's average. + ''' + + valid_methods = ['monthly_mean', 'direction'] + if method not in valid_methods: + raise ValueError(f"Invalid method '{method}'. Available methods are: {', '.join(valid_methods)}.") + + # Standardize the 2D wave spectrum dataset to match the WINDSURFER/NORA3 format. + data = spec_funcs.standardize_wave_dataset(data) + + # Filters the dataset to the specified time period. + if isinstance(period, str): + raise ValueError("Period must include both start and end times, e.g. ('2024-01-01T00', '2024-12-31T23').") + filtered_data, period_label = spec_funcs.filter_period(data=data, period=period) + + datspec=filtered_data[var] + frequencies = datspec.coords['freq'] + directions = datspec.coords['direction'] + + # Normalize and wrap directional spectral data to 0-360°. + spec, dirc = spec_funcs.wrap_directional_spectrum(data, var=var) + + cmap = plt.cm.hot_r + + if method == 'monthly_mean': + max_spec = np.round(spec.groupby('time.month').mean(dim='time').max().item(),1) + norm = plt.Normalize(vmin=0, vmax=max_spec) + + months=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'] + fig, axs = plt.subplots(nrows=4, ncols=3, figsize=(15, 16), subplot_kw=dict(projection="polar")) + mon=1 + for r in range(4): + for c in range(3): + data_spec=spec.sel(time=spec.time.dt.month.isin([mon])).mean(dim='time') + axs[r,c].set_theta_zero_location('N') + axs[r,c].set_theta_direction(-1) + + if radius=='period': + rad_data=1/frequencies.values + elif radius=='frequency': + rad_data=frequencies.values + + if plot_type == 'contour': + cmap = cm.hot_r(np.linspace(0,1,10)) + cmap = mcol.ListedColormap(cmap) + cp=axs[r,c].contourf(np.radians(dirc.values), rad_data, data_spec.values, cmap=cmap, vmin=0, vmax=max_spec if cbar=='single' else None) + elif plot_type == 'pcolormesh': + cp = axs[r,c].pcolormesh(np.radians(dirc.values), rad_data, data_spec.values, cmap=cmap, shading='auto', vmin=0, vmax=max_spec if cbar=='single' else None) + + if dir_letters: + ticks_loc = axs[r,c].get_xticks().tolist() + axs[r,c].xaxis.set_major_locator(mticker.FixedLocator(ticks_loc)) + axs[r,c].set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW']) + + if cbar == 'multiple': + plt.colorbar(cp, label=r'[' + re.sub(r'\*\*(\d+)', r'$^\1$', datspec.units) + ']', pad=0.15) + + axs[r,c].text(0.015,0.98,months[mon-1],transform=axs[r,c].transAxes,fontsize=14,weight='bold') + mon=mon+1 + axs[r,c].grid(True) + + + elif method == 'direction': + idir=data[var][:,:,:].argmax(dim=['freq','direction'])['direction'].values + dir_mean_spec=np.zeros((len(directions),len(frequencies),len(directions)+1)) + count=np.zeros((len(directions))) + + for i in range(np.shape(datspec)[0]): + id=int(idir[i]) + dir_mean_spec[id,:,:]=dir_mean_spec[id,:,:]+spec[i,:,:].values + count[id]=count[id]+1 + for i in range(len(directions)): + if count[i]!=0: + dir_mean_spec[i,:,:]=dir_mean_spec[i,:,:]/count[i] + else: + dir_mean_spec[i,:,:]=np.full((len(frequencies),len(directions)+1),np.nan) + + # Combines the directions into 30 deg sectors + dir_mean_spec_1=np.zeros((int(len(directions)/2),len(frequencies),len(directions)+1)) + count1=np.zeros((int(len(directions)/2))) + i1=0 + for i in range(-1,len(directions)-1,2): + dir_mean_spec_1[i1,:,:]=(dir_mean_spec[i,:,:]+dir_mean_spec[i+1,:,:])/2 + count1[i1]=count[i]+count[i+1] + i1=i1+1 + + labels=['0$^{\circ}$','30$^{\circ}$','60$^{\circ}$','90$^{\circ}$','120$^{\circ}$','150$^{\circ}$','180$^{\circ}$','210$^{\circ}$', + '240$^{\circ}$','270$^{\circ}$','300$^{\circ}$','330$^{\circ}$'] + + fig, axs = plt.subplots(nrows=4, ncols=3, figsize=(15, 16), subplot_kw=dict(projection="polar")) + i1=0 + + max_spec = np.round(np.nanmax(dir_mean_spec_1),1) + norm = plt.Normalize(vmin=0, vmax=max_spec) + for r in range(4): + for c in range(3): + axs[r,c].set_theta_zero_location('N') + axs[r,c].set_theta_direction(-1) + + if radius=='period': + rad_data=1/frequencies.values + elif radius=='frequency': + rad_data=frequencies.values + + if plot_type == 'contour': + cmap = cm.hot_r(np.linspace(0,1,10)) + cmap = mcol.ListedColormap(cmap) + cp=axs[r,c].contourf(np.radians(dirc.values), rad_data, dir_mean_spec_1[i1,:,:], cmap=cmap, vmin=0, vmax=max_spec if cbar=='single' else None) + elif plot_type == 'pcolormesh': + cp = axs[r,c].pcolormesh(np.radians(dirc.values), rad_data, dir_mean_spec_1[i1,:,:], cmap=cmap, shading='auto', vmin=0, vmax=max_spec if cbar=='single' else None) + + if dir_letters: + ticks_loc = axs[r,c].get_xticks().tolist() + axs[r,c].xaxis.set_major_locator(mticker.FixedLocator(ticks_loc)) + axs[r,c].set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW']) + + if cbar == 'multiple': + plt.colorbar(cp, label=r'[' + re.sub(r'\*\*(\d+)', r'$^\1$', datspec.units) + ']', pad=0.15) + + axs[r,c].text(0,1,labels[i1],transform=axs[r,c].transAxes,fontsize=14,weight='bold') + axs[r,c].text(1,1,str(int(np.rint(count1[i1]/np.sum(count1)*100)))+'%',transform=axs[r,c].transAxes,fontsize=14) + i1=i1+1 + axs[r,c].grid(True) + + + if cbar == 'single': + cbar = fig.add_axes([0.92, 0.15, 0.02, 0.7]) # [left, bottom, width, height] + cbar = fig.colorbar(cp, cax=cbar) if plot_type=='pcolormesh' else fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap=cmap), cax=cbar, ticks=np.linspace(0, max_spec, 11)) + cbar.ax.tick_params(labelsize=14) + cbar.ax.set_title(r'[' + re.sub(r'\*\*(\d+)', r'$^\1$', datspec.units) + ']', pad=10, fontsize=14) + + fig.text(0.1,0.07,rf'$\mathbf{{Period}}$: {period_label}',transform=fig.transFigure,fontsize=16, ha='left') + + lat_str = ", ".join( + f"{abs(lat):.1f}°{'N' if lat >= 0 else 'S'}" # Format latitudes with N/S + for lat in np.unique(data['latitude'].round(1).values)) + + lon_str = ", ".join( + f"{abs(lon):.1f}°{'E' if lon >= 0 else 'W'}" # Format longitudes with E/W + for lon in np.unique(data['longitude'].round(1).values)) + + label_position = ( + rf"$\mathbf{{Position:}}$ Lon: {lon_str}, Lat: {lat_str}" # Position label + + "\u2007" * max(0, 27 - len(f"Lon: {lon_str}, Lat: {lat_str}")) + "\n" + ) + + fig.text(0.1,0.035,label_position,transform=fig.transFigure,fontsize=16, ha='left') + + if output_file is not None: + plt.savefig(output_file, dpi=200, bbox_inches='tight') + plt.close() + return fig + + +def plot_1d_spectrum_on_ax(data, ax, var='SPEC', month=None, color='#0463d7ff', alpha=0.22, method='mean'): + ''' + Plots a 1D spectrum for a specified month on the given Matplotlib axis. Data is aggregated based on method. + If month=None the full period is selected. + + Parameters: + - data : xarray.Dataset + Dataset containing 2D wave spectra over time. + - ax : matplotlib.axes.Axes + Axes object on which to plot the spectrum. + - var : str, optional, defaul = 'SPEC' + Name of the spectral variable to plot. + - month : int, optional, default = None + Month number to filter (1 = January, ..., 12 = December). If None, the full dataset is used. + - color : str, optional, default = '#0463d7ff' + Color of the plotted line (any Matplotlib-compatible color). + - alpha : float, optional, default = 0.22 + Transparency level of the line (between 0 and 1). An offset of +0.3 is applied but capped at 1. + - method : str, optional, method = 'mean' + Aggregation method: + - 'mean' : Average over time. + - 'top_1_percent_mean' : Average over times where Hm0 ≥ 99th percentile. + - 'hm0_max' : Use time step with maximum Hm0. + - 'hm0_top3_mean' : Average top 3 time steps with highest Hm0. + + Returns + - df : pandas.DataFrame + DataFrame containing the computed 1D spectra for each month, including Hm0 values. + ''' + # DataFrame containing one aggregated 1D spectrum per month and a final row for the overall mean spectrum based on the selected method. + df = spectra.table_monthly_freq_1dspectrum(data, var=var, method=method, average_over='whole_dataset') + + if month is None: + row = df[df['Month'] == 'Average'] + else: + row = df[df['Month_no'] == month] + + # Ensure the row was found + if row.empty: + raise ValueError(f"No data found for month: {month}") + + freq = [col for col in df.columns if col not in ['Month', 'Month_no', 'Hm0']] # Select all columns except 'Month', 'Month_no', and 'Hm0' + spec_values = row[freq].values.flatten() # Extract variance density values as numpy array + freq = [float(f) for f in freq] + + ax.plot(freq,spec_values,color=color,linewidth=1.5, alpha=min(alpha + 0.55, 1.0)) + + for spine in ax.spines.values(): # Remove all spines to start clean + spine.set_visible(False) + + ax.set_xlim(left=0) # Set limits so origin is visible + ax.set_ylim(bottom=0) + + x_max = np.max(freq) + y_max = np.max(spec_values) + x_ext = max(0.05 * x_max, 0.05) # Add extension (5% of data range or a minimum margin) + y_ext = max(0.05 * y_max, 0.05) + ax.set_xlim(0, x_max + x_ext) + ax.set_ylim(0, y_max + y_ext) + + freq_step = 0.1 if x_max <= 0.6 else 0.5 # Set x-axis tick spacing: 0.1 if x_max ≤ 0.6, otherwise 0.5 to keep ticks clear + + ax.set_xticks(np.arange(0, x_max + x_ext, freq_step)) + + list_yticks=ax.get_yticks(minor=False) + ax.set_yticks(list_yticks[0:-2]) + + arrow_props = dict(arrowstyle='->', linewidth=1, color='black') + ax.annotate('', xy=(x_max + x_ext, 0), xytext=(0, 0), arrowprops=arrow_props) # Draw arrows slightly beyond max data values + ax.annotate('', xy=(0, y_max + y_ext), xytext=(0, 0), arrowprops=arrow_props) + + ax.text(x_max + x_ext * 1.05, 0, 'Frequency [Hz]', va='center', ha='left', fontsize=10) # Axis labels just beyond arrows + ax.text(0.1, y_max + y_ext * 1.05, r'Variance Density [m$^2$/Hz]', va='bottom', ha='center', fontsize=10, color=None) + + return df + +def plot_direction_arrow(ax, direction, hm0, method, month, x0, y0, color='black', label='', draw_ticks=False, + speed_data=None, arrow_length=0.25, style='arrow'): + ''' + Aggregate direction and plot it as an arrow or line on the specified axes. + + Parameters: + - ax : matplotlib.axes.Axes + Axis to plot on. + - direction : xarray.DataArray or ndarray + Direction data in oceanographic degrees (0° = North, increasing clockwise). + - hm0 : xarray.DataArray + Wave height data used as weights for aggregation. + - method : str + Aggregation method ('mean', 'top_1_percent_mean', 'hm0_max', 'hm0_top3_mean'). + - month : int or None + Month for aggregation or None for all data. + - x0, y0 : float + Base position of arrow in axes fraction coordinates (0 to 1). + - color : str, optional, default = 'balck' + Color of arrow or line. + - label : str, optional, default = '' + Label for printing mean direction. + - draw_ticks : bool, optional, default = False + Whether to draw wind speed ticks (barbs) along the arrow. + - speed_data : xarray.DataArray, optional, default = None + Wind speed data in m/s, required if draw_ticks is True. + - arrow_length : float, optional, default = 0.25 + Length of the arrow or line in axes fraction units. + - style : str, optional, default = 'arrow' + 'arrow' to plot arrow with head, or 'line' for plain line. + ''' + + # Compute mean_rad if not already aggregated: + if direction.max() > 2 * np.pi: + radians = np.deg2rad((450 - direction) % 360) # Convert meteorological degrees (0°=N, clockwise) to mathematical radians (0°=E, counter clockwise) + u = spec_funcs.aggregate_spectrum(xr.ufuncs.sin(radians), hm0=hm0, method=method, month=month) # Aggregate based on method + v = spec_funcs.aggregate_spectrum(xr.ufuncs.cos(radians), hm0=hm0, method=method, month=month) + + mean_rad = np.arctan2(u, v) + else: mean_rad = direction + + dx, dy = np.cos(mean_rad) * arrow_length, np.sin(mean_rad) * arrow_length # Calculate arrow vector components + start = (x0 + dx, y0 + dy) + end = (x0, y0) + + if style == 'arrow': + ax.annotate('', xy=start, xytext=end, + arrowprops=dict(facecolor=color, edgecolor='none', width=1, headwidth=3), + xycoords='axes fraction', zorder=5) + elif style == 'line': + ax.plot([start[0], end[0]], [start[1], end[1]], color=color, + linewidth=1, transform=ax.transAxes, zorder=5) + + # Draw wind speed ticks along the arrow if requested + if draw_ticks and speed_data is not None: + speed_knots = spec_funcs.aggregate_spectrum(speed_data, hm0=hm0, method=method, month=month).values * 1.94384 # Aggregate wind speed and convert to knots + vec = -np.array([dx, dy]) / np.linalg.norm([dx, dy]) # Unit vector along arrow direction (from arrow tip back) + perp = np.array([-vec[1], vec[0]]) # Perpendicular unit vector for ticks + + full_ticks = int(speed_knots // 10) # Number of full 10-knot ticks + half_tick = (speed_knots % 10) >= 5 # Half tick if remainder ≥ 5 knots + tick_len = 0.1 # Tick length in axes fraction units + spacing = 0.025 # Distance between ticks along arrow + + # Draw full ticks + for i in range(full_ticks): + s = np.array(start) + vec * i * spacing + e = s + perp * tick_len + ax.plot([s[0], e[0]], [s[1], e[1]], color=color, linewidth=1, transform=ax.transAxes) + + # Draw half tick if applicable + if half_tick: + s = np.array(start) + vec * full_ticks * spacing + e = s + perp * tick_len / 2 + ax.plot([s[0], e[0]], [s[1], e[1]], color=color, linewidth=1, transform=ax.transAxes) + + # Uncomment the following lines to print the mean direction + # mean_deg = (450 - np.rad2deg(mean_rad)) % 360 + # print(f"Mean {label} direction:", mean_deg.values if hasattr(mean_deg, 'values') else mean_deg) diff --git a/metocean_stats/stats/__pycache__/__init__.cpython-312.pyc b/metocean_stats/stats/__pycache__/__init__.cpython-312.pyc deleted file mode 100644 index ee8b781e..00000000 Binary files a/metocean_stats/stats/__pycache__/__init__.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/stats/__pycache__/aux_funcs.cpython-312.pyc b/metocean_stats/stats/__pycache__/aux_funcs.cpython-312.pyc deleted file mode 100644 index c3697cf8..00000000 Binary files a/metocean_stats/stats/__pycache__/aux_funcs.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/stats/__pycache__/extreme.cpython-312.pyc b/metocean_stats/stats/__pycache__/extreme.cpython-312.pyc deleted file mode 100644 index a75c19c4..00000000 Binary files a/metocean_stats/stats/__pycache__/extreme.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/stats/__pycache__/general.cpython-312.pyc b/metocean_stats/stats/__pycache__/general.cpython-312.pyc deleted file mode 100644 index da06f87e..00000000 Binary files a/metocean_stats/stats/__pycache__/general.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/stats/__pycache__/spec_funcs.cpython-312.pyc b/metocean_stats/stats/__pycache__/spec_funcs.cpython-312.pyc deleted file mode 100644 index ec8d4da5..00000000 Binary files a/metocean_stats/stats/__pycache__/spec_funcs.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/stats/spec_funcs.py b/metocean_stats/stats/spec_funcs.py index 64342ffe..3c257468 100644 --- a/metocean_stats/stats/spec_funcs.py +++ b/metocean_stats/stats/spec_funcs.py @@ -1,6 +1,7 @@ import numpy as np import xarray as xr import scipy +import pandas as pd def jonswap(f,hs,tp,gamma='fit', sigma_low=.07, sigma_high=.09): """ @@ -276,7 +277,7 @@ def interpolate_dataarray_spec( spec: xr.DataArray, new_coordinates[dir_var] = new_directions return xr.DataArray(new_spec,new_coordinates) -def integrated_parameters( +def integrated_parameters_dict( spec: np.ndarray|xr.DataArray, frequencies:np.ndarray|xr.DataArray, directions: np.ndarray|xr.DataArray) -> dict: @@ -352,4 +353,504 @@ def integrated_parameters( "peak_dir":peak_dir } - return spec_parameters \ No newline at end of file + return spec_parameters + + +def integrated_parameters(data, var='SPEC', params=None): + ''' + Calculate significant wave height (Hm0), peak frequency, and peak direction + from a directional wave spectrum. You can select which parameters to compute. + + Parameters + - data : xarray.Dataset or xarray.DataArray + Wave spectrum with dimensions including 'freq' and 'direction'. If a Dataset, + the spectral variable is specified by `var`. + - var : str, optional, default = 'SPEC' + Name of the spectral variable in `data` if `data` is a Dataset. + - params : list of str or None, optional, default = None + List of parameters to compute. Possible values: 'hm0', 'peak_freq', 'peak_dir'. + If None, all parameters are computed. + + Returns + - tuple or single xarray.DataArray + The requested parameters, in order ('hm0', 'peak_freq', 'peak_dir'). If only + one parameter is requested, returns it directly. + ''' + if params is None: + params = ['hm0', 'peak_freq', 'peak_dir'] + else: + # Validate params list + valid_params = {'hm0', 'peak_freq', 'peak_dir'} + if not set(params).issubset(valid_params): + raise ValueError(f"Invalid params {params}. Allowed values are {valid_params}") + + try: + data = data[var] + except KeyError: + pass + + directions = data['direction'] + frequencies = data['freq'] + + # Convert directions to radians if in degrees + if directions.max().item() > 2 * np.pi: + directions = np.deg2rad(directions) + + # Assign radian direction coordinates to spectrum (for correct integration) + data = data.assign_coords(direction=directions) + + results = {} + + # Compute peak_freq and/or peak_dir only if requested + if 'peak_freq' in params or 'peak_dir' in params: + spec_flat = data.stack(freq_dir=("freq", "direction")) # Combine frequency and direction dimensions into a single dimension to find global peak + peak_flat_idx = spec_flat.argmax(dim="freq_dir") # Find index of maximum spectral energy along the combined freq-direction dimension + freq_idx, dir_idx = np.unravel_index(peak_flat_idx, (len(frequencies), len(directions))) # Convert the flat index back into separate frequency and direction indices + + if 'peak_freq' in params: # Select frequencies and directions at the peak indices + results['peak_freq'] = frequencies.isel(freq=freq_idx) + if 'peak_dir' in params: + results['peak_dir'] = directions.isel(direction=dir_idx) + + # Compute hm0 only if requested + if 'hm0' in params: + spec_1d = data.integrate(coord='direction') + m0 = spec_1d.integrate(coord='freq') + results['hm0'] = 4 * np.sqrt(m0) + + # Return in fixed order, but only requested params + output_order = ['hm0', 'peak_freq', 'peak_dir'] + output = tuple(results[param] for param in output_order if param in params) + + # If only one param requested, return it directly (not tuple) + if len(output) == 1: + return output[0] + else: + return output + + +def combine_spec_wind(dataspec, datawind): + ''' + Merge spectral and wind datasets after dropping conflicting spatial coordinates ('x' and 'y'). + Both input datasets must correspond to a single geographic location. + + Parameters + - dataspec : xarray.Dataset + Spectral data. + - datawind : xarray.Dataset + Wind data. + + Returns + - xarray.Dataset + Combined dataset with variables from both inputs. + ''' + + datawind = datawind.drop_vars(['x', 'y'], errors='ignore') + dataspec = dataspec.drop_vars('x', errors='ignore') + return xr.merge([dataspec, datawind]) + + +def from_2dspec_to_1dspec(data, var='SPEC', dataframe=True, hm0=False): + ''' + Converts a 2D directional wave spectrum to a 1D frequency spectrum by integrating over direction. + Optionally includes significant wave height (Hm0) if present in the input dataset. + + Parameters + - data : xarray.Dataset or xarray.DataArray + Input containing 2D spectral data with dimensions: + - (time, frequency, direction), or + - (frequency, direction) for single time steps. + - var : str, optional, default='SPEC' + Name of the spectral variable in the dataset. Ignored if data is already a DataArray. + - dataframe : bool, optional, default=True + If True, return a pandas DataFrame; otherwise, return an xarray.DataArray. + - hm0 : bool, optional, default=False + If True, include the Hm0 column in the output DataFrame when present. + + Returns + - pandas.DataFrame or xarray.DataArray + 1D frequency spectrum. Hm0 included if requested and present. + ''' + + try: # Try to extract the 2D spectral variable from the dataset; if not found, assume input is already the DataArray + spec2d = data[var] + except KeyError: + spec2d = data + + if spec2d['direction'].max() > 2 * np.pi: # Convert direction coordinates from degrees to radians if needed (assuming max direction value > 2π) + spec2d = spec2d.assign_coords(direction=np.deg2rad(spec2d['direction'])) + + spec_1d = spec2d.integrate(coord='direction') # Integrate over the 'direction' coordinate to collapse 2D spectrum into 1D frequency spectrum + + if dataframe: + df = spec_1d.to_pandas().reset_index() # Convert xarray.DataArray to pandas DataFrame and reset index so 'time' becomes a column + df.columns.name = None + + if hm0 and isinstance(data, xr.Dataset) and 'hm0' in data: # Optionally add 'Hm0' column if requested and present in the input dataset + df['Hm0'] = data['hm0'].values + + return df + + if hm0 and isinstance(data, xr.Dataset) and 'hm0' in data: # If returning xarray.DataArray and hm0 requested, assign Hm0 as a coordinate if present + spec_1d = spec_1d.assign_coords(hm0=data['hm0']) + + return spec_1d + + +def standardize_wave_dataset(data): + ''' + Standardize a 2D wave spectrum dataset to match the WINDSURFER/NORA3 format. + + Parameters: + - data : xarray.Dataset + Input dataset to standardize. + + Returns: + - xarray.Dataset + The standardized dataset. + ''' + + # Detect NORAC product + if 'product_name' in data.attrs and data.attrs['product_name'].startswith("ww3"): + # Rename dims and vars + data = data.rename({ + 'frequency': 'freq', + 'efth': 'SPEC'}) + + return data + + +def filter_period(data, period): + ''' + Filters the dataset to a specified time period. + + Parameters + - data : xarray.Dataset or xarray.DataArray + The input dataset containing a 'time' dimension. + period : tuple or None + A tuple specifying the desired time range. + - (start_time, end_time): Filters between start_time and end_time. + - (start_time,): Filters to a single timestamp. + - None: Uses the full time range available in data. + Both start_time and end_time may be strings or datetime-like objects. + + Returns + filtered_data : xarray.Dataset or xarray.DataArray + Subset of the data within the specified time range. + period_label : str + A string label describing the filtered time period. + If start and end times are equal, only that time is returned. + ''' + + # Finds data time range + data_start = pd.to_datetime(data.time.min().values) + data_end = pd.to_datetime(data.time.max().values) + + # Uses the full dataset if period is None + if period is None: + start_time, end_time = data_start, data_end + + # Uses the choosen period range + elif isinstance(period,tuple): + start_time, end_time = pd.to_datetime(period[0]), pd.to_datetime(period[1]) + + if start_time < data_start or end_time > data_end: + raise ValueError(f"Period {start_time} to {end_time} is outside data range {data_start} to {data_end}.") + + # Uses one timestamp if choosen + else: + start_time = end_time = pd.to_datetime(period) + if start_time not in data.time.values: + raise ValueError(f"({start_time}) is outside data range {data_start} to {data_end}.") + + filtered_data = data.sel(time=slice(start_time, end_time)) + + if start_time == end_time: + period_label = f"{start_time.strftime('%Y-%m-%dT%H')}Z" + else: + period_label = f"{start_time.strftime('%Y-%m-%dT%H')}Z to {end_time.strftime('%Y-%m-%dT%H')}Z" + + return filtered_data, period_label + +def wrap_directional_spectrum(data, var='SPEC'): + ''' + Normalize and wrap directional spectral data to 0-360°. + + Parameters: + - data : xarray.Dataset or DataArray + Input spectral dataset. + - var : str, optional, default='SPEC' + Variable name. + + Returns: + - spec : xarray.DataArray + Direction-sorted and wrapped spectral data. + - dirc : xarray.DataArray + Corresponding direction coordinates. + ''' + try: + datspec=data[var] + except KeyError: + datspec=data + + directions = datspec.coords['direction'] + + # Ensure directions are 0–360 and sorted + directions = directions % 360 # Normalize direction values to [0, 360) + sorted_idx = np.argsort(directions.values) # Sort direction values in increasing order + directions = directions[sorted_idx] + datspec = datspec.isel({'direction': sorted_idx}) # Reorder spectrum data accordingly + + # Ensure continuity in the directional spectrum by appending the initial directional value (e.g., 0°) at the 360° position + if directions[0] != directions[-1]: + if directions[0] < 360: + dirc = xr.concat([directions, directions[0:1] + 360], dim='direction') # Add 360° copy of first value + else: + dirc = xr.concat([directions, directions[0:1]], dim='direction') # Fallback: just repeat first value + spec = xr.concat([datspec, datspec.isel({'direction': 0})], dim='direction') # Add same data point at the end + else: + dirc = directions + spec = datspec + return spec, dirc + + +def aggregate_spectrum(data, hm0, var='SPEC', method='mean', month=None): + ''' + Aggregate a variable over time using the specified method, with optional filtering by month. + + Parameters + - data : xarray.Dataset or xarray.DataArray + Input data to aggregate. + Can be: + - A 2D directional-frequency variable (e.g., wave spectrum with 'freq' and 'direction'), or + - A simpler variable (e.g., wind speed) without those dimensions. + If a Dataset, `var` specifies which variable to use. + If a DataArray, it's used directly. + - hm0 : xarray.DataArray + Significant wave height used for percentile-based and max-based methods. + - var : str, optional, default='SPEC' + Variable name to aggregate (ignored if `data` is a DataArray). + - method : str, optional, default='mean' + Aggregation method: + - 'mean' : Average over time. + - 'top_1_percent_mean' : Average over times where Hm0 ≥ 99th percentile. + - 'hm0_max' : Use time step with maximum Hm0. + - 'hm0_top3_mean' : Average of Hm0 over the three time steps with the highest values. + - month : int, optional, default=None + If set (1-12), filters data to that month. Otherwise, uses all available data. + + Returns + - data_aggregated : xarray.DataArray + The aggregated result with time collapsed. + ''' + + try: + data=data[var] + except KeyError: + data = data + + # Apply monthly filtering if needed + if month is not None: + time_mask = data['time'].dt.month == month + data = data.sel(time=time_mask) + hm0 = hm0.sel(time=time_mask) + + if method == 'mean': + data_aggregated = data.mean(dim='time') + + elif method == 'top_1_percent_mean': + # Get time steps where hm0 >= 99th percentile + p99 = hm0.quantile(0.99) + high_hm0_times = hm0['time'][hm0 >= p99] + + # Select those time steps from data + spec_p99 = data.sel(time=high_hm0_times) + data_aggregated = spec_p99.mean(dim='time') + + elif method == 'hm0_max': + # Find time coordinate where hm0 is maximum + hm0_max_time = hm0.idxmax(dim='time') + + # Select data at that max hm0 time step + spec_hm0_max = data.sel(time=hm0_max_time) + data_aggregated = spec_hm0_max + + elif method == 'hm0_top3_mean': + # Sort hm0 descending and get top 3 time indices + top3_times = hm0.sortby(hm0, ascending=False)['time'][:3] + + # Select those times from data + spec_top3 = data.sel(time=top3_times) + data_aggregated = spec_top3.mean(dim='time') + + # Store top 3 times as metadata + data_aggregated.attrs['hm0_top3_timesteps'] = [str(t) for t in top3_times.values] + + return data_aggregated + + +def compute_mean_wave_direction(data, var='SPEC', method='mean', month=None, hm0=None): + ''' + Compute the mean wave direction from a directional wave energy spectrum. + + This function calculates the mean wave direction dir_mean based on the + discrete approximation of the integrals: + + a = ∫∫ cos(dir) * F(freq, dir) dfreq ddir + b = ∫∫ sin(dir) * F(freq, dir) dfreq ddir + dir_mean = arctan2(b, a) + + where: + - F(freq, dir) is the spectral energy density as a function of frequency and direction. + - The integrals are approximated by summations over the frequency and + direction bins weighted by the bin widths. + + Parameters: + - data : xarray.Dataset or xarray.DataArray + Wave spectrum with dimensions including 'freq' and 'direction'. If a Dataset, + the spectral variable is specified by `var`. + - var : str, optional, default = 'SPEC' + Name of the spectral variable in `data` if `data` is a Dataset. + - method : str, optional, default = 'mean' + Aggregation method to apply over time ('mean', 'top_1_percent_mean', 'hm0_max', 'hm0_top3_mean') after computing a and b. + month : int or None, optional, default = None + If specified, selects data only for that month (1-12) before computing direction. + hm0 : xarray.DataArray or None, optional, default = None + Significant wave height time series. This is required for weighted aggregation. + It must either be passed explicitly or be present as `data[var]['hm0']`. + + + Returns: + - mean_dir_rad : xarray.DataArray + Mean wave direction in radians using the mathematical convention: + 0 = East, positive counter-clockwise (CCW). + + Notes: + - Direction angles in the input are assumed to be oceanographic (0° = North, increasing clockwise). + They are converted to mathematical radians (0 = East, increasing CCW). + - Frequency and direction bin widths are computed using gradients and used to weight the integration. + - The final direction is based on vector summation (a, b) and converted using arctangent. + - Based on the method in the WAVEWATCH III User Manual (v6.07, NOAA/NCEP, 2019). + ''' + + try: + spectrum=data[var] + except KeyError: + spectrum = data # Extract the spectrum variable (dimensions expected: ['time', 'freq', 'direction']) + + directions_rad = np.deg2rad((450 - spectrum['direction']) % 360) # Convert oceanographic directions (degrees, coming from North clockwise) + # to mathematical directions (radians, pointing to East counterclockwise) + + delta_freq = np.gradient(spectrum.freq.values) # Calculate frequency bin widths + + delta_dir = np.gradient(spectrum['direction']) # Calculate direction bin widths in degrees and convert to radians + delta_dir_rad = np.deg2rad(delta_dir) + + dfreq_2d = xr.DataArray(delta_freq, dims=['freq']) # Create DataArrays for bin widths to broadcast over spectrum dims + ddir_2d = xr.DataArray(delta_dir_rad, dims=['direction']) + + area_element = dfreq_2d.broadcast_like(spectrum) * ddir_2d.broadcast_like(spectrum) # Compute the area element dfreq ddir for each frequency-direction bin by outer product + + a = (xr.ufuncs.cos(directions_rad) * spectrum * area_element).sum(dim=['freq', 'direction']) # Compute weighted sums a and b over freq and direction dimensions + b = (xr.ufuncs.sin(directions_rad) * spectrum * area_element).sum(dim=['freq', 'direction']) + + if hm0 is None: + try: + hm0 = data['hm0'] + except KeyError: + hm0 = integrated_parameters(data=data, var=var, params=['hm0']) + + a = aggregate_spectrum(data=a, hm0=hm0, method=method, month=month) # Aggregate based on method + b = aggregate_spectrum(data=b, hm0=hm0, method=method, month=month) + + mean_dir_rad = np.arctan2(b, a) # Compute mean wave direction as arctangent of vector sum components + + return mean_dir_rad + + + + +def Spectral_Partition_wind(data, beta =1.3, method='mean', month=None): + ''' + Partition wave spectrum into wind sea and swell using the dimensionless parameter A. + + Based on the formulation by Komen et al. (1984). + A threshold of beta ≤ 1.3 is typically used to identify pure wind sea (Hasselmann et al. 1996; Voorrips et al. 1997; Bidlot 2001). + + Parameters + - data : xarray.Dataset + Must contain 'SPEC', 'freq', 'direction', 'wind_speed', 'wind_direction' (all with time). + - beta : float, optional, default = 1.3 + Partition threshold (typically ≤ 1.3). + - method : str, default 'mean' + Aggregation method for computing mean direction. + - month : int or None + If set, compute mean direction only for that month. + + Returns + - data : xarray.Dataset + Dataset with added variables for each partition: + - 'SPEC_swell', 'SPEC_windsea' + - 'Hm0_*', 'Tp_*', 'pdir_*', 'mean_dir_*' + ''' + + # Calculate the phase speed (cp) using the dispersion relation for deep water waves + data['cp'] = 9.81/(2*np.pi*data['freq']) + + # Initialize an array to store the directional difference between wave direction and wind direction + data['diff_dir'] = xr.DataArray(np.zeros((data.time.size, data.direction.size)), + dims=['time', 'direction'], + coords={'time': data.time, 'direction': data.direction}) + + # Loop over each time step to compute the angular difference between wave direction and wind direction + for i in range (len(data['diff_dir'].time)): + data['diff_dir'][i,:] = angular_difference((data.direction), ((data['wind_direction'].sel(height=10)[i].item() - 180))) + + # Calculate the dimensionless parameter A, which is used to distinguish between swell and windsea + data['A'] = beta*(data['wind_speed'].sel(height=10)/data['cp'])*np.cos(np.deg2rad(data['diff_dir'])) + + # Partition the wave spectrum into swell and windsea components based on the value of A + data['SPEC_swell'] = data['SPEC'].where(data['A']<=1,0) + data['SPEC_windsea'] = data['SPEC'].where(data['A']>1,0) + + # Estimate integrated parameters + part = ['swell','windsea'] + for k in range(len(part)): + data['Hm0_'+part[k]] = (4*(data['SPEC_'+part[k]].integrate("freq").integrate("direction"))**0.5).assign_attrs(units='m', standard_name = 'significant_wave_height_from_spectrum_'+part[k]) + data['Tp_'+part[k]] = (1/data['SPEC_'+part[k]].integrate('direction').idxmax(dim='freq')).assign_attrs(units='s', standard_name = 'peak_wave_period_'+part[k]) + data['pdir_'+part[k]] = (data['SPEC_'+part[k]].integrate('freq').idxmax(dim='direction')) + data['mean_dir_'+part[k]+'_rad'] = compute_mean_wave_direction(data = data['SPEC_' + part[k]], hm0 = data['hm0'], method = method, month=month) + data['mean_dir_'+part[k]+'_deg'] = (450 - np.rad2deg(data['mean_dir_'+part[k]+'_rad'].values)) % 360 # Meteorological convention (0° = North, clockwise) + + # Uncomment the following lines to print the mean direction + # print('Mean wave direction swell:', data['mean_dir_swell_deg'].values) + # print('Mean wave direction sea', data['mean_dir_windsea_deg'].values) + + return data + +def angular_difference(deg1, deg2): + """ + Calculate the smallest difference between two angles in degrees. + + Parameters: + deg1 (float or array-like): The first angle(s) in degrees. + deg2 (float or array-like): The second angle(s) in degrees. + + Returns: + float or np.ndarray: The smallest difference(s) between the two angles in degrees. + """ + # Convert inputs to numpy arrays for element-wise operations + deg1 = np.asarray(deg1) + deg2 = np.asarray(deg2) + + # Calculate the absolute difference + diff = np.abs(deg1 - deg2) + + # Normalize the difference to the range [0, 360) + diff = diff % 360 + + # Adjust if the difference is greater than 180 to get the shortest path + diff = np.where(diff > 180, 360 - diff, diff) + + return diff + diff --git a/metocean_stats/tables/__init__.py b/metocean_stats/tables/__init__.py index 6a57bad2..283bbeec 100644 --- a/metocean_stats/tables/__init__.py +++ b/metocean_stats/tables/__init__.py @@ -8,3 +8,4 @@ from .dir import * from .extreme import * from .general import * +from .spectra import * diff --git a/metocean_stats/tables/__pycache__/__init__.cpython-312.pyc b/metocean_stats/tables/__pycache__/__init__.cpython-312.pyc deleted file mode 100644 index 90a8157c..00000000 Binary files a/metocean_stats/tables/__pycache__/__init__.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/tables/__pycache__/dir.cpython-312.pyc b/metocean_stats/tables/__pycache__/dir.cpython-312.pyc deleted file mode 100644 index c6654227..00000000 Binary files a/metocean_stats/tables/__pycache__/dir.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/tables/__pycache__/extreme.cpython-312.pyc b/metocean_stats/tables/__pycache__/extreme.cpython-312.pyc deleted file mode 100644 index b4109779..00000000 Binary files a/metocean_stats/tables/__pycache__/extreme.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/tables/__pycache__/general.cpython-312.pyc b/metocean_stats/tables/__pycache__/general.cpython-312.pyc deleted file mode 100644 index b6656995..00000000 Binary files a/metocean_stats/tables/__pycache__/general.cpython-312.pyc and /dev/null differ diff --git a/metocean_stats/tables/spectra.py b/metocean_stats/tables/spectra.py new file mode 100644 index 00000000..3ed6f35e --- /dev/null +++ b/metocean_stats/tables/spectra.py @@ -0,0 +1,184 @@ +import pandas as pd +import numpy as np +from ..stats import spec_funcs + +def table_monthly_freq_1dspectrum(data, var='SPEC', method='mean', month=None, average_over='month', output_file='monthly_table_spectrum.csv'): + ''' + Aggregate 1D wave spectra by month or by year for a specific month, and optionally save to CSV. + + Parameters + - data : xarray.Dataset + Dataset containing 2D wave spectra over time. + - var : str, optional, default='SPEC' + Name of the spectral variable in the input data. + - month : int or None, optional, default=None + - If None: aggregate by calendar month across all years (e.g. all Januaries averaged together). + - If 1-12: aggregate by year for the given calendar month (e.g. January for each year). + - method : str, optional, default='mean' + Aggregation method: + - 'mean' : Average over time. + - 'top_1_percent_mean' : Average over times where Hm0 ≥ 99th percentile. + - 'hm0_max' : Use time step with maximum Hm0. + - 'hm0_top3_mean' : Average of Hm0 over the three time steps with the highest values. + - average_over : str, optional, default='month' + How to compute the final 'Average' row: + - 'month': takes the mean of the previously computed monthly (or specified-month) aggregates. + - 'whole_dataset': aggregates the full dataset directly using the selected method. + - output_file : str or None, default='monthly_table_spectrum.csv' + Path to save the output CSV file. If None, no file is saved. + + Returns + - pandas.DataFrame + Aggregated spectra with frequency columns and Hm0 values. + Includes one row per month/year and a final 'Average' row. + ''' + spec = spec_funcs.standardize_wave_dataset(data) # Standardize the 2D wave spectrum dataset to match the WINDSURFER/NORA3 format. # Standardize dataset + + if not isinstance(spec, pd.DataFrame): # If not a preprocessed 1D spectrum (pandas.DataFrame) + spec = spec_funcs.from_2dspec_to_1dspec(spec, var=var, dataframe=True, hm0=True) # Calculate 1D frequency spectrum by integrating over directions from a 2D directional spectrum. + + freq_cols = [c for c in spec.columns if c not in ['time', 'Hm0']] # Define frequency columns (exclude 'time' and 'Hm0') + + if 'Hm0' not in spec.columns: + hm0 = spec_funcs.integrated_parameters(data=data, var=var, params=['hm0']) + spec['Hm0'] = hm0.values + + month_map = { # Define month names, with 'Average' as a label for the overall mean + 'Jan':1,'Feb':2,'Mar':3,'Apr':4,'May':5,'Jun':6, + 'Jul':7,'Aug':8,'Sep':9,'Oct':10,'Nov':11,'Dec':12, + 'Average':0} + + + ######## Monthly aggregation across all years ######## + if month is None: + months_idx = np.unique(spec['time'].dt.month) # Extract unique month numbers from the time dimension (e.g., [2, 3, 4, ...]) + months = [k for k,v in month_map.items() if v in months_idx] # Get corresponding month names for the months present in the data + months.append('Average') + + df = pd.DataFrame(columns=freq_cols) # Create empty DataFrame to store aggregated spectra + list_hm0 = [] + + # Loop over each month, aggregate data, and append to df + for i, m in enumerate(months_idx): + df_month = spec[spec['time'].dt.month == m] # Select the subset of the data corresponding to the current month 'm' + spec_mean, hm0_mean = aggregate_group(df_month, freq_cols, method) # Aggregate the spectral data and Hm0 for this month using the chosen method + df.loc[i] = spec_mean.tolist() + list_hm0.append(hm0_mean) + + # Compute the final 'Average' row depending on the average_over method + if average_over == 'whole_dataset': # Aggregate over the entire dataset directly using the chosen method + spec_mean, hm0_mean = aggregate_group(spec, freq_cols, method=method) + avg_vals = spec_mean + avg_hm0 = hm0_mean + else: # Take the average of the previously computed yearly aggregates + if len(df) > 1: + avg_vals = df[freq_cols].mean() + avg_hm0 = np.mean(list_hm0) + else: + avg_vals = df[freq_cols].iloc[0] # If only one year present, use its values directly + avg_hm0 = list_hm0[0] + + # Create the last row + df.loc[len(df)] = pd.Series(data=list(avg_vals) + [avg_hm0], index=freq_cols + ['Hm0']) + list_hm0.append(avg_hm0) + + # Add Hm0 values and Month labels to the DataFrame + df['Hm0'] = list_hm0 + df.insert(0, 'Month', months) + df.insert(1, 'Month_no', df['Month'].map(month_map).fillna(0).astype(int)) + + + ######## Yearly aggregation for a specific month ######## + else: + if not (1 <= month <= 12): + raise ValueError("month must be between 1 and 12") + + spec_month = spec[spec['time'].dt.month == month].copy() # Filter data for the given month only + spec_month[f'Year'] = spec_month['time'].dt.year # Extract the years + years = sorted(spec_month['Year'].unique()) # Sorted list of unique years in data for that month + + df = pd.DataFrame(columns=freq_cols) # Create empty DataFrame to store aggregated spectra per year for the specified month + list_hm0 = [] # Create empty list to store aggregated hm0 per year for the specified month + + # Loop over each year, aggregate data, and append to df and list_hm0 + for i, y in enumerate(years): + df_year = spec_month[spec_month['Year'] == y] # Select data for the current year 'y' + spec_mean, hm0_mean = aggregate_group(df_year, freq_cols, method) # Aggregate the spec data and Hm0 for this year using the specified method + df.loc[i] = spec_mean.tolist() + list_hm0.append(hm0_mean) + + # Compute the final 'Average' row depending on the average_over method + if average_over == 'whole_dataset': + spec_mean, hm0_mean = aggregate_group(spec, freq_cols, method=method) # Aggregate over the entire dataset directly using the chosen method + avg_vals = spec_mean + avg_hm0 = hm0_mean + else: # Take the average of the previously computed yearly aggregates + if len(df) > 1: + avg_vals = df[freq_cols].mean() + avg_hm0 = np.mean(list_hm0) + else: + avg_vals = df[freq_cols].iloc[0] # If only one year present, use its values directly + avg_hm0 = list_hm0[0] + + # Create the last row + df.loc[len(df)] = pd.Series(data=list(avg_vals) + [avg_hm0], index=freq_cols + ['Hm0']) + list_hm0.append(avg_hm0) + + # Add Hm0 values and Year labels (including 'Average') to the DataFrame + df['Hm0'] = list_hm0 + df.insert(0, 'Year', years + ['Average']) + + # Set the index name to the specified month’s abbreviation + df.index.name = next((key for key, value in month_map.items() if value == month), None) + + # Write output CSV if requested + if output_file is not None: + if not output_file.endswith('.csv'): + output_file += '.csv' + with open(output_file, 'w') as f: + f.write('# Frequencies (in Hz) are given in the header\n') + f.write('# Spectra in m**2 s\n') + df.to_csv(f, index=False) + + return df + +def aggregate_group(df, freq_cols, method): + ''' + Aggregate frequency columns and the 'Hm0' column from a 1D frequency spectrum DataFrame based on the specified method. + + Parameters: + - df : pandas.DataFrame + Input DataFrame containing frequency columns and an 'Hm0' column. + - freq_cols : list of str + List of column names in df that contain frequency data to aggregate. + - method : str + Aggregation method to apply. Supported methods are: + - 'mean': average all rows. + - 'top_1_percent_mean': average only rows where 'Hm0' is in the top 1% quantile. + - 'hm0_max': return the row with the maximum 'Hm0' value. + - 'hm0_top3_mean': average the top 3 rows with highest 'Hm0' values. + + Returns: + - tuple of (aggregated_freq, aggregated_hm0) + - aggregated_freq: pandas Series with aggregated frequency columns. + - aggregated_hm0: float, aggregated 'Hm0' value. + ''' + + if method == 'mean': + return df[freq_cols].mean(), df['Hm0'].mean() # Average all rows for frequency columns and 'Hm0' + + elif method == 'top_1_percent_mean': + p99 = df['Hm0'].quantile(0.99) # Calculate 99th percentile of 'Hm0' + top = df[df['Hm0'] >= p99] # Filter rows where 'Hm0' >= 99th percentile + return top[freq_cols].mean(), top['Hm0'].mean() # Average frequency columns and 'Hm0' over top 1% rows + + elif method in ('hm0_max'): + max_row = df.loc[df['Hm0'].idxmax()] # Find row with maximum 'Hm0' value + return max_row[freq_cols], max_row['Hm0'] # Return frequency columns and 'Hm0' from that row + + elif method == 'hm0_top3_mean': + top3 = df.nlargest(3, 'Hm0') # Select top 3 rows with highest 'Hm0' values + return top3[freq_cols].mean(), top3['Hm0'].mean() # Average frequency columns and 'Hm0' over top 3 rows + + else: + raise ValueError(f"Unknown method: {method}") \ No newline at end of file diff --git a/tests/data/synthetic_dataset.py b/tests/data/synthetic_dataset.py new file mode 100644 index 00000000..03a33d1c --- /dev/null +++ b/tests/data/synthetic_dataset.py @@ -0,0 +1,80 @@ +import numpy as np +import pandas as pd +import xarray as xr + +def synthetic_dataset_spectra(): + ''' + Create a synthetic xarray.Dataset representing a 2D wave spectrum time series + at a single geographical point. + + The dataset contains: + - Time: Hourly timestamps over 2 years (2021-2022) + - Frequency: 30 values from 0.035 to 0.55 Hz + - Direction: 24 values from 7.5° to 352.5° + - SPEC: Randomly generated 2D wave spectral data (time x freq x direction) + - Latitude and Longitude: Constant single-point coordinates (59°N, 4°E) + + Returns + - xr.Dataset + Synthetic wave spectrum dataset. + ''' + # Define dimensions + n_time = 24 * 365 * 2 # Number of time steps: 2 years of hourly data + n_freq = 30 + n_dir = 24 + + # Create coordinate arrays + time = pd.date_range('2021-01-01', periods=n_time, freq='h') # Hourly time steps from Jan 1, 2021 + freq = np.linspace(0.03, 0.55, n_freq).astype(np.float32) # Frequency values from 0.03 to 0.55 Hz + direction = np.linspace(7.5, 352.5, n_dir).astype(np.float32) # Direction values evenly spaced in 360° range + + # Define a base 2D spectrum shape: Gaussian frequency profile × uniform direction profile + freq_profile = np.exp(-((freq - 0.15) / 0.05) ** 2).astype(np.float32) # Gaussian centered at 0.15 Hz + dir_profile = np.ones(n_dir, dtype=np.float32) # Flat directional distribution + base_spectrum = np.outer(freq_profile, dir_profile) # Create base 2D spectrum (freq × dir) + + # Create a seasonal modulation factor that varies smoothly over time (between 0.7 and 1.3) + time_hours = np.arange(n_time) + hours_per_year = 24 * 365 + seasonal_cycle = 1 + 0.3 * np.sin(2 * np.pi * time_hours / hours_per_year) + + # Initialize 3D SPEC array (time × freq × dir) with seasonal variations and random noise + np.random.seed(42) + SPEC = np.empty((n_time, n_freq, n_dir), dtype=np.float32) + + for t in range(n_time): + noise = 0.1 * np.random.randn(n_freq, n_dir).astype(np.float32) # Add small Gaussian noise + SPEC[t] = base_spectrum * seasonal_cycle[t] + noise # Apply seasonal modulation + + SPEC = np.clip(SPEC, 0, None) # Ensure all values are non-negative + + # Wrap SPEC array in a DataArray with units and metadata + spec = xr.DataArray( + SPEC, + dims=("time", "freq", "direction"), + coords={"time": time, "freq": freq, "direction": direction}, + attrs={ + "units": "m**2 s", + "long_name": "Variance density spectrum" + } + ) + + # Construct the full Dataset with SPEC and metadata + ds = xr.Dataset( + { + "SPEC": spec + }, + coords={ + "time": time, + "freq": freq, + "direction": direction, + "x": 4, # Dummy x-dimension for compatibility + "longitude": 4.0, # Constant longitude (as coordinate) + "latitude": 60.0 # Constant latitude (as coordinate) + }, + attrs={ + "title": "Simplified synthetic wave dataset, single point, 2 years", + } + ) + + return ds diff --git a/tests/test_plots.py b/tests/test_plots.py index f73f90cb..906b5eb0 100644 --- a/tests/test_plots.py +++ b/tests/test_plots.py @@ -5,12 +5,14 @@ from metocean_stats import plots from metocean_stats.stats.aux_funcs import readNora10File +from .data import synthetic_dataset # Define TimeSeries-object for NORA3 dirname = os.path.dirname(__file__) ds = readNora10File(os.path.join(dirname, 'data/NORA_test.txt')) ds_ocean = pd.read_csv(os.path.join(dirname, 'data/NorkystDA_test.csv'),comment='#',index_col=0, parse_dates=True) depth = ['0m', '1m', '2.5m', '5m', '10m', '15m', '20m', '25m', '30m', '40m', '50m', '75m', '100m', '150m', '200m'] +ds_synthetic_spectra = synthetic_dataset.synthetic_dataset_spectra() def test_plot_scatter(ds=ds): @@ -379,10 +381,12 @@ def test_plot_monthly_stats_month_xticks(): assert isinstance(fig, plt.Figure), "The output is not a Matplotlib Figure." -def test_plot_taylor_diagram(): +def test_plot_taylor_diagram(monkeypatch): + monkeypatch.setattr(plt, "show", lambda: None) # override plt.show to no-op fig = plots.taylor_diagram(ds,var_ref=['HS'],var_comp=['HS.1','HS.2'],norm_std=True,output_file="") # Check that the output is a Matplotlib Figure assert isinstance(fig, plt.Figure), "The output is not a Matplotlib Figure." + plt.close(fig) # def test_plot_environmental_contours_hs_tp(): # figures = plot_environmental_contours(ds,'HS','TP',config='DNVGL_hs_tp',save_path='total_sea_hs_tp_') @@ -397,3 +401,27 @@ def test_plot_cca_profile(): fig = plots.plot_cca_profiles(ds_ocean,var='current_speed_',month=None,return_period=10,output_file="") # Check that the output is a Matplotlib Figure assert isinstance(fig, plt.Figure), "The output is not a Matplotlib Figure." + +def test_plot_spectra_1d(): + output_file = 'test_plot_monthly_spectra_1d.png' + fig = plots.plot_spectra_1d(data=ds_synthetic_spectra, var='SPEC', period=None, month=None, method='mean', output_file=output_file) + if os.path.exists(output_file): + os.remove(output_file) + +def test_plot_spectrum_2d(): + output_file = 'test_plot_spectrum_2d.png' + fig = plots.plot_spectrum_2d(data=ds_synthetic_spectra, var='SPEC', period=None, month=None, method='mean', output_file=output_file) + if os.path.exists(output_file): + os.remove(output_file) + +def test_plot_diana_spectrum(): + output_file = 'test_plot_diana_spectrum.png' + fig = plots.plot_diana_spectrum(data=ds_synthetic_spectra, var='SPEC', period=None, month=None, method='mean', partition=False, output_file=output_file) + if os.path.exists(output_file): + os.remove(output_file) + +def test_plot_spectra_2d(): + output_file = 'test_plot_dir_mean_2dspectrum.png' + fig = plots.plot_spectra_2d(data=ds_synthetic_spectra, var='SPEC', method='monthly_mean', output_file=output_file) + if os.path.exists(output_file): + os.remove(output_file) \ No newline at end of file diff --git a/tests/test_tables.py b/tests/test_tables.py index c3ee6852..b0f4c085 100644 --- a/tests/test_tables.py +++ b/tests/test_tables.py @@ -4,12 +4,13 @@ from metocean_stats import tables from metocean_stats.stats.aux_funcs import readNora10File +from .data import synthetic_dataset # Define TimeSeries-object for NORA3 ds = readNora10File('tests/data/NORA_test.txt') ds_ocean = pd.read_csv('tests/data/NorkystDA_test.csv',comment='#',index_col=0, parse_dates=True) depth = ['0m', '1m', '2.5m', '5m', '10m', '15m', '20m', '25m', '30m', '40m', '50m', '75m', '100m', '150m', '200m'] - +ds_synthetic_spectra = synthetic_dataset.synthetic_dataset_spectra() def test_scatter_diagram(ds=ds): output_file = 'test_scatter_diagram.csv' @@ -586,4 +587,8 @@ def test_table_cca_profile(): else: raise ValueError("Shape is not correct") - +def test_table_monthly_freq_1dspectrum(): + output_file='test_table_monthly_freq_1dspectrum.csv' + df = tables.table_monthly_freq_1dspectrum(data=ds_synthetic_spectra,var='SPEC',output_file=output_file) + if os.path.exists(output_file): + os.remove(output_file) diff --git a/wave_spectrum_diana_partition.png b/wave_spectrum_diana_partition.png new file mode 100644 index 00000000..9abc5a6a Binary files /dev/null and b/wave_spectrum_diana_partition.png differ