Merge pull request #732 from aodn/qc_tests_adjustments_part3

Qc tests adjustments part3

AutomaticQC/imosEchoIntensitySetQC.m

@@ -0,0 +1,184 @@
+function [sample_data, varChecked, paramsLog] = imosEchoIntensitySetQC(sample_data, ~)
+%function [sample_data, varChecked, paramsLog] = imosEchoIntensitySetQC(sample_data, ~)
+%
+% A Echo intensity test for ADCPs using a threshold.
+%
+% The detection is done by inspecting the ADCP bin profile
+% echo amplitude first-order gradient. Everything above the
+% `echo_intensity_threshold` option is marked as bad.
+%
+% Everything beyond the first threshold marking may also be
+% marked as bad, if the `propagate` option is true. This is useful
+% for surface detection.
+%
+% As well, if `bound_by_depth` or `bound_by_index` is True, than
+% bad markings are bounded by the `bound_value` option.
+% This is useful to filter the echo amplitude QC markings to only
+% further specific depths or bin indexes.
+% 
+% See imosEchoIntensitySetQC.txt options.
+%
+%
+% Example:
+% % see test_imosEchoIntensitySetQC.m
+%
+% author: [email protected] [refactored from multiple versions of imosSurfaceDetectionQC,imosEchoIntensityVelocitySetQC, and others].
+%
+%
+narginchk(1, 2);
+varChecked = {};
+paramsLog = [];
+currentQCtest = mfilename;
+if ~isstruct(sample_data), error('sample_data must be a struct'); end
+
+[valid, reason] = IMOS.validate_dataset(sample_data, currentQCtest);
+
+if ~valid
+    %TODO: we may need to include a global verbose flag to avoid pollution here.
+    for k = 1:numel(reason)
+        dispmsg('Skipping %s. Reason: \n', sample_data.toolbox_input_file, reason{k})
+    end
+
+    return
+end
+
+avail_variables = IMOS.get(sample_data.variables, 'name');
+absic_counter = sum(contains(avail_variables, 'ABSIC'));
+absic_vars = cell(1, absic_counter);
+
+for k = 1:absic_counter
+    absic_vars{k} = ['ABSIC' num2str(k)];
+end
+
+if isfield(sample_data, 'history') && ~contains(sample_data.history, 'adcpBinMappingPP.m')
+    dispmsg('%s is not using bin-mapped variables.', sample_data.toolbox_input_file)
+end
+
+options = IMOS.resolve.function_parameters();
+threshold = options('echo_intensity_threshold');
+propagate = options('propagate');
+bound_by_depth = options('bound_by_depth');
+bound_by_index = options('bound_by_index');
+bound_value = options('bound_value');
+
+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');
+    dims_tz = {'TIME', 'DIST_ALONG_BEAMS'};
+    nbeams = numel(absic_vars);
+else
+    bin_dist = IMOS.get_data(sample_data.dimensions, 'HEIGHT_ABOVE_SENSOR');
+    dims_tz = {'TIME', 'HEIGHT_ABOVE_SENSOR'};
+    nbeams = numel(absic_vars);
+end
+
+flag_vars = IMOS.variables_with_dimensions(sample_data.dimensions, sample_data.variables, dims_tz);
+
+switch nbeams
+    case 3
+        non_flag_vars = {'PERG1', 'PERG2', 'PERG3', 'CMAG1', 'CMAG2', 'CMAG3'};
+        flag_vars = setdiff(flag_vars, non_flag_vars);
+    case 4
+        non_flag_vars = {'PERG1', 'PERG2', 'PERG3', 'PERG4', 'CMAG1', 'CMAG2', 'CMAG3', 'CMAG4'};
+        flag_vars = setdiff(flag_vars, non_flag_vars);
+end
+
+% Now compute threshold and apply mask.
+% The order is:
+% 1.Compute threshold, 2. Bound detections by index or depth, 3. propagate detections further away in the beam.
+
+nz = numel(bin_dist);
+absic_extended = zeros(nt, nz + 1);
+beyond_threshold = zeros(nt, nz, 'logical');
+first_order_diff = 1;
+along_z = 2;
+
+for k = 1:absic_counter
+    varname = absic_vars{k};
+    absic = IMOS.get_data(sample_data.variables, varname);
+    absic_extended(:, 1) = absic(:, 1);
+    absic_extended(:, 2:end) = absic;
+    echodiff = abs(diff(absic_extended, first_order_diff, along_z));
+    beyond_threshold = beyond_threshold | echodiff > threshold;
+end
+
+if bound_by_index || bound_by_depth
+    do_index_bound = bound_by_index && bound_value > 0 && bound_value <= nz;
+    do_depth_bound = ~do_index_bound && bound_value > 0 && bound_by_depth && ~isempty(IMOS.get_data(sample_data.variables, 'DEPTH'));
+
+    invalid_options = ~do_index_bound && ~do_depth_bound;
+
+    if invalid_options
+        dispmsg('Bound option invalid in parameter file %s.', [mfilename '.txt']);
+    end
+
+else
+    do_index_bound = false;
+    do_depth_bound = false;
+end
+
+if do_index_bound
+    beyond_threshold(:,1:bound_value) = 0; %keep further bin as detected.
+elseif do_depth_bound
+    idepth = IMOS.get_data(sample_data.variables, 'DEPTH');
+    [t, z] = find(beyond_threshold);
+    bin_depths_at_invalid_points = idepth(t) - bin_dist(z); %TxZ arr
+    upward_looking = all(bin_dist > 0);
+
+    if upward_looking
+        bins_before_depth = find(bin_depths_at_invalid_points >= bound_value);
+    else
+        bins_before_depth = find(bin_depths_at_invalid_points <= bound_value);
+    end
+
+    rind = sub2ind(size(beyond_threshold), t(bins_before_depth), z(bins_before_depth));
+    beyond_threshold(rind) = 0;
+end
+
+if propagate
+    [t, z] = find(beyond_threshold);
+
+    if ~isempty(t)
+        beyond_threshold(t(1), z(1):end) = 1;
+        prev_t = t(1);
+
+        for k = 2:numel(t)
+
+            if t(k) == prev_t
+                continue
+            else
+                prev_t = t(k);
+                %because array is NTxNZ, z is ordered
+                beyond_threshold(t(k), z(k):end) = 1;
+            end
+
+        end
+
+    end
+
+end
+
+qcSet = str2double(readProperty('toolbox.qc_set'));
+badFlag = imosQCFlag('bad', qcSet, 'flag');
+goodFlag = imosQCFlag('good', qcSet, 'flag');
+
+flags = ones(size(beyond_threshold), 'int8') * goodFlag;
+flags(beyond_threshold) = 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;
+nbadbins = sum(beyond_threshold, 1);
+near_bad_bin = find(nbadbins, 1, 'first');
+further_bad_bin = find(nbadbins, 1, 'last');
+
+paramsLog = ['echo_intensity_threshold=' threshold, 'propagate=', propagate, 'near_bad_bin=' near_bad_bin, 'further_bad_bin=' further_bad_bin];
+
+writeDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'echo_intensity_threshold', threshold);
+
+end

