Skip to content

Commit

Permalink
Merge pull request #777 from aodn/bug_binmapping
Browse files Browse the repository at this point in the history 
Bug binmapping
  • Loading branch information
@lbesnard
lbesnard authored Apr 13, 2022
2 parents ff58bf2 + 4a14925 commit 48fcf66
Show file tree
Hide file tree
Showing 3 changed files with 153 additions and 32 deletions.
60 changes: 60 additions & 0 deletions Preprocessing/+NortekADCP/adcp_adjusted_distances_nortek3beam.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
function nonMappedHeightAboveSensorBeam = adcp_adjusted_distances_nortek3beam(roll, pitch, distAlongBeams, beamAngle, number_of_beams)
% function nonMappedHeightAboveSensorBeam = adcp_adjusted_distances_nortek3beam(roll, pitch, distAlongBeams, beamAngle, number_of_beams)
%
% Adjusts height above sensor for each beam in adcpBinMapping.m
% Can read 3 beam Nortek instruments.
% Uses the beam angle and the roll/pitch of the instrument.
% No adjustement when tilt is zero,
% nonMappedHeightAboveSensorBeam will be the same than distAlongBeams.
%
% Inputs:
%
% roll - 1D single precision array (Ntime), the internal roll angles in degrees
% pitch - 1D single precision array (Ntime), the internal pitch angles in degrees
% distAlongBeams - 1D array (Nbins), the internal distance along beam for each bin
% beamAngle - scalar, beam angle from the instrument
% number_of_beams - interger, number of beams on the ADCP
%
% Outputs:
%
% nonMappedHeightAboveSensorBeam - 3D single precision array (Ntime, Nbins, Nbeams),
% tilt adjusted height above sensor for
% each of the 3 beams
%
% Example:
%
% roll_ini = [-180:10:180]';
% pitch_ini = [-180:10:180]';
% distAlongBeams = [20:20:600];
% beamAngle = 0.3;
% number_of_beams = 3;
%
% x = NortekADCP.adcp_adjusted_distances_nortek3beam(roll_ini,pitch_ini,distAlongBeams,beamAngle,number_of_beams);
%
% idx0 = find(roll_ini == 0);
% for i=1:number_of_beams; assert(isequal(x(idx0,:,i), distAlongBeams)); assert(~isequal(x(~idx0,:,i), distAlongBeams)); end
%
%

narginchk(5,5)

number_of_beams = 3;
nBins = length(distAlongBeams);

% set to single precision
nonMappedHeightAboveSensorBeam = nan(length(pitch), length(distAlongBeams), number_of_beams, 'single');

nonMappedHeightAboveSensorBeam(:,:,1) = (cos(beamAngle - pitch) / cos(beamAngle)) * distAlongBeams;
nonMappedHeightAboveSensorBeam(:,:,1) = repmat(cos(roll), 1, nBins) .* nonMappedHeightAboveSensorBeam(:,:,1);

beamAngleX = atan(tan(beamAngle) * cos(60 * pi / 180)); % beams 2 and 3 angle projected on the X axis
beamAngleY = atan(tan(beamAngle) * cos(30 * pi / 180)); % beams 2 and 3 angle projected on the Y axis

nonMappedHeightAboveSensorBeam(:,:,2) = (cos(beamAngleX + pitch) / cos(beamAngleX)) * distAlongBeams;
nonMappedHeightAboveSensorBeam(:,:,2) = repmat(cos(beamAngleY + roll) / cos(beamAngleY), 1, nBins) .* nonMappedHeightAboveSensorBeam(:,:,2);

nonMappedHeightAboveSensorBeam(:,:,3) = (cos(beamAngleX + pitch) / cos(beamAngleX)) * distAlongBeams;
nonMappedHeightAboveSensorBeam(:,:,3) = repmat(cos(beamAngleY - roll) / cos(beamAngleY), 1, nBins) .* nonMappedHeightAboveSensorBeam(:,:,3);

end

