Skip to content

Commit

Permalink
Dolfyn: Add horizontal speed and velocity calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
@simmsa
simmsa committed Feb 13, 2025
1 parent c98a1d5 commit 8ad973b
Show file tree
Hide file tree
Showing 2 changed files with 340 additions and 0 deletions.
143 changes: 143 additions & 0 deletions mhkit/dolfyn/tools/calculate_horizontal_speed_and_direction.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
function ds_out = calculate_horizontal_speed_and_direction(ds)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate horizontal speed and direction from velocity components
%
% Parameters
% ------------
% ds: structure
% Dolfyn data structure containing velocity data. Must include:
% - vel: structure with fields:
% - data: array with velocity components
% - dims: cell array of dimension names
% - coords: structure with coordinate information including:
% - dir: cell array of direction labels ('E', 'N', 'U1', 'U2')
% - attrs: structure with metadata
% - coord_sys: string
% Coordinate system identifier (e.g., "earth", "beam", "inst", "principal")
%
% Returns
% ------------
% ds_out: structure
% Input structure with additional fields:
% - U_mag: structure
% - data: array of velocity magnitudes (m/s)
% - dims: cell array of dimension names (excluding dir)
% - coords: coordinate information
% - units: "m s-1"
% - standard_name: "sea_water_speed"
% - long_name: "Water Speed"
% - U_dir: structure
% - data: array of flow directions
% - dims: cell array of dimension names (excluding dir)
% - coords: coordinate information
% - units: "degrees_CW_from_[REF]" where [REF] depends on coord_sys:
% - earth: "degrees_CW_from_N"
% - others: "degrees_CW_from_E" (or relevant streamwise component)
% - standard_name: "sea_water_to_direction"
% - long_name: "Water Direction"
% Note: Direction is NaN where velocity magnitude is zero
% Note: For earth coordinates, angles are CW from North [0, 360]
% For other coordinates, angles are CW from East/streamwise [0, 360]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Input validation - check for required fields
validateattributes(ds, {'struct'}, {'nonempty'}, mfilename, 'ds');
assert(isfield(ds, 'attrs'), 'MATLAB:assert:failed', 'Field attrs missing from ds');
assert(isfield(ds, 'vel'), 'MATLAB:assert:failed', 'Field vel missing from ds');
validateattributes(ds.vel, {'struct'}, {'nonempty'}, mfilename, 'ds.vel');

% Check for required vel subfields
required_fields = {'data', 'dims', 'coords'};
for i = 1:length(required_fields)
assert(isfield(ds.vel, required_fields{i}), 'MATLAB:assert:failed', ...
['Field ' required_fields{i} ' missing from ds.vel']);
end

% Check for coords.dir
assert(isfield(ds.vel.coords, 'dir'), 'MATLAB:assert:failed', ...
'Field dir missing from ds.vel.coords');
validateattributes(ds.vel.coords.dir, {'cell'}, {'nonempty'}, mfilename, 'ds.vel.coords.dir');

% Find the index of the 'dir' dimension
dir_idx = find(strcmp(ds.vel.dims, 'dir'));
assert(~isempty(dir_idx), 'MATLAB:assert:failed', ...
"'dir' dimension not found in vel.dims");

% Get the size of the data array
data_size = size(ds.vel.data);

% Ensure data array has enough dimensions
assert(length(data_size) >= dir_idx, 'MATLAB:assert:failed', ...
'Data array has fewer dimensions than the dir index');

% Find indices for East and North components in the direction coordinates
dir_coords = ds.vel.coords.dir;

E_idx = find(strcmp(dir_coords, 'E'));
N_idx = find(strcmp(dir_coords, 'N'));

% Verify both components exist
assert(~isempty(E_idx), 'MATLAB:assert:failed', ...
'East component not found in vel.coords.dir');
assert(~isempty(N_idx), 'MATLAB:assert:failed', ...
'North component not found in vel.coords.dir');

% Create indexing cell array for the entire data array
idx_cell = repmat({':'}, 1, length(data_size));

% Extract East and North velocities using dynamic indexing
idx_cell{dir_idx} = E_idx;

% Remove singleton dimensions
vel = squeeze(ds.vel.data);
u = squeeze(vel(idx_cell{:})); % East component

idx_cell{dir_idx} = N_idx;
v = squeeze(vel(idx_cell{:})); % North component

% Calculate complex velocity
U = u + v * 1j;
U_mag = abs(U);

% Prepare output structure for magnitude
ds_out = ds;
ds_out.U_mag = struct();
ds_out.U_mag.data = U_mag;
ds_out.U_mag.dims = ds.vel.dims(1:end ~= dir_idx);
ds_out.U_mag.coords = ds.vel.coords;
ds_out.U_mag.units = "m s-1";
ds_out.U_mag.standard_name = "sea_water_speed";
ds_out.U_mag.long_name = "Water Speed";

% Calculate direction
angle_rad = angle(U);
angle_deg = angle_rad * (180/pi);

