Adcp updates

AutomaticQC/imosRingingBinSetQC.m

@@ -0,0 +1,57 @@
+function [sample_data, varChecked, paramsLog] = imosRingingBinSetQC(sample_data, ~)
+%function [sample_data, varChecked, paramsLog] = imosRingingBinSetQC(sample_data, ~)
+%
+% A ringing bin removal for ADCPs.
+%
+% The user enters the bin number(s) to be flagged bad in
+% imosRingingBinSetQC.txt.
+%
+% Every datapoint in the bin(s) is marked as bad.
+%
+% See imosRingingBinSetQC.txt to enter bins for flagging.
+%
+%
+%
+% author: Rebecca.Cowley@csiro.au
+%
+%
+narginchk(1, 2);
+varChecked = {};
+paramsLog = [];
+currentQCtest = mfilename;
+if ~isstruct(sample_data), error('sample_data must be a struct'); end
+
+%get the bins to be flagged
+options = IMOS.resolve.function_parameters();
+%probably a better way to do this.
+bins = options('bin');
+if ischar(bins)
+    bins = str2num(bins);
+end
+
+nt = numel(IMOS.get_data(sample_data.dimensions, 'TIME'));
+if IMOS.adcp.is_along_beam(sample_data)
+    bin_dist = IMOS.get_data(sample_data.dimensions, 'DIST_ALONG_BEAMS');
+else
+    bin_dist = IMOS.get_data(sample_data.dimensions, 'HEIGHT_ABOVE_SENSOR');
+end
+
+
+flag_vars = IMOS.adcp.current_variables(sample_data);
+qcSet = str2double(readProperty('toolbox.qc_set'));
+badFlag = imosQCFlag('bad', qcSet, 'flag');
+goodFlag = imosQCFlag('good', qcSet, 'flag');
+
+flags = ones(nt,length(bin_dist), 'int8') * goodFlag;
+flags(:,bins) = badFlag;
+flag_vars_inds = IMOS.find(sample_data.variables, flag_vars);
+for k = 1:numel(flag_vars_inds)
+    vind = flag_vars_inds(k);
+    sample_data.variables{vind}.flags = flags;
+end
+
+varChecked = flag_vars;
+
+paramsLog = ['ringing_bin_removed=' num2str(bins)];
+writeDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'ringing_bin_removed', bins);
+end

AutomaticQC/imosRingingBinSetQC.txt

@@ -0,0 +1,4 @@
+% list of bins to be flagged as bad, separated by a comma and a space
+% eg: bin = 1 or bin = 1, 2, 3
+bin= 2, 4
+

IMOS/imosSites.txt

@@ -68,6 +68,7 @@ SAM7DS, 	135.8439, 		-36.1809, 		0.025, 	0.025, 	2.5
 SAM8SG, 	136.69, 		-35.25, 		0.025, 	0.025, 	2.5
 TAN100, 	113.9065833, 	-21.8493, 		0.025, 	0.025, 	2.5
 EAC0500,	153.89,		-27.32,			0.05, 	0.05, 	2.5
+EAC1520,	153.96,		-27.31,			0.05, 	0.05, 	2.5
 EAC2000,	153.99,		-27.31,			0.05, 	0.05, 	2.5
 EAC3200,	154.13,		-27.28,			0.05, 	0.05, 	2.5
 EAC4200,	154.29,		-27.24,			0.05, 	0.05, 	2.5

Preprocessing/+TeledyneADCP/workhorse_3beamsolution.m

