Merge pull request #173 from simmsa/feat_dolfyn_turbulence

Feature: Add Dolfyn Turbulence Functionality

Author: simmsa

examples/adcp_example.html

@@ -1,28 +1,33 @@
 function out = water_depth_from_pressure(ds, options)
 
-%%%%%%%%%%%%%%%%%%%%
-% Calculates the distance to the water surface. Temperature and salinity
-% are used to calculate seawater density, which is in turn used with the
-% pressure data to calculate depth.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Calculate water depth from pressure using seawater density
+%
+% This function calculates the distance to the water surface using pressure
+% ADCP sensor data. Temperature and salinity are used to calculate seawater density,
+% which is then used with pressure data to calculate depth%
 %
 % Parameters
-% ----------
-% ds: Dataset
-%   The full adcp dataset
-% salinity: numeric
-%   Water salinity in psu
+% ------------
+%   ds : structure
+%       ADCP dataset structure containing pressure and temperature data
+%       ds.pressure.data : Pressure measurements [dbar]
+%       ds.temp.data : Temperature measurements [°C]  
+%       ds.attrs.h_deploy : Optional deployment height above seafloor [m]
+%   salinity : double, optional (name-value)
+%       Water salinity [psu] (default: 35)
 %
 % Returns
-% -------
-% out: Dataset
-%   adds the variables "water_density" and "depth" to the input dataset.
-%
-% Notes
-% -----
-% Requires that the instrument's pressure sensor was calibrated/zeroed
-% before deployment to remove atmospheric pressure.
+% ---------
+%   out : structure
+%       Input dataset with added depth and density fields:
+%       out.depth.data : Water depth measurements [m]
+%       out.water_density.data : Seawater density [kg/m³]
+%       out.depth.long_name : Descriptive name
+%       out.depth.units : "m"
 %
-%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     arguments
         ds

examples/adcp_example.m

@@ -1,28 +1,48 @@
 function ds=dolfyn_read(filename,options)
 
-%%%%%%%%%%%%%%%%%%%%
-%     Read a binary Nortek (e.g., .VEC, .wpr, .ad2cp, etc.) or RDI
-%    (.000, .PD0, .ENX, etc.) data file.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Read binary ADCP data files from Nortek or RDI instruments
+%
+% This function reads binary ADCP data files in .VEC, .wpr, .ad2cp (Nortek)-
+% and .000, .PD0, .ENX, (RDI) formats into MATLAB structures that can be saved
+% as .nc files and used with MHKiT-MATLAB dolfyn ADCP analysis functions.
 %
 % Parameters
 % ------------
-%     filename: string
-%         Filename of instrument file to read.
-%     userdata: bool or string (optional)
-%         true, false, or string of userdata.json filename (default true)
-%         Whether to read the '<base-filename>.userdata.json' file.
-%     nens: nan, int, or 2-element array (optional)
-%         nan (default: read entire file), int, or 2-element tuple
-%         (start, stop) Number of pings to read from the file.
-%
-%     call with options -> dolfyn_read(filename,'userdata',false,'nens',12)
+%   filename : string
+%       Path to instrument file to read
+%   userdata : logical or string, optional (name-value)
+%       Whether to read '<base-filename>.userdata.json' file
+%       - true: read userdata.json file (default)
+%       - false: skip userdata.json file  
+%       - string: path to specific userdata.json file
+%   nens : double or array, optional (name-value)
+%       Number of pings/ensembles to read
+%       - nan: read entire file (default)
+%       - scalar: read first N pings
+%       - [start, stop]: read pings from start to stop
 %
 % Returns
 % ---------
-%     ds: structure
-%         Structure from the binary instrument data
+%   ds : structure
+%       ADCP dataset structure containing:
+%       ds.vel.data : Velocity data [range x time x beam] [m/s]
+%       ds.coords : Coordinate information (time, range, beam)
+%       ds.attrs : Instrument metadata and deployment information
+%
+% Examples
+% --------
+% Read entire ADCP file, may be slow for large (>1GB) files
+% ds = dolfyn_read('deployment.ad2cp');
+%
+% Read first 1000 ensembles without userdata
+% ds = dolfyn_read('deployment.ad2cp', 'userdata', false, 'nens', 1000);
+%
+% Read specific ensemble range
+% ds = dolfyn_read('deployment.ad2cp', 'nens', [500, 1500]);
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     arguments
         filename
         options.userdata = true;

examples/adcp_example.mlx