AutomaticQC/imosEchoIntensitySetQC.txt

@@ -0,0 +1,10 @@
+% the first-order difference ABSIC count threshold
+echo_intensity_threshold = 50
+% a boolean to bound the threshold detection by depth. Turn this on to limit bad markings to bins further away from a specific depth (inclusive). For upward (downward) looking adcps, the Bins deeper (shallower) than a certain depth bound value will always be good.
+bound_by_depth = 0
+% a boolean to bound the threshold detection by index. Turn this on limit bad markings to bins further away from a bin index (inclusive). Bins with index value between 1:bound_value will always be good.
+bound_by_index = 0
+%bound value is positive in meters (bound_by_depth) or in bin index (bound_by_index).
+bound_value = 99999
+% propagate will expand the markings for bins further away. So if bin N from Nbins was marked, N:Nbin will be marked as bad.
+propagate = 0

AutomaticQC/imosEchoRangeSetQC.m

@@ -22,24 +22,6 @@
 % Author:       Guillaume Galibert <[email protected]>
 %               Rebecca Cowley <[email protected]>
 %
-
-%
-% 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(1, 2);
 if ~isstruct(sample_data), error('sample_data must be a struct'); end
 
@@ -110,7 +92,7 @@
 
 paramsLog = ['ea_fishthresh=' num2str(ea_fishthresh)];
 
-
+%TODO: refactor from below
 % Run QC
 % Following code is adapted from the UWA 'adcpfishdetection.m' code
 