% Convert angle based on coordinate system
if strcmp(ds.coord_sys, "earth")
% Convert "deg CCW from East" to "deg CW from North" [0, 360]
angle_compass = -(angle_deg - 90);
relative_to = ds.vel.coords.dir{N_idx}; % North component
else
% Switch to clockwise and from [-180, 180] to [0, 360]
angle_compass = -angle_deg;
angle_compass(angle_compass < 0) = angle_compass(angle_compass < 0) + 360;
relative_to = ds.vel.coords.dir{E_idx}; % East/streamwise component
end

% Ensure angles are in [0, 360] range
angle_compass = mod(angle_compass, 360);
% Set direction to NaN where velocity is zero
angle_compass(U_mag == 0) = NaN;

ds_out.U_dir = struct();
% Create direction output
ds_out.U_dir.data = angle_compass;
ds_out.U_dir.dims = ds.vel.dims(1:end ~= dir_idx);
ds_out.U_dir.coords = ds.vel.coords;
ds_out.U_dir.units = "degrees clockwise from " + relative_to;
ds_out.U_dir.standard_name = "sea_water_to_direction";
ds_out.U_dir.long_name = "Water To Direction";
end
197 changes: 197 additions & 0 deletions mhkit/tests/Dolfyn_Test_VAP.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,197 @@
classdef Dolfyn_Test_VAP < matlab.unittest.TestCase
methods(Test)
function test_basic_velocity_magnitude(testCase)
% Test basic velocity magnitude calculation with simple values
ds.vel.data = zeros(2,2,1,4); % time, range, dir, components
ds.vel.coords.dir = {'E', 'N', 'U1', 'U2'};
ds.vel.data(:,:,:,1) = 3 * ones(2,2,1); % East component
ds.vel.data(:,:,:,2) = 4 * ones(2,2,1); % North component
ds.vel.dims = {'time', 'range', 'dir'};
ds.vel.coords.time = 1:2;
ds.vel.coords.range = 1:2;
ds.attrs = struct();
ds.coord_sys = "earth";

result = calculate_horizontal_speed_and_direction(ds);

% Magnitude should be 5 (3-4-5 triangle)
expected_magnitude = 5 * ones(2,2,1);
testCase.verifyEqual(result.U_mag.data, expected_magnitude, 'AbsTol', 1e-10);
testCase.verifyEqual(result.U_mag.units, "m s-1");
end

function test_different_dir_index(testCase)
% Test with dir in different dimension positions
ds.vel.data = zeros(2,4,2); % time, dir, range
ds.vel.dims = {'time', 'dir', 'range'};
ds.vel.coords.dir = {'E', 'N', 'U1', 'U2'};
ds.vel.coords.time = 1:2;
ds.vel.coords.range = 1:2;
ds.attrs = struct();
ds.coord_sys = "earth";

% Set velocity components
idx_east = repmat({':'}, 1, ndims(ds.vel.data));
idx_north = repmat({':'}, 1, ndims(ds.vel.data));
idx_east{2} = 1; % E component in dir dimension
idx_north{2} = 2; % N component in dir dimension

ds.vel.data(idx_east{:}) = 3; % East component
ds.vel.data(idx_north{:}) = 4; % North component

result = calculate_horizontal_speed_and_direction(ds);

% Verify dimensions are correct (dir dimension removed)
testCase.verifyEqual(result.U_mag.dims, {'time', 'range'});
% Verify magnitude calculation
expected_magnitude = 5 * ones(2,2); % 3-4-5 triangle
testCase.verifyEqual(result.U_mag.data, expected_magnitude, 'AbsTol', 1e-10);
end

function test_squeezed_dimensions(testCase)
% Test with singleton dimensions that need squeezing
ds.vel.data = zeros(2,1,4,1,2); % time, singleton, dir, singleton, range
ds.vel.dims = {'time', 'dir', 'range'};
ds.vel.coords.dir = {'E', 'N', 'U1', 'U2'};
ds.vel.coords.time = 1:2;
ds.vel.coords.range = 1:2;
ds.attrs = struct();
ds.coord_sys = "earth";

% Set velocity components
idx_east = repmat({':'}, 1, ndims(ds.vel.data));
idx_north = repmat({':'}, 1, ndims(ds.vel.data));
idx_east{3} = 1; % E component in dir dimension
idx_north{3} = 2; % N component in dir dimension

ds.vel.data(idx_east{:}) = 3; % East component
ds.vel.data(idx_north{:}) = 4; % North component

result = calculate_horizontal_speed_and_direction(ds);

% Verify dimensions are correct after squeezing
testCase.verifyEqual(result.U_mag.dims, {'time', 'range'});
% Verify magnitude calculation
expected_magnitude = 5 * ones(2,2); % 3-4-5 triangle
testCase.verifyEqual(result.U_mag.data, expected_magnitude, 'AbsTol', 1e-10);
end

function test_zero_velocity(testCase)
% Test behavior with zero velocity
ds.vel.data = zeros(2,2,4);
ds.vel.coords.dir = {'E', 'N', 'U1', 'U2'};
ds.coord_sys = "earth";
ds.vel.dims = {'time', 'range', 'dir'};
ds.vel.coords.time = 1:2;
ds.vel.coords.range = 1:2;
ds.attrs = struct();