@@ -0,0 +1,50 @@
+function [b1,b2,b3,b4] = workhorse_3beamsolution(T, beamdat)
+%function btim = workhorse_3beamsolution(T, beamdat)
+%
+% Return the 3 beam solutions for instrument coordinate data with screened
+% data indicating one bad beam
+%
+% 
+% Reference: page 14-15 of ADCP Coordinate Transformation,
+% Formulas and Calculations, Teledyne RD Instruments.
+% P/N 951-6079-00 (January 2008).
+%
+%
+% Inputs:
+%
+%   T=transformation matrix from beam to instrument axis as obtained from
+%   workhorse_beam2inst
+%   beamdat: mx4 matrix with m depth cells, and 4 beams (one profile)
+%
+% Outputs:
+%
+% btim [double] - A 4x4 beam data matrix with 3-beam solutions included.
+%
+%
+% author: rebecca.cowley@csiro.au
+%
+%
+narginchk(2, 2)
+
+if nargin < 2
+    errormsg('Not enough arguments to continue')
+end
+
+%code snippet from any3beam.m from UWA code set.
+three_beam = find(isfinite(beamdat)*ones(4, 1) == 3);
+btim=beamdat;
+%disp(length(three_beam))
+if any(three_beam)
+	data = beamdat(three_beam, :); % CB selecting depth cells  with exactly 3 "good" beams 
+	mask = isnan(data); % CB finding which of the 4 beam is bad
+	data(mask) = 0; %CB forcing "bad beam"  to zero instead of NaN. Now the error vel is 0.
+	tran = ones(length(three_beam), 1) * T(4, :); % CB matrix with rows representing the depth cells with 3 beams, each row represents the error vel transformation
+	err = (data .* tran) * ones(4, 1);% mult the data with the error vel trans matrix
+	data(mask) = -err ./ tran(mask); % CB Assigning the bad beam to its value dependent on other 3 beams
+	btim(three_beam, :) = data;
+end
+b1 = btim(:,1)';
+b2 = btim(:,2)';
+b3 = btim(:,3)';
+b4 = btim(:,4)';
+end

Preprocessing/+TeledyneADCP/workhorse_beam2earth.m

@@ -84,6 +84,9 @@
 end
 
 for t = 1:size(v1, 1)