@@ -0,0 +1,115 @@
+function tilt = calc_tilt(pitch, roll, options)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Calculate instrument tilt from pitch and roll angles.
+%
+% This function calculates the total tilt angle of an ADCP instrument from
+% pitch and roll measurements. This function warns the user if the tilt is above 5°
+% as tilts above this threshold are likely to have a negative affect on the accuracy 
+% of flow and turbulence calculations.
+%
+% Parameters
+% ------------
+%   pitch: array
+%       time series of pitch angle (forward/backward tilt)
+%   roll: array  
+%       time series of roll angle (side-to-side tilt)
+%   units: string
+%       Units of input angles: 'degrees', 'deg', 'radians' or 'rad'
+%       Default: 'degrees'
+%   output_units: string
+%       Units for output tilt: 'degrees', 'deg', 'radians' or 'rad'
+%       Default: same as input units
+%
+% Returns
+% ---------
+%   tilt: array
+%       tilt angle in specified output units
+%
+% Algorithm
+% ---------
+% tilt_rad = atan( √( tan(roll_rad)² + tan(pitch_rad)² ) )
+%
+% Example
+% -------
+% % Calculate a time series of tilt from pitch and roll time series in degrees
+% tilt = calc_tilt(pitch_data, roll_data, 'units', 'degrees');
+%
+% Notes
+% -----
+% - Large tilts (> 5°) can affect turbulence measurements
+% - This function issues warnings for tilts > 5°
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    arguments
+        pitch
+        roll
+        options.units = 'degrees'
+        options.output_units = ''
+    end
+    
+    % Validate inputs
+    if ~isnumeric(pitch) || ~isnumeric(roll)
+        error('mhkit:dolfyn:calc_tilt: Pitch and roll must be numeric arrays');
+    end
+    
+    if ~isequal(size(pitch), size(roll))
+        error('mhkit:dolfyn:calc_tilt: Pitch and roll arrays must have the same size');
+    end
+    
+    % Validate units
+    valid_units = {'degrees', 'radians', 'deg', 'rad'};
+    if ~ismember(lower(options.units), valid_units)
+        error('mhkit:dolfyn:calc_tilt: Units must be ''degrees'' or ''radians''');
+    end
+    
+    % Set default output units
+    if isempty(options.output_units)
+        options.output_units = options.units;
+    end
+    
+    if ~ismember(lower(options.output_units), valid_units)
+        error('mhkit:dolfyn:calc_tilt: Output units must be ''degrees'' or ''radians''');
+    end
+    
+    % Normalize unit names
+    input_is_degrees = ismember(lower(options.units), {'degrees', 'deg'});
+    output_is_degrees = ismember(lower(options.output_units), {'degrees', 'deg'});
+    
+    % Convert to radians if needed for calculation
+    if input_is_degrees
+        pitch_rad = deg2rad(pitch);
+        roll_rad = deg2rad(roll);
+    else
+        pitch_rad = pitch;
+        roll_rad = roll;
+    end
+    
+    % Calculate tilt using trigonometric relationship
+    % Source: mhkit_python:dolfyn.rotate.base.calc_tilt, author @jmcvey3
+    tilt_rad = atan(sqrt(tan(roll_rad).^2 + tan(pitch_rad).^2));
+
+    % Quality assessment and warnings
+    % Convert tilt to degrees for quality assessment regardless of output units
+    tilt_deg = rad2deg(tilt_rad);
+    max_tilt_deg = max(tilt_deg(:));
+    mean_tilt_deg = mean(tilt_deg(:));
+    
+    % Issue warnings for large tilts
+    if max_tilt_deg > 5
+        warning('mhkit:dolfyn: Maximum tilt %.1f° exceeds recommended 5° limit for accurate turbulence measurements', max_tilt_deg);
+        fprintf('Tilt Analysis Summary:\n');
+        fprintf('  Mean tilt: %.2f°\n', mean_tilt_deg);
+        fprintf('  Max tilt:  %.2f°\n', max_tilt_deg);
+        fprintf('  Std tilt:  %.2f°\n', std(tilt_deg(:)));
+    end
+    
+    % Convert to desired output units
+    if output_is_degrees
+        tilt = rad2deg(tilt_rad);
+    else
+        tilt = tilt_rad;
+    end
+end

mhkit/dolfyn/adp/water_depth_from_pressure.m

@@ -45,14 +45,19 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-   if nargin < 3
-        dim_to_find = 'time';
-    end
-    if nargin < 2
-        error('n_samples is required');
-    end
+arguments
+    ds (1,1) struct
+    n_samples (1,1) double {mustBePositive, mustBeInteger}
+    dim_to_find {mustBeTextScalar} = 'time'
+end
 
     ds_out = ds;  % Start with a copy of the input
+    
+    % Store n_bin in attrs for compatibility with turbulence functions
+    if ~isfield(ds_out, 'attrs')
+        ds_out.attrs = struct();
+    end
+    ds_out.attrs.n_bin = n_samples;
 
     % Validate dimension name exists in at least one field
     dim_exists = false;
@@ -124,14 +129,30 @@
                     inv_perm_order(perm_order) = 1:length(perm_order);
                     ds_out.(current_field).data = permute(reshape(mean_data, [n_bins sz(perm_order(2:end))]), inv_perm_order);
 
