From 057d2034bbb590818ab7eb20e962879760e959e9 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 15:27:53 +0000 Subject: [PATCH 1/7] Initial plan From 45e0cbf1b3ed4fd07f1bc059ebb3b09d1fb9652c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 15:31:21 +0000 Subject: [PATCH 2/7] Enhance docstrings for main API functions with comprehensive documentation Co-authored-by: SBFRF <8375832+SBFRF@users.noreply.github.com> --- pydiwasp/dirspec.py | 182 ++++++++++++++++++++++++++++++++++------- pydiwasp/infospec.py | 146 +++++++++++++++++++++++++++++---- pydiwasp/interpspec.py | 121 +++++++++++++++++++++++---- pydiwasp/plotspec.py | 127 +++++++++++++++++++++++----- pydiwasp/writespec.py | 106 +++++++++++++++++++++--- 5 files changed, 589 insertions(+), 93 deletions(-) diff --git a/pydiwasp/dirspec.py b/pydiwasp/dirspec.py index d9bad60..ca4f97f 100644 --- a/pydiwasp/dirspec.py +++ b/pydiwasp/dirspec.py @@ -21,37 +21,157 @@ def dirspec(ID, SM, EP, Options_=None): """ - DIWASP V1.4 function - dirspec: main spectral estimation routine - - [SMout,EPout]=dirspec(ID,SM,EP,{options}) - - Outputs: - SMout A spectral matrix structure containing the results - Epout The estimation parameters structure with the values actually used for the computation including any default settings. - - Inputs: - ID An instrument data structure containing the measured data - SM A spectral matrix structure; data in field SM.S is ignored. - EP The estimation parameters structure. For default values enter EP as [] - [options] options entered as cell array with parameter/value pairs: e.g.{'MESSAGE',1,'PLOTTYPE',2}; - Available options with default values: - 'MESSAGE',1, Level of screen display: 0,1,2 (increasing output) - 'PLOTTYPE',1, Plot type: 0 none, 1 3d surface, 2 polar type plot, 3 3d surface(compass angles), 4 polar plot(compass angles) - 'FILEOUT','' Filename for output file: empty string means no file output - - Input structures ID and SM are required. Either [EP] or [options] can be included but must be in order if both are included. - "help data_structures" for information on the DIWASP data structures - - All of the implemented calculation algorithms are as described by: - Hashimoto,N. 1997 "Analysis of the directional wave spectrum from field data" - In: Advances in Coastal Engineering Vol.3. Ed:Liu,P.L-F. Pub:World Scientific,Singapore - - - Original copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth - - Translated by Chuan Li and Spicer Bak, - Field Research Facility, US Army Corps of Engineers + Estimate directional wave spectrum from instrument array measurements. + + This is the main function for directional wave spectrum analysis. It takes + measurements from an array of wave instruments and estimates the distribution + of wave energy as a function of frequency and direction. + + Parameters + ---------- + ID : dict + Instrument data structure containing measured wave data. Required fields: + + * 'layout' : ndarray, shape (3, N) + Instrument positions as [x; y; z] coordinates in meters. + x, y are horizontal positions, z is vertical (typically 0 or depth). + * 'datatypes' : list of str + Type of measurement for each instrument. Valid types: + 'pres' (pressure), 'elev' (surface elevation), + 'velx' (x-velocity), 'vely' (y-velocity), + 'vels' (horizontal velocity magnitude), 'accs' (acceleration). + * 'depth' : float + Water depth in meters. + * 'fs' : float + Sampling frequency in Hz. + * 'data' : ndarray, shape (nsamples, ninstruments) + Time series measurements from each instrument. + + SM : dict + Spectral matrix structure defining the output grid. Required fields: + + * 'freqs' : ndarray + Frequency values for the output spectrum (Hz or rad/s). + * 'dirs' : ndarray + Direction values for the output spectrum (radians or degrees). + + Optional fields: + + * 'funit' : str, optional + Frequency units: 'Hz' (default) or 'rad/s'. + * 'dunit' : str, optional + Direction units: 'rad' (default) or 'deg'. + * 'xaxisdir' : float, optional + Direction of the x-axis in compass degrees (default: 90 = East). + + EP : dict + Estimation parameters. Use empty dict {} for default values. Fields: + + * 'method' : str, optional + Estimation method: 'IMLM' (default) or 'EMEP'. + * 'nfft' : int, optional + FFT length for spectral estimation. Auto-calculated if not provided. + * 'dres' : int, optional + Directional resolution (number of direction bins, default: 180). + * 'iter' : int, optional + Number of iterations for iterative methods (default: 100). + * 'smooth' : str, optional + Spectral smoothing: 'ON' (default) or 'OFF'. + + Options_ : list, optional + Optional parameters as alternating name/value pairs. + Example: ['MESSAGE', 1, 'PLOTTYPE', 2, 'FILEOUT', 'output.txt'] + + Available options: + + * 'MESSAGE' : int, default 1 + Level of console output (0=none, 1=normal, 2=verbose). + * 'PLOTTYPE' : int, default 1 + Plot type: 0=none, 1=3D surface, 2=polar, + 3=3D compass, 4=polar compass. + * 'FILEOUT' : str, default '' + Output filename (empty string = no file output). + + Returns + ------- + SMout : dict + Output spectral matrix with computed directional spectrum. Contains all + input SM fields plus: + + * 'S' : ndarray, shape (nfreqs, ndirs) + Spectral density values [m²/Hz/degree or m²s/rad]. + + EPout : dict + Estimation parameters actually used, including any default values that + were applied. + + Examples + -------- + Basic example with three pressure sensors in a triangular array: + + >>> import numpy as np + >>> from pydiwasp import dirspec + >>> + >>> # Define instrument array + >>> ID = { + ... 'layout': np.array([[0, 10, 5], [0, 0, 8.66], [0, 0, 0]]), + ... 'datatypes': ['pres', 'pres', 'pres'], + ... 'depth': 10.0, + ... 'fs': 2.0, + ... 'data': wave_measurements # shape: (nsamples, 3) + ... } + >>> + >>> # Define output frequency-direction grid + >>> SM = { + ... 'freqs': np.linspace(0.05, 0.5, 50), # 0.05 to 0.5 Hz + ... 'dirs': np.linspace(-180, 180, 36) # -180 to 180 degrees + ... } + >>> + >>> # Use default estimation parameters (IMLM method) + >>> EP = {'method': 'IMLM', 'iter': 100} + >>> + >>> # Compute directional spectrum + >>> SMout, EPout = dirspec(ID, SM, EP) + >>> + >>> # Access the directional spectrum + >>> spectrum = SMout['S'] # shape: (50, 36) + + With custom options: + + >>> # Suppress output and skip plotting + >>> options = ['MESSAGE', 0, 'PLOTTYPE', 0] + >>> SMout, EPout = dirspec(ID, SM, EP, options) + + Notes + ----- + The implemented estimation algorithms are described in: + + Hashimoto, N. (1997). "Analysis of the directional wave spectrum from field data". + In: Advances in Coastal Engineering Vol. 3. Ed: Liu, P.L-F. + Pub: World Scientific, Singapore. + + The function automatically: + - Detects and removes trends from input data + - Calculates cross-spectral densities between all instrument pairs + - Computes wave numbers using linear wave theory + - Calculates instrument transfer functions + - Applies the selected estimation method + - Interpolates results onto the requested output grid + - Optionally smooths the spectrum + + See Also + -------- + infospec : Calculate wave statistics from directional spectrum + plotspec : Plot directional spectrum + interpspec : Interpolate spectrum to different grid + + References + ---------- + Original MATLAB version: Copyright (C) 2002 Coastal Oceanography Group, + CWR, UWA, Perth + + Python translation: Chuan Li and Spicer Bak, Field Research Facility, + US Army Corps of Engineers """ Options = {'MESSAGE':1, 'PLOTTYPE':1, 'FILEOUT':''} diff --git a/pydiwasp/infospec.py b/pydiwasp/infospec.py index 6cdd515..ad00648 100644 --- a/pydiwasp/infospec.py +++ b/pydiwasp/infospec.py @@ -3,27 +3,96 @@ def infospec(SM): """ - DIWASP V1.4 function - infospec: calculates and displays information about a directional spectrum + Calculate and display wave statistics from a directional spectrum. - [Hsig,Tp,DTp,Dp]=infospec(SM) + This function computes key wave parameters from a directional spectrum, + including significant wave height, peak period, and dominant directions. + Results are printed to the console. - Outputs: - Hsig Signficant wave height - Tp Peak period - DTp Direction of spectral peak - Dp Dominant direction + Parameters + ---------- + SM : dict + Spectral matrix structure containing the directional spectrum. + Required fields: + + * 'freqs' : ndarray + Frequency values (Hz or rad/s). + * 'dirs' : ndarray + Direction values (radians or degrees). + * 'S' : ndarray, shape (nfreqs, ndirs) + Spectral density values. + + Optional fields: + + * 'xaxisdir' : float, optional + Direction of x-axis in compass degrees (default: 90 = East). + Used for converting to compass bearings. - Inputs: - SM A spectral matrix structure containing the file data + Returns + ------- + Hsig : float + Significant wave height in meters. Defined as 4 times the standard + deviation of the sea surface elevation, calculated from the integrated + spectrum: Hsig = 4 * sqrt(m0), where m0 is the zeroth moment. + + Tp : float + Peak period in seconds. The period corresponding to the frequency with + the maximum energy in the frequency spectrum (1D spectrum integrated + over all directions). + + DTp : float + Direction of the spectral peak in the units specified by SM['dunit']. + The direction at which the 2D spectrum has its maximum value. + This is the direction of waves at the peak frequency. + + Dp : float + Dominant direction in the units specified by SM['dunit']. + The direction with the highest integrated energy across all frequencies. + This represents the overall dominant wave propagation direction. - Hsig is the significant wave height. Tp is the peak frequency, the highest point in the one dimensional spectrum. - DTp is the main direction of the peak period (i.e the highest point in the two-dimensional directional spectrum). - Dp is the dominant direction defined as the direction with the highest energy integrated over all frequencies. + Examples + -------- + Calculate wave statistics from a computed spectrum: - "help data_structures" for information on the DIWASP data structures - - Copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth + >>> from pydiwasp import dirspec, infospec + >>> + >>> # After computing spectrum with dirspec + >>> SMout, EPout = dirspec(ID, SM, EP) + >>> + >>> # Get wave statistics + >>> Hsig, Tp, DTp, Dp = infospec(SMout) + Infospec:: + Significant wave height: 2.5 + Peak period: 10.0 + Direction of peak period: 45.0 axis angle / 45.0 compass bearing + Dominant direction: 50.0 axis angle / 40.0 compass bearing + >>> + >>> print(f"Wave height: {Hsig:.2f} m") + Wave height: 2.50 m + >>> print(f"Peak period: {Tp:.1f} s") + Peak period: 10.0 s + + Notes + ----- + The function prints results to stdout with both axis angles and compass + bearings. The conversion between these depends on the 'xaxisdir' field + in the SM structure. + + Axis angles are measured counter-clockwise from the x-axis. + Compass bearings are measured clockwise from North (0°). + + The significant wave height Hsig is approximately equal to the average + height of the highest one-third of waves in a wave record. + + See Also + -------- + dirspec : Compute directional spectrum + compangle : Convert between axis angles and compass bearings + + References + ---------- + Original MATLAB version: Copyright (C) 2002 Coastal Oceanography Group, + CWR, UWA, Perth """ H = hsig(SM) @@ -48,4 +117,49 @@ def infospec(SM): return H, Tp, DTp, Dp def compangle(dirs, xaxisdir): + """ + Convert between axis angles and compass bearings. + + Converts directional angles from a mathematical coordinate system + (counter-clockwise from x-axis) to nautical compass bearings + (clockwise from North). + + Parameters + ---------- + dirs : float or ndarray + Direction(s) in degrees, measured counter-clockwise from the x-axis. + + xaxisdir : float + Compass bearing of the x-axis in degrees. Typically 90 (East). + + Returns + ------- + bearings : float or ndarray + Compass bearing(s) in degrees (0-360), measured clockwise from North. + + Examples + -------- + >>> from pydiwasp import compangle + >>> + >>> # Convert 45° axis angle to compass bearing (x-axis points East) + >>> bearing = compangle(45, 90) + >>> print(bearing) + 45.0 + >>> + >>> # Convert multiple directions + >>> import numpy as np + >>> dirs = np.array([0, 45, 90, 180, 270]) + >>> bearings = compangle(dirs, 90) + >>> print(bearings) + [90. 45. 0. 270. 180.] + + Notes + ----- + The conversion formula is: bearing = (180 + xaxisdir - dirs) mod 360 + + This accounts for: + - The 180° difference between "direction to" and "direction from" + - The rotation from mathematical (CCW from x) to nautical (CW from N) + - The orientation of the x-axis + """ return (180 + xaxisdir * np.ones(np.shape(dirs)) - dirs) % 360 diff --git a/pydiwasp/interpspec.py b/pydiwasp/interpspec.py index b229b58..9ea0eb1 100644 --- a/pydiwasp/interpspec.py +++ b/pydiwasp/interpspec.py @@ -6,22 +6,111 @@ def interpspec(SMin, SMout, method='linear'): """ - DIWASP V1.4 function - interpspec: interpolates between spectral matrix bases - - SMout=interpspec(SMin,SMout) - - Outputs: - SMout Output spectral matrix structure with interpolated power density - - Inputs: - SMin A spectral matrix structure containing the original spectra - SMout A spectral matrix defining the new spectral matrix - - SMout only needs to have the frequency and directional axes defined - - spectral density ignored - - "help data_structures" for information on the DIWASP data structures + Interpolate directional spectrum onto a different frequency/direction grid. + + This function resamples a directional spectrum from one frequency-direction + grid to another. It uses 2D interpolation in frequency-direction space and + preserves the total wave energy (significant wave height) to within 2%. + + Parameters + ---------- + SMin : dict + Input spectral matrix containing the original spectrum. Required fields: + + * 'freqs' : ndarray + Frequency values (Hz or rad/s). + * 'dirs' : ndarray + Direction values (radians or degrees). + * 'S' : ndarray, shape (nfreqs, ndirs) + Spectral density values to be interpolated. + + SMout : dict + Output spectral matrix defining the target grid. Required fields: + + * 'freqs' : ndarray + Target frequency values (Hz or rad/s). + * 'dirs' : ndarray + Target direction values (radians or degrees). + + Note: The 'S' field in SMout is ignored if present; it will be + filled with interpolated values. + + method : str, optional + Interpolation method passed to scipy.interpolate.griddata. + Options are: + + * 'linear' : Linear interpolation (default, recommended) + * 'nearest' : Nearest-neighbor interpolation + * 'cubic' : Cubic interpolation (smoother but slower) + + Returns + ------- + SMout : dict + Output spectral matrix with interpolated spectrum. Contains all input + SMout fields with 'S' field filled with interpolated spectral density. + + Examples + -------- + Interpolate spectrum to a coarser grid: + + >>> from pydiwasp import dirspec, interpspec + >>> import numpy as np + >>> + >>> # Compute high-resolution spectrum + >>> SM_hires = { + ... 'freqs': np.linspace(0.05, 0.5, 100), + ... 'dirs': np.linspace(-180, 180, 72) + ... } + >>> SMout_hires, EPout = dirspec(ID, SM_hires, EP) + >>> + >>> # Define coarser output grid + >>> SM_coarse = { + ... 'freqs': np.linspace(0.05, 0.5, 25), + ... 'dirs': np.linspace(-180, 180, 36) + ... } + >>> + >>> # Interpolate to coarse grid + >>> SMout_coarse = interpspec(SMout_hires, SM_coarse, method='linear') + + Interpolate to different frequency range: + + >>> # Zoom in on specific frequency range + >>> SM_zoom = { + ... 'freqs': np.linspace(0.1, 0.2, 50), # Focus on 0.1-0.2 Hz + ... 'dirs': np.linspace(-180, 180, 36) + ... } + >>> SMout_zoom = interpspec(SMout_hires, SM_zoom) + + Notes + ----- + The interpolation is performed in a transformed coordinate system where + each point is represented by (f*sin(θ), f*cos(θ)) to properly handle + the periodic nature of directions and the coupling between frequency + and direction in wave spectra. + + The function automatically checks that the interpolated spectrum preserves + the significant wave height to within 2%. If the error is larger, a + warning is issued suggesting that the output grid may be too coarse. + + NaN values in the interpolated result (typically at edges) are set to zero. + + If the input and output grids are identical, no interpolation is performed + and a warning is issued. + + Warnings + -------- + If the significant wave height changes by more than 2% during interpolation, + the function warns that the output grid may be too coarse. Consider using + a finer frequency or direction resolution in SMout. + + See Also + -------- + dirspec : Compute directional spectrum + + References + ---------- + Original MATLAB version: Copyright (C) 2002 Coastal Oceanography Group, + CWR, UWA, Perth """ Hs1 = hsig(SMin) diff --git a/pydiwasp/plotspec.py b/pydiwasp/plotspec.py index 1adb393..7dde03e 100644 --- a/pydiwasp/plotspec.py +++ b/pydiwasp/plotspec.py @@ -4,30 +4,117 @@ def plotspec(SM, ptype): """ - DIWASP V1.4 function - plotspec: plots the spectral matrix in 3D or polar form + Plot directional wave spectrum in 3D or polar form. - plotspec(SM,ptype) + Creates visualizations of the directional spectrum showing how wave energy + is distributed across frequencies and directions. Supports multiple plot + types including 3D surface plots and polar contour plots. - Inputs: - SM A spectral matrix structure - ptype plot type: - 1 3D surface plot - 2 polar type plot - 3 3D surface plot (compass bearing angles direction from) - 4 polar type plot (compass bearing angles direction from) + Parameters + ---------- + SM : dict + Spectral matrix structure containing the directional spectrum. + Required fields: + + * 'freqs' : ndarray + Frequency values (Hz or rad/s). + * 'dirs' : ndarray + Direction values (radians or degrees). + * 'S' : ndarray, shape (nfreqs, ndirs) + Spectral density values. + + Optional fields: + + * 'xaxisdir' : float, optional + Direction of x-axis in compass degrees (default: 90 = East). + Used for plot types 3 and 4. + * 'funit' : str, optional + Frequency units: 'Hz' or 'rad/s'. + * 'dunit' : str, optional + Direction units: 'rad' or 'deg'. + + ptype : int + Plot type selection: + + * 1 : 3D surface plot with mathematical angles + - X-axis: frequency (Hz) + - Y-axis: direction (degrees, -180 to 180) + - Z-axis: spectral density + - Directions measured counter-clockwise from x-axis + + * 2 : Polar contour plot with mathematical angles + - Angle: direction (radians) + - Radius: frequency (Hz) + - Contours: spectral density levels + - Directions measured counter-clockwise from x-axis + + * 3 : 3D surface plot with compass bearings + - X-axis: frequency (Hz) + - Y-axis: compass bearing (0-360°) + - Z-axis: spectral density + - Bearings are "direction from" (oceanographic convention) + + * 4 : Polar contour plot with compass bearings + - Angle: compass bearing (radians) + - Radius: frequency (Hz) + - Contours: spectral density levels + - Bearings are "direction from" (oceanographic convention) - The 3D surface plot type is a MATLAB surface plot with SM.freqs on the x axis, SM.dirs on the y axis and the spectral density, SM.S as the z value. - The polar type plot is a MATLAB polar plot with the direction showing values in SM.dirs, the radius showing values in SM.freqs - and contours representing the spectral density, SM.S. An example of the polar type plot is shown on the front cover of the manual. - For plot types 1 and 2, the direction is the direction of propagation relative to the Cartesian axis. - For options 3 and 4 the direction is coming from as a true compass bearing (this has changed from previous versions). - Directions are corrected internally from the SM.xaxisdir and SM.dunit - fields that define the orientation of the axes and directional units in the spectral matrix. + Returns + ------- + None + The function displays a matplotlib figure. Use plt.show() to display + or plt.savefig() to save the figure. - "help data_structures" for information on the DIWASP data structures - - Copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth + Examples + -------- + Create different visualizations of a directional spectrum: + + >>> from pydiwasp import dirspec, plotspec + >>> import matplotlib.pyplot as plt + >>> + >>> # Compute spectrum + >>> SMout, EPout = dirspec(ID, SM, EP) + >>> + >>> # 3D surface plot (mathematical angles) + >>> plotspec(SMout, 1) + >>> plt.show() + >>> + >>> # Polar contour plot (compass bearings) + >>> plotspec(SMout, 4) + >>> plt.show() + >>> + >>> # Create multiple subplots + >>> fig = plt.figure(figsize=(12, 5)) + >>> plt.subplot(121) + >>> plotspec(SMout, 1) + >>> plt.subplot(122) + >>> plotspec(SMout, 2) + >>> plt.tight_layout() + >>> plt.show() + + Notes + ----- + The 3D surface plots provide a clear view of the spectral shape and are + good for presentations. The polar plots are more traditional in wave + analysis and are compact for showing directional spreading. + + For plot types 1 and 2, directions are in mathematical convention + (counter-clockwise from x-axis). For types 3 and 4, directions are + nautical compass bearings representing "direction from" (where waves + are coming from), which is the oceanographic convention. + + The spectral density units shown are m²s/deg for the 3D plots. + + See Also + -------- + dirspec : Compute directional spectrum + infospec : Calculate wave statistics + + References + ---------- + Original MATLAB version: Copyright (C) 2002 Coastal Oceanography Group, + CWR, UWA, Perth """ fig = plt.figure(tight_layout=True) diff --git a/pydiwasp/writespec.py b/pydiwasp/writespec.py index c573f32..87603e3 100644 --- a/pydiwasp/writespec.py +++ b/pydiwasp/writespec.py @@ -2,20 +2,106 @@ def writespec(SM, filename): """ - DIWASP V1.4 function - writespec: writes spectrum matrix to file using DIWASP format + Write directional spectrum to file in DIWASP format. - writespec(SM,filename) + Exports a spectral matrix structure to a text file using the standard + DIWASP format. This format can be read by other DIWASP tools and is + useful for archiving results or sharing data. - Inputs: - SM A spectral matrix structure - filename String containing the filename including file extension if required + Parameters + ---------- + SM : dict + Spectral matrix structure to write. Required fields: + + * 'freqs' : ndarray + Frequency values (Hz or rad/s). + * 'dirs' : ndarray + Direction values (radians or degrees). + * 'S' : ndarray, shape (nfreqs, ndirs) + Spectral density values. + * 'xaxisdir' : float + Direction of x-axis in compass degrees (typically 90 = East). + + filename : str + Output filename including path and extension if desired. + Common extensions are .txt or .spc. - All inputs required + Returns + ------- + None + The function writes to a file and returns nothing. - "help data_structures" for information on the DIWASP data structures - - Copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth + File Format + ----------- + The DIWASP file format is a simple text format with the following structure:: + + xaxisdir + nfreqs + ndirs + freq1 + freq2 + ... + freqN + dir1 + dir2 + ... + dirM + 999 (separator) + S[0,0] + S[0,1] + ... + S[N-1,M-1] + + All values are written as floating point numbers, one per line. + The spectral density values S are written in row-major order (frequency varies fastest). + + Examples + -------- + Save a computed spectrum to file: + + >>> from pydiwasp import dirspec, writespec + >>> + >>> # Compute spectrum + >>> SMout, EPout = dirspec(ID, SM, EP) + >>> + >>> # Write to file + >>> writespec(SMout, 'output_spectrum.txt') + >>> + >>> # Or with full path + >>> import os + >>> output_dir = '/path/to/results' + >>> filename = os.path.join(output_dir, 'spectrum_20240101.spc') + >>> writespec(SMout, filename) + + Save multiple spectra from a time series: + + >>> import numpy as np + >>> + >>> # Process multiple time windows + >>> for i, data_window in enumerate(time_windows): + ... ID['data'] = data_window + ... SMout, EPout = dirspec(ID, SM, EP) + ... writespec(SMout, f'spectrum_{i:03d}.txt') + + Notes + ----- + The file format is plain text and can be inspected with any text editor. + The separator value 999 marks the transition from metadata to spectral data. + + Only the real part of the spectral density is written. Any imaginary + components are discarded. + + The saved file can be read back using custom parsing or by other DIWASP- + compatible software. + + See Also + -------- + dirspec : Compute directional spectrum + + References + ---------- + Original MATLAB version: Copyright (C) 2002 Coastal Oceanography Group, + CWR, UWA, Perth """ nf = np.max(SM['freqs'].shape) From 8d0e1f775c2a2918d2422f3dc373fcd91dcadb1a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 15:38:47 +0000 Subject: [PATCH 3/7] Add comprehensive Sphinx documentation structure Co-authored-by: SBFRF <8375832+SBFRF@users.noreply.github.com> --- .gitignore | 1 + docs/Makefile | 20 ++ docs/api/dirspec.rst | 4 + docs/api/index.rst | 147 ++++++++ docs/api/infospec.rst | 4 + docs/api/interpspec.rst | 4 + docs/api/plotspec.rst | 4 + docs/api/writespec.rst | 4 + docs/conf.py | 69 ++++ docs/developer/contributing.rst | 303 ++++++++++++++++ docs/developer/structure.rst | 373 ++++++++++++++++++++ docs/examples.rst | 326 ++++++++++++++++++ docs/index.rst | 103 ++++++ docs/installation.rst | 82 +++++ docs/license.rst | 142 ++++++++ docs/make.bat | 35 ++ docs/quickstart.rst | 157 +++++++++ docs/references.rst | 231 +++++++++++++ docs/requirements.txt | 3 + docs/user_guide/estimation_methods.rst | 282 +++++++++++++++ docs/user_guide/index.rst | 12 + docs/user_guide/instrument_types.rst | 264 ++++++++++++++ docs/user_guide/troubleshooting.rst | 402 ++++++++++++++++++++++ docs/user_guide/understanding_spectra.rst | 167 +++++++++ 24 files changed, 3139 insertions(+) create mode 100644 docs/Makefile create mode 100644 docs/api/dirspec.rst create mode 100644 docs/api/index.rst create mode 100644 docs/api/infospec.rst create mode 100644 docs/api/interpspec.rst create mode 100644 docs/api/plotspec.rst create mode 100644 docs/api/writespec.rst create mode 100644 docs/conf.py create mode 100644 docs/developer/contributing.rst create mode 100644 docs/developer/structure.rst create mode 100644 docs/examples.rst create mode 100644 docs/index.rst create mode 100644 docs/installation.rst create mode 100644 docs/license.rst create mode 100644 docs/make.bat create mode 100644 docs/quickstart.rst create mode 100644 docs/references.rst create mode 100644 docs/requirements.txt create mode 100644 docs/user_guide/estimation_methods.rst create mode 100644 docs/user_guide/index.rst create mode 100644 docs/user_guide/instrument_types.rst create mode 100644 docs/user_guide/troubleshooting.rst create mode 100644 docs/user_guide/understanding_spectra.rst diff --git a/.gitignore b/.gitignore index f7dd7f8..42a173f 100644 --- a/.gitignore +++ b/.gitignore @@ -62,3 +62,4 @@ coverage.xml Thumbs.db +_build/ diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/api/dirspec.rst b/docs/api/dirspec.rst new file mode 100644 index 0000000..9b3d312 --- /dev/null +++ b/docs/api/dirspec.rst @@ -0,0 +1,4 @@ +dirspec +======= + +.. autofunction:: pydiwasp.dirspec.dirspec diff --git a/docs/api/index.rst b/docs/api/index.rst new file mode 100644 index 0000000..3799f02 --- /dev/null +++ b/docs/api/index.rst @@ -0,0 +1,147 @@ +API Reference +============= + +This section provides detailed documentation for all public functions in pyDIWASP. + +Main Functions +-------------- + +.. toctree:: + :maxdepth: 2 + + dirspec + infospec + plotspec + interpspec + writespec + +Utility Functions +----------------- + +.. autofunction:: pydiwasp.infospec.compangle + +Data Structures +--------------- + +pyDIWASP uses Python dictionaries to represent data structures. This section describes +the expected structure and fields for each type. + +Instrument Data Structure (ID) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Dictionary containing information about the instrument array and measurements. + +**Required Fields:** + +* ``layout`` : ndarray, shape (3, N) + Instrument positions as [x; y; z] in meters. Each column represents one instrument. + +* ``datatypes`` : list of str, length N + Type of measurement for each instrument. Valid values: 'pres', 'elev', 'velx', 'vely', 'vels', 'accs'. + +* ``depth`` : float + Water depth in meters at the array location. + +* ``fs`` : float + Sampling frequency in Hz. + +* ``data`` : ndarray, shape (nsamples, N) + Time series measurements. Each column corresponds to one instrument. + +**Example:** + +.. code-block:: python + + ID = { + 'layout': np.array([[0, 10, 20], [0, 0, 0], [0, 0, 0]]), + 'datatypes': ['pres', 'pres', 'pres'], + 'depth': 10.0, + 'fs': 2.0, + 'data': measurements # shape: (nsamples, 3) + } + +Spectral Matrix Structure (SM) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Dictionary defining the frequency-direction grid for the spectrum. + +**Required Fields (Input):** + +* ``freqs`` : ndarray + Frequency values for the output grid (Hz or rad/s). + +* ``dirs`` : ndarray + Direction values for the output grid (radians or degrees). + +**Optional Fields (Input):** + +* ``funit`` : str, default 'Hz' + Frequency units: 'Hz' or 'rad/s'. + +* ``dunit`` : str, default 'rad' + Direction units: 'rad' or 'deg'. + +* ``xaxisdir`` : float, default 90 + Direction of positive x-axis in compass degrees (90 = East, 0 = North). + +**Output Fields (Added by dirspec):** + +* ``S`` : ndarray, shape (nfreqs, ndirs) + Computed spectral density values. + +**Example:** + +.. code-block:: python + + # Input + SM = { + 'freqs': np.linspace(0.05, 0.5, 50), + 'dirs': np.linspace(-180, 180, 36), + 'funit': 'Hz', + 'dunit': 'deg' + } + + # After dirspec() + SMout = { + 'freqs': ..., + 'dirs': ..., + 'S': ... # shape (50, 36) + } + +Estimation Parameters Structure (EP) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Dictionary specifying parameters for the spectral estimation algorithm. + +**Optional Fields (all have defaults):** + +* ``method`` : str, default 'IMLM' + Estimation method: 'IMLM' or 'EMEP'. + +* ``nfft`` : int, default: auto-calculated + FFT length for spectral estimation. Should be power of 2. + +* ``dres`` : int, default 180 + Directional resolution (number of direction bins for internal calculation). + +* ``iter`` : int, default 100 + Number of iterations for iterative methods. + +* ``smooth`` : str, default 'ON' + Spectral smoothing: 'ON' or 'OFF'. + +**Example:** + +.. code-block:: python + + # Use defaults + EP = {} + + # Or specify parameters + EP = { + 'method': 'IMLM', + 'iter': 100, + 'nfft': 512, + 'dres': 180, + 'smooth': 'ON' + } diff --git a/docs/api/infospec.rst b/docs/api/infospec.rst new file mode 100644 index 0000000..0698d4c --- /dev/null +++ b/docs/api/infospec.rst @@ -0,0 +1,4 @@ +infospec +======== + +.. autofunction:: pydiwasp.infospec.infospec diff --git a/docs/api/interpspec.rst b/docs/api/interpspec.rst new file mode 100644 index 0000000..f4dbff8 --- /dev/null +++ b/docs/api/interpspec.rst @@ -0,0 +1,4 @@ +interpspec +========== + +.. autofunction:: pydiwasp.interpspec.interpspec diff --git a/docs/api/plotspec.rst b/docs/api/plotspec.rst new file mode 100644 index 0000000..20def06 --- /dev/null +++ b/docs/api/plotspec.rst @@ -0,0 +1,4 @@ +plotspec +======== + +.. autofunction:: pydiwasp.plotspec.plotspec diff --git a/docs/api/writespec.rst b/docs/api/writespec.rst new file mode 100644 index 0000000..daab27e --- /dev/null +++ b/docs/api/writespec.rst @@ -0,0 +1,4 @@ +writespec +========= + +.. autofunction:: pydiwasp.writespec.writespec diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..2225108 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,69 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +import os +import sys +sys.path.insert(0, os.path.abspath('..')) + +# -- Project information ----------------------------------------------------- +project = 'pyDIWASP' +copyright = '2024, SBFRF' +author = 'SBFRF' +release = '0.1.0' + +# -- General configuration --------------------------------------------------- +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.napoleon', + 'sphinx.ext.viewcode', + 'sphinx.ext.intersphinx', + 'sphinx.ext.mathjax', + 'sphinx_rtd_theme', +] + +templates_path = ['_templates'] +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] + +# -- Options for HTML output ------------------------------------------------- +html_theme = 'sphinx_rtd_theme' +html_static_path = ['_static'] +html_theme_options = { + 'navigation_depth': 4, + 'titles_only': False, +} + +# -- Extension configuration ------------------------------------------------- +# Napoleon settings +napoleon_google_docstring = True +napoleon_numpy_docstring = True +napoleon_include_init_with_doc = False +napoleon_include_private_with_doc = False +napoleon_include_special_with_doc = True +napoleon_use_admonition_for_examples = False +napoleon_use_admonition_for_notes = False +napoleon_use_admonition_for_references = False +napoleon_use_ivar = False +napoleon_use_param = True +napoleon_use_rtype = True +napoleon_preprocess_types = False +napoleon_type_aliases = None +napoleon_attr_annotations = True + +# Intersphinx mapping +intersphinx_mapping = { + 'python': ('https://docs.python.org/3', None), + 'numpy': ('https://numpy.org/doc/stable/', None), + 'scipy': ('https://docs.scipy.org/doc/scipy/', None), + 'matplotlib': ('https://matplotlib.org/stable/', None), +} + +# Autodoc settings +autodoc_default_options = { + 'members': True, + 'member-order': 'bysource', + 'special-members': '__init__', + 'undoc-members': True, + 'exclude-members': '__weakref__' +} diff --git a/docs/developer/contributing.rst b/docs/developer/contributing.rst new file mode 100644 index 0000000..33875d9 --- /dev/null +++ b/docs/developer/contributing.rst @@ -0,0 +1,303 @@ +Contributing to pyDIWASP +========================= + +We welcome contributions to pyDIWASP! This guide explains how to contribute. + +Types of Contributions +----------------------- + +We appreciate all types of contributions: + +* **Bug reports**: Help us identify and fix issues +* **Feature requests**: Suggest new functionality +* **Documentation**: Improve or expand documentation +* **Code contributions**: Bug fixes, new features, optimizations +* **Examples**: Add new example notebooks or scripts + +Getting Started +--------------- + +1. **Fork the repository** on GitHub + +2. **Clone your fork**:: + + git clone https://github.com/YOUR-USERNAME/pyDIWASP.git + cd pyDIWASP + +3. **Create a development environment**:: + + python -m venv venv + source venv/bin/activate # On Windows: venv\\Scripts\\activate + pip install -e . + pip install pytest pytest-cov flake8 + +4. **Create a branch** for your changes:: + + git checkout -b feature/my-new-feature + # or + git checkout -b fix/issue-123 + +Development Workflow +-------------------- + +1. **Make your changes** in your feature branch + +2. **Write tests** for new functionality:: + + # Add tests to tests/test_*.py + def test_my_feature(): + # Test code here + assert result == expected + +3. **Run tests** to ensure everything works:: + + pytest tests/ + +4. **Check code style**:: + + flake8 pydiwasp tests + +5. **Commit your changes**:: + + git add . + git commit -m "Add: Brief description of changes" + +6. **Push to your fork**:: + + git push origin feature/my-new-feature + +7. **Create a Pull Request** on GitHub + +Code Guidelines +--------------- + +Style +~~~~~ + +* Follow PEP 8 style guidelines +* Use meaningful variable names +* Add docstrings to all public functions +* Keep functions focused and modular + +Docstrings +~~~~~~~~~~ + +Use NumPy-style docstrings: + +.. code-block:: python + + def my_function(param1, param2): + """ + Brief description of function. + + Longer description if needed, explaining what the function does + and any important details. + + Parameters + ---------- + param1 : type + Description of param1. + param2 : type + Description of param2. + + Returns + ------- + result : type + Description of return value. + + Examples + -------- + >>> result = my_function(1, 2) + >>> print(result) + 3 + """ + return param1 + param2 + +Testing +~~~~~~~ + +* Write tests for all new functions +* Aim for high test coverage (>80%) +* Test edge cases and error conditions +* Use descriptive test names + +.. code-block:: python + + def test_dirspec_with_pressure_sensors(): + """Test dirspec with a simple pressure sensor array.""" + # Setup + ID = {...} + SM = {...} + EP = {...} + + # Execute + SMout, EPout = dirspec(ID, SM, EP) + + # Assert + assert SMout is not None + assert 'S' in SMout + assert SMout['S'].shape == (len(SM['freqs']), len(SM['dirs'])) + +Documentation +~~~~~~~~~~~~~ + +* Update documentation for any API changes +* Add examples for new features +* Keep README.md up to date + +Reporting Bugs +-------------- + +When reporting bugs, please include: + +1. **Description**: Clear description of the problem + +2. **Minimal example**: Code that reproduces the issue:: + + import numpy as np + from pydiwasp import dirspec + + # Minimal code that shows the bug + ID = {...} + SMout, EPout = dirspec(ID, SM, EP) + # Error occurs here + +3. **Error message**: Full traceback if applicable + +4. **Environment**: + * Python version + * pyDIWASP version + * Operating system + * Package versions: ``pip list`` + +5. **Expected behavior**: What you expected to happen + +6. **Actual behavior**: What actually happened + +Requesting Features +------------------- + +For feature requests, please: + +1. **Check existing issues** to avoid duplicates + +2. **Describe the feature** clearly + +3. **Explain the use case**: Why is this feature needed? + +4. **Provide examples**: Show how you'd use the feature + +5. **Consider implementation**: If possible, suggest how it might work + +Adding New Estimation Methods +------------------------------ + +To add a new spectrum estimation method: + +1. **Create the method file**:: + + # pydiwasp/private/MYMETHOD.py + + def MYMETHOD(xps, trm, kx, Ss, pidirs, miter, displ): + """ + My new estimation method. + + Parameters + ---------- + xps : ndarray + Cross-spectral density matrix + trm : ndarray + Transfer function matrix + kx : ndarray + Wavenumber separation matrix + Ss : ndarray + Auto-spectra + pidirs : ndarray + Direction array + miter : int + Number of iterations + displ : int + Display level + + Returns + ------- + S : ndarray + Estimated directional spectrum + """ + # Implementation here + return S + +2. **Import in dirspec.py**:: + + from .private.MYMETHOD import MYMETHOD + +3. **Add to method call**: + + The method is called via ``eval(EP['method'])`` in dirspec.py + +4. **Write tests**:: + + # tests/test_methods.py + + def test_MYMETHOD(): + # Test your method + pass + +5. **Update documentation**: + + * Add to README.md + * Add to docs/user_guide/estimation_methods.rst + * Add docstring examples + +6. **Add citation**: Include reference paper/algorithm description + +Release Process +--------------- + +For maintainers: + +1. Update version in ``pyproject.toml`` and ``setup.py`` + +2. Update CHANGES_SUMMARY.md + +3. Create git tag:: + + git tag -a v0.2.0 -m "Release version 0.2.0" + git push origin v0.2.0 + +4. Build distribution:: + + python -m build + +5. Upload to PyPI:: + + twine upload dist/* + +Code of Conduct +--------------- + +* Be respectful and inclusive +* Welcome newcomers +* Focus on constructive feedback +* Assume good intentions + +License +------- + +By contributing, you agree that your contributions will be licensed under the +GNU General Public License v3.0 with the same addendum as the original DIWASP: + +* Use for educational purposes or scientific research without financial gain +* Written consent required for commercial use +* No warranties provided + +Questions? +---------- + +If you have questions about contributing: + +* Open a GitHub issue with the "question" label +* Check existing documentation +* Contact the maintainers + +Thank you for contributing to pyDIWASP! diff --git a/docs/developer/structure.rst b/docs/developer/structure.rst new file mode 100644 index 0000000..0e7f02d --- /dev/null +++ b/docs/developer/structure.rst @@ -0,0 +1,373 @@ +Code Structure +============== + +This document describes the organization of the pyDIWASP codebase. + +Directory Structure +------------------- + +:: + + pyDIWASP/ + ├── pydiwasp/ # Main package directory + │ ├── __init__.py # Package initialization, public API + │ ├── dirspec.py # Main spectral estimation function + │ ├── infospec.py # Wave statistics calculation + │ ├── plotspec.py # Spectrum visualization + │ ├── interpspec.py # Spectrum interpolation + │ ├── writespec.py # File output + │ └── private/ # Internal functions + │ ├── __init__.py + │ ├── IMLM.py # IMLM estimation method + │ ├── EMEP.py # EMEP estimation method + │ ├── wavenumber.py # Dispersion relation + │ ├── pres.py # Pressure transfer function + │ ├── elev.py # Elevation transfer function + │ ├── velx.py # X-velocity transfer function + │ ├── vely.py # Y-velocity transfer function + │ ├── vels.py # Velocity magnitude transfer + │ ├── accs.py # Acceleration transfer + │ ├── diwasp_csd.py # Cross-spectral density + │ ├── smoothspec.py # Spectral smoothing + │ ├── hsig.py # Significant wave height + │ ├── check_data.py # Input validation + │ └── spectobasis.py # Coordinate transformations + ├── tests/ # Test suite + │ ├── test_api.py # API tests + │ ├── test_core.py # Core functionality tests + │ └── test_integration.py # Integration tests + ├── examples/ # Example notebooks and scripts + │ └── pyDIWASP_example.ipynb + ├── docs/ # Documentation + │ ├── conf.py # Sphinx configuration + │ ├── index.rst # Documentation index + │ └── ... + ├── setup.py # Package setup script + ├── pyproject.toml # Modern Python package metadata + ├── README.md # Project README + └── requirements.txt # Dependencies + +Module Overview +--------------- + +Public API (pydiwasp/) +~~~~~~~~~~~~~~~~~~~~~~ + +These modules form the user-facing API: + +**dirspec.py** + Main function for directional spectrum estimation. Orchestrates the entire + analysis workflow: + + 1. Validates input data structures + 2. Computes cross-spectral densities + 3. Calculates wavenumbers and transfer functions + 4. Calls the selected estimation method + 5. Interpolates to user-specified grid + 6. Optionally smooths and plots results + +**infospec.py** + Computes wave statistics from a directional spectrum: + + * Significant wave height (Hsig) + * Peak period (Tp) + * Direction at peak (DTp) + * Dominant direction (Dp) + + Also includes ``compangle()`` for converting between angle conventions. + +**plotspec.py** + Creates visualizations of directional spectra: + + * 3D surface plots + * Polar contour plots + * Mathematical or nautical angle conventions + +**interpspec.py** + Interpolates spectra between different frequency-direction grids. + Uses 2D interpolation in transformed coordinates to handle periodicity. + +**writespec.py** + Exports spectra to DIWASP text format for archiving or compatibility. + +Private Modules (pydiwasp/private/) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +These are internal implementation details, not part of the public API: + +Estimation Methods +^^^^^^^^^^^^^^^^^^ + +**IMLM.py** + Implements the Iterated Maximum Likelihood Method. + + * Iteratively refines directional estimate + * Maximizes likelihood of observed cross-spectra + * Good for complex, multi-modal distributions + +**EMEP.py** + Implements the Extended Maximum Entropy Principle. + + * Based on maximum entropy principle + * Good for narrow directional spreads + * Non-iterative formulation + +Transfer Functions +^^^^^^^^^^^^^^^^^^ + +Each instrument type has a transfer function relating measurements to surface elevation: + +**pres.py** + Pressure sensor transfer function. Accounts for pressure attenuation with depth + according to linear wave theory. + +**elev.py** + Surface elevation transfer function (identity, no transformation needed). + +**velx.py, vely.py** + Horizontal velocity transfer functions. Convert velocity measurements to + equivalent surface elevation using linear wave theory. + +**vels.py** + Horizontal velocity magnitude transfer function. + +**accs.py** + Acceleration transfer function. + +Core Utilities +^^^^^^^^^^^^^^ + +**wavenumber.py** + Solves the dispersion relation to compute wave numbers from frequency and depth. + Uses iterative method based on hyperbolic tangent. + +**diwasp_csd.py** + Computes cross-spectral density between two time series using Welch's method. + +**smoothspec.py** + Applies smoothing to directional spectra to reduce noise. + +**hsig.py** + Computes significant wave height from integrated spectrum energy. + +**check_data.py** + Validates input data structures (ID, SM, EP) and fills in default values. + +**spectobasis.py** + Handles coordinate system transformations for spectral interpolation. + +Data Flow +--------- + +The typical data flow through pyDIWASP: + +1. **User Input**: + + * ID (instrument data) + * SM (spectral matrix grid) + * EP (estimation parameters) + +2. **dirspec() Main Function**: + + a. **Validation** (check_data.py): + - Verify data structures + - Fill in defaults + + b. **Preprocessing**: + - Detrend data + - Determine FFT parameters + + c. **Cross-Spectral Density** (diwasp_csd.py): + - Compute CSD for all instrument pairs + - Uses Welch's method with windowing + + d. **Wave Numbers** (wavenumber.py): + - Compute k from frequency and depth + - Uses dispersion relation + + e. **Transfer Functions** (pres.py, etc.): + - Compute transfer matrix for each instrument + - Accounts for instrument type and position + + f. **Estimation** (IMLM.py or EMEP.py): + - Apply selected method + - Estimate S(f, θ) + + g. **Interpolation** (interpspec.py): + - Map to user-specified grid + - Preserve energy + + h. **Post-processing**: + - Optional smoothing (smoothspec.py) + - Optional plotting (plotspec.py) + - Optional file output (writespec.py) + +3. **Output**: + + * SMout (spectrum with estimated S field) + * EPout (parameters actually used) + +Key Algorithms +-------------- + +Dispersion Relation (wavenumber.py) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Solves: + +.. math:: + + \\omega^2 = gk \\tanh(kh) + +where ω is angular frequency, g is gravity, k is wavenumber, h is depth. + +Uses Newton-Raphson iteration starting from an initial guess based on the +deep-water approximation. + +Cross-Spectral Density (diwasp_csd.py) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Computes: + +.. math:: + + C_{ij}(f) = \\langle X_i(f) X_j^*(f) \\rangle + +where X_i(f) is the Fourier transform of signal i, * denotes complex conjugate, +and ⟨⟩ denotes averaging over segments. + +IMLM Algorithm (IMLM.py) +~~~~~~~~~~~~~~~~~~~~~~~~~ + +Iteratively updates the directional distribution: + +.. math:: + + E^{(n+1)}(\\theta) = E^{(n)}(\\theta) \\cdot f(E^{(n)}) + +where f depends on the comparison between observed and predicted cross-spectra. + +Converges to maximum likelihood estimate. + +Extension Points +---------------- + +Adding New Features +~~~~~~~~~~~~~~~~~~~ + +**New Instrument Type:** + +1. Create transfer function in ``private/newtype.py`` +2. Follow same signature as existing transfer functions +3. Add to valid datatypes in ``check_data.py`` +4. Update documentation + +**New Estimation Method:** + +1. Create method in ``private/NEWMETHOD.py`` +2. Follow same signature as IMLM.py +3. Method will be called via eval() in dirspec.py +4. Add tests and documentation + +**New Output Format:** + +1. Add function to main package or private/ +2. Follow similar pattern to writespec.py +3. Document file format + +Testing Strategy +---------------- + +The test suite uses pytest and is organized by scope: + +**test_api.py** + Tests the public API functions with various inputs. Ensures backward compatibility. + +**test_core.py** + Tests core algorithms and private functions. Verifies mathematical correctness. + +**test_integration.py** + End-to-end tests with realistic scenarios. Ensures components work together. + +Run tests with:: + + pytest tests/ + +For coverage report:: + + pytest --cov=pydiwasp tests/ + +Dependencies +------------ + +Core Dependencies +~~~~~~~~~~~~~~~~~ + +* **NumPy**: Array operations, FFT +* **SciPy**: Interpolation, signal processing +* **Matplotlib**: Plotting + +Development Dependencies +~~~~~~~~~~~~~~~~~~~~~~~~ + +* **pytest**: Testing framework +* **pytest-cov**: Coverage reporting +* **flake8**: Code style checking +* **sphinx**: Documentation generation + +Performance Considerations +-------------------------- + +Computational Bottlenecks +~~~~~~~~~~~~~~~~~~~~~~~~~ + +1. **Cross-spectral density computation**: O(N² × nfft × log(nfft)) + + * Scales quadratically with number of instruments + * Consider parallel computation for large arrays + +2. **Estimation methods**: O(nfreqs × ndirs × iters) + + * IMLM: Linear in iterations + * Most time in matrix operations + +3. **Interpolation**: O(nfreqs_in × ndirs_in × log(npoints)) + + * Uses scipy's griddata (QHull algorithm) + +Memory Usage +~~~~~~~~~~~~ + +Main memory consumers: + +* Cross-spectral matrix: N × N × nfreqs (complex) +* Transfer matrices: N × nfreqs × ndirs +* Wavenumber matrix: N × N × nfreqs × ndirs + +For typical arrays (3-5 instruments), memory is not a concern. + +Optimization Opportunities +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Potential improvements (contributions welcome): + +* Parallel cross-spectral density computation +* Cython/Numba acceleration for IMLM inner loops +* Sparse matrix operations where applicable +* GPU acceleration for large arrays + +Coding Conventions +------------------ + +* Variables follow MATLAB naming from original DIWASP for easier comparison +* Function names are lowercase with underscores +* Private functions start with underscore (though in private/ dir) +* Dictionary keys for data structures match MATLAB field names +* Comments explain *why*, not *what* (code should be self-documenting) + +See Also +-------- + +* :doc:`contributing` for contribution guidelines +* Original DIWASP documentation for algorithm details +* Hashimoto (1997) for mathematical foundations diff --git a/docs/examples.rst b/docs/examples.rst new file mode 100644 index 0000000..c7b707a --- /dev/null +++ b/docs/examples.rst @@ -0,0 +1,326 @@ +Examples +======== + +This section provides example usage of pyDIWASP for various scenarios. + +Example Notebook +---------------- + +The main example is provided as a Jupyter notebook in the repository: + +* `pyDIWASP_example.ipynb `_ + +This notebook demonstrates: + +* Setting up instrument data structures +* Computing directional spectra +* Analyzing wave statistics +* Creating visualizations +* Comparing estimation methods + +Basic Usage Examples +-------------------- + +Example 1: Three Pressure Sensors +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec, infospec, plotspec + import matplotlib.pyplot as plt + + # Generate synthetic wave data (for demonstration) + t = np.arange(0, 1024) / 2.0 # 512 seconds at 2 Hz + + # Simulate waves from 45° with period 8s + freq = 1/8 # 0.125 Hz + wave = np.sin(2*np.pi*freq*t) + + # Three sensors in triangular array + sensor1 = wave + 0.1*np.random.randn(len(t)) + sensor2 = wave + 0.1*np.random.randn(len(t)) + sensor3 = wave + 0.1*np.random.randn(len(t)) + + data = np.column_stack([sensor1, sensor2, sensor3]) + + # Define instrument data + ID = { + 'layout': np.array([[0, 10, 5], [0, 0, 8.66], [0, 0, 0]]), + 'datatypes': ['pres', 'pres', 'pres'], + 'depth': 10.0, + 'fs': 2.0, + 'data': data + } + + # Define spectral matrix + SM = { + 'freqs': np.linspace(0.05, 0.5, 50), + 'dirs': np.linspace(-180, 180, 36) + } + + # Estimation parameters + EP = {'method': 'IMLM', 'iter': 100} + + # Compute spectrum + SMout, EPout = dirspec(ID, SM, EP, ['PLOTTYPE', 0]) + + # Get statistics + Hsig, Tp, DTp, Dp = infospec(SMout) + + # Plot + plotspec(SMout, 4) # polar plot with compass bearings + plt.show() + +Example 2: Mixed Instrument Types +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec + + # Load your data + pressure_data = np.loadtxt('pressure.csv', delimiter=',') + velocity_x = np.loadtxt('velocity_x.csv', delimiter=',') + velocity_y = np.loadtxt('velocity_y.csv', delimiter=',') + elevation = np.loadtxt('elevation.csv', delimiter=',') + + # Combine into single array + data = np.column_stack([ + pressure_data[:, 0], + pressure_data[:, 1], + velocity_x, + velocity_y, + elevation + ]) + + # Mixed array configuration + ID = { + 'layout': np.array([ + [0, 10, 15, 15, 5], # x positions + [0, 0, 5, 5, 10], # y positions + [0, 0, -2, -2, 0] # z positions + ]), + 'datatypes': ['pres', 'pres', 'velx', 'vely', 'elev'], + 'depth': 12.0, + 'fs': 4.0, + 'data': data + } + + SM = { + 'freqs': np.linspace(0.05, 0.5, 50), + 'dirs': np.linspace(-np.pi, np.pi, 36) + } + + EP = {'method': 'IMLM'} + + SMout, EPout = dirspec(ID, SM, EP) + +Example 3: Time Series Analysis +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec, infospec + import matplotlib.pyplot as plt + + # Load long time series + full_data = np.loadtxt('wave_timeseries.csv', delimiter=',') + + # Analysis parameters + fs = 2.0 # sampling frequency + window_length = 20 * 60 * fs # 20 minutes + overlap = 0.5 + step = int(window_length * (1 - overlap)) + + # Fixed structures + ID = { + 'layout': np.array([[0, 10, 20], [0, 0, 0], [0, 0, 0]]), + 'datatypes': ['pres', 'pres', 'pres'], + 'depth': 10.0, + 'fs': fs, + 'data': None # will be filled for each window + } + + SM = { + 'freqs': np.linspace(0.05, 0.5, 50), + 'dirs': np.linspace(-180, 180, 36) + } + + EP = {'method': 'IMLM', 'iter': 100} + + # Process each window + times = [] + Hsigs = [] + Tps = [] + Dps = [] + + for i in range(0, len(full_data) - window_length, step): + window = full_data[i:i+int(window_length), :] + ID['data'] = window + + # Compute spectrum + SMout, EPout = dirspec(ID, SM, EP, ['MESSAGE', 0, 'PLOTTYPE', 0]) + + # Get statistics + Hsig, Tp, DTp, Dp = infospec(SMout) + + times.append(i / fs / 3600) # time in hours + Hsigs.append(Hsig) + Tps.append(Tp) + Dps.append(Dp) + + # Plot time series of parameters + fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True) + + axes[0].plot(times, Hsigs) + axes[0].set_ylabel('Hsig (m)') + axes[0].grid(True) + + axes[1].plot(times, Tps) + axes[1].set_ylabel('Tp (s)') + axes[1].grid(True) + + axes[2].plot(times, Dps) + axes[2].set_ylabel('Dp (deg)') + axes[2].set_xlabel('Time (hours)') + axes[2].grid(True) + + plt.tight_layout() + plt.show() + +Example 4: Comparing Estimation Methods +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec, infospec, plotspec + import matplotlib.pyplot as plt + + # Your data setup + ID = {...} + SM = {...} + + # Compare IMLM and EMEP + methods = ['IMLM', 'EMEP'] + results = {} + + for method in methods: + EP = {'method': method, 'iter': 100} + SMout, EPout = dirspec(ID, SM, EP, ['PLOTTYPE', 0]) + + Hsig, Tp, DTp, Dp = infospec(SMout) + + results[method] = { + 'spectrum': SMout, + 'Hsig': Hsig, + 'Tp': Tp, + 'DTp': DTp, + 'Dp': Dp + } + + # Plot comparison + fig = plt.figure(figsize=(14, 5)) + + for i, method in enumerate(methods): + plt.subplot(1, 2, i+1) + plotspec(results[method]['spectrum'], 1) + plt.title(f"{method}: Hsig={results[method]['Hsig']:.2f}m, " + f"Tp={results[method]['Tp']:.1f}s") + + plt.tight_layout() + plt.show() + + # Print statistics comparison + print("Method Comparison:") + print("-" * 60) + for method in methods: + r = results[method] + print(f"{method}:") + print(f" Hsig = {r['Hsig']:.3f} m") + print(f" Tp = {r['Tp']:.2f} s") + print(f" DTp = {r['DTp']:.1f}°") + print(f" Dp = {r['Dp']:.1f}°") + print() + +Example 5: Exporting Results +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec, infospec, writespec + + # Compute spectrum + SMout, EPout = dirspec(ID, SM, EP) + Hsig, Tp, DTp, Dp = infospec(SMout) + + # Write spectrum to DIWASP format + writespec(SMout, 'output_spectrum.txt') + + # Export statistics to CSV + stats = { + 'Hsig': Hsig, + 'Tp': Tp, + 'DTp': DTp, + 'Dp': Dp + } + + with open('wave_statistics.csv', 'w') as f: + f.write('Parameter,Value\\n') + for key, value in stats.items(): + f.write(f'{key},{value}\\n') + + # Export spectrum to numpy format + np.savez('spectrum.npz', + freqs=SMout['freqs'], + dirs=SMout['dirs'], + S=SMout['S'], + Hsig=Hsig, + Tp=Tp, + DTp=DTp, + Dp=Dp) + + # Load back later + loaded = np.load('spectrum.npz') + freqs = loaded['freqs'] + dirs = loaded['dirs'] + S = loaded['S'] + +Example 6: Custom Frequency Range +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec + + # Focus on swell (long period waves) + SM_swell = { + 'freqs': np.linspace(0.04, 0.15, 30), # 6.7 to 25 seconds + 'dirs': np.linspace(-180, 180, 36) + } + + EP = {'method': 'EMEP'} # EMEP good for swell + SMout_swell, _ = dirspec(ID, SM_swell, EP) + + # Focus on wind waves (short period) + SM_wind = { + 'freqs': np.linspace(0.15, 0.5, 40), # 2 to 6.7 seconds + 'dirs': np.linspace(-180, 180, 36) + } + + EP = {'method': 'IMLM'} # IMLM good for complex seas + SMout_wind, _ = dirspec(ID, SM_wind, EP) + +More Examples +------------- + +For more examples, see: + +* The `examples directory `_ in the repository +* The test files in `tests/ `_ directory +* The original MATLAB DIWASP manual for conceptual guidance diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..02fbf33 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,103 @@ +.. pyDIWASP documentation master file + +pyDIWASP: Directional Wave Spectrum Analysis in Python +======================================================= + +.. image:: https://img.shields.io/pypi/v/pyDIWASP.svg + :target: https://pypi.org/project/pyDIWASP/ + :alt: PyPI version + +.. image:: https://img.shields.io/pypi/pyversions/pyDIWASP.svg + :target: https://pypi.org/project/pyDIWASP/ + :alt: Python versions + +**pyDIWASP** is a Python implementation of the DIWASP (DIrectional WAve SPectrum) toolbox for estimating directional wave spectra from wave measurement data. + +Key Features +------------ + +* **Directional Wave Analysis**: Estimate directional wave spectra from instrument array data +* **Multiple Estimation Methods**: Supports IMLM (Iterated Maximum Likelihood Method) and EMEP (Extended Maximum Entropy Principle) +* **Flexible Input**: Works with various instrument types (pressure, velocity, elevation sensors) +* **Visualization**: Built-in plotting functions for 3D and polar spectral plots +* **Wave Statistics**: Calculate significant wave height, peak period, and dominant direction + +Quick Start +----------- + +Installation:: + + pip install pyDIWASP + +Basic usage: + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec, infospec, plotspec + + # Define instrument data + ID = { + 'layout': np.array([[0, 10, 20], [0, 0, 0], [0, 0, 0]]), + 'datatypes': ['pres', 'pres', 'pres'], + 'depth': 10.0, + 'fs': 2.0, + 'data': wave_data # your measurement data + } + + # Define spectral matrix + SM = { + 'freqs': np.linspace(0.05, 0.5, 50), + 'dirs': np.linspace(-np.pi, np.pi, 36) + } + + # Define estimation parameters + EP = {'method': 'IMLM', 'iter': 100} + + # Compute directional spectrum + SMout, EPout = dirspec(ID, SM, EP) + + # Display wave statistics + Hsig, Tp, DTp, Dp = infospec(SMout) + + # Plot the spectrum + plotspec(SMout, 1) + +Documentation Contents +---------------------- + +.. toctree:: + :maxdepth: 2 + :caption: User Guide + + installation + quickstart + user_guide/index + examples + +.. toctree:: + :maxdepth: 2 + :caption: API Reference + + api/index + +.. toctree:: + :maxdepth: 1 + :caption: Developer Guide + + developer/contributing + developer/structure + +.. toctree:: + :maxdepth: 1 + :caption: Additional Information + + references + license + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/installation.rst b/docs/installation.rst new file mode 100644 index 0000000..5a8a5f2 --- /dev/null +++ b/docs/installation.rst @@ -0,0 +1,82 @@ +Installation +============ + +pyDIWASP can be installed from source or via pip (once published to PyPI). + +Requirements +------------ + +pyDIWASP requires: + +* Python 3.8 or later +* NumPy +* SciPy +* Matplotlib + +From PyPI +--------- + +Once published, you can install pyDIWASP using pip:: + + pip install pyDIWASP + +From Source +----------- + +To install from source: + +1. Clone the repository:: + + git clone https://github.com/SBFRF/pyDIWASP.git + cd pyDIWASP + +2. Install in development mode:: + + pip install -e . + +This installs the package in "editable" mode, so changes to the source code are immediately reflected without reinstalling. + +For Contributors +---------------- + +If you plan to contribute to pyDIWASP, install the development dependencies:: + + pip install -e . + pip install pytest pytest-cov flake8 + +This includes testing and linting tools. + +Verifying Installation +---------------------- + +To verify that pyDIWASP is correctly installed, run: + +.. code-block:: python + + import pydiwasp + print(pydiwasp.__version__) + +You should see the version number printed without any errors. + +You can also run the test suite:: + + pytest tests/ + +Troubleshooting +--------------- + +**Import errors** + +If you get import errors, make sure all dependencies are installed:: + + pip install numpy scipy matplotlib + +**Module not found** + +If Python cannot find the pydiwasp module, ensure you're running Python from the correct environment where you installed the package. + +Check your Python path:: + + python -c "import sys; print(sys.path)" + +The directory containing pydiwasp should be in this list. diff --git a/docs/license.rst b/docs/license.rst new file mode 100644 index 0000000..7a8e5d0 --- /dev/null +++ b/docs/license.rst @@ -0,0 +1,142 @@ +License +======= + +pyDIWASP License +---------------- + +pyDIWASP is a Python implementation of the DIWASP toolbox and inherits the same +license terms as the original MATLAB version. + +GNU General Public License v3.0 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +pyDIWASP is free software; you can redistribute it and/or modify it under the +terms of the GNU General Public License as published by the Free Software Foundation, +either version 3 of the License, or (at your option) any later version. + +Important License Addendum +--------------------------- + +The DIWASP license includes the following important addendum concerning its usage: + + **This software and any derivatives of it shall only be used for educational + purposes or scientific research without the intention of any financial gain.** + + **Use of this software or derivatives for any purpose that results in financial + gain for a person or organization without written consent from the author is a + breach of the license agreement.** + +This means: + +Permitted Uses +~~~~~~~~~~~~~~ + +✓ **Academic Research**: Use in university research projects + +✓ **Education**: Teaching and learning about wave analysis + +✓ **Non-commercial Studies**: Scientific investigations without financial benefit + +✓ **Open Source Projects**: Contributions to free, open-source software + +Requires Permission +~~~~~~~~~~~~~~~~~~~ + +✗ **Commercial Applications**: Products or services sold for profit + +✗ **Consulting Work**: Paid consulting using this software + +✗ **Industrial Use**: Use by for-profit companies without agreement + +✗ **Proprietary Software**: Incorporation into closed-source commercial software + +**If you wish to use pyDIWASP for commercial purposes, you must obtain written +consent from the original DIWASP author.** + +Warranty Disclaimer +------------------- + +This software is distributed in the hope that it will be useful, but **WITHOUT ANY +WARRANTY**; without even the implied warranty of **MERCHANTABILITY** or **FITNESS +FOR A PARTICULAR PURPOSE**. + +In addition, the author is **not liable in any way** for consequences arising from +the application of software output for any design or decision-making process. + +Full License Text +----------------- + +The full GNU General Public License v3.0 can be found at: +https://www.gnu.org/licenses/gpl-3.0.en.html + +Or see the LICENSE.md file in the repository root. + +Copyright Information +--------------------- + +Original DIWASP +~~~~~~~~~~~~~~~ + +Copyright (C) 2002 David Johnson + +Coastal Oceanography Group, Centre for Water Research, University of Western Australia, Perth + +Python Implementation +~~~~~~~~~~~~~~~~~~~~~ + +Translated by Chuan Li and Spicer Bak + +Field Research Facility, US Army Corps of Engineers + +Current Maintenance +~~~~~~~~~~~~~~~~~~~ + +Maintained by SBFRF (Systems, Biology, and Behavior Research Foundation) + +GitHub: https://github.com/SBFRF + +Contact +------- + +For license inquiries or permission requests: + +* **GitHub Issues**: https://github.com/SBFRF/pyDIWASP/issues +* **Email**: support@sbfrf.org + +Contributing +------------ + +By contributing to pyDIWASP, you agree that your contributions will be licensed +under the same GNU General Public License v3.0 with the non-commercial use addendum. + +See CONTRIBUTING.md for more details on contributing to the project. + +Third-Party Licenses +-------------------- + +pyDIWASP depends on several open-source libraries: + +NumPy +~~~~~ + +Copyright (c) 2005-2023, NumPy Developers. + +Licensed under the BSD 3-Clause License. + +SciPy +~~~~~ + +Copyright (c) 2001-2023 SciPy Developers. + +Licensed under the BSD 3-Clause License. + +Matplotlib +~~~~~~~~~~ + +Copyright (c) 2002-2023 John D. Hunter, Michael Droettboom, and the Matplotlib +development team. + +Licensed under the PSF-based license. + +These libraries are not affected by pyDIWASP's commercial use restrictions and +are subject to their own licenses. diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..922152e --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/quickstart.rst b/docs/quickstart.rst new file mode 100644 index 0000000..c4c2156 --- /dev/null +++ b/docs/quickstart.rst @@ -0,0 +1,157 @@ +Quick Start Guide +================= + +This guide will get you started with pyDIWASP in just a few minutes. + +What is Directional Wave Spectrum Analysis? +-------------------------------------------- + +A directional wave spectrum describes how wave energy is distributed across: + +* **Frequencies** (or wave periods): Different wave oscillation rates +* **Directions**: Where waves are coming from + +This is more informative than a simple frequency spectrum because it tells you not just how much energy exists at each frequency, but also the direction that energy is traveling. + +Basic Workflow +-------------- + +The typical pyDIWASP workflow consists of three steps: + +1. **Prepare data**: Define instrument array geometry and load measurements +2. **Compute spectrum**: Use ``dirspec()`` to estimate the directional spectrum +3. **Analyze results**: Use ``infospec()`` for statistics and ``plotspec()`` for visualization + +Minimal Example +--------------- + +Here's a complete minimal example: + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec, infospec, plotspec + + # Step 1: Define instrument data + ID = { + 'layout': np.array([ + [0, 10, 5], # x positions (m) + [0, 0, 8.66], # y positions (m) + [0, 0, 0] # z positions (m) + ]), + 'datatypes': ['pres', 'pres', 'pres'], # pressure sensors + 'depth': 10.0, # water depth (m) + 'fs': 2.0, # sampling frequency (Hz) + 'data': wave_data # time series data (nsamples x 3) + } + + # Step 2: Define output grid + SM = { + 'freqs': np.linspace(0.05, 0.5, 50), # frequencies (Hz) + 'dirs': np.linspace(-180, 180, 36) # directions (degrees) + } + + # Step 3: Set estimation parameters + EP = { + 'method': 'IMLM', # estimation method + 'iter': 100 # iterations + } + + # Step 4: Compute spectrum + SMout, EPout = dirspec(ID, SM, EP) + + # Step 5: Get wave statistics + Hsig, Tp, DTp, Dp = infospec(SMout) + + # Step 6: Visualize + plotspec(SMout, 1) # 3D surface plot + +Understanding the Data Structures +---------------------------------- + +Instrument Data (ID) +~~~~~~~~~~~~~~~~~~~~ + +The ``ID`` dictionary describes your instrument array: + +* **layout**: 3×N array with x, y, z positions of N instruments +* **datatypes**: List specifying what each instrument measures +* **depth**: Water depth at the array location +* **fs**: Data sampling rate +* **data**: The actual time series measurements (nsamples × N) + +Valid instrument types: + +* ``'pres'``: Pressure sensor (most common) +* ``'elev'``: Surface elevation (e.g., wave staff) +* ``'velx'``, ``'vely'``: Horizontal velocities +* ``'vels'``: Horizontal velocity magnitude + +Spectral Matrix (SM) +~~~~~~~~~~~~~~~~~~~~~ + +The ``SM`` dictionary defines the frequency-direction grid for the output: + +* **freqs**: Array of frequencies where you want spectrum values +* **dirs**: Array of directions where you want spectrum values + +After computation, ``SMout`` also contains: + +* **S**: The computed spectral density (nfreqs × ndirs) + +Estimation Parameters (EP) +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The ``EP`` dictionary controls the analysis: + +* **method**: ``'IMLM'`` (default, recommended) or ``'EMEP'`` +* **iter**: Number of iterations (default: 100) +* **nfft**: FFT length (auto-computed if omitted) +* **dres**: Directional resolution (default: 180) + +Loading Your Data +----------------- + +pyDIWASP works with NumPy arrays. Here's how to load data from common formats: + +From CSV: + +.. code-block:: python + + import numpy as np + + # Load from CSV file + wave_data = np.loadtxt('wave_measurements.csv', delimiter=',') + + ID['data'] = wave_data + +From MATLAB: + +.. code-block:: python + + from scipy.io import loadmat + + # Load from .mat file + mat_data = loadmat('wave_measurements.mat') + wave_data = mat_data['wave_data'] + + ID['data'] = wave_data + +From netCDF: + +.. code-block:: python + + import netCDF4 as nc + + # Load from netCDF file + dataset = nc.Dataset('wave_measurements.nc') + wave_data = dataset.variables['wave_data'][:] + + ID['data'] = wave_data + +Next Steps +---------- + +* See :doc:`user_guide/index` for detailed information on all features +* Check :doc:`examples` for more complex use cases +* Read :doc:`api/index` for complete function documentation diff --git a/docs/references.rst b/docs/references.rst new file mode 100644 index 0000000..6ec011b --- /dev/null +++ b/docs/references.rst @@ -0,0 +1,231 @@ +References +========== + +This page lists key references for pyDIWASP and directional wave spectrum analysis. + +Primary Reference +----------------- + +The estimation algorithms implemented in pyDIWASP are described in: + + Hashimoto, N. (1997). "Analysis of the directional wave spectrum from field data". + In: *Advances in Coastal Engineering Vol. 3*. Ed: Liu, P.L-F. + Pub: World Scientific, Singapore. + +This comprehensive chapter covers: + +* Theoretical foundations of directional spectrum estimation +* Description of multiple estimation methods (MLM, IMLM, EMLM, EMEP, BDM, DFTM) +* Array design considerations +* Validation and comparison studies + +Wave Theory +----------- + +Linear Wave Theory +~~~~~~~~~~~~~~~~~~ + +Dean, R. G., & Dalrymple, R. A. (1991). *Water wave mechanics for engineers and scientists*. +World Scientific. + +Comprehensive treatment of linear wave theory, dispersion relations, and wave kinematics. + +Mei, C. C., Stiassnie, M., & Yue, D. K. P. (2005). *Theory and applications of ocean surface waves*. +World Scientific. + +Advanced treatment of water wave theory including directional spectra. + +Directional Spectra +------------------- + +Theoretical Background +~~~~~~~~~~~~~~~~~~~~~~ + +Longuet-Higgins, M. S., Cartwright, D. E., & Smith, N. D. (1963). +"Observations of the directional spectrum of sea waves using the motions of a floating buoy". +In: *Ocean Wave Spectra*, Prentice-Hall, 111-136. + +Early foundational work on estimating directional spectra from buoy measurements. + +Capon, J. (1969). "High-resolution frequency-wavenumber spectrum analysis". +*Proceedings of the IEEE*, 57(8), 1408-1418. + +Introduction of Maximum Likelihood Method for spectral estimation. + +Maximum Entropy Methods +~~~~~~~~~~~~~~~~~~~~~~~~ + +Lygre, A., & Krogstad, H. E. (1986). +"Maximum entropy estimation of the directional distribution in ocean wave spectra". +*Journal of Physical Oceanography*, 16(12), 2052-2060. + +Development of the Maximum Entropy Principle for directional wave estimation. + +Estimation Methods Comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Benoit, M., Frigaard, P., & Schäffer, H. A. (1997). +"Analyzing multidirectional wave spectra: a tentative classification of available methods". +In: *IAHR seminar multidirectional waves and their interaction with structures*. + +Comprehensive comparison of various directional spectrum estimation methods. + +Array Design +------------ + +Barber, N. F., & Ursell, F. (1948). +"The generation and propagation of ocean waves and swell". +*Philosophical Transactions of the Royal Society of London. Series A*, 240(824), 527-560. + +Early work on array-based wave measurement. + +Krogstad, H. E. (1988). +"Maximum likelihood estimation of ocean wave spectra from general arrays of wave gauges". +*Modeling, Identification and Control*, 9(2), 81-97. + +Theory for array design and optimal sensor placement. + +Applications +------------ + +Coastal Engineering +~~~~~~~~~~~~~~~~~~~ + +Goda, Y. (2010). *Random seas and design of maritime structures* (3rd ed.). +World Scientific. + +Application of directional wave spectra in coastal structure design. + +Wave Modeling +~~~~~~~~~~~~~ + +Komen, G. J., Cavaleri, L., Donelan, M., Hasselmann, K., Hasselmann, S., & Janssen, P. A. E. M. (1994). +*Dynamics and modelling of ocean waves*. Cambridge University Press. + +Use of directional spectra in numerical wave models. + +DIWASP Software +--------------- + +Original MATLAB Version +~~~~~~~~~~~~~~~~~~~~~~~ + +Johnson, D. (2002). *DIWASP, a directional wave spectrum analysis package*. +Research Report WP-1601-DJ, Centre for Water Research, University of Western Australia. + +Documentation for the original MATLAB implementation that pyDIWASP is based on. + +Available at: https://github.com/metocean/diwasp + +Standards and Guidelines +------------------------ + +Measurement Standards +~~~~~~~~~~~~~~~~~~~~~ + +ISO 19901-1:2015. *Petroleum and natural gas industries — Specific requirements for offshore structures — Part 1: Metocean design and operating considerations*. + +Includes guidelines for wave measurement and spectrum estimation. + +IAHR (1989). *List of sea-state parameters*. +Journal of Waterway, Port, Coastal, and Ocean Engineering, 115(6), 793-808. + +Standardized definitions of wave parameters. + +Related Software +---------------- + +Python Wave Analysis +~~~~~~~~~~~~~~~~~~~~ + +* **WAVEWATCH III**: Spectral wave model (FORTRAN, with Python interfaces) +* **SWAN**: Simulating WAves Nearshore model +* **oceanwaves**: Python package for wave data analysis +* **wafo**: Wave Analysis for Fatigue and Oceanography (Python port of MATLAB toolbox) + +MATLAB Tools +~~~~~~~~~~~~ + +* **DIWASP**: Original MATLAB version +* **WAFO**: Wave Analysis for Fatigue and Oceanography +* **CDIP toolbox**: Coastal Data Information Program tools + +Historical Context +------------------ + +The field of directional wave spectrum estimation has evolved significantly since the 1950s: + +* **1940s-1950s**: Initial theoretical work (Barber & Ursell, Longuet-Higgins) +* **1960s-1970s**: Development of estimation methods (Capon, Oltman-Shay) +* **1980s**: Maximum Entropy methods (Lygre & Krogstad) +* **1990s**: Iterative methods (IMLM), comparison studies +* **2000s**: Software packages (DIWASP), operational systems +* **2010s**: Advanced methods, machine learning approaches + +Online Resources +---------------- + +Organizations and Data +~~~~~~~~~~~~~~~~~~~~~~ + +* **CDIP** (Coastal Data Information Program): https://cdip.ucsd.edu/ + + Provides wave data and analysis tools + +* **NDBC** (National Data Buoy Center): https://www.ndbc.noaa.gov/ + + Real-time wave measurements from buoys + +* **Copernicus Marine Service**: https://marine.copernicus.eu/ + + European wave forecasts and reanalysis + +* **WAVEWATCH III**: https://github.com/NOAA-EMC/WW3 + + Global wave model + +Educational Resources +~~~~~~~~~~~~~~~~~~~~~ + +* **Coastal Engineering Manual**: https://www.publications.usace.army.mil/USACE-Publications/Engineer-Manuals/ + + Comprehensive coastal engineering reference + +* **Ocean Wave Spectra**: Prentice-Hall symposium proceedings (1963) + + Classic collection of papers on wave spectra + +Software Documentation +~~~~~~~~~~~~~~~~~~~~~~ + +* **NumPy**: https://numpy.org/doc/ +* **SciPy**: https://docs.scipy.org/ +* **Matplotlib**: https://matplotlib.org/stable/contents.html + +Citation +-------- + +If you use pyDIWASP in your research, please cite both pyDIWASP and the original DIWASP: + +**pyDIWASP**:: + + SBFRF. (2024). pyDIWASP: Directional Wave Spectrum Analysis in Python. + GitHub repository: https://github.com/SBFRF/pyDIWASP + +**Original DIWASP**:: + + Johnson, D. (2002). DIWASP, a directional wave spectrum analysis package. + Research Report WP-1601-DJ, Centre for Water Research, University of Western Australia. + +**Algorithms**:: + + Hashimoto, N. (1997). Analysis of the directional wave spectrum from field data. + In: Advances in Coastal Engineering Vol. 3. Ed: Liu, P.L-F. + World Scientific, Singapore. + +Contact and Support +------------------- + +* **GitHub Issues**: https://github.com/SBFRF/pyDIWASP/issues +* **Email**: support@sbfrf.org +* **Documentation**: https://pydiwasp.readthedocs.io/ (when published) diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..10478f1 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,3 @@ +sphinx>=5.0 +sphinx-rtd-theme>=1.0 +sphinx-autodoc-typehints>=1.19 diff --git a/docs/user_guide/estimation_methods.rst b/docs/user_guide/estimation_methods.rst new file mode 100644 index 0000000..c3bfd01 --- /dev/null +++ b/docs/user_guide/estimation_methods.rst @@ -0,0 +1,282 @@ +Estimation Methods +================== + +pyDIWASP implements two methods for estimating directional wave spectra from +instrument array measurements. This guide explains when to use each method. + +Available Methods +----------------- + +IMLM: Iterated Maximum Likelihood Method +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Default and recommended method for most applications.** + +The IMLM method iteratively refines the directional spectrum estimate to maximize +the likelihood that the observed cross-spectra would result from the estimated +directional spectrum. + +**Characteristics:** + +* Good balance of accuracy and computational efficiency +* Handles multi-modal directional distributions well +* Robust to noise +* Recommended for general ocean wave conditions + +**Usage:** + +.. code-block:: python + + EP = { + 'method': 'IMLM', + 'iter': 100 # number of iterations + } + + SMout, EPout = dirspec(ID, SM, EP) + +**Parameters:** + +* ``iter``: Number of iterations (default: 100) + - More iterations → better convergence but slower computation + - 50-200 iterations is typical + - Monitor convergence by checking if results stabilize + +EMEP: Extended Maximum Entropy Principle +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The EMEP method is based on the principle of maximum entropy. It finds the directional +spectrum that is "smoothest" (maximum entropy) while still consistent with the +measurements. + +**Characteristics:** + +* Works well for narrow directional spreads +* Good for swell-dominated conditions +* May over-smooth in complex sea states +* Faster than IMLM (non-iterative in original form) + +**Usage:** + +.. code-block:: python + + EP = { + 'method': 'EMEP', + 'iter': 100 # still used for consistency + } + + SMout, EPout = dirspec(ID, SM, EP) + +**Best For:** + +* Long-period swell with narrow directional spread +* Conditions with dominant single wave direction +* Quick preliminary analyses + +Comparing Methods +----------------- + +The choice between methods depends on your wave conditions: + +Wind Waves (Complex, Multi-Directional) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +* **Preferred**: IMLM +* **Reason**: Better handles multiple wave components from different directions + +Swell (Long-Period, Narrow Spread) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +* **Preferred**: EMEP or IMLM (both work well) +* **Reason**: EMEP's entropy principle matches the physics of narrow-spread swell + +Mixed Seas (Wind Waves + Swell) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +* **Preferred**: IMLM +* **Reason**: Better separation of multiple wave systems + +Example Comparison +------------------ + +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec, infospec, plotspec + import matplotlib.pyplot as plt + + # Same data, same grid + ID = {...} # your instrument data + SM = {...} # your spectral grid + + # IMLM estimation + EP_imlm = {'method': 'IMLM', 'iter': 100} + SM_imlm, EP_imlm_out = dirspec(ID, SM, EP_imlm, ['PLOTTYPE', 0]) + + # EMEP estimation + EP_emep = {'method': 'EMEP', 'iter': 100} + SM_emep, EP_emep_out = dirspec(ID, SM, EP_emep, ['PLOTTYPE', 0]) + + # Compare visually + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) + + plt.subplot(121) + plotspec(SM_imlm, 1) + plt.title('IMLM Method') + + plt.subplot(122) + plotspec(SM_emep, 1) + plt.title('EMEP Method') + + plt.tight_layout() + plt.show() + + # Compare statistics + print("IMLM:") + H_imlm, Tp_imlm, DTp_imlm, Dp_imlm = infospec(SM_imlm) + + print("\\nEMEP:") + H_emep, Tp_emep, DTp_emep, Dp_emep = infospec(SM_emep) + +Other Parameters +---------------- + +Beyond the method choice, several other parameters affect the estimation: + +Directional Resolution (dres) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Controls the number of directions used in internal calculations: + +.. code-block:: python + + EP = { + 'method': 'IMLM', + 'dres': 180 # default + } + +* Higher ``dres`` → finer directional resolution but slower computation +* Typical range: 36-360 +* Default (180) works well for most cases + +FFT Length (nfft) +~~~~~~~~~~~~~~~~~ + +Controls frequency resolution: + +.. code-block:: python + + EP = { + 'method': 'IMLM', + 'nfft': 512 # auto-calculated if omitted + } + +* Must be power of 2 (256, 512, 1024, ...) +* Larger ``nfft`` → finer frequency resolution +* Auto-calculation usually optimal +* Constrained by data length: nfft ≤ length of data + +Spectral Smoothing (smooth) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Controls whether to apply spectral smoothing: + +.. code-block:: python + + EP = { + 'method': 'IMLM', + 'smooth': 'ON' # default + } + +* ``'ON'``: Applies smoothing (recommended for noisy data) +* ``'OFF'``: No smoothing (use if data is already clean) + +Smoothing reduces noise but may blur sharp features. + +Number of Iterations (iter) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For IMLM, controls how many refinement iterations to perform: + +.. code-block:: python + + EP = { + 'method': 'IMLM', + 'iter': 100 # default + } + +* More iterations generally improve accuracy +* Returns diminish after ~50-100 iterations +* Check convergence if concerned about accuracy + +Performance Considerations +-------------------------- + +Computational Cost +~~~~~~~~~~~~~~~~~~ + +Typical computation times (3 sensors, 50 frequencies, 36 directions): + +* IMLM with 100 iterations: ~2-5 seconds +* EMEP: ~1-3 seconds + +Factors affecting speed: + +* Number of instruments: O(N²) scaling +* Number of frequencies: Linear scaling +* Number of directions: Linear scaling +* Number of iterations (IMLM only): Linear scaling + +Memory Usage +~~~~~~~~~~~~ + +Memory requirements scale with: + +* Number of cross-spectral pairs: N×N×nfreqs +* Directional resolution: O(dres) + +For typical arrays (3-10 sensors), memory is not usually a concern. + +Validation and Quality Control +------------------------------- + +After computing a spectrum, always: + +1. **Visual inspection**: Plot the spectrum and check if it looks physically reasonable + +2. **Check statistics**: Verify Hsig, Tp match expectations + +3. **Compare methods**: If results differ significantly between IMLM and EMEP, + investigate your data quality + +4. **Energy conservation**: The integrated spectrum should match the variance of + measurements + +5. **Frequency limits**: Check that spectral energy is concentrated in expected + frequency range + +When Results Are Poor +~~~~~~~~~~~~~~~~~~~~~ + +If estimated spectra look unrealistic: + +* Check array geometry (spacing, configuration) +* Verify instrument synchronization +* Inspect raw data for quality issues +* Try different frequency/direction grids +* Adjust smoothing or iteration count + +See Also +-------- + +* :doc:`understanding_spectra` for background on directional spectra +* :doc:`instrument_types` for array design considerations +* :doc:`troubleshooting` for debugging problems + +References +---------- + +Hashimoto, N. (1997). "Analysis of the directional wave spectrum from field data". +In: Advances in Coastal Engineering Vol. 3. Ed: Liu, P.L-F. +Pub: World Scientific, Singapore. + +This reference describes the mathematical basis for both IMLM and EMEP methods. diff --git a/docs/user_guide/index.rst b/docs/user_guide/index.rst new file mode 100644 index 0000000..968fbdd --- /dev/null +++ b/docs/user_guide/index.rst @@ -0,0 +1,12 @@ +User Guide +========== + +This guide provides detailed information on using pyDIWASP for directional wave spectrum analysis. + +.. toctree:: + :maxdepth: 2 + + understanding_spectra + instrument_types + estimation_methods + troubleshooting diff --git a/docs/user_guide/instrument_types.rst b/docs/user_guide/instrument_types.rst new file mode 100644 index 0000000..7abc382 --- /dev/null +++ b/docs/user_guide/instrument_types.rst @@ -0,0 +1,264 @@ +Working with Different Instrument Types +======================================== + +pyDIWASP supports various types of wave measurement instruments. This guide explains +how to configure your instrument array for each type. + +Supported Instrument Types +--------------------------- + +pyDIWASP supports the following instrument types (specified in the ``datatypes`` field): + +* ``'pres'``: Pressure sensors +* ``'elev'``: Surface elevation (e.g., wave staff, laser, ultrasonic) +* ``'velx'``: Horizontal velocity in x-direction +* ``'vely'``: Horizontal velocity in y-direction +* ``'vels'``: Horizontal velocity magnitude +* ``'accs'``: Acceleration (vertical) + +Each type has a different transfer function relating the measurement to surface elevation. + +Pressure Sensors +---------------- + +Pressure sensors are the most commonly used instruments for directional wave analysis. + +**Principle:** + +Pressure sensors measure the pressure variation due to passing waves. The relationship +between surface elevation and pressure depends on water depth and wave frequency due +to pressure attenuation with depth. + +**Configuration:** + +.. code-block:: python + + ID = { + 'layout': np.array([ + [0, 10, 5], # x positions + [0, 0, 8.66], # y positions + [0, 0, 0] # z positions (seabed = 0) + ]), + 'datatypes': ['pres', 'pres', 'pres'], + 'depth': 10.0, # water depth + 'fs': 2.0, + 'data': pressure_data # in Pascals or meters of water + } + +**Considerations:** + +* Pressure attenuation limits the maximum frequency that can be resolved +* In shallow water (depth/wavelength < 0.5), attenuation is minimal +* In deep water, high-frequency waves may not be detectable near the bottom +* Rule of thumb: kh > 1 for good measurement (k = wavenumber, h = depth) + +**Frequency Limits:** + +For depth h and frequency f, the pressure attenuation factor is: + +.. math:: + + A(f, h) = \\frac{\\cosh(kh)}{\\cosh(k(h-z))} + +where z is the sensor elevation (z=0 at seabed). As A approaches zero, the signal +becomes too weak. + +Surface Elevation Sensors +-------------------------- + +Surface elevation sensors directly measure the water surface position. + +**Types:** + +* Resistance wave staffs +* Capacitance wave staffs +* Ultrasonic sensors (downward-looking) +* Radar sensors +* Laser sensors + +**Configuration:** + +.. code-block:: python + + ID = { + 'layout': np.array([ + [0, 10, 5], # x positions + [0, 0, 8.66], # y positions + [0, 0, 0] # z positions (not critical for elevation) + ]), + 'datatypes': ['elev', 'elev', 'elev'], + 'depth': 10.0, + 'fs': 2.0, + 'data': elevation_data # in meters + } + +**Considerations:** + +* No frequency-dependent attenuation (ideal for all frequencies) +* Wave staffs work in shallow to moderate depths +* Ultrasonic/radar sensors may have limitations in extreme wave conditions +* Best for laboratory or nearshore applications + +Velocity Sensors +---------------- + +Velocity sensors (e.g., Acoustic Doppler Velocimeters, ADVs) measure water particle +velocity. + +**Configuration (separate x, y components):** + +.. code-block:: python + + ID = { + 'layout': np.array([ + [0, 0, 10, 10], # x positions + [0, 0, 0, 0], # y positions + [-2, -2, -2, -2] # z positions (depth below surface) + ]), + 'datatypes': ['velx', 'vely', 'velx', 'vely'], + 'depth': 10.0, + 'fs': 2.0, + 'data': velocity_data # in m/s + } + +**Configuration (velocity magnitude):** + +.. code-block:: python + + ID = { + 'layout': np.array([ + [0, 10, 5], + [0, 0, 8.66], + [-2, -2, -2] + ]), + 'datatypes': ['vels', 'vels', 'vels'], + 'depth': 10.0, + 'fs': 2.0, + 'data': velocity_magnitude_data + } + +**Considerations:** + +* Velocity measurements also have depth-dependent attenuation +* Horizontal velocities attenuate less rapidly than pressure +* ADVs work well in moderate-depth applications +* The measurement volume should be small compared to wavelength + +Mixed Arrays +------------ + +You can combine different instrument types in the same array: + +.. code-block:: python + + ID = { + 'layout': np.array([ + [0, 10, 20, 30], + [0, 0, 0, 0], + [0, 0, -1, -1] + ]), + 'datatypes': ['pres', 'pres', 'velx', 'vely'], + 'depth': 15.0, + 'fs': 2.0, + 'data': mixed_data + } + +This allows you to leverage the strengths of different sensor types. + +Array Geometry Guidelines +-------------------------- + +Regardless of instrument type, follow these guidelines for array design: + +Minimum Requirements +~~~~~~~~~~~~~~~~~~~~ + +* At least 3 instruments (for unique directional resolution) +* Non-collinear arrangement (not all in a line) +* Known, precise locations + +Optimal Configurations +~~~~~~~~~~~~~~~~~~~~~~ + +**Triangular Array:** + +.. code-block:: python + + # Equilateral triangle, 10m sides + layout = np.array([ + [0, 10, 5], + [0, 0, 8.66], + [0, 0, 0] + ]) + +Good for omnidirectional resolution. + +**Cross Array:** + +.. code-block:: python + + # Cross pattern, 10m arms + layout = np.array([ + [0, 10, -10, 0, 0], + [0, 0, 0, 10, -10], + [0, 0, 0, 0, 0] + ]) + +Excellent directional resolution, especially for dominant directions aligned with axes. + +Spacing Considerations +~~~~~~~~~~~~~~~~~~~~~~ + +The instrument spacing should satisfy: + +.. math:: + + d < \\frac{\\lambda_{min}}{2} + +where d is spacing and λₘᵢₙ is the shortest wavelength of interest. + +For deep water: λ ≈ 1.56 T² (T in seconds, λ in meters) + +Example: For waves up to 0.5 Hz (T=2s), λ ≈ 6.2m, so use spacing < 3m. + +Data Quality Considerations +---------------------------- + +Regardless of instrument type: + +* **Sampling rate**: Should satisfy Nyquist criterion (fs > 2*fₘₐₓ) + Typically 1-4 Hz is sufficient for ocean waves + +* **Record length**: Longer records give better statistical estimates + Typical: 20-30 minutes for ocean waves + +* **Synchronization**: All instruments must be precisely time-synchronized + +* **Calibration**: Ensure all instruments are properly calibrated to physical units + +* **Quality control**: Remove spikes, check for sensor failures, verify data ranges + +Common Issues and Solutions +---------------------------- + +**Issue**: High-frequency energy appears unrealistic + +**Solution**: Check pressure sensor depth; may be too deep for high frequencies + +**Issue**: Directional resolution is poor + +**Solution**: Check array aperture; may be too small relative to wavelength + +**Issue**: Spurious peaks in spectrum + +**Solution**: Check for electronic noise, sensor aliasing, or motion artifacts + +**Issue**: Missing low-frequency energy + +**Solution**: Record length may be too short; use longer time series + +See Also +-------- + +* :doc:`estimation_methods` for choosing the right algorithm +* :doc:`troubleshooting` for debugging common problems diff --git a/docs/user_guide/troubleshooting.rst b/docs/user_guide/troubleshooting.rst new file mode 100644 index 0000000..48b3ed0 --- /dev/null +++ b/docs/user_guide/troubleshooting.rst @@ -0,0 +1,402 @@ +Troubleshooting +=============== + +This guide helps you diagnose and fix common problems when using pyDIWASP. + +Installation Issues +------------------- + +Module Not Found +~~~~~~~~~~~~~~~~ + +**Error**: ``ModuleNotFoundError: No module named 'pydiwasp'`` + +**Solutions**: + +1. Verify installation:: + + pip list | grep pyDIWASP + +2. If not installed, install it:: + + pip install -e . # from source + # or + pip install pyDIWASP # from PyPI + +3. Check you're using the correct Python environment:: + + which python + python -c "import sys; print(sys.path)" + +Dependency Issues +~~~~~~~~~~~~~~~~~ + +**Error**: Import errors for numpy, scipy, or matplotlib + +**Solution**:: + + pip install numpy scipy matplotlib + +Data Structure Issues +--------------------- + +Invalid Data Structure +~~~~~~~~~~~~~~~~~~~~~~ + +**Error**: ``check_data`` warnings or errors + +**Common causes**: + +1. **Missing required fields**:: + + # Wrong: missing 'datatypes' + ID = { + 'layout': ..., + 'depth': 10.0, + 'fs': 2.0, + 'data': ... + } + + # Correct: + ID = { + 'layout': ..., + 'datatypes': ['pres', 'pres', 'pres'], + 'depth': 10.0, + 'fs': 2.0, + 'data': ... + } + +2. **Wrong array shapes**:: + + # Wrong: layout should be (3, N) not (N, 3) + layout = np.array([[0, 0, 0], [10, 0, 0], [20, 0, 0]]) + + # Correct: + layout = np.array([[0, 10, 20], [0, 0, 0], [0, 0, 0]]) + +3. **Mismatched dimensions**:: + + # Wrong: 3 instruments but data has 4 columns + ID['layout'] = np.array([[0, 10, 20], [0, 0, 0], [0, 0, 0]]) + ID['datatypes'] = ['pres', 'pres', 'pres'] + ID['data'] = np.random.randn(1000, 4) # 4 columns! + + # Correct: match number of columns to number of instruments + ID['data'] = np.random.randn(1000, 3) + +Computation Issues +------------------ + +Data Length Too Small +~~~~~~~~~~~~~~~~~~~~~ + +**Error**: ``Data length of N too small`` + +**Cause**: The specified ``nfft`` is larger than the data length. + +**Solutions**: + +1. Reduce ``nfft``:: + + EP = { + 'method': 'IMLM', + 'nfft': 256 # reduce from default + } + +2. Use longer data records (recommended) + +3. Let pyDIWASP auto-calculate ``nfft`` (omit the parameter) + +NaN or Inf Values +~~~~~~~~~~~~~~~~~ + +**Issue**: Spectrum contains NaN or infinite values + +**Causes and solutions**: + +1. **Constant or near-constant data**: + - Check that your data contains actual wave variations + - Remove DC offset if needed + +2. **Very small variance**: + - Check units of input data + - Verify instruments are working + +3. **Numerical instability**: + - Try different estimation method + - Reduce frequency range + - Check array geometry + +Poor Spectral Estimates +----------------------- + +Unrealistic Peak Directions +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Issue**: Estimated directions don't match expectations + +**Checks**: + +1. Verify array orientation: + - Check that layout coordinates match actual positions + - Verify x-axis direction (``xaxisdir``) + +2. Check instrument order: + - Ensure ``data`` columns match ``layout`` columns + - Verify ``datatypes`` order matches instruments + +3. Coordinate system: + - Confirm mathematical vs. nautical conventions + - Use ``compangle`` to convert between conventions + +Missing Energy at High Frequencies +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Issue**: Spectrum shows little energy at high frequencies + +**Causes**: + +1. **Pressure sensor too deep**: + - Pressure attenuates with depth + - Solution: Use sensors closer to surface or use velocimeters + +2. **Sampling rate too low**: + - Check Nyquist criterion: fs > 2*fₘₐₓ + - Solution: Increase sampling rate + +3. **Natural wave conditions**: + - May be realistic if no high-frequency waves present + +Over-Smoothed Spectrum +~~~~~~~~~~~~~~~~~~~~~~~ + +**Issue**: Spectrum looks too smooth, missing details + +**Solutions**: + +1. Disable smoothing:: + + EP = {'method': 'IMLM', 'smooth': 'OFF'} + +2. Increase frequency resolution:: + + SM = { + 'freqs': np.linspace(0.05, 0.5, 100), # more points + 'dirs': np.linspace(-180, 180, 72) + } + +3. Try different estimation method (IMLM vs EMEP) + +Noisy Spectrum +~~~~~~~~~~~~~~ + +**Issue**: Spectrum has many small spurious peaks + +**Solutions**: + +1. Enable smoothing:: + + EP = {'method': 'IMLM', 'smooth': 'ON'} + +2. Check data quality: + - Remove spikes or outliers + - Verify instrument calibration + - Check for electronic noise + +3. Reduce directional resolution:: + + EP = {'dres': 90} # coarser resolution + +4. Average multiple spectra: + - Compute spectra from overlapping windows + - Average the results + +Poor Directional Resolution +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Issue**: Cannot distinguish between nearby directions + +**Causes**: + +1. **Array too small**: + - Aperture limits resolution + - Solution: Use larger array if possible + +2. **Insufficient angular resolution**: + - Increase ``dres``:: + + EP = {'dres': 360} + +3. **Waves too short relative to array**: + - Array spacing must be appropriate for wavelength + - Solution: Adjust focus to lower frequencies + +Interpolation Warnings +~~~~~~~~~~~~~~~~~~~~~~ + +**Warning**: ``User defined grid may be too coarse`` + +**Cause**: The output grid cannot adequately represent the spectrum without +significant loss of energy. + +**Solutions**: + +1. Increase frequency resolution:: + + SM = { + 'freqs': np.linspace(0.05, 0.5, 100), # more points + 'dirs': ... + } + +2. Increase directional resolution:: + + SM = { + 'freqs': ..., + 'dirs': np.linspace(-180, 180, 72) # more points + } + +3. If the warning persists but Hsig error is small (<2%), it may be acceptable + +Visualization Issues +-------------------- + +Plot Appears Empty +~~~~~~~~~~~~~~~~~~ + +**Issue**: ``plotspec`` creates figure but shows no data + +**Checks**: + +1. Verify spectrum was computed:: + + print(SMout['S'].shape) + print(np.max(SMout['S'])) + +2. Check for NaN values:: + + print(np.sum(np.isnan(SMout['S']))) + +3. Ensure matplotlib backend is working:: + + import matplotlib.pyplot as plt + plt.plot([1, 2, 3]) + plt.show() # should display a simple plot + +Wrong Direction Convention +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Issue**: Directions appear rotated or flipped + +**Solution**: Adjust ``xaxisdir`` and use correct plot type: + +.. code-block:: python + + # For mathematical angles + plotspec(SMout, 1) # or 2 + + # For compass bearings + plotspec(SMout, 3) # or 4 + + # Adjust x-axis direction if needed + SMout['xaxisdir'] = 90 # for East = 0° in math coordinates + +Performance Issues +------------------ + +Slow Computation +~~~~~~~~~~~~~~~~ + +**Issue**: Analysis takes too long + +**Solutions**: + +1. Reduce number of iterations:: + + EP = {'iter': 50} # instead of 100 + +2. Reduce directional resolution:: + + EP = {'dres': 90} # instead of 180 + +3. Reduce number of frequencies:: + + SM = { + 'freqs': np.linspace(0.05, 0.5, 25), # fewer points + ... + } + +4. Use EMEP instead of IMLM (faster but different results) + +5. Process fewer instruments (if some are redundant) + +Memory Errors +~~~~~~~~~~~~~ + +**Issue**: Out of memory errors + +**Solutions**: + +1. Reduce ``nfft``:: + + EP = {'nfft': 256} + +2. Process data in segments + +3. Reduce number of directions:: + + EP = {'dres': 90} + +Diagnostic Checklist +-------------------- + +When encountering problems, work through this checklist: + +Input Data +~~~~~~~~~~ + +- [ ] Check data shape: ``print(ID['data'].shape)`` +- [ ] Verify no NaN: ``print(np.sum(np.isnan(ID['data'])))`` +- [ ] Check data range: ``print(np.min(ID['data']), np.max(ID['data']))`` +- [ ] Verify variance: ``print(np.var(ID['data'], axis=0))`` + +Array Configuration +~~~~~~~~~~~~~~~~~~~ + +- [ ] Check layout shape: ``print(ID['layout'].shape)`` → should be (3, N) +- [ ] Verify datatypes length: ``print(len(ID['datatypes']))`` → should match N +- [ ] Check spacing: instruments too close or too far? +- [ ] Verify depth: reasonable for your location? + +Parameters +~~~~~~~~~~ + +- [ ] Frequency range: appropriate for your waves? +- [ ] Direction range: covers all possible directions? +- [ ] Sampling frequency: satisfies Nyquist? +- [ ] Record length: sufficient for lowest frequency? + +Output +~~~~~~ + +- [ ] Check spectrum shape: ``print(SMout['S'].shape)`` +- [ ] Verify no NaN: ``print(np.sum(np.isnan(SMout['S'])))`` +- [ ] Check energy: ``print(np.max(SMout['S']))`` → should be > 0 +- [ ] Compare Hsig: does it match expectations? + +Getting Help +------------ + +If you're still stuck: + +1. **Check examples**: See the example notebook for working code + +2. **Review documentation**: Ensure you're using the API correctly + +3. **Search issues**: Check GitHub issues for similar problems + +4. **Create an issue**: If you've found a bug, report it on GitHub with: + - Minimal reproducible example + - Error message (full traceback) + - System info (Python version, OS, package versions) + +5. **Contact maintainers**: See CONTRIBUTING.md for contact information diff --git a/docs/user_guide/understanding_spectra.rst b/docs/user_guide/understanding_spectra.rst new file mode 100644 index 0000000..f8e78ae --- /dev/null +++ b/docs/user_guide/understanding_spectra.rst @@ -0,0 +1,167 @@ +Understanding Directional Wave Spectra +======================================= + +This section explains the concepts behind directional wave spectrum analysis. + +What is a Wave Spectrum? +------------------------- + +A wave spectrum describes how wave energy is distributed across different frequencies. +In nature, ocean waves are not simple sinusoids but a superposition of many different +wave components with different frequencies, amplitudes, and phases. + +A **frequency spectrum** S(f) shows the distribution of wave energy as a function of +frequency f. The units are typically m²/Hz or m²s (energy per unit frequency). + +Why Directional Spectra? +------------------------- + +While a frequency spectrum tells you *how much* energy exists at each frequency, it +doesn't tell you *where* that energy is coming from. Ocean waves can arrive from +multiple directions simultaneously: + +* Locally generated wind waves from one direction +* Swell from distant storms from another direction +* Reflected waves from coastlines or structures + +A **directional spectrum** S(f, θ) extends the frequency spectrum to include direction θ. +It shows the distribution of wave energy as a function of both frequency and direction. + +Mathematical Definition +----------------------- + +The directional spectrum S(f, θ) satisfies: + +.. math:: + + \\int_{f_1}^{f_2} \\int_{\\theta_1}^{\\theta_2} S(f, \\theta) \\, d\\theta \\, df + +This integral gives the energy in the frequency band [f₁, f₂] coming from the +directional sector [θ₁, θ₂]. + +The frequency spectrum (1D) can be recovered by integrating over all directions: + +.. math:: + + S(f) = \\int_{-\\pi}^{\\pi} S(f, \\theta) \\, d\\theta + +Key Wave Parameters +-------------------- + +From the directional spectrum, we can compute important wave statistics: + +Significant Wave Height (Hₛᵢ��) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The significant wave height is defined as: + +.. math:: + + H_{sig} = 4\\sqrt{m_0} + +where m₀ is the zeroth moment (total energy): + +.. math:: + + m_0 = \\int_{0}^{\\infty} \\int_{-\\pi}^{\\pi} S(f, \\theta) \\, d\\theta \\, df + +Hₛᵢ𝗀 approximates the average height of the highest one-third of waves. + +Peak Period (Tₚ) +~~~~~~~~~~~~~~~~~ + +The peak period is the period (reciprocal of frequency) at which the frequency +spectrum has its maximum: + +.. math:: + + T_p = 1/f_p, \\quad \\text{where} \\quad f_p = \\arg\\max_f S(f) + +Dominant Direction (Dₚ) +~~~~~~~~~~~~~~~~~~~~~~~~ + +The dominant direction is the direction from which the most energy arrives, +integrated over all frequencies: + +.. math:: + + D_p = \\arg\\max_{\\theta} \\int_{0}^{\\infty} S(f, \\theta) \\, df + +Direction at Peak (DTₚ) +~~~~~~~~~~~~~~~~~~~~~~~~ + +The direction at the spectral peak is the direction θ at the point (fₚ, θ) where +the 2D spectrum has its maximum value: + +.. math:: + + DT_p = \\arg\\max_{\\theta} S(f_p, \\theta) + +Directional Conventions +----------------------- + +Directions in oceanography can be specified in different ways: + +Mathematical Convention +~~~~~~~~~~~~~~~~~~~~~~~ + +* Measured counter-clockwise from the positive x-axis +* Range: -180° to +180° or 0° to 360° +* Used in most mathematical and computational contexts + +Nautical/Oceanographic Convention +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +* Measured clockwise from North +* Called "compass bearing" or "azimuth" +* Range: 0° to 360° (North = 0°, East = 90°, South = 180°, West = 270°) + +Direction To vs. Direction From +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +* **Direction to**: Direction waves are traveling toward (mathematical convention) +* **Direction from**: Direction waves are coming from (oceanographic convention) +* These differ by 180° + +pyDIWASP can work with either convention using the ``xaxisdir`` parameter in the +spectral matrix structure. + +Array-Based Estimation +---------------------- + +pyDIWASP estimates directional spectra from measurements by an array of instruments. +The basic principle is: + +1. Measure waves at multiple spatial locations simultaneously +2. Compute cross-spectral densities between all instrument pairs +3. The phase relationships between instruments contain directional information +4. Apply an estimation algorithm (IMLM or EMEP) to reconstruct S(f, θ) + +The array geometry is critical: + +* **Aperture**: The array size determines the shortest wavelength (highest frequency) + that can be resolved +* **Spacing**: Instrument spacing should be less than half the shortest wavelength + of interest +* **Configuration**: Triangular or cross arrays work well; linear arrays have + directional ambiguity + +For more details on array design, see the original DIWASP manual or Hashimoto (1997). + +Applications +------------ + +Directional wave spectra are used in: + +* **Coastal engineering**: Design of breakwaters, jetties, and coastal structures +* **Ship design**: Assessing vessel performance in different sea states +* **Wave modeling**: Validation of numerical wave models +* **Renewable energy**: Optimizing wave energy converter designs +* **Ocean forecasting**: Understanding wave climate and extremes + +References +---------- + +Hashimoto, N. (1997). "Analysis of the directional wave spectrum from field data". +In: Advances in Coastal Engineering Vol. 3. Ed: Liu, P.L-F. +Pub: World Scientific, Singapore. From a6af3c1042032197e6b470549c314e13aad9f6d9 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 15:40:42 +0000 Subject: [PATCH 4/7] Add comprehensive examples to README and docs README Co-authored-by: SBFRF <8375832+SBFRF@users.noreply.github.com> --- README.md | 172 ++++++++++++++++++++++++++++++++ docs/README.md | 260 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 432 insertions(+) create mode 100644 docs/README.md diff --git a/README.md b/README.md index d132d0a..620824a 100644 --- a/README.md +++ b/README.md @@ -65,6 +65,51 @@ EP = { SMout, EPout = dirspec(ID, SM, EP) ``` +### Complete Example with Output + +```python +import numpy as np +from pydiwasp import dirspec, infospec, plotspec +import matplotlib.pyplot as plt + +# Define instrument data structure +ID = { + 'layout': np.array([[0, 10, 20], [0, 0, 0], [0, 0, 0]]), # x, y, z positions + 'datatypes': ['pres', 'pres', 'pres'], # instrument types + 'depth': 10.0, # water depth in meters + 'fs': 2.0, # sampling frequency in Hz + 'data': wave_data # nsamples x ninstruments array +} + +# Define spectral matrix structure +SM = { + 'freqs': np.linspace(0.05, 0.5, 50), # frequency bins in Hz + 'dirs': np.linspace(-180, 180, 36) # direction bins in degrees +} + +# Define estimation parameters +EP = { + 'method': 'IMLM', # estimation method + 'iter': 100 # number of iterations +} + +# Compute directional spectrum +SMout, EPout = dirspec(ID, SM, EP) + +# Get wave statistics +Hsig, Tp, DTp, Dp = infospec(SMout) +# Output: +# Infospec:: +# Significant wave height: 2.5 +# Peak period: 10.0 +# Direction of peak period: 45.0 axis angle / 45.0 compass bearing +# Dominant direction: 50.0 axis angle / 40.0 compass bearing + +# Visualize the spectrum +plotspec(SMout, 1) # 3D surface plot +plt.show() +``` + For more detailed examples, see the [example notebook](examples/pyDIWASP_example.ipynb). ## Main Functions @@ -180,6 +225,109 @@ See the [examples directory](examples/) for Jupyter notebooks demonstrating: - Comparing estimation methods - Visualization options +## Examples + +### Basic Usage + +```python +import numpy as np +from pydiwasp import dirspec, infospec, plotspec + +# Load your wave data (example) +wave_data = np.loadtxt('wave_measurements.csv', delimiter=',') + +# Configure instrument array (3 pressure sensors in triangular pattern) +ID = { + 'layout': np.array([ + [0, 10, 5], # x positions (m) + [0, 0, 8.66], # y positions (m) + [0, 0, 0] # z positions (m, seabed = 0) + ]), + 'datatypes': ['pres', 'pres', 'pres'], + 'depth': 10.0, # water depth (m) + 'fs': 2.0, # sampling frequency (Hz) + 'data': wave_data +} + +# Define output grid +SM = { + 'freqs': np.linspace(0.05, 0.5, 50), # 0.05-0.5 Hz + 'dirs': np.linspace(-180, 180, 36) # full circle +} + +# Estimation parameters (use defaults) +EP = {'method': 'IMLM', 'iter': 100} + +# Compute spectrum +SMout, EPout = dirspec(ID, SM, EP) +``` + +### Analyzing Results + +```python +# Get wave statistics +Hsig, Tp, DTp, Dp = infospec(SMout) + +print(f"Significant wave height: {Hsig:.2f} m") +print(f"Peak period: {Tp:.1f} s") +print(f"Peak direction: {DTp:.0f}°") +print(f"Dominant direction: {Dp:.0f}°") +``` + +### Visualization + +```python +import matplotlib.pyplot as plt + +# Create a figure with multiple plots +fig = plt.figure(figsize=(14, 5)) + +# 3D surface plot +plt.subplot(121) +plotspec(SMout, 1) # mathematical angles + +# Polar plot with compass bearings +plt.subplot(122) +plotspec(SMout, 4) # nautical bearings + +plt.tight_layout() +plt.show() +``` + +### Comparing Methods + +```python +# Compare IMLM and EMEP methods +for method in ['IMLM', 'EMEP']: + EP = {'method': method, 'iter': 100} + SMout, EPout = dirspec(ID, SM, EP, ['PLOTTYPE', 0]) + + Hsig, Tp, DTp, Dp = infospec(SMout) + print(f"\n{method} Results:") + print(f" Hsig = {Hsig:.3f} m") + print(f" Tp = {Tp:.2f} s") +``` + +### Exporting Results + +```python +from pydiwasp import writespec + +# Write spectrum to file (DIWASP format) +writespec(SMout, 'output_spectrum.txt') + +# Or save as numpy archive +np.savez('spectrum.npz', + freqs=SMout['freqs'], + dirs=SMout['dirs'], + S=SMout['S']) +``` + +For more examples, see: +- [Example Notebook](examples/pyDIWASP_example.ipynb): Interactive Jupyter notebook +- [Examples Documentation](docs/examples.rst): Additional usage patterns +- [Test Suite](tests/): Example usage in tests + ## Contributing Contributions are welcome! Please feel free to submit pull requests or open issues for: @@ -212,6 +360,30 @@ The GNU General Public License forms the main part of the license agreement incl Copyright (C) 2002 David Johnson Coastal Oceanography Group, CWR, UWA, Perth +## Documentation + +Comprehensive documentation is available covering: + +- **[Installation Guide](docs/installation.rst)**: Detailed installation instructions +- **[Quick Start](docs/quickstart.rst)**: Get started in minutes +- **[API Reference](docs/api/index.rst)**: Complete function documentation +- **[User Guide](docs/user_guide/index.rst)**: In-depth tutorials and guides + - Understanding directional wave spectra + - Working with different instrument types + - Choosing estimation methods + - Troubleshooting common issues +- **[Examples](docs/examples.rst)**: Practical usage examples +- **[Developer Guide](docs/developer/contributing.rst)**: Contributing to pyDIWASP + +To build the documentation locally: + +```bash +cd docs +pip install -r requirements.txt +make html +# Open _build/html/index.html in your browser +``` + ## CI/CD Pipeline This repository includes GitHub Actions workflows for: diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 0000000..53a9f9d --- /dev/null +++ b/docs/README.md @@ -0,0 +1,260 @@ +# pyDIWASP Documentation + +This directory contains the source files for pyDIWASP's documentation, built using [Sphinx](https://www.sphinx-doc.org/). + +## Building the Documentation + +### Prerequisites + +Install the documentation dependencies: + +```bash +pip install -r requirements.txt +``` + +Or install with the main package: + +```bash +cd .. +pip install -e . +cd docs +pip install -r requirements.txt +``` + +### Building HTML + +To build the HTML documentation: + +```bash +make html +``` + +The generated documentation will be in `_build/html/`. Open `_build/html/index.html` in your browser to view it. + +### Other Formats + +You can build the documentation in other formats: + +```bash +make latexpdf # PDF (requires LaTeX) +make epub # ePub +make man # Man pages +``` + +### Cleaning + +To remove all built documentation: + +```bash +make clean +``` + +## Documentation Structure + +``` +docs/ +├── conf.py # Sphinx configuration +├── index.rst # Main documentation index +├── installation.rst # Installation guide +├── quickstart.rst # Quick start guide +├── examples.rst # Usage examples +├── api/ # API reference +│ ├── index.rst # API index +│ ├── dirspec.rst # dirspec function +│ ├── infospec.rst # infospec function +│ ├── plotspec.rst # plotspec function +│ ├── interpspec.rst # interpspec function +│ └── writespec.rst # writespec function +├── user_guide/ # User guide +│ ├── index.rst # User guide index +│ ├── understanding_spectra.rst # Theory and concepts +│ ├── instrument_types.rst # Instrument configuration +│ ├── estimation_methods.rst # Methods comparison +│ └── troubleshooting.rst # Common issues +├── developer/ # Developer documentation +│ ├── contributing.rst # Contribution guide +│ └── structure.rst # Code structure +├── references.rst # References and citations +├── license.rst # License information +├── requirements.txt # Documentation dependencies +├── Makefile # Build script (Unix) +├── make.bat # Build script (Windows) +└── README.md # This file +``` + +## Contributing to Documentation + +We welcome contributions to improve the documentation! Here's how: + +### Adding New Pages + +1. Create a new `.rst` file in the appropriate directory +2. Add it to the `toctree` in the parent `index.rst` +3. Build and check for errors +4. Submit a pull request + +### Editing Existing Pages + +1. Edit the `.rst` file directly +2. Build locally to verify changes +3. Check for formatting issues +4. Submit a pull request + +### Writing reStructuredText + +The documentation uses reStructuredText (reST) format. Here are some basics: + +#### Headers + +```rst +Main Title +========== + +Section +------- + +Subsection +~~~~~~~~~~ + +Subsubsection +^^^^^^^^^^^^^ +``` + +#### Code Blocks + +```rst +.. code-block:: python + + import numpy as np + from pydiwasp import dirspec + + # Your code here +``` + +#### Links + +```rst +`Link text `_ +:doc:`other_page` +:ref:`section-label` +``` + +#### Lists + +```rst +* Bullet point 1 +* Bullet point 2 + +1. Numbered item 1 +2. Numbered item 2 +``` + +#### Notes and Warnings + +```rst +.. note:: + This is a note. + +.. warning:: + This is a warning. + +.. seealso:: + See also this related topic. +``` + +#### Math + +```rst +.. math:: + + E = mc^2 + +Inline math: :math:`\\alpha = \\beta + \\gamma` +``` + +### Docstring Documentation + +API documentation is automatically generated from docstrings in the code. Follow NumPy docstring style: + +```python +def function_name(param1, param2): + """ + Brief description. + + Longer description explaining what the function does. + + Parameters + ---------- + param1 : type + Description of param1. + param2 : type + Description of param2. + + Returns + ------- + result : type + Description of return value. + + Examples + -------- + >>> result = function_name(1, 2) + >>> print(result) + 3 + + See Also + -------- + related_function : Related functionality + + Notes + ----- + Any additional notes or caveats. + + References + ---------- + .. [1] Author. "Title". Journal. Year. + """ + return param1 + param2 +``` + +### Building and Checking + +Before submitting documentation changes: + +1. **Build locally**: `make html` +2. **Check for warnings**: Sphinx will report any issues +3. **View in browser**: Open `_build/html/index.html` +4. **Check links**: Verify all internal and external links work +5. **Proofread**: Check for typos and formatting + +### Style Guidelines + +- Use clear, concise language +- Provide practical examples +- Include expected output when relevant +- Use consistent terminology +- Link to related documentation +- Keep code examples simple and focused +- Test all code examples to ensure they work + +## Documentation TODO + +Areas where documentation could be expanded: + +- [ ] More detailed mathematical background +- [ ] Additional examples for edge cases +- [ ] Performance optimization guide +- [ ] Comparison with other wave analysis tools +- [ ] Video tutorials or animated examples +- [ ] FAQ section +- [ ] Glossary of terms +- [ ] Case studies from real deployments + +## Questions? + +If you have questions about the documentation: + +- Open an issue on GitHub +- Check the [Contributing Guide](developer/contributing.rst) +- Contact the maintainers + +Thank you for helping improve pyDIWASP's documentation! From 7b8cdb9d131908aa8ee352ad041ed4e3e7d6dd67 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 15:42:11 +0000 Subject: [PATCH 5/7] Fix character encoding issues in documentation Co-authored-by: SBFRF <8375832+SBFRF@users.noreply.github.com> --- docs/license.rst | 2 +- docs/user_guide/understanding_spectra.rst | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/license.rst b/docs/license.rst index 7a8e5d0..0f6a1e3 100644 --- a/docs/license.rst +++ b/docs/license.rst @@ -91,7 +91,7 @@ Field Research Facility, US Army Corps of Engineers Current Maintenance ~~~~~~~~~~~~~~~~~~~ -Maintained by SBFRF (Systems, Biology, and Behavior Research Foundation) +Maintained by SBFRF GitHub: https://github.com/SBFRF diff --git a/docs/user_guide/understanding_spectra.rst b/docs/user_guide/understanding_spectra.rst index f8e78ae..a8dd36d 100644 --- a/docs/user_guide/understanding_spectra.rst +++ b/docs/user_guide/understanding_spectra.rst @@ -50,8 +50,8 @@ Key Wave Parameters From the directional spectrum, we can compute important wave statistics: -Significant Wave Height (Hₛᵢ��) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Significant Wave Height (Hsig) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The significant wave height is defined as: @@ -65,9 +65,9 @@ where m₀ is the zeroth moment (total energy): m_0 = \\int_{0}^{\\infty} \\int_{-\\pi}^{\\pi} S(f, \\theta) \\, d\\theta \\, df -Hₛᵢ𝗀 approximates the average height of the highest one-third of waves. +Hsig approximates the average height of the highest one-third of waves. -Peak Period (Tₚ) +Peak Period (Tp) ~~~~~~~~~~~~~~~~~ The peak period is the period (reciprocal of frequency) at which the frequency @@ -77,7 +77,7 @@ spectrum has its maximum: T_p = 1/f_p, \\quad \\text{where} \\quad f_p = \\arg\\max_f S(f) -Dominant Direction (Dₚ) +Dominant Direction (Dp) ~~~~~~~~~~~~~~~~~~~~~~~~ The dominant direction is the direction from which the most energy arrives, @@ -87,10 +87,10 @@ integrated over all frequencies: D_p = \\arg\\max_{\\theta} \\int_{0}^{\\infty} S(f, \\theta) \\, df -Direction at Peak (DTₚ) +Direction at Peak (DTp) ~~~~~~~~~~~~~~~~~~~~~~~~ -The direction at the spectral peak is the direction θ at the point (fₚ, θ) where +The direction at the spectral peak is the direction θ at the point (fp, θ) where the 2D spectrum has its maximum value: .. math:: From 31b9985ae203113d0ec5992c046274e0e51b4aa2 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 15:43:43 +0000 Subject: [PATCH 6/7] Extract version from package instead of hardcoding in docs conf Co-authored-by: SBFRF <8375832+SBFRF@users.noreply.github.com> --- docs/conf.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 2225108..13ada64 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -7,11 +7,14 @@ import sys sys.path.insert(0, os.path.abspath('..')) +# Import version from package +import pydiwasp + # -- Project information ----------------------------------------------------- project = 'pyDIWASP' copyright = '2024, SBFRF' author = 'SBFRF' -release = '0.1.0' +release = pydiwasp.__version__ # -- General configuration --------------------------------------------------- extensions = [ From 9a6883480c9838bff97714c3256350941ed3f14f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 10 Feb 2026 15:48:39 +0000 Subject: [PATCH 7/7] Fix documentation issues based on code review feedback - Fix compangle example output to match actual function behavior (225.0 instead of 45.0) - Update infospec example to show realistic compass bearing conversion - Correct writespec docstring about flattening order (direction varies fastest, not frequency) - Clarify interpspec warning behavior (only warns on Hsig increase >2%) - Document actual accepted values for funit ('hz' vs other) and dunit ('cart' vs 'naut') - Update infospec to specify it expects Hz and degrees (cartesian) - Add explicit funit/dunit to interpspec examples for clarity Co-authored-by: SBFRF <8375832+SBFRF@users.noreply.github.com> --- pydiwasp/dirspec.py | 25 ++++++++++++++++++++----- pydiwasp/infospec.py | 15 ++++++++------- pydiwasp/interpspec.py | 20 ++++++++++++++------ pydiwasp/plotspec.py | 17 +++++++++++++---- pydiwasp/writespec.py | 4 +++- 5 files changed, 58 insertions(+), 23 deletions(-) diff --git a/pydiwasp/dirspec.py b/pydiwasp/dirspec.py index ca4f97f..8a38716 100644 --- a/pydiwasp/dirspec.py +++ b/pydiwasp/dirspec.py @@ -51,18 +51,33 @@ def dirspec(ID, SM, EP, Options_=None): Spectral matrix structure defining the output grid. Required fields: * 'freqs' : ndarray - Frequency values for the output spectrum (Hz or rad/s). + Frequency values for the output spectrum. Values should match + the units specified by 'funit' (see below). * 'dirs' : ndarray - Direction values for the output spectrum (radians or degrees). + Direction values for the output spectrum. Values should match + the convention specified by 'dunit' (see below). Optional fields: * 'funit' : str, optional - Frequency units: 'Hz' (default) or 'rad/s'. + Frequency units (case-insensitive). Accepted values: + + - 'hz' (default): Frequencies in cycles per second (Hz) + - Any other value: Treated as angular frequency (rad/s) + * 'dunit' : str, optional - Direction units: 'rad' (default) or 'deg'. + Direction convention (case-insensitive). Accepted values: + + - 'cart' (default): Cartesian/mathematical directions in degrees, + measured counter-clockwise from x-axis + - 'naut': Nautical/meteorological directions in degrees, + measured clockwise from North (compass bearings) + + Note: Both conventions expect directions in degrees, not radians. + * 'xaxisdir' : float, optional - Direction of the x-axis in compass degrees (default: 90 = East). + Direction of the positive x-axis in compass degrees (default: 90 = East). + Used for coordinate transformations. EP : dict Estimation parameters. Use empty dict {} for default values. Fields: diff --git a/pydiwasp/infospec.py b/pydiwasp/infospec.py index ad00648..5e0b567 100644 --- a/pydiwasp/infospec.py +++ b/pydiwasp/infospec.py @@ -16,9 +16,10 @@ def infospec(SM): Required fields: * 'freqs' : ndarray - Frequency values (Hz or rad/s). + Frequency values. Must be in Hz (not rad/s). * 'dirs' : ndarray - Direction values (radians or degrees). + Direction values. Must be in degrees (axis angles, counter-clockwise + from x-axis). The function assumes cartesian convention. * 'S' : ndarray, shape (nfreqs, ndirs) Spectral density values. @@ -26,7 +27,7 @@ def infospec(SM): * 'xaxisdir' : float, optional Direction of x-axis in compass degrees (default: 90 = East). - Used for converting to compass bearings. + Used for converting to compass bearings in output. Returns ------- @@ -64,8 +65,8 @@ def infospec(SM): Infospec:: Significant wave height: 2.5 Peak period: 10.0 - Direction of peak period: 45.0 axis angle / 45.0 compass bearing - Dominant direction: 50.0 axis angle / 40.0 compass bearing + Direction of peak period: 45.0 axis angle / 225.0 compass bearing + Dominant direction: 50.0 axis angle / 220.0 compass bearing >>> >>> print(f"Wave height: {Hsig:.2f} m") Wave height: 2.50 m @@ -144,14 +145,14 @@ def compangle(dirs, xaxisdir): >>> # Convert 45° axis angle to compass bearing (x-axis points East) >>> bearing = compangle(45, 90) >>> print(bearing) - 45.0 + 225.0 >>> >>> # Convert multiple directions >>> import numpy as np >>> dirs = np.array([0, 45, 90, 180, 270]) >>> bearings = compangle(dirs, 90) >>> print(bearings) - [90. 45. 0. 270. 180.] + [270. 225. 180. 90. 0.] Notes ----- diff --git a/pydiwasp/interpspec.py b/pydiwasp/interpspec.py index 9ea0eb1..a432d47 100644 --- a/pydiwasp/interpspec.py +++ b/pydiwasp/interpspec.py @@ -10,7 +10,8 @@ def interpspec(SMin, SMout, method='linear'): This function resamples a directional spectrum from one frequency-direction grid to another. It uses 2D interpolation in frequency-direction space and - preserves the total wave energy (significant wave height) to within 2%. + attempts to preserve the total wave energy (significant wave height). A warning + is issued if Hsig increases by more than 2%. Parameters ---------- @@ -59,14 +60,18 @@ def interpspec(SMin, SMout, method='linear'): >>> # Compute high-resolution spectrum >>> SM_hires = { ... 'freqs': np.linspace(0.05, 0.5, 100), - ... 'dirs': np.linspace(-180, 180, 72) + ... 'dirs': np.linspace(-180, 180, 72), + ... 'funit': 'hz', + ... 'dunit': 'cart' ... } >>> SMout_hires, EPout = dirspec(ID, SM_hires, EP) >>> >>> # Define coarser output grid >>> SM_coarse = { ... 'freqs': np.linspace(0.05, 0.5, 25), - ... 'dirs': np.linspace(-180, 180, 36) + ... 'dirs': np.linspace(-180, 180, 36), + ... 'funit': 'hz', + ... 'dunit': 'cart' ... } >>> >>> # Interpolate to coarse grid @@ -77,7 +82,9 @@ def interpspec(SMin, SMout, method='linear'): >>> # Zoom in on specific frequency range >>> SM_zoom = { ... 'freqs': np.linspace(0.1, 0.2, 50), # Focus on 0.1-0.2 Hz - ... 'dirs': np.linspace(-180, 180, 36) + ... 'dirs': np.linspace(-180, 180, 36), + ... 'funit': 'hz', + ... 'dunit': 'cart' ... } >>> SMout_zoom = interpspec(SMout_hires, SM_zoom) @@ -99,9 +106,10 @@ def interpspec(SMin, SMout, method='linear'): Warnings -------- - If the significant wave height changes by more than 2% during interpolation, + If the significant wave height increases by more than 2% during interpolation, the function warns that the output grid may be too coarse. Consider using - a finer frequency or direction resolution in SMout. + a finer frequency or direction resolution in SMout. Note that the function + only warns on increases, not decreases in Hsig. See Also -------- diff --git a/pydiwasp/plotspec.py b/pydiwasp/plotspec.py index 7dde03e..079d76f 100644 --- a/pydiwasp/plotspec.py +++ b/pydiwasp/plotspec.py @@ -17,9 +17,9 @@ def plotspec(SM, ptype): Required fields: * 'freqs' : ndarray - Frequency values (Hz or rad/s). + Frequency values. Should match units specified by 'funit'. * 'dirs' : ndarray - Direction values (radians or degrees). + Direction values. Should match convention specified by 'dunit'. * 'S' : ndarray, shape (nfreqs, ndirs) Spectral density values. @@ -29,9 +29,18 @@ def plotspec(SM, ptype): Direction of x-axis in compass degrees (default: 90 = East). Used for plot types 3 and 4. * 'funit' : str, optional - Frequency units: 'Hz' or 'rad/s'. + Frequency units (case-insensitive): + + - 'hz' (default): Frequencies in Hz + - Other: Treated as rad/s + * 'dunit' : str, optional - Direction units: 'rad' or 'deg'. + Direction convention (case-insensitive): + + - 'cart' (default): Cartesian degrees (counter-clockwise from x-axis) + - 'naut': Nautical degrees (clockwise from North) + + Both conventions expect directions in degrees. ptype : int Plot type selection: diff --git a/pydiwasp/writespec.py b/pydiwasp/writespec.py index 87603e3..f4129db 100644 --- a/pydiwasp/writespec.py +++ b/pydiwasp/writespec.py @@ -53,7 +53,9 @@ def writespec(SM, filename): S[N-1,M-1] All values are written as floating point numbers, one per line. - The spectral density values S are written in row-major order (frequency varies fastest). + The spectral density values S are flattened using NumPy's default C-order, + where the last axis (direction) varies fastest and the first axis (frequency) + varies slowest. Examples --------