bluemath_tk.waves package

Submodules

bluemath_tk.waves.binwaves module

bluemath_tk.waves.binwaves.generate_swan_cases(frequencies_array: ndarray = None, directions_array: ndarray = None, direction_range: tuple = (0, 360), direction_divisions: int = 24, direction_sector: tuple = None, frequency_range: tuple = (0.035, 0.5), frequency_divisions: int = 29, gamma: float = 50, spr: float = 2) → Dataset[source]

Generate the SWAN cases monocromatic wave parameters.

Parameters:

  • frequencies_array (np.ndarray, optional) – The frequencies array. If None, it is generated using frequency_range and frequency_divisions.

  • directions_array (np.ndarray, optional) – The directions array. If None, it is generated using direction_range and direction_divisions.

  • direction_range (tuple) – (min, max) range for directions in degrees.

  • direction_divisions (int) – Number of directional divisions.

  • frequency_range (tuple) – (min, max) range for frequencies in Hz.

  • frequency_divisions (int) – Number of frequency divisions.

Returns:

The SWAN monocromatic cases Dataset with coordinates freq and dir.

Return type:

xr.Dataset

bluemath_tk.waves.binwaves.plot_selected_cases_grid(frequencies: ndarray, directions: ndarray, figsize: Tuple[int, int] = (8, 8), **kwargs)[source]

Plot the selected subset parameters.

Parameters:

  • frequencies (np.ndarray) – The frequencies array.

  • directions (np.ndarray) – The directions array.

  • figsize (tuple, optional) – The figure size. Default is (8, 8).

  • **kwargs (dict, optional) – Additional keyword arguments to be passed to the pcolormesh function.

bluemath_tk.waves.binwaves.plot_selected_subset_parameters(selected_subset: DataFrame, color: str = 'blue', **kwargs) → Tuple[figure, axes][source]

Plot the selected subset parameters.

Parameters:

  • selected_subset (pd.DataFrame) – The selected subset parameters.

  • color (str, optional) – The color to use in the plot. Default is “blue”.

  • **kwargs (dict, optional) – Additional keyword arguments to be passed to the scatter plot function.

Returns:

  • plt.figure – The figure object containing the plot.

  • plt.axes – Array of axes objects for the subplots.

bluemath_tk.waves.binwaves.process_kp_coefficients(list_of_input_spectra: List[str], list_of_output_spectra: List[str]) → Dataset[source]

Process the kp coefficients from the output and input spectra.

Parameters:

  • list_of_input_spectra (List[str]) – The list of input spectra files.

  • list_of_output_spectra (List[str]) – The list of output spectra files.

Returns:

The kp coefficients Dataset in frequency and direction.

Return type:

xr.Dataset

bluemath_tk.waves.binwaves.reconstruc_spectra(offshore_spectra: Dataset, kp_coeffs: Dataset, num_workers: int = None, memory_limit: float = 0.5, chunk_sizes: dict = {'time': 24}, verbose: bool = False)[source]

Reconstruct the onshore spectra using offshore spectra and kp coefficients.

Parameters:

  • offshore_spectra (xr.Dataset) – The offshore spectra dataset.

  • kp_coeffs (xr.Dataset) – The kp coefficients dataset.

  • num_workers (int, optional) – The number of workers to use. Default is None.

  • memory_limit (float, optional) – The memory limit to use. Default is 0.5.

  • chunk_sizes (dict, optional) – The chunk sizes to use. Default is {“time”: 24}.

  • verbose (bool, optional) – Whether to print verbose output. Default is False. If False, Dask logs are suppressed. If True, Dask logs are shown.

Returns:

The reconstructed onshore spectra dataset.

Return type:

xr.Dataset

bluemath_tk.waves.calibration module

class bluemath_tk.waves.calibration.CalVal[source]

Bases: BlueMathModel

Calibrates wave data using reference data.

This class provides a framework for calibrating wave model outputs (e.g., hindcast or reanalysis) using reference data (e.g., satellite or buoy observations). It supports directionally-dependent calibration for both sea and swell components.