+                    % Calculate standard deviation for velocity fields (following Python DOLfYN)
+                    if strcmp(current_field, 'vel') || contains(current_field, 'vel')
+                        % Calculate standard deviation along first dimension
+                        std_data = squeeze(std(reshaped, 0, 1, 'omitnan'));
+                        
+                        % Create vel_std field
+                        std_field_name = [current_field '_std'];
+                        ds_out.(std_field_name) = ds_out.(current_field); % Copy structure
+                        ds_out.(std_field_name).data = permute(reshape(std_data, [n_bins sz(perm_order(2:end))]), inv_perm_order);
+                        ds_out.(std_field_name).long_name = 'Velocity Standard Deviation';
+                        ds_out.(std_field_name).description = 'Standard deviation calculated during ensemble averaging';
+                    end
+
                     % Update coordinates for this dimension if they exist
                     if isfield(ds.(current_field), 'coords')
                         coord_fields = fieldnames(ds.(current_field).coords);
                         for k = 1:length(coord_fields)
                             coord_field = coord_fields{k};
                             if contains(lower(coord_field), lower(dim_to_find))
                                 coord_data = ds.(current_field).coords.(coord_field);
-                                ds_out.(current_field).coords.(coord_field) = coord_data(1:n_samples:usable_samples);
+                                % For numeric data (like Unix timestamps), take arithmetic mean of each bin
+                                % This matches Python's sequential chunking with mean per bin
+                                coord_reshaped = reshape(coord_data(1:usable_samples), n_samples, n_bins);
+                                ds_out.(current_field).coords.(coord_field) = mean(coord_reshaped, 1, 'omitnan')';
                             end
                         end
                     end
@@ -147,8 +168,41 @@
             if contains(lower(coord_fields{i}), lower(dim_to_find))
                 coord_data = ds.coords.(coord_fields{i});
                 usable_samples = floor(length(coord_data)/n_samples) * n_samples;
-                ds_out.coords.(coord_fields{i}) = coord_data(1:n_samples:usable_samples);
+                n_bins = usable_samples / n_samples;
+                % For numeric data (like Unix timestamps), take arithmetic mean of each bin
+                % This matches Python's sequential chunking with mean per bin
+                coord_reshaped = reshape(coord_data(1:usable_samples), n_samples, n_bins);
+                ds_out.coords.(coord_fields{i}) = mean(coord_reshaped, 1, 'omitnan')';
             end
         end
     end
+    
+    % Calculate U_std from horizontal velocity components (following Python DOLfYN)
+    if isfield(ds_out, 'vel') && isfield(ds_out, 'vel_std')
+        % Calculate horizontal velocity magnitude standard deviation
+        if size(ds_out.vel.data, 3) >= 2  % Check we have at least u and v components
+            u_mean = ds_out.vel.data(:, :, 1);
+            v_mean = ds_out.vel.data(:, :, 2);
+            u_std = ds_out.vel_std.data(:, :, 1);
+            v_std = ds_out.vel_std.data(:, :, 2);
+            
+            % Calculate U_mag standard deviation using error propagation formula
+            % U_std = sqrt((u*u_std)^2 + (v*v_std)^2) / U_mag
+            u_mag = sqrt(u_mean.^2 + v_mean.^2);
+            u_std_mag = sqrt((u_mean .* u_std).^2 + (v_mean .* v_std).^2) ./ u_mag;
+            
+            % Handle division by zero
+            u_std_mag(u_mag == 0) = 0;
+            
+            % Create U_std field
+            ds_out.U_std = struct();
+            ds_out.U_std.data = single(u_std_mag);
+            ds_out.U_std.dims = ds_out.vel.dims(1:2);  % Remove direction dimension
+            ds_out.U_std.coords = ds_out.vel.coords;
+            ds_out.U_std.coords = rmfield(ds_out.U_std.coords, 'dir');  % Remove dir coordinate
+            ds_out.U_std.units = "m s-1";
+            ds_out.U_std.long_name = "Water Velocity Standard Deviation";
+            ds_out.U_std.description = 'Horizontal velocity magnitude standard deviation from ensemble averaging';
+        end
+    end
 end

mhkit/dolfyn/io/dolfyn_read.m

@@ -288,11 +288,12 @@ function plot_3d_data(current_dim)
         % Set colormap based on field type
         switch field
             case 'vel'
-                colormap(bluewhitered_colormap(256))
+                colormap(cmocean('balance', 256))
             case 'corr'
-                colormap(viridis_colormap(256))
+                colormap(cmocean('haline', 256))
             otherwise
-                colormap('default')  % Use MATLAB's default colormap
+                % Default colormap
+                colormap(cmocean('thermal', 256))
         end
 
         % Set colorbar limits if provided