75 changes: 75 additions & 0 deletions Preprocessing/+TeledyneADCP/adcp_adjusted_distances_rdi4beam.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
function nonMappedHeightAboveSensorBeam = adcp_adjusted_distances_rdi4beam(roll, pitch, pitchSign, distAlongBeams, beamAngle, number_of_beams)
% function nonMappedHeightAboveSensorBeam = adcp_adjusted_distances_rdi4beam(roll, pitch, pitchSign, distAlongBeams, beamAngle, number_of_beams)
%
% Adjusts height above sensor for each beam in adcpBinMapping.m
% Can read 4 beam RDI instruments.
% Uses the beam angle and the roll/pitch of the instrument.
% No adjustement when tilt is zero,
% nonMappedHeightAboveSensorBeam will be the same than distAlongBeams.
% Deals with up and down facing ADCPs
%
% Inputs:
%
% roll - 1D single precision array (Ntime), the internal roll angles in degrees
% pitch - 1D single precision array (Ntime), the internal pitch angles in degrees
% pitchSign - integer (-1 or 1), pitch sign based on the ADCP's orientation
% distAlongBeams - 1D array (Nbins), the internal distance along beam for each bin
% beamAngle - scalar, beam angle from the instrument
% number_of_beams - interger, number of beams on the ADCP
%
% Outputs:
%
% nonMappedHeightAboveSensorBeam - 3D single precision array (Ntime, Nbins, Nbeams),
% tilt adjusted height above sensor for
% each of the 4 beams
%
% Example:
%
% %testing for up and down facing ADCPs - no change for beam 1 and 2 (roll)
% %changes only for beam 3 and 4 (pitch)
%
% roll_ini = [-180:10:180]';
% pitch_ini = [-180:10:180]';
% pitchSign_up = 1;
% pitchSign_down = -1;
% distAlongBeams = [20:20:600];
% beamAngle = 0.3;
% number_of_beams = 4;
%
% x_up = TeledyneADCP.adcp_adjusted_distances_rdi4beam(roll_ini,pitch_ini,pitchSign_up,distAlongBeams,beamAngle,number_of_beams);
% x_down = TeledyneADCP.adcp_adjusted_distances_rdi4beam(roll_ini,pitch_ini,pitchSign_down,distAlongBeams,beamAngle,number_of_beams);
%
% assert(isequal(x_up(:,:,1), x_down(:,:,1)));
% assert(isequal(x_up(:,:,2), x_down(:,:,2)));
% assert(~isequal(x_up(:,:,3), x_down(:,:,3)));
% assert(~isequal(x_up(:,:,4), x_down(:,:,4)));
%
% idx0 = find(roll_ini == 0);
% for i=1:number_of_beams; assert(isequal(x_up(idx0,:,i), distAlongBeams)); assert(isequal(x_down(idx0,:,i), distAlongBeams)); end
% for i=1:number_of_beams; assert(~isequal(x_up(~idx0,:,i), distAlongBeams)); assert(~isequal(x_down(~idx0,:,i), distAlongBeams)); end
%
%

narginchk(6,6)

nmh=1;
if pitchSign < 0
% adjust nonMappedHeightAboveSensorBeam sign when downfacing
nmh = -1;
% adjust roll when downfacing
roll = roll-pi;
end

CP = cos(pitchSign * pitch);
CR = cos(roll);

% set to single precision
nonMappedHeightAboveSensorBeam = nan(length(pitch), length(distAlongBeams), number_of_beams, 'single');

nonMappedHeightAboveSensorBeam(:,:,1) = nmh * (CP .* (cos(beamAngle + roll) / cos(beamAngle) .* distAlongBeams));
nonMappedHeightAboveSensorBeam(:,:,2) = nmh * (CP .* (cos(beamAngle - roll) / cos(beamAngle) .* distAlongBeams));
nonMappedHeightAboveSensorBeam(:,:,3) = nmh * (CR .* (cos(beamAngle - pitchSign * pitch) / cos(beamAngle) .* distAlongBeams));
nonMappedHeightAboveSensorBeam(:,:,4) = nmh * (CR .* (cos(beamAngle + pitchSign * pitch) / cos(beamAngle) .* distAlongBeams));

end

50 changes: 18 additions & 32 deletions Preprocessing/adcpBinMappingPP.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,9 @@
% RDI 4 beams ADCPs:
% It is assumed that the beams are in a convex configuration such as beam 1
% and 2 (respectively 3 and 4) are aligned on the pitch (respectively roll)
% axis. When pitch is positive beam 3 is closer to the surface while beam 4
% gets further away. When roll is positive beam 2 is closer to the surface
% while beam 1 gets further away.
% axis. For an upward looking ADCP, when pitch is positive beam 3 is closer
% to the surface while beam 4 gets further away. When roll is positive beam
% 2 is closer to the surface while beam 1 gets further away.
%
% Nortek 3 beams ADCPs:
% It is assumed that the beams are in a convex configuration such as beam 1
Expand All @@ -81,47 +81,33 @@
% while beams 2 and 3 get further away. When roll is positive beam 3 is
% closer to the surface while beam 2 gets further away.
%