direction_bin_size

Size of directional bins in degrees.

Type:

int

direction_bins

Array of bin edges for directions.

Type:

np.ndarray

calibration_model

The calibration model, more details in statsmodels.api.OLS.

Type:

sm.OLS

calibrated_data

DataFrame with columns [‘Hs’, ‘Hs_CORR’, ‘Hs_CAL’] after calibration. The time domain is the same as the model data.

Type:

pd.DataFrame

calibration_params

Dictionary with ‘sea_correction’ and ‘swell_correction’ correction coefficients.

Type:

dict

property calibrated_data: DataFrame

Returns the calibrated data.

property calibration_model: OLS

Returns the calibration model.

property calibration_params: Series

Returns the calibration parameters.

correct(data: DataFrame | Dataset) → DataFrame | Dataset[source]

Apply the calibration correction to new data.

Parameters:

data (pd.DataFrame or xr.Dataset) –

Data to correct. If DataFrame, must contain columns:

  • ’Hs’, ‘Hsea’, ‘Dirsea’, ‘Hswell1’, ‘Dirswell1’, …

If xarray.Dataset, must have variable ‘efth’ and dimension ‘part’.

Returns:

Corrected data. For DataFrame, returns columns [‘Hs’, ‘Hs_CORR’] (original and corrected SWH). For Dataset, adds variables ‘corr_coeffs’ and ‘corr_efth’.

Return type:

pd.DataFrame or xr.Dataset

Notes

  • The correction is directionally dependent and uses the coefficients obtained from fit.

direction_bin_size: int = 22.5
direction_bins: ndarray = array([ 22.5,  45. ,  67.5,  90. , 112.5, 135. , 157.5, 180. , 202.5,        225. , 247.5, 270. , 292.5, 315. , 337.5, 360. ])
fit(data: DataFrame, data_longitude: float, data_latitude: float, data_to_calibrate: DataFrame, max_time_diff: int = 2) → None[source]

Calibrate the model data using reference (calibration) data.

This method matches the model data and calibration data in time, constructs directionally-binned sea and swell matrices, and fits a linear regression to obtain correction coefficients for each direction bin.

Parameters:

  • data (pd.DataFrame) –

    Model data to calibrate. Must contain columns:

    • ’Hs’ (float): Significant wave height

    • ’Hsea’ (float): Sea component significant wave height

    • ’Dirsea’ (float): Sea component mean direction (degrees)

    • ’Hswell1’, ‘Dirswell1’, … (float): Swell components (at least one required)

    The index must be datetime-like.

  • data_longitude (float) – Longitude of the model location (used for plotting and filtering).

  • data_latitude (float) – Latitude of the model location (used for plotting and filtering).

  • data_to_calibrate (pd.DataFrame) –

    Reference data for calibration. Must contain column:

    • ’Hs_CAL’ (float): Calibrated significant wave height (e.g., from satellite)

    The index must be datetime-like.

  • max_time_diff (int, optional) – Maximum time difference (in hours) allowed when matching model and calibration data. Default is 2.

Notes

  • After calling this method, the calibration parameters are stored in self.calibration_params

and the calibrated data is available in self.calibrated_data. - The calibration is directionally dependent, meaning it uses different correction coefficients for different wave directions. - The coefficients with p-values greater than 0.05 or negative values are set to 1.0, indicating no correction is applied for those directions.

plot_calibration_results() → Tuple[Figure, list][source]

Plot the calibration results, including: - Pie charts of correction coefficients for sea and swell - Scatter plots of model vs. reference (before and after correction) - Polar density plots of sea and swell wave climate

Returns:

The matplotlib Figure and a list of Axes objects for all subplots.

Return type:

Tuple[Figure, list]

validate_calibration(data_to_validate: DataFrame) → Tuple[Figure, list][source]

Validate the calibration using independent validation data.

This method compares the original and corrected model data to the validation data, both as time series and with scatter plots.

Parameters:

data_to_validate (pd.DataFrame) –