+    %get values for 3-beam solutions
+    [b1(t, :), b2(t, :), b3(t, :), b4(t, :)] = TeledyneADCP.workhorse_3beamsolution(I, [b1(t, :); b2(t, :); b3(t, :); b4(t, :)]');
+
     h_t = heading(t);
     p_t = gpitch(t);
     r_t = adjusted_roll(t);

Preprocessing/adcpWorkhorseCorrMagPP.m

@@ -0,0 +1,132 @@
+function sample_data = adcpWorkhorseCorrMagPP( sample_data, qcLevel, ~ )
+%ADCPWORKHORSECORRMAGPP Screening procedure for Teledyne Workhorse (and similar)
+% ADCP instrument data, using the correlation magnitude velocity diagnostic variable.
+%
+%
+% Inputs:
+%   sample_data - struct containing the entire data set and dimension data.
+%   auto - logical, run QC in batch mode
+%
+% Outputs:
+%   sample_data - same as input, with QC flags added for variable/dimension
+%                 data.
+% Author:       Guillaume Galibert <guillaume.galibert@utas.edu.au>
+%               Rebecca Cowley <rebecca.cowley@csiro.au>
+%
+
+%
+% Copyright (C) 2017, Australian Ocean Data Network (AODN) and Integrated
+% Marine Observing System (IMOS).
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation version 3 of the License.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+
+% You should have received a copy of the GNU General Public License
+% along with this program.
+% If not, see <https://www.gnu.org/licenses/gpl-3.0.en.html>.
+%
+narginchk(2, 3);
+if ~iscell(sample_data), error('sample_data must be a cell array'); end
+
+% auto logical in input to enable running under batch processing
+if nargin<2, auto=false; end
+
+% no modification of data is performed on the raw FV00 dataset except
+% local time to UTC conversion
+if strcmpi(qcLevel, 'raw'), return; end
+
+for k = 1:length(sample_data)
+    % do not process if not RDI nor Nortek
+    isRDI = false;
+    if strcmpi(sample_data{k}.meta.instrument_make, 'Teledyne RDI'), isRDI = true; end
+    if ~isRDI, continue; end
+    % get all necessary dimensions and variables id in sample_data struct
+    idVel1 = 0;
+    idVel2 = 0;
+    idVel3 = 0;
+    idVel4 = 0;
+    idCMAG = cell(4, 1);
+    for j=1:4
+        idCMAG{j}  = 0;
+    end
+    lenVar = length(sample_data{k}.variables);
+    for i=1:lenVar
+        paramName = sample_data{k}.variables{i}.name;
+
+        if strncmpi(paramName, 'VEL1', 4),  idVel1 = i; end
+        if strncmpi(paramName, 'VEL2', 4),  idVel2 = i; end
+        if strncmpi(paramName, 'VEL3',4),      idVel3 = i; end
+        if strncmpi(paramName, 'VEL4',4),      idVel4 = i; end
+        for j=1:4
+            cc = int2str(j);
+            if strcmpi(paramName, ['CMAG' cc]), idCMAG{j} = i; end
+        end
+    end
+
+    % check if the data is compatible with the QC algorithm
+    idMandatory = (idVel1 | idVel2 | idVel3 | idVel4 );
+    for j=1:4
+        idMandatory = idMandatory & idCMAG{j};
+    end
+    if ~idMandatory, return; end
+
+    % let's get the associated vertical dimension
+    idVertDim = sample_data{k}.variables{idCMAG{1}}.dimensions(2);
+    if strcmpi(sample_data{k}.dimensions{idVertDim}.name, 'DIST_ALONG_BEAMS')
+        disp(['Warning : imosCorrMagVelocitySetQC applied with a non tilt-corrected CMAGn (no bin mapping) on dataset ' sample_data{k}.toolbox_input_file]);
+    end
+
+    qcSet           = str2double(readProperty('toolbox.qc_set'));
+    badFlag         = imosQCFlag('bad',             qcSet, 'flag');
+    goodFlag        = imosQCFlag('good',            qcSet, 'flag');
+    rawFlag         = imosQCFlag('raw',             qcSet, 'flag');
+
+    %Pull out correlation magnitude
+    sizeData = size(sample_data{k}.variables{idCMAG{1}}.data);
+
+    % read in filter parameters
+    propFile = fullfile('Preprocessing', 'adcpWorkhorseCorrMagPP.txt');
+    cmag     = str2double(readProperty('cmagPP',   propFile));
+
+    % read dataset QC parameters if exist and override previous
+    % parameters file
+    currentQCtest = mfilename;
+    cmag = readDatasetParameter(sample_data{k}.toolbox_input_file, currentQCtest, 'cmag', cmag);
+
+    paramsLog = ['cmag=' num2str(cmag)];
+
+    sizeCur = size(sample_data{k}.variables{idVel1}.flags);
+
+    % same flags are given to any variable
+    for a = 1:4
+        flags = ones(sizeCur, 'int8')*rawFlag;
+
+        % Run QC. For this screening test, each beam is tested independently
+        eval(['iPass = sample_data{k}.variables{idCMAG{' num2str(a) '}}.data > cmag;'])
+        % Run QC filter (iFail) on velocity data
+        flags(~iPass) = badFlag;
+        flags(iPass) = goodFlag;
+
+        eval(['sample_data{k}.variables{idVel' num2str(a) '}.flags = flags;'])
+
+    end
+
+
+    % write/update dataset QC parameters
+    writeDatasetParameter(sample_data{k}.toolbox_input_file, currentQCtest, 'cmagPP', cmag);
+    cmagcomment = ['adcpWorkhorseCorrMagPP.m: Correlation Magnitude preprocessing screening applied to beam velocities. Threshold = ' num2str(cmag)];
+
+
+    if isfield(sample_data{k}, 'history') && ~isempty(sample_data{k}.history)
+        sample_data{k}.history = sprintf('%s\n%s - %s', sample_data{k}.history, datestr(now_utc, readProperty('exportNetCDF.dateFormat')), cmagcomment);
+    else
+        sample_data{k}.history = sprintf('%s - %s', datestr(now_utc, readProperty('exportNetCDF.dateFormat')), cmagcomment);
+    end
+end
+end

Preprocessing/adcpWorkhorseCorrMagPP.txt

@@ -0,0 +1 @@
+cmagPP      = 64