% Set pitch sign to deal with downward facing ADCPs
pitchSign = 1;

distAlongBeams = sample_data{k}.dimensions{distAlongBeamsIdx}.data';
if all(diff(distAlongBeams)<0)
%invert distAlongBeams so we have increasing values
% this is required for the interpolation function.
distAlongBeams = distAlongBeams*-1;
pitchSign = -1;
end
pitch = sample_data{k}.variables{pitchIdx}.data * pi / 180;
roll = sample_data{k}.variables{rollIdx}.data * pi / 180;
beamAngle = sample_data{k}.meta.beam_angle * pi / 180;

%TODO: the adjusted distances should be exported
% as variables to enable further diagnostic plots and debugging.
%TODO: the adjusted distances should be a Nx4 array or Nx3 array.
if isRDI% RDI 4 beams
if isRDI && absic4Idx ~= 0 % RDI 4 beams
number_of_beams = 4;
%TODO: block function.
CP = cos(pitch);
% H[TxB] = P[T] x (+-R[T] x B[B])
nonMappedHeightAboveSensorBeam1 = CP .* (cos(beamAngle + roll) / cos(beamAngle) .* distAlongBeams);
nonMappedHeightAboveSensorBeam2 = CP .* (cos(beamAngle - roll) / cos(beamAngle) .* distAlongBeams);

% H[TxB] = R[T] x (-+P[T] x B[B])
CR = cos(roll);
nonMappedHeightAboveSensorBeam3 = CR .* (cos(beamAngle - pitch) / cos(beamAngle) .* distAlongBeams);
nonMappedHeightAboveSensorBeam4 = CR .* (cos(beamAngle + pitch) / cos(beamAngle) .* distAlongBeams);
else
number_of_beams = 3;
nBins = length(distAlongBeams);
%TODO: block function, include tests.
% Nortek 3 beams
nonMappedHeightAboveSensorBeam1 = (cos(beamAngle - pitch) / cos(beamAngle)) * distAlongBeams;
nonMappedHeightAboveSensorBeam1 = repmat(cos(roll), 1, nBins) .* nonMappedHeightAboveSensorBeam1;

beamAngleX = atan(tan(beamAngle) * cos(60 * pi / 180)); % beams 2 and 3 angle projected on the X axis
beamAngleY = atan(tan(beamAngle) * cos(30 * pi / 180)); % beams 2 and 3 angle projected on the Y axis
nonMappedHeightAboveSensorBeam = TeledyneADCP.adcp_adjusted_distances_rdi4beam(roll, pitch, pitchSign, distAlongBeams, beamAngle, number_of_beams);

elseif isNortek && absic4Idx == 0 % Nortek 3 beams
number_of_beams = 3;

nonMappedHeightAboveSensorBeam2 = (cos(beamAngleX + pitch) / cos(beamAngleX)) * distAlongBeams;
nonMappedHeightAboveSensorBeam2 = repmat(cos(beamAngleY + roll) / cos(beamAngleY), 1, nBins) .* nonMappedHeightAboveSensorBeam2;
nonMappedHeightAboveSensorBeam = NortekADCP.adcp_adjusted_distances_nortek3beam(roll, pitch, distAlongBeams, beamAngle, number_of_beams);

nonMappedHeightAboveSensorBeam3 = (cos(beamAngleX + pitch) / cos(beamAngleX)) * distAlongBeams;
nonMappedHeightAboveSensorBeam3 = repmat(cos(beamAngleY - roll) / cos(beamAngleY), 1, nBins) .* nonMappedHeightAboveSensorBeam3;
end

nSamples = length(pitch);
Expand All @@ -136,13 +122,13 @@

switch n
case 1
dvar = nonMappedHeightAboveSensorBeam1;
dvar = nonMappedHeightAboveSensorBeam(:,:,n);
case 2
dvar = nonMappedHeightAboveSensorBeam2;
dvar = nonMappedHeightAboveSensorBeam(:,:,n);
case 3
dvar = nonMappedHeightAboveSensorBeam3;
dvar = nonMappedHeightAboveSensorBeam(:,:,n);
case 4
dvar = nonMappedHeightAboveSensorBeam4;
dvar = nonMappedHeightAboveSensorBeam(:,:,n);
otherwise
errormsg('Beam number %s not supported', num2str(n))
end
Expand Down

0 comments on commit 48fcf66

Please sign in to comment.