Validation data. Must contain column:

  • ’Hs_VAL’ (float): Validation significant wave height (e.g., from buoy)

The index must be datetime-like.

Returns:

The matplotlib Figure and a list of Axes objects: [time series axis, scatter (no correction), scatter (corrected)].

Return type:

Tuple[Figure, list]

bluemath_tk.waves.calibration.get_matching_times_between_arrays(times1: ndarray, times2: ndarray, max_time_diff: int) → Tuple[ndarray, ndarray][source]

Finds matching time indices between two arrays of timestamps.

For each time in times1, finds the closest time in times2 that is within max_time_diff hours. Returns the indices of matching times in both arrays.

Parameters:

  • times1 (np.ndarray) – First array of timestamps (reference times, e.g., from model data).

  • times2 (np.ndarray) – Second array of timestamps (e.g., from satellite or validation data).

  • max_time_diff (int) – Maximum time difference in hours for considering times as matching.

Returns:

Two arrays containing the indices of matching times: - First array: indices in times1 that have matches - Second array: corresponding indices in times2 that match

Return type:

Tuple[np.ndarray, np.ndarray]

Example

>>> idx1, idx2 = get_matching_times_between_arrays(
...     model_df.index.values,
...     sat_df.index.values,
...     max_time_diff=2,
... )
bluemath_tk.waves.calibration.process_imos_satellite_data(satellite_df: DataFrame, ini_lat: float, end_lat: float, ini_lon: float, end_lon: float, depth_threshold: float = -200) → DataFrame[source]

Processes IMOS satellite data for calibration.

This function filters and processes IMOS satellite altimeter data to be used as reference data for calibration (e.g., as data_to_calibrate in CalVal.fit).

Parameters:

  • satellite_df (pd.DataFrame) –

    IMOS satellite data. Must contain columns:

    • ’LATITUDE’ (float): Latitude in decimal degrees

    • ’LONGITUDE’ (float): Longitude in decimal degrees

    • ’SWH_KU_quality_control’ (float): Quality control flag for Ku-band

    • ’SWH_KA_quality_control’ (float): Quality control flag for Ka-band

    • ’SWH_KU_CAL’ (float): Calibrated significant wave height (Ku-band)

    • ’SWH_KA_CAL’ (float): Calibrated significant wave height (Ka-band)

    • ’BOT_DEPTH’ (float): Bathymetry (negative values for ocean)

  • ini_lat (float) – Minimum latitude (southern boundary) for filtering.

  • end_lat (float) – Maximum latitude (northern boundary) for filtering.

  • ini_lon (float) – Minimum longitude (western boundary) for filtering.

  • end_lon (float) – Maximum longitude (eastern boundary) for filtering.

  • depth_threshold (float, optional) – Only include points with BOT_DEPTH < depth_threshold. Default is -200.

Returns:

Filtered and processed satellite data, suitable for use as data_to_calibrate in CalVal.fit. Includes a new column ‘Hs_CAL’ (combination of Ku-band and Ka-band calibrated significant wave heights).

Return type:

pd.DataFrame

Notes

  • The returned DataFrame can be used directly as the data_to_calibrate argument in CalVal.fit.

bluemath_tk.waves.climate module

bluemath_tk.waves.greenswell module

bluemath_tk.waves.partitioning module

bluemath_tk.waves.series module

bluemath_tk.waves.series.energy_spectrum(hs: float, tp: float, gamma: float, duration: int) → List[float][source]

Calculate JONSWAP energy spectrum for given wave parameters.

This function computes the energy spectrum S(f) using the JONSWAP spectral model, which is commonly used for wind-generated waves in deep water.

Parameters:

  • hs (float) – Significant wave height in meters.

  • tp (float) – Peak period in seconds.

  • gamma (float) – JONSWAP peak enhancement factor (typically 1-7, default 3.3).

  • duration (int) – Number of frequency points for the spectrum.

Returns:

Energy spectrum values S(f) for each frequency.

Return type:

List[float]

Notes