result = calculate_horizontal_speed_and_direction(ds);

% Verify zero magnitude
testCase.verifyEqual(sum(result.U_mag.data, 'all'), 0);
% Direction should be NaN for zero velocity
testCase.verifyTrue(isnan(sum(result.U_dir.data, 'all')));
end

function test_cardinal_directions(testCase)
% Create base dataset
ds.vel.data = zeros(2,2,1,4);
ds.vel.coords.dir = {'E', 'N', 'U1', 'U2'};
ds.vel.dims = {'time', 'range', 'dir'};
ds.vel.coords.time = 1:2;
ds.vel.coords.range = 1:2;
ds.attrs = struct();
ds.coord_sys = "earth";

% Test East direction (0 degrees mathematical, 90 degrees compass)
ds.vel.data(:,:,:,1) = 1; % East = 1
ds.vel.data(:,:,:,2) = 0; % North = 0
result = calculate_horizontal_speed_and_direction(ds);
testCase.verifyEqual(result.U_dir.data, 90 * ones(2,2,1), 'AbsTol', 1e-10);

% Test North direction (90 degrees mathematical, 0 degrees compass)
ds.vel.data(:,:,:,1) = 0; % East = 0
ds.vel.data(:,:,:,2) = 1; % North = 1
result = calculate_horizontal_speed_and_direction(ds);
testCase.verifyEqual(result.U_dir.data, 0 * ones(2,2,1), 'AbsTol', 1e-10);

% Test West direction (180 degrees mathematical, 270 degrees compass)
ds.vel.data(:,:,:,1) = -1; % East = -1
ds.vel.data(:,:,:,2) = 0; % North = 0
result = calculate_horizontal_speed_and_direction(ds);
testCase.verifyEqual(result.U_dir.data, 270 * ones(2,2,1), 'AbsTol', 1e-10);

% Test South direction (270 degrees mathematical, 180 degrees compass)
ds.vel.data(:,:,:,1) = 0; % East = 0
ds.vel.data(:,:,:,2) = -1; % North = -1
result = calculate_horizontal_speed_and_direction(ds);
testCase.verifyEqual(result.U_dir.data, 180 * ones(2,2,1), 'AbsTol', 1e-10);
end

function test_intercardinal_directions(testCase)
% Create base dataset
ds.vel.data = zeros(2,2,1,4);
ds.vel.coords.dir = {'E', 'N', 'U1', 'U2'};
ds.vel.dims = {'time', 'range', 'dir'};
ds.vel.coords.time = 1:2;
ds.vel.coords.range = 1:2;
ds.attrs = struct();
ds.coord_sys = "earth";

% Test Northeast (45 degrees mathematical, 45 degrees compass)
ds.vel.data(:,:,:,1) = 1; % East = 1
ds.vel.data(:,:,:,2) = 1; % North = 1
result = calculate_horizontal_speed_and_direction(ds);
testCase.verifyEqual(result.U_dir.data, 45 * ones(2,2,1), 'AbsTol', 1e-10);

% Test Northwest (135 degrees mathematical, 315 degrees compass)
ds.vel.data(:,:,:,1) = -1; % East = -1
ds.vel.data(:,:,:,2) = 1; % North = 1
result = calculate_horizontal_speed_and_direction(ds);
testCase.verifyEqual(result.U_dir.data, 315 * ones(2,2,1), 'AbsTol', 1e-10);
end

function test_non_earth_coordinate_system(testCase)
% Create base dataset
ds.vel.data = zeros(2,2,1,4);
ds.vel.coords.dir = {'E', 'N', 'U1', 'U2'};
ds.vel.dims = {'time', 'range', 'dir'};
ds.vel.coords.time = 1:2;
ds.vel.coords.range = 1:2;
ds.attrs = struct();
ds.coord_sys = "principal"; % Change to non-earth coordinate system

result = calculate_horizontal_speed_and_direction(ds);

% Verify that U_dir field is present
testCase.verifyTrue(isfield(result, 'U_dir'));
testCase.verifyEqual(result.U_dir.units, "degrees clockwise from E");

% Verify that magnitude is still calculated
testCase.verifyTrue(isfield(result, 'U_mag'));
end

function test_error_conditions(testCase)
% Test missing vel field
ds = struct('attrs', struct(), 'coord_sys', "earth");
testCase.verifyError(@() calculate_horizontal_speed_and_direction(ds), ...
'MATLAB:assert:failed');

% Test missing coords.dir
ds.vel.data = zeros(2,2,1,4);
ds.vel.coords = struct();
ds.vel.dims = {'time', 'range', 'dir'};
testCase.verifyError(@() calculate_horizontal_speed_and_direction(ds), ...
'MATLAB:assert:failed');

% Test missing E component
ds.vel.data = zeros(2,2,1,4);
ds.vel.coords.dir = {'N', 'U1', 'U2'}; % Remove E
ds.vel.dims = {'time', 'range', 'dir'};
testCase.verifyError(@() calculate_horizontal_speed_and_direction(ds), ...
'MATLAB:assert:failed');
end

end
end

0 comments on commit 8ad973b

Please sign in to comment.