From 351a41b8d2c0d3e09978cbde2f98865a71f99875 Mon Sep 17 00:00:00 2001 From: BecCowley Date: Tue, 22 Mar 2022 15:55:00 +1100 Subject: [PATCH 1/2] fix: allowing down-facing instrument bin mapping for RDI 4 beams --- Preprocessing/adcpBinMappingPP.m | 40 ++++++++++++++++++++++---------- 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/Preprocessing/adcpBinMappingPP.m b/Preprocessing/adcpBinMappingPP.m index feef80220..486ca8a05 100644 --- a/Preprocessing/adcpBinMappingPP.m +++ b/Preprocessing/adcpBinMappingPP.m @@ -62,7 +62,7 @@ end end - + % We apply tilt corrections to project DIST_ALONG_BEAMS onto the vertical % axis HEIGHT_ABOVE_SENSOR. % @@ -80,32 +80,48 @@ % each other. When pitch is positive beam 1 is closer to the surface % 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. - % + + isUpwardLooking = true; + 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; + isUpwardLooking = false; 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; + % Signs for upward and downfacing instruments + if ~isUpwardLooking %for down facing adcps + sg1 = 1; sg3 = 1; + nmh = -1; + CP = cos(-pitch); + roll = roll-pi; + CR = cos(roll); + else %for up looking + sg1 = 1; sg3 = -1; + nmh = 1; + CP = cos(pitch); + CR = cos(roll); + end + %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 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); + nonMappedHeightAboveSensorBeam1 = CP .* (cos(beamAngle + sg1*roll) / cos(beamAngle) .* distAlongBeams); + nonMappedHeightAboveSensorBeam2 = CP .* (cos(beamAngle - sg1*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); + nonMappedHeightAboveSensorBeam3 = CR .* (cos(beamAngle + sg3*pitch) / cos(beamAngle) .* distAlongBeams); + nonMappedHeightAboveSensorBeam4 = CR .* (cos(beamAngle - sg3*pitch) / cos(beamAngle) .* distAlongBeams); else number_of_beams = 3; nBins = length(distAlongBeams); @@ -136,13 +152,13 @@ switch n case 1 - dvar = nonMappedHeightAboveSensorBeam1; + dvar = nonMappedHeightAboveSensorBeam1 * nmh; case 2 - dvar = nonMappedHeightAboveSensorBeam2; + dvar = nonMappedHeightAboveSensorBeam2 * nmh; case 3 - dvar = nonMappedHeightAboveSensorBeam3; + dvar = nonMappedHeightAboveSensorBeam3 * nmh; case 4 - dvar = nonMappedHeightAboveSensorBeam4; + dvar = nonMappedHeightAboveSensorBeam4 * nmh; otherwise errormsg('Beam number %s not supported', num2str(n)) end From 4a149257c3e8e1eb39b53dfc2463e5d5d9f446b5 Mon Sep 17 00:00:00 2001 From: evacougnon Date: Fri, 8 Apr 2022 12:16:59 +1000 Subject: [PATCH 2/2] feat(tests): refactor and add test for adcpBinMapping adcpBinMapping refactoring a small section on the adjusted height above sensor for each beam which is now set as a block function for each instruments (RDI and Nortek) Also added examples for testing --- .../adcp_adjusted_distances_nortek3beam.m | 60 +++++++++++++++ .../adcp_adjusted_distances_rdi4beam.m | 75 +++++++++++++++++++ Preprocessing/adcpBinMappingPP.m | 70 +++++------------ 3 files changed, 155 insertions(+), 50 deletions(-) create mode 100644 Preprocessing/+NortekADCP/adcp_adjusted_distances_nortek3beam.m create mode 100644 Preprocessing/+TeledyneADCP/adcp_adjusted_distances_rdi4beam.m diff --git a/Preprocessing/+NortekADCP/adcp_adjusted_distances_nortek3beam.m b/Preprocessing/+NortekADCP/adcp_adjusted_distances_nortek3beam.m new file mode 100644 index 000000000..f53065f94 --- /dev/null +++ b/Preprocessing/+NortekADCP/adcp_adjusted_distances_nortek3beam.m @@ -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 + diff --git a/Preprocessing/+TeledyneADCP/adcp_adjusted_distances_rdi4beam.m b/Preprocessing/+TeledyneADCP/adcp_adjusted_distances_rdi4beam.m new file mode 100644 index 000000000..6cf726231 --- /dev/null +++ b/Preprocessing/+TeledyneADCP/adcp_adjusted_distances_rdi4beam.m @@ -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 + diff --git a/Preprocessing/adcpBinMappingPP.m b/Preprocessing/adcpBinMappingPP.m index 486ca8a05..6016cae1a 100644 --- a/Preprocessing/adcpBinMappingPP.m +++ b/Preprocessing/adcpBinMappingPP.m @@ -62,16 +62,16 @@ end end - + % We apply tilt corrections to project DIST_ALONG_BEAMS onto the vertical % axis HEIGHT_ABOVE_SENSOR. % % 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 @@ -80,64 +80,34 @@ % each other. When pitch is positive beam 1 is closer to the surface % 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. - - isUpwardLooking = true; - + % + + % 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; - isUpwardLooking = false; + 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; - % Signs for upward and downfacing instruments - if ~isUpwardLooking %for down facing adcps - sg1 = 1; sg3 = 1; - nmh = -1; - CP = cos(-pitch); - roll = roll-pi; - CR = cos(roll); - else %for up looking - sg1 = 1; sg3 = -1; - nmh = 1; - CP = cos(pitch); - CR = cos(roll); - end - %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. - % H[TxB] = P[T] x (+-R[T] x B[B]) - nonMappedHeightAboveSensorBeam1 = CP .* (cos(beamAngle + sg1*roll) / cos(beamAngle) .* distAlongBeams); - nonMappedHeightAboveSensorBeam2 = CP .* (cos(beamAngle - sg1*roll) / cos(beamAngle) .* distAlongBeams); - - % H[TxB] = R[T] x (-+P[T] x B[B]) - nonMappedHeightAboveSensorBeam3 = CR .* (cos(beamAngle + sg3*pitch) / cos(beamAngle) .* distAlongBeams); - nonMappedHeightAboveSensorBeam4 = CR .* (cos(beamAngle - sg3*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); @@ -152,13 +122,13 @@ switch n case 1 - dvar = nonMappedHeightAboveSensorBeam1 * nmh; + dvar = nonMappedHeightAboveSensorBeam(:,:,n); case 2 - dvar = nonMappedHeightAboveSensorBeam2 * nmh; + dvar = nonMappedHeightAboveSensorBeam(:,:,n); case 3 - dvar = nonMappedHeightAboveSensorBeam3 * nmh; + dvar = nonMappedHeightAboveSensorBeam(:,:,n); case 4 - dvar = nonMappedHeightAboveSensorBeam4 * nmh; + dvar = nonMappedHeightAboveSensorBeam(:,:,n); otherwise errormsg('Beam number %s not supported', num2str(n)) end