The JONSWAP spectrum combines the Pierson-Moskowitz spectrum with a peak enhancement factor. The spectrum is defined as:

S(f) = Beta * (Hs^2) * (Tp^-4) * (f^-5) * exp(-1.25 * (Tp*f)^-4) * gamma^exp(-((Tp*f-1)^2)/(2*sigma^2))

where sigma varies based on frequency relative to peak frequency.

bluemath_tk.waves.series.series_Jonswap(waves: Dict[str, float | int]) → ndarray[source]

Generate surface elevation time series from JONSWAP spectrum.

This function generates a time series of surface elevation using the JONSWAP spectral model. The elevation is computed by summing harmonic components with random phases.

Parameters:

waves (Dict[str, Union[float, int]]) – Dictionary containing wave parameters: - ‘H’: Significant wave height (m) - ‘T’: Peak period (s) - ‘gamma’: JONSWAP peak enhancement factor - ‘warmup’: Spin-up time (s) - ‘deltat’: Time step (s) - ‘tendc’: Simulation period (s)

Returns:

2D array with shape (n_time_steps, 2) containing: - Column 0: Time values (s) - Column 1: Surface elevation values (m)

Return type:

np.ndarray

Notes

The surface elevation is generated by: 1. Computing the JONSWAP energy spectrum 2. Converting spectrum to amplitude components 3. Summing harmonic components with random phases

bluemath_tk.waves.series.series_Jonswap_bimodal(waves: Dict[str, float | int]) → ndarray[source]

Generate surface elevation time series from bimodal JONSWAP spectrum.

This function generates a time series of surface elevation using a bimodal JONSWAP spectral model, combining two wave systems with different characteristics.

Parameters:

waves (Dict[str, Union[float, int]]) – Dictionary containing wave parameters for two wave systems: - ‘Hs1’: Significant wave height of system 1 (m) - ‘Tp1’: Peak period of system 1 (s) - ‘gamma1’: JONSWAP peak enhancement factor for system 1 - ‘Hs2’: Significant wave height of system 2 (m) - ‘Tp2’: Peak period of system 2 (s) - ‘gamma2’: JONSWAP peak enhancement factor for system 2 - ‘warmup’: Spin-up time (s) - ‘deltat’: Time step (s) - ‘tendc’: Simulation period (s)

Returns:

2D array with shape (n_time_steps, 2) containing: - Column 0: Time values (s) - Column 1: Surface elevation values (m)

Return type:

np.ndarray

Notes

The bimodal spectrum is created by summing two JONSWAP spectra, each representing a different wave system (e.g., wind sea and swell).

bluemath_tk.waves.series.series_TMA(waves: Dict[str, float | int], depth: float) → ndarray[source]

Generate surface elevation time series from TMA spectrum.

This function generates a time series of surface elevation using the TMA spectral model, which extends the JONSWAP spectrum to include depth effects.

Parameters:

  • waves (Dict[str, Union[float, int]]) – Dictionary containing wave parameters: - ‘H’: Significant wave height (m) - ‘T’: Peak period (s) - ‘gamma’: JONSWAP peak enhancement factor - ‘warmup’: Spin-up time (s) - ‘deltat’: Time step (s) - ‘comptime’: Simulation period (s)

  • depth (float) – Water depth in meters.

Returns:

2D array with shape (n_time_steps, 2) containing: - Column 0: Time values (s) - Column 1: Surface elevation values (m)

Return type:

np.ndarray

Notes

The TMA spectrum modifies the JONSWAP spectrum to account for depth effects using a depth function fi that depends on wave number and depth.

bluemath_tk.waves.series.series_TMA_bimodal(waves: Dict[str, float | int], depth: float) → ndarray[source]

Generate surface elevation time series from bimodal TMA spectrum.

This function generates a time series of surface elevation using a bimodal TMA spectral model, combining two wave systems with different characteristics.

