Rotation of Velocities to East-North-Up Co-ordinate System
Velocity data acquired from RDI and Nortek ADCPs can take on several forms. The original, unaltered velocity data, read from the raw files is normally presented in the MAT file product as adcp.velocity (RDI) or data.velocity (Nortek). (Please note that the data structure names for RDI is adcp and for Nortek we use data, subsitute data for adcp when looking at Nortek data). Also in the MAT file product is the config struct that contains information on how the data was collected. Based on the configuration, the original data is processed to present the velocities relative to East-North-Up: this is the final data in the netCDF, MAT files and in the daily current plot PNG and PDF (see adcp.u, adcp.v, adcp.w in the MAT file and in the MAT file description below). u,v,w velocities are always relative to geographic North and local gravity (for Up).
...
| Expand |
|---|
| title | Nortek ADCP Beam co-ordinates transform to True Earth |
|---|
|
| Code Block |
|---|
velocityTrue = nan(3, n*p);
maxNumPosData = max([length(heading) length(pitch) length(roll)]); % either 1 or p, the number of time stamps
for i = 1:maxNumPosData
% heading/pitch/roll matrices
% this rotation is from Nortek: http://www.nortek-as.com/lib/forum-attachments/coordinate-transformation/view
% I considered various ways to get the index - pointers, use an eval, etc, but using 'min(length(heading),i)' is actually faster
hThis = (heading(min(length(heading),i)) - 90) * pi/180; % convert heading to yaw, the Nortek code doesn't have a '-' here, but the sine term is negated below
pThis = pitch(min(length(pitch),i)) * pi/180;
rThis = roll(min(length(roll),i)) * pi/180;
H = [cos(hThis) sin(hThis) 0; -sin(hThis) cos(hThis) 0; 0 0 1];
PR = [cos(pThis) -sin(pThis)*sin(rThis) -cos(rThis)*sin(pThis);...
0 cos(rThis) -sin(rThis); ...
sin(pThis) sin(rThis)*cos(pThis) cos(pThis)*cos(rThis)];
% apply heading/pitch/roll rotations
if maxNumPosData == 1
velocityTrue = H * PR * T * beamReshape;
else
velocityTrue(:, (i-1)*n+1:1:i*n) = H * PR * T * beamReshape(:, (i-1)*n+1:i*n);
end
end
% extract w values (u,v below) - use permute to ensure we don't lose singleton dimensions
velocityTrue = reshape(velocityTrue, [3 p n]);
ADCP.u = permute(velocityTrue(1,:,:), [3 2 1]);
ADCP.v = permute(velocityTrue(2,:,:), [3 2 1]);
ADCP.w = permute(velocityTrue(3,:,:), [3 2 1]); |
|
Correction for Tilt: Bin-mapping
In addition to coordinate system rotations, the orientation of the instrument affects the vertical range on which the measurement bins are spaced (depth cell and measurement bin are synonymous). If the instrument is significantly tilted, the bins used to determine the velocities at each range will be a different water depths. Here is Figure 21 from Acoustic Doppler Current Profiler: Principles of Operation - A Practical Primer, notice in (B) that the highlighted depth cells are not the same as in (A):
...
| Expand |
|---|
| title | Ott method bin-mapping |
|---|
|
| Code Block |
|---|
%% do the Ott method bin mapping, including Ott's filling/smoothing method
binSpacing = mean(diff(range));
minRangeExtrap = min(range) - binSpacing;
maxRangeExtrap = max(range) + binSpacing;
for i = 1:m % loop over beams
% do this calc here since it doesn't change in time
matrixBeamElevationAzimuth = [-cosd(beamElevation(i))*sind(beamAzimuth(i)); cosd(beamElevation(i))*cosd(beamAzimuth(i)); -sind(beamElevation(i))] ./ sind(abs(beamElevation(i)));
for k = 1:p % loop over time stamps
iP = min(length(pitch),k);
iR = min(length(roll),k);
% NaN any pings that don't have good pitch/roll
pThis = pitch(iP);
rThis = roll(iR);
if any([pitchOverLim(iP) rollOverLim(iR)])
beam(i,:,k) = NaN;
doMaxPitchRollWarning = true;
continue
end
if any([isnan(pThis) isnan(rThis)])
beam(i,:,k) = NaN;
doNaNPitchRollWarning = true;
continue
end
% skip any pings that are all NaN
thisBeamTime = squeeze(beam(i,:,k));
iNaNthisBeamTime = isnan(thisBeamTime);
numNaNthisBeamTime = sum(iNaNthisBeamTime);
if numNaNthisBeamTime >= min(n * maxFracNaN, n - 1) % skip this beam time if too much data is NaN (need at least 2 valid points)
beam(i,:,k) = NaN;
continue
end
% for this time and beam, fill the NaNs with linear interpolation, except where there are too many NaNs in a row
% this won't fill on outside of valid data - linear interpolation doesn't extrapolate by default
if any(iNaNthisBeamTime)
% find consecutive NaNs for this beam and time, only if there are enough of them to flag
iTooManyConsecNaNs = false(1, n);
if numNaNthisBeamTime >= maxConsecutiveNaNs
numConsec = 0;
for j = find(iNaNthisBeamTime, 1, 'first') : find(iNaNthisBeamTime, 1, 'last') % loop over depth bins
if iNaNthisBeamTime(j)
numConsec = numConsec + 1;
if numConsec > maxConsecutiveNaNs
iTooManyConsecNaNs((j-numConsec+1):j) = true;
end
else
numConsec = 0;
end
end
end
iNaNsInterp = iNaNthisBeamTime & ~iTooManyConsecNaNs;
thisBeamTime(iNaNsInterp) = interp1(binIndex(~iNaNthisBeamTime), thisBeamTime(~iNaNthisBeamTime), binIndex(iNaNsInterp));
end
% this is simplicaton of the commented code below
Mp = [1 0 0;0 cos(pThis) -sin(pThis);0 sin(pThis) cos(pThis)];
Mr = [cos(rThis) 0 sin(rThis);0 1 0;-sin(rThis) 0 cos(rThis)];
beamXYZ = Mr * Mp * matrixBeamElevationAzimuth;
actualVerticalRange = range * abs(beamXYZ(3));
% % range calculation code from beam2uvw.m
% Mp = [1 0 0;0 cos(pThis) -sin(pThis);0 sin(pThis) cos(pThis)];
% Mr = [cos(rThis) 0 sin(rThis);0 1 0;-sin(rThis) 0 cos(rThis)];
% beamXYZ = Mr * Mp * [-cosd(beamElevation(i))*sind(beamAzimuth(i)); cosd(beamElevation(i))*cosd(beamAzimuth(i)); -sind(beamElevation(i))];
% actualVerticalRange = range * abs(beamXYZ(3) ./ sind(abs(beamElevation(i))));
% map the beam bins to the vertical range - allow some extrapolation - won't extrap through NaNs though (good).
% Extrapolation limited to +/- one bin spacing.
beam(i,:,k) = interp1( actualVerticalRange, thisBeamTime, range, 'linear', 'extrap' );
beam(i, actualVerticalRange > maxRangeExtrap | actualVerticalRange < minRangeExtrap, k) = NaN;
end
end |
|
Three-Beam Solutions and Screening (For RDI ADCPs Only)
In systems with four beams, the transformation from radial velocites to Earth co-ordinates velocities has redundant information; only three beams are needed to resolve currents in the three orthogonal directions, while the 4th beam effectively provides error estimation. In the case where exactly one bin out of of the four at any level is flagged and NaN'ed, then currents can still be calculated from the remaining three radial velocities, using the three-beam transformation solution. The screening steps that can flag the data 'NaN' prior to transformation include fish detection and correlation threshold, as well as any bin-mapping (particularly the RDI nearest vertical bin which can set many bins to NaN). The three-beam solution works by assigning the error velocity to be zero and then solving the transformation for the missing bin (for more details see Three Beam Solutions in adcp coordinate transformation_Jan08.pdf). With the missing data filled in, the normal transform from Beam to Instrument co-ordinate data is applied.
...
| Expand |
|---|
| title | Three beam solution and transform |
|---|
|
| Code Block |
|---|
%% apply the three beam solution
% Note, for instr data, this may be a bit redundant unless some data has been screened or bin-mapped out
if Config.coord_EX(4) == '1' % only 3-beam if the EX parameter is set, this replicates winADCP's behaviour
isnanBM = isnan(beamMatrix);
is3beam = squeeze(sum(isnanBM, 1)) == 1; % index of all bins over time that have 1 NaN only (3 beam)
transformMatrixRow4 = transformMatrix(4,:).';
for j = find(any(is3beam,2)).' % loop over the bin indices that have at least one 3 beam instance
for k = find(is3beam(j,:)) % loop over the times when the current bin has a 3 beam instance
isnanThis = isnan(beamMatrix(:,j,k));
temp = beamMatrix(:,j,k) .* transformMatrixRow4;
beamMatrix(isnanThis,j,k) = -sum(temp(~isnanThis)) / transformMatrix(4, isnanThis);
end
end
ADCP.processingComments(end+1,1) = {['Three beam solution applied to ' num2str(sum(is3beam(:))) ' individual bins (of ' num2str(numel(beamMatrix)) ' total). ']};
else
ADCP.processingComments(end+1,1) = {'Three beam solution not applied as the EX command was set to xxx0x. (This matches the behaviour of winADCP which responds to EX(4). We can override this on request, contact us. As a side note, winADCP does not respond to EX(5), the command for bin-mapping; winADCP always does bin-mapping on beam co-ord data.'};
end
%% convert bin-mapped, screened and 3-beamed data to XYZ
vTemp = reshape(beamMatrix, 4, n * p);
velocityInstrOrEarth = reshape(transformMatrix * vTemp, [m n p]); |
|
Overall Processing Flow:
From Instrument or Beam data:
...
The steps above with option available are tuneable by the user as detailed in the data product options below. The ADCP.processingComments field in both Nortek and RDI MAT file data products provides the user with a log of the decisions made in this processing flow. The MAT files (RDI only) also has a structure called ADCP.processingOptions that details what options the user choose and what options were actually applied, as some options, such as the error velocity screening threshold are not always applicable. See the data product options and MAT file format sections below for more information. This is complicated, contact us for help.
Data Verification
As noted earlier, the ultimate requirement for the ADCP data products is that they replicate the results of the manufacturer's software. We test against the manufacturers' software to verify the data at every software release (regression testing). The testing suite includes manual and automated testing where the output of the manufacturer's software is compared to the data products and where data products are compared to the captured files from the previous release. The test suite is regularly updated for any new scenarios. Improvements in the test procedure are also made as needed. In addition to testing the software, there has a been a significant effort to insure the metadata, heading in particular, is correct. These steps include careful measurement and documentation at deployment time, verifying and vetting ROV heading sensors, comparing data to nearby ADCPs, comparing results between various processing methods, and perhaps the best check is to compare the currents with the known/modelled tides. Contact us for information on the testing procedures and metadata verification.
...