@@ -122,7 +104,18 @@
 [B, Ix]=sort(ea,1); %sort echo from highest to lowest along each bin
     %step one - really only useful for data in beam coordinates. If one
     %beam fails, then can do 3-beam solutions
-if unique(sample_data.meta.fixedLeader.coordinateTransform) ~=7 %7 is the value for ENU setup. Not sure what it is for beam coords. Can be refined here.
+
+try
+    frame_of_reference = sample_data.meta.adcp_info.coords.frame_of_reference;
+    is_enu = strcmpi(frame_of_reference,'enu');
+catch
+    try
+        is_enu = unique(sample_data.meta.fixedLeader.coordinateTransform) ~=7;
+    catch
+        is_enu = false;
+    end
+end
+if is_enu
     df = B(4,:,:)-B(1,:,:); %get the difference from highest value to lowest
     ind=df>ea_fishthresh; % problematic depth cells
     bad_ea(ind) = true; %flag these as bad (individual cells)

AutomaticQC/imosSurfaceDetectionByDepthSetQC.m

@@ -0,0 +1,71 @@
+function [sample_data, varChecked, paramsLog] = imosSurfaceDetectionByDepthSetQC(sample_data, ~)
+%function [sample_data, varChecked, paramsLog] = imosSurfaceDetectionByDepthSetQC(sample_data, ~)
+%
+% A Surface detection test for ADCPs using Depth information.
+%
+% The detection is done by inspecting the ADCP bins and the
+% observable water column height. Any bin beyond this height
+% is marked as bad.
+%
+% This test only works for datasets with available TIME dimension, 
+% HEIGHT_ABOVE_SENSOR or DIST_ALONG_BEAMS dimensions, DEPTH variable,
+% and adcp_site_nominal_depth or site_depth_at_deployment metadata.
+%
+% author: [email protected] [refactored from older versions of imosSurfaceDetectSetQC].
+%
+narginchk(1, 2);
+varChecked = {};
+paramsLog = [];
+currentQCtest = mfilename;
+
+if ~isstruct(sample_data), error('sample_data must be a struct'); end
+
+[valid,reason] = IMOS.validate_dataset(sample_data,currentQCtest);
+if ~valid
+    %TODO: we may need to include a global verbose flag to avoid pollution here.
+    for k=1:numel(reason)
+        dispmsg('Skipping %s. Reason: \n', sample_data.toolbox_input_file, reason{k})
+    end
+    return
+end
+
+idepth = IMOS.get_data(sample_data.variables,'DEPTH');
+bathy = IMOS.meta.site_bathymetry(sample_data);
+beam_face_config = IMOS.meta.beam_face_config(sample_data);
+
+avail_dimensions = IMOS.get(sample_data.dimensions, 'name');
+along_height = inCell(avail_dimensions, 'HEIGHT_ABOVE_SENSOR');
+along_beams = inCell(avail_dimensions, 'DIST_ALONG_BEAMS');
+if along_height
+    bin_dist = IMOS.get_data(sample_data.dimensions,'HEIGHT_ABOVE_SENSOR');
+elseif along_beams
+    bin_dist = IMOS.get_data(sample_data.dimensions,'DIST_ALONG_BEAMS');
+end
+
+[wbins, last_water_bin] = IMOS.adcp.bin_in_water(idepth,bin_dist,beam_face_config,bathy);
+
+
+qcSet = str2double(readProperty('toolbox.qc_set'));
+badFlag = imosQCFlag('bad', qcSet, 'flag');
+goodFlag = imosQCFlag('good', qcSet, 'flag');
+
+flags = ones(size(wbins), 'int8') * badFlag;
+flags(wbins) = goodFlag;
+
+if along_beams
+    dims_tz = {'TIME','DIST_ALONG_BEAMS'};
+    dispmsg('No bin-mapping performed. Surface Detections will be contaminated by missing tilt corrections.')
+else
+    dims_tz = {'TIME','HEIGHT_ABOVE_SENSOR'};
+end
+
+vars_tz = IMOS.variables_with_dimensions(sample_data.dimensions,sample_data.variables,dims_tz);
+vars_tz_inds = IMOS.find(sample_data.variables,vars_tz);
+
+for k=1:numel(vars_tz_inds)
+    sample_data.variables{k}.flags = flags;
+end
+varChecked = vars_tz;
+paramsLog = ['min_surface_bin=' min(last_water_bin), 'max_surface_bin=' max(last_water_bin)];
+
+end