Parameters:

  • waves (Dict[str, Union[float, int]]) – Dictionary containing wave parameters for two wave systems: - ‘Hs1’: Significant wave height of system 1 (m) - ‘Tp1’: Peak period of system 1 (s) - ‘gamma1’: JONSWAP peak enhancement factor for system 1 - ‘Hs2’: Significant wave height of system 2 (m) - ‘Tp2’: Peak period of system 2 (s) - ‘gamma2’: JONSWAP peak enhancement factor for system 2 - ‘warmup’: Spin-up time (s) - ‘deltat’: Time step (s) - ‘tendc’: Simulation period (s)

  • depth (float) – Water depth in meters.

Returns:

2D array with shape (n_time_steps, 2) containing: - Column 0: Time values (s) - Column 1: Surface elevation values (m)

Return type:

np.ndarray

Notes

The bimodal spectrum is created by summing two TMA spectra, each representing a different wave system (e.g., wind sea and swell).

bluemath_tk.waves.series.series_regular_bichromatic(waves: Dict[str, float | int]) → ndarray[source]

Generate bichromatic regular wave time series.

This function generates a time series of surface elevation for bichromatic regular waves, combining two monochromatic components.

Parameters:

waves (Dict[str, Union[float, int]]) – Dictionary containing wave parameters: - ‘T1’: Period of component 1 (s) - ‘T2’: Period of component 2 (s) - ‘H’: Wave height (m) - ‘WL’: Water level (m) - ‘warmup’: Spin-up time (s) - ‘deltat’: Time step (s) - ‘tendc’: Simulation period (s)

Returns:

2D array with shape (n_time_steps, 2) containing: - Column 0: Time values (s) - Column 1: Surface elevation values (m)

Return type:

np.ndarray

Notes

The surface elevation is the sum of two cosine components: η(t) = (H/2) * cos(2πt/T1) + (H/2) * cos(2πt/T2)

bluemath_tk.waves.series.series_regular_monochromatic(waves: Dict[str, float | int]) → ndarray[source]

Generate monochromatic regular wave time series.

This function generates a time series of surface elevation for monochromatic regular waves with constant amplitude and frequency.

Parameters:

waves (Dict[str, Union[float, int]]) – Dictionary containing wave parameters: - ‘T’: Wave period (s) - ‘H’: Wave height (m) - ‘warmup’: Spin-up time (s) - ‘deltat’: Time step (s) - ‘comptime’: Simulation period (s)

Returns:

2D array with shape (n_time_steps, 2) containing: - Column 0: Time values (s) - Column 1: Surface elevation values (m)

Return type:

np.ndarray

Notes

The surface elevation follows a simple cosine function: η(t) = (H/2) * cos(2πt/T)

bluemath_tk.waves.series.waves_dispersion(T: float, h: float) → Tuple[float, float, float][source]

Solve the wave dispersion relation to calculate wavelength, wave number, and phase speed.

This function solves the linear wave dispersion relation to determine wave properties in water of finite depth.

Parameters:

  • T (float) – Wave period in seconds.

  • h (float) – Water depth in meters.

Returns:

A tuple containing: - L: Wavelength in meters - k: Wave number in radians per meter - c: Phase speed in meters per second

Return type:

Tuple[float, float, float]

Notes

The dispersion relation is solved iteratively using: L = (g*T²)/(2π) * tanh(2πh/L)

where g is gravitational acceleration (9.81 m/s²).

The wave number k = 2π/L and phase speed c = √(g*tanh(kh)/k).

bluemath_tk.waves.snakes module

bluemath_tk.waves.spectra module

bluemath_tk.waves.spectra.spectral_analysis(water_level: ndarray, delta_time: float) → Tuple[float, float, float, float][source]

Performs spectral analysis of water level variable to separate wave components.

This function uses Welch’s method to estimate the power spectral density of the water level time series and then separates the energy into different frequency bands: incident waves (IC), infragravity waves (IG), and very low frequency waves (VLF). The significant wave heights are calculated for each component.

Parameters:

  • water_level (np.ndarray) – Water level time series in meters. Should be a 1D array of measurements.

  • delta_time (float) – Time interval between consecutive samples in seconds.

Returns:

A tuple containing the significant wave heights in meters: - Hsi: Total significant wave height - HIC: Incident wave significant height (f > 0.04 Hz) - HIG: Infragravity wave significant height (0.004 Hz < f < 0.04 Hz) - HVLF: Very low frequency wave significant height (0.001 Hz < f < 0.004 Hz)

Return type:

Tuple[float, float, float, float]

Notes

The function uses the following frequency bands for wave component separation: - Incident waves (IC): f > 0.04 Hz - Infragravity waves (IG): 0.004 Hz < f < 0.04 Hz - Very low frequency waves (VLF): 0.001 Hz < f < 0.004 Hz

The significant wave height is calculated as H_s = 4 * sqrt(m0), where m0 is the zeroth moment of the spectrum.

Examples

>>> import numpy as np
>>> water_level = np.random.randn(1000)  # 1000 samples
>>> delta_time = 1.0  # 1 second intervals
>>> Hsi, HIC, HIG, HVLF = spectral_analysis(water_level, delta_time)

bluemath_tk.waves.statistics module

bluemath_tk.waves.statistics.upcrossing(series: ndarray) → Tuple[List[float], ndarray][source]

Calculate wave heights and periods using the zero-crossing method.

This function performs wave analysis by detecting zero-crossings in a surface elevation time series. It removes the mean from the series and identifies upcrossings (where the signal crosses zero from negative to positive). For each wave cycle between upcrossings, it calculates the wave height (difference between maximum and minimum) and wave period (time between crossings).

Parameters:

series (np.ndarray) – A 2D array with shape (n_samples, 2) where: - series[:, 0] contains the time values - series[:, 1] contains the surface elevation values

Returns:

A tuple containing: - Hi : List of wave heights (amplitudes) for each wave cycle - Ti : Array of wave periods (time intervals between zero crossings)

Return type:

Tuple[List[float], np.ndarray]

Notes

  • The function removes the mean from the surface elevation series before analysis

  • Wave heights are calculated as the difference between maximum and minimum values within each wave cycle

  • Wave periods are the time intervals between consecutive zero upcrossings

  • This method is commonly used in oceanography for wave analysis

Examples

>>> import numpy as np
>>> time = np.linspace(0, 10, 1000)
>>> elevation = np.sin(2 * np.pi * 0.5 * time) + 0.1 * np.random.randn(1000)
>>> series = np.column_stack([time, elevation])
>>> heights, periods = upcrossing(series)

bluemath_tk.waves.superpoint module

bluemath_tk.waves.superpoint.superpoint_calculation(stations_data: DataArray, stations_dimension_name: str, sectors_for_each_station: Dict[str, Tuple[float, float]], overlap_angle: float = 0.0) → DataArray[source]

Join multiple station spectral data for each directional sector.

Parameters:

  • stations_data (xr.DataArray) – DataArray containing spectral data for multiple stations.

  • stations_dimension_name (str) – Name of the dimension representing different stations in the DataArray.

  • sectors_for_each_station (Dict[str, Tuple[float, float]]) – Dictionary mapping each station ID to a tuple of (min_direction, max_direction) representing the directional sector for that station.

  • overlap_angle (float) – Angle in degrees that defines the overlap between two sectors.

Returns:

A new DataArray where each point is the sum of spectral data from all stations for the specified directional sector.

Return type:

xr.DataArray

Notes

If your stations_data is saved in different files, you can load them all and then concatenate them using xr.open_mfdataset function. Example below:

```python files = [

“path/to/station1.nc”, “path/to/station2.nc”, “path/to/station3.nc”

]

def load_station_data(ds: xr.Dataset) -> xr.DataArray:

return ds.efth.expand_dims(“station”)

stations_data = xr.open_mfdataset(

files, concat_dim=”station”, preprocess=load_station_data,

)

Module contents

Project: BlueMath_tk Sub-Module: waves Author: GeoOcean Research Group, Universidad de Cantabria Repository: https://github.com/GeoOcean/BlueMath_tk.git Status: Under development (Working)