diff --git a/AutomaticQC/imosEchoIntensitySetQC.m b/AutomaticQC/imosEchoIntensitySetQC.m new file mode 100644 index 000000000..ae0234953 --- /dev/null +++ b/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: hugo.oliveira@utas.edu.au [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 diff --git a/AutomaticQC/imosEchoIntensitySetQC.txt b/AutomaticQC/imosEchoIntensitySetQC.txt new file mode 100644 index 000000000..de137fe1d --- /dev/null +++ b/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 diff --git a/AutomaticQC/imosEchoRangeSetQC.m b/AutomaticQC/imosEchoRangeSetQC.m index 4b86ea0cc..2c8cfab32 100644 --- a/AutomaticQC/imosEchoRangeSetQC.m +++ b/AutomaticQC/imosEchoRangeSetQC.m @@ -22,24 +22,6 @@ % Author: Guillaume Galibert % Rebecca Cowley % - -% -% 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 . -% 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) diff --git a/AutomaticQC/imosSurfaceDetectionByDepthSetQC.m b/AutomaticQC/imosSurfaceDetectionByDepthSetQC.m new file mode 100644 index 000000000..ffb071294 --- /dev/null +++ b/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: hugo.oliveira@utas.edu.au [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 diff --git a/AutomaticQC/imosSurfaceDetectionSetQC.m b/AutomaticQC/imosSurfaceDetectionSetQC.m deleted file mode 100644 index 0d186778c..000000000 --- a/AutomaticQC/imosSurfaceDetectionSetQC.m +++ /dev/null @@ -1,270 +0,0 @@ -function [sample_data, varChecked, paramsLog] = imosSurfaceDetectionSetQC( sample_data, auto ) -%IMOSSURFACEDETECTIONSETQC Quality control procedure for Teledyne Workhorse (and similar) -% ADCP instrument data, using the echo intensity velocity diagnostic variable, -% in conjunction with a side-lobe test and enhancements. -% -% Echo Amplitude test : -% this test looks at the difference between consecutive vertical bin values of ea and -% if the value exceeds the threshold, then the bin fails, and all bins -% above/below it are also considered to have failed. This test is designed to get -% rid of above/below surface/bottom bins. -% Enhancements to this test (originally imosEchoIntensityVelocitySetQC): -% 1. only bins within 3 * bin size of the surface/bottom are examined -% 2. NaN values in bins are counted as a pass -% 3. Only one beam has to fail the threshold test (was 3) -% 4. The imosSideLobeVelocitySetQC test is applied for bins that still -% have 'good' data that is above the surface + 1 bin (below the bottom - 1 bin) -% -% Assumptions: -% 1. The echo data has been bin mapped (but it will work without the bin -% mapping) -% 2. The depth is known for each time stamp -% -% 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. -% varChecked - cell array of variables' name which have been checked -% paramsLog - string containing details about params' procedure to include in QC log -% -% Author: Guillaume Galibert -% - -% -% 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 . -% -narginchk(1, 2); -if ~isstruct(sample_data), error('sample_data must be a struct'); end - -% auto logical in input to enable running under batch processing -if nargin<2, auto=false; end - -varChecked = {}; -paramsLog = []; - -% get all necessary dimensions and variables id in sample_data struct -idUcur = 0; -idVcur = 0; -idWcur = 0; -idCspd = 0; -idCdir = 0; -idABSIC = cell(4, 1); -for j=1:4 - idABSIC{j} = 0; -end -lenVar = length(sample_data.variables); -for i=1:lenVar - paramName = sample_data.variables{i}.name; - - if strncmpi(paramName, 'UCUR', 4), idUcur = i; end - if strncmpi(paramName, 'VCUR', 4), idVcur = i; end - if strcmpi(paramName, 'WCUR'), idWcur = i; end - if strcmpi(paramName, 'CSPD'), idCspd = i; end - if strncmpi(paramName, 'CDIR', 4), idCdir = i; end - for j=1:4 - cc = int2str(j); - if strcmpi(paramName, ['ABSIC' cc]), idABSIC{j} = i; end - end -end -% check if the data is compatible with the QC algorithm, otherwise quit -% silently -idHeight = getVar(sample_data.dimensions, 'HEIGHT_ABOVE_SENSOR'); -% check if the data is compatible with the QC algorithm -idMandatory = idHeight; -if ~idMandatory, return; end - -Bins = sample_data.dimensions{idHeight}.data'; -idDepth = getVar(sample_data.variables, 'DEPTH'); -depth = sample_data.variables{idDepth}.data; - -% let's get the associated vertical dimension -idVertDim = sample_data.variables{idABSIC{1}}.dimensions(2); -if strcmpi(sample_data.dimensions{idVertDim}.name, 'DIST_ALONG_BEAMS') - disp(['Warning : imosEchoIntensityVelocitySetQC applied with a non tilt-corrected ABSICn (no bin mapping) on dataset ' sample_data.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 echo intensity -sizeData = size(sample_data.variables{idABSIC{1}}.data); -ea = nan(4, sizeData(1), sizeData(2)); -for j=1:4 - ea(j, :, :) = sample_data.variables{idABSIC{j}}.data; -end - -% read in filter parameters -propFile = fullfile('AutomaticQC', 'imosSurfaceDetectionSetQC.txt'); -ea_thresh = str2double(readProperty('ea_thresh', propFile)); -nBinSize = str2double(readProperty('nBinSize', propFile)); - -% read dataset QC parameters if exist and override previous -% parameters file -currentQCtest = mfilename; -ea_thresh = readDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'ea_thresh', ea_thresh); -nBinSize = readDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'nBinSize', nBinSize); - -paramsLog = ['ea_thresh=' num2str(ea_thresh) ', nBinSize=' num2str(nBinSize)]; - -sizeCur = size(sample_data.variables{idUcur}.flags); - -% same flags are given to any variable -flags = ones(sizeCur, 'int8')*rawFlag; - -isUpwardLooking = true; -if all(Bins <= 0), isUpwardLooking = false; end -% we handle the case of a downward looking ADCP -if ~isUpwardLooking - if isempty(sample_data.site_nominal_depth) && isempty(sample_data.site_depth_at_deployment) - error(['Downward looking ADCP in file ' sample_data.toolbox_input_file ' => Fill site_nominal_depth or site_depth_at_deployment!']); - else - % the distance between transducer and obstacle is not depth anymore but - % (site_nominal_depth - depth) - if ~isempty(sample_data.site_nominal_depth) - site_nominal_depth = sample_data.site_nominal_depth; - end - if ~isempty(sample_data.site_depth_at_deployment) - site_nominal_depth = sample_data.site_depth_at_deployment; - end - end -end - -% Run QC -% this test looks at the difference between consecutive vertical bin values of ea and -% if the value exceeds the threshold, then the bin fails, and all bins -% above it are also considered to have failed. This test is designed to get -% rid of above/below surface/bottom bins. -lenTime = sizeCur(1); -lenBin = sizeCur(2); -%Let's refine this to only look at bins within +/- 3 bins of the surface/bottom -% and profiles where the last bin reaches the surface/bottom -% currently assumes that there is a depth calcuated. -binDepth = depth - repmat(Bins,lenTime,1); -binRange = Bins(2) - Bins(1); -if isUpwardLooking - isurface = binDepth <= 5*binRange ; -else - isurface = binDepth >= site_nominal_depth - 5*binRange; -end - -% if the following test is successfull, the bin gets good -ib = uint8(abs(diff(squeeze(ea(1, :,:)),1,2)) <= ea_thresh) + ... - uint8(abs(diff(squeeze(ea(2, :,:)),1,2)) <= ea_thresh) + ... - uint8(abs(diff(squeeze(ea(3, :,:)),1,2)) <= ea_thresh) + ... - uint8(abs(diff(squeeze(ea(4, :,:)),1,2)) <= ea_thresh); - - -% allow for NaNs in the ea - false fails where more than 1 beam has NaN -inan = isnan(abs(diff(squeeze(ea(1, :,:)),1,2))) + ... - isnan(abs(diff(squeeze(ea(2, :,:)),1,2))) + ... - isnan(abs(diff(squeeze(ea(3, :,:)),1,2))) + ... - isnan(abs(diff(squeeze(ea(4, :,:)),1,2))); - -% we look for the bins where all bins pass the tests, as we know the depth -% and we know where the surface is, so any beam that sees is indicates -% where bad data is. -ib = ib + uint8(inan) == 4; %these are 'good', no beam sees the surface/bottom -ib = [true(lenTime, 1), ib]; %add one row at the start to account for diff -ib(~isurface) = 4; % bins not within 5*binRange of surface or bottom do not fail - -% any good bin further than a bad one should be set bad -jkf = repmat(single(1:1:lenBin), [lenTime, 1]); - -iii = single(~ib).*jkf; -iii(iii == 0) = NaN; -iif = min(iii, [], 2); -iifNotNan = ~isnan(iif); - -iPass = true(lenTime, lenBin); -if any(iifNotNan) - % all bins further than the first bad one is reset to bad - iPass(jkf >= repmat(iif, [1, lenBin])) = false; -end -clear iifNotNan iif jkf ib iii - -% finally, tidy up bins that still have good data above the surface with -% the side-lobe test: -% calculate contaminated depth -% -% http://www.nortekusa.com/usa/knowledge-center/table-of-contents/doppler-velocity#Sidelobes -% -% by default substraction of 1/2*BinSize to the non-contaminated height in order to be -% conservative and be sure that the first bin below the contaminated depth -% hasn't been computed from any contaminated signal. -% read in filter parameters -BinSize = sample_data.meta.binSize; -ipass = single(iPass).*binDepth; %bin depths that pass - -if isUpwardLooking - distanceTransducerToObstacle = depth; - cDepth = distanceTransducerToObstacle - (distanceTransducerToObstacle * cos(sample_data.meta.beam_angle*pi/180) - nBinSize*BinSize); -else - distanceTransducerToObstacle=site_nominal_depth - depth; - cDepth = site_nominal_depth - (distanceTransducerToObstacle - (distanceTransducerToObstacle * cos(sample_data.meta.beam_angle*pi/180) - nBinSize*BinSize)); -end -cd = repmat(cDepth, [1, length(Bins)]); - -% test bins depths against contaminated depth -% where ibad, apply the cDepth cutoff: -if isUpwardLooking - % upward looking : all bins above the contaminated depth are flagged, - % except where the maximum bin depth is > 0 - ipass(ipass < cd & ipass <= 0 + binRange) = NaN; -else - % downward looking : all bins below the contaminated depth are flagged, - % except where the maximum bin depth is < bottom - ipass(ipass > cd & ipass >= site_nominal_depth - binRange) = NaN; -end - -%ipass is a matrix of depths that passes -%fix the pass/fail -ipass(ipass==0) = NaN; -iPass = ~isnan(ipass); -iFail = ~iPass; - -% Run QC filter (iFail) on velocity data -flags(iFail) = badFlag; -flags(iPass) = goodFlag; - -sample_data.variables{idUcur}.flags = flags; -sample_data.variables{idVcur}.flags = flags; -sample_data.variables{idWcur}.flags = flags; - -varChecked = {sample_data.variables{idUcur}.name, ... - sample_data.variables{idVcur}.name, ... - sample_data.variables{idWcur}.name}; - -if idCdir - sample_data.variables{idCdir}.flags = flags; - varChecked = [varChecked, {sample_data.variables{idCdir}.name}]; -end - -if idCspd - sample_data.variables{idCspd}.flags = flags; - varChecked = [varChecked, {sample_data.variables{idCspd}.name}]; -end - -% write/update dataset QC parameters -writeDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'ea_thresh', ea_thresh); -writeDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'nBinSize', nBinSize); - -end diff --git a/Geomag/sample_coords.txt b/Geomag/sample_coords.txt deleted file mode 100644 index d5dd73e72..000000000 --- a/Geomag/sample_coords.txt +++ /dev/null @@ -1,6 +0,0 @@ -2018,12,17 D M-108.000000 -27.284100 154.130100 -2019,01,02 D M-108.000000 -27.284100 154.130100 -2019,01,02 D M-638.000000 -27.284100 154.130100 -2019,01,02 D M-1000.000000 -27.284100 154.130100 -2019,01,02 D M-1000.000000 -27.284100 154.130100 -2019,01,02 D M-1000.000000 -27.284100 154.130100 diff --git a/Geomag/sample_out_IGRF12.txt b/Geomag/sample_out_IGRF12.txt deleted file mode 100644 index 197b45d6c..000000000 --- a/Geomag/sample_out_IGRF12.txt +++ /dev/null @@ -1,7 +0,0 @@ -Date Coord-System Altitude Latitude Longitude D_deg D_min I_deg I_min H_nT X_nT Y_nT Z_nT F_nT dD_min dI_min dH_nT dX_nT dY_nT dZ_nT dF_nT -2018,12,17 D M-108.000000 -27.284100 154.130100 11d 15m -57d 10m 28649.3 28099.0 5588.2 -44389.4 52831.8 0.6 -1.0 -18.2 -18.8 1.6 -0.7 -9.2 -2019,01,02 D M-108.000000 -27.284100 154.130100 11d 15m -57d 10m 28648.5 28098.1 5588.3 -44389.4 52831.4 0.6 -1.0 -18.2 -18.8 1.6 -0.7 -9.2 -2019,01,02 D M-638.000000 -27.284100 154.130100 11d 15m -57d 10m 28656.0 28105.5 5590.0 -44401.4 52845.6 0.6 -1.0 -18.2 -18.8 1.6 -0.7 -9.2 -2019,01,02 D M-1000.000000 -27.284100 154.130100 11d 15m -57d 10m 28661.2 28110.6 5591.1 -44409.6 52855.3 0.6 -1.0 -18.2 -18.8 1.6 -0.7 -9.2 -2019,01,02 D M-1000.000000 -27.284100 154.130100 11d 15m -57d 10m 28661.2 28110.6 5591.1 -44409.6 52855.3 0.6 -1.0 -18.2 -18.8 1.6 -0.7 -9.2 -2019,01,02 D M-1000.000000 -27.284100 154.130100 11d 15m -57d 10m 28661.2 28110.6 5591.1 -44409.6 52855.3 0.6 -1.0 -18.2 -18.8 1.6 -0.7 -9.2 diff --git a/Parser/readWorkhorseEnsembles.m b/Parser/readWorkhorseEnsembles.m index 153461080..54891aa53 100644 --- a/Parser/readWorkhorseEnsembles.m +++ b/Parser/readWorkhorseEnsembles.m @@ -71,10 +71,11 @@ % read in the whole file into 'data' fid = -1; data = []; -% try -% m = memmapfile(filename,'Format','uint8','Writable',false); -% data = m.data; -% catch e +try + %reduce a lot of peak memory consumption. + m = memmapfile(filename,'Format','uint8','Writable',false); + data = m.data; +catch try fid = fopen(filename, 'rb'); data = fread(fid, inf, '*uint8'); @@ -83,7 +84,7 @@ if fid ~= -1, fclose(fid); end rethrow e; end -% end +end % get ensemble start key headerID = 127; % hexadecimal '7F' @@ -362,7 +363,7 @@ % - - - - - 0 0 1 150-kHz SYSTEM % - - - - - 0 1 0 300-kHz SYSTEM % - - - - - 0 1 1 600-kHz SYSTEM -% - - - -data - 1 0 0 1200-kHz SYSTEM +% - - - - - 1 0 0 1200-kHz SYSTEM % - - - - - 1 0 1 2400-kHz SYSTEM % - - - - 0 - - - CONCAVE BEAM PAT. % - - - - 1 - - - CONVEX BEAM PAT. @@ -382,8 +383,17 @@ % 0 1 0 0 - - - - 4-BEAM JANUS CONFIG % 0 1 0 1 - - - - 5-BM JANUS CFIG DEMOD) % 1 1 1 1 - - - - 5-BM JANUS CFIG.(2 DEMD) - LSB = logical(de2bi(data(idx+4),8,'left-msb')); - MSB = logical(de2bi(data(idx+5),8,'left-msb')); +% +LSB = dec2bin(double(data(idx+4))); +MSB = dec2bin(double(data(idx+5))); +while(size(LSB, 2) < 8), LSB = [repmat('0', size(LSB, 1), 1) LSB]; end +while(size(MSB, 2) < 8), MSB = [repmat('0', size(MSB, 1), 1) MSB]; end +LSB = bin2logical(LSB); +MSB = bin2logical(MSB); +% compatibility layer, code above is equivalent to: +%LSB = logical(de2bi(data(idx+4),8,'left-msb')); +%MSB = logical(de2bi(data(idx+5),8,'left-msb')); +% sect.systemConfiguration = [LSB MSB]; sect.realSimFlag = double(data(idx+6)); sect.lagLength = double(data(idx+7)); @@ -409,12 +419,29 @@ % xxxxx1xx = TILTS (PITCH AND ROLL) USED IN SHIP OR EARTH TRANSFORMATION % xxxxxx1x = 3-BEAM SOLUTION USED IF ONE BEAM IS BELOW THE CORRELATION THRESHOLD SET BY THE WC command % xxxxxxx1 = BIN MAPPING USED - sect.coordinateTransform = logical(de2bi(data(idx+25),8,'left-msb')); + % + ct = dec2bin(double(data(idx+25))); + while(size(ct, 2) < 8), ct = [repmat('0', size(ct,1), 1) ct]; end + ct = bin2logical(ct); + sect.coordinateTransform = ct; + %compatibility layer, code above is equivalent to: + %sect.coordinateTransform = logical(de2bi(data(idx+25),8,'left-msb')); block = indexData(data, idx+26, idx+29, 'int16', cpuEndianness)'; sect.headingAlignment = block(:,1); sect.headingBias = block(:,2); - sect.sensorSource = logical(de2bi(data(idx+30),8,'left-msb')); - sect.sensorsAvailable = logical(de2bi(data(idx+31),8,'left-msb')); + % + ss = dec2bin(double(data(idx+30))); + while(size(ss, 2) < 8), ss = [repmat('0', size(ss,1), 1) ss]; end + ss = bin2logical(ss); + sect.sensorSource = ss; + %compatibility layer, code above is equivalent to: + %sect.sensorSource = logical(de2bi(data(idx+30),8,'left-msb')); + sa = dec2bin(double(data(idx+31))); + while(size(sa, 2) < 8), sa = [repmat('0', size(sa,1), 1) sa]; end + sa = bin2logical(sa); + sect.sensorsAvailable = sa; + %compatibility layer, code above is equivalent to: + %sect.sensorsAvailable = logical(de2bi(data(idx+31),8,'left-msb')); block = indexData(data, idx+32, idx+35, 'uint16', cpuEndianness)'; sect.bin1Distance = block(:,1); sect.xmitPulseLength = block(:,2); diff --git a/Preprocessing/adcpBinMappingPP.m b/Preprocessing/adcpBinMappingPP.m index 3d62eb9ba..a41cac2a6 100644 --- a/Preprocessing/adcpBinMappingPP.m +++ b/Preprocessing/adcpBinMappingPP.m @@ -82,6 +82,11 @@ % closer to the surface while beam 2 gets further away. % 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; + 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; @@ -125,8 +130,6 @@ %TODO: remove nested logics to outer loop, by pre-computing the vars/ % interpolation conditions. This would allow a progress bar here. - Tlen = size(IMOS.get_data(sample_data{k}.dimensions, 'TIME'), 1); - for n = 1:number_of_beams beam_vars = get_beam_vars(sample_data{k}, n); [all_beam_vars] = IMOS.concatenate_variables(sample_data{k}.variables, beam_vars); diff --git a/Preprocessing/adcpWorkhorseVelocityBeam2EnuPP.m b/Preprocessing/adcpWorkhorseVelocityBeam2EnuPP.m index 0e24c7c39..c642fc77d 100644 --- a/Preprocessing/adcpWorkhorseVelocityBeam2EnuPP.m +++ b/Preprocessing/adcpWorkhorseVelocityBeam2EnuPP.m @@ -18,7 +18,7 @@ % sample_data - the same data sets, with updated velocity variables in ENU coordinates. % % -% author: hugo.oliveira@utas.edu.au +% author: hugo.oliveira@utas.edu.au [based on several draft and BecCowley 675 branch] % narginchk(2, 3) @@ -60,7 +60,7 @@ if info.coords.used_binmapping %compat: The toolbox is silent about this. - %warnmsg('Instrument frame of reference is beam but bin_mapping was selected at instrument settings. Performing rotation to earth coordinates anyway...') + %dispmsg('Instrument frame of reference is beam but bin_mapping was selected at instrument settings. Performing rotation to earth coordinates anyway...') % Even if the EX bit is set for bin mapping, % if the data is collected in beam coordinates, % bin mapping does not occur on board. diff --git a/Util/+IMOS/+adcp/beam_height.m b/Util/+IMOS/+adcp/beam_height.m new file mode 100644 index 000000000..81ea808ab --- /dev/null +++ b/Util/+IMOS/+adcp/beam_height.m @@ -0,0 +1,32 @@ +function [bheight] = beam_height(height_above_sensor) +% function [bheight] = beam_height(height_above_sensor) +% +% Compute the beam height from an ADCP. +% +% Inputs: +% +% height_above_sensor[numeric] - the bin centre distance from +% the transducer face. +% +% Outputs: +% +% bheight - the total beam height +% +% Example: +% +% %basic usage +% x = IMOS.adcp.beam_height([-10,-20,-30]); +% assert(x==35) +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(1,1) + +if isempty(height_above_sensor) + bheight = []; + return +end + +bheight = abs(height_above_sensor(end)) + 0.5 .* ( abs(height_above_sensor(end)) - abs(height_above_sensor(end-1))); + +end diff --git a/Util/+IMOS/+adcp/bin_in_water.m b/Util/+IMOS/+adcp/bin_in_water.m new file mode 100644 index 000000000..828a6f458 --- /dev/null +++ b/Util/+IMOS/+adcp/bin_in_water.m @@ -0,0 +1,115 @@ +function [wb, fbin] = bin_in_water(idepth,bin_dist,beam_face_config,bathy) +% function [wb,fbin] = bin_in_water(idepth,bin_dist,beam_face_config,bathy) +% +% Find the ADCP water bins and the +% further bin that is within the water envelope: +% +% Dowward out > x +% Looking out > x +% ------------- --------------zeta(t,x,y) +% Fbin +% vvvv <-- idepth wb +% wb wb +% wb wb +% wb ^^^^ <-- idepth +% Fbin +% ------------- ------------- h(x,y) +% out > x Upward +% out > x Looking +% +% Inputs: +% +% idepth [double] [1x1 or Nx1] - the instrument depth. +% bin_dist [double] [Mx1] - the adcp bin centre distances. +% beam_face_config [str] - the adcp face config ["up" or "down"] +% bathy [double] [1x1] - the local bathymetry (required for "down" measurements). +% +% Outputs: +% +% wb [logical] TxZ - the bins in the water envelope. +% fbin [array] 1xT - the last/further bin inside the envelope. +% +% Example: +% +% %basic usage - no bin in water +% idepth = 190*ones(100,1); +% beam_face_config = 'down'; +% bin_dist = -1*(10:10:500)'; +% bathy = 200; +% [wb,fbin] = IMOS.adcp.bin_in_water(idepth,bin_dist,beam_face_config,bathy); +% assert(~any(wb,'all')); +% assert(all(fbin==0)); +% +% % one bin in water +% bathy = 210; +% [wb,fbin] = IMOS.adcp.bin_in_water(idepth,bin_dist,beam_face_config,bathy); +% assert(sum(wb,'all')==100); +% assert(all(fbin==1)); +% +% % 10 bins in water +% bathy = 300; +% [wb,fbin] = IMOS.adcp.bin_in_water(idepth,bin_dist,beam_face_config,bathy); +% assert(sum(wb,'all')==1000); % 100*10 +% assert(all(fbin==10)); +% +% %upward +% beam_face_config = 'up'; +% idepth = 10*ones(100,1); +% [wb,fbin] = IMOS.adcp.bin_in_water(idepth,bin_dist,beam_face_config); +% assert(sum(wb,'all')==0); +% assert(all(fbin==0)); +% +% idepth = 20*ones(100,1); +% [wb,fbin] = IMOS.adcp.bin_in_water(idepth,bin_dist,beam_face_config); +% assert(sum(wb,'all')==100); +% assert(all(fbin==1)); +% +% idepth = 110*ones(100,1); +% [wb,fbin] = IMOS.adcp.bin_in_water(idepth,bin_dist,beam_face_config); +% assert(sum(wb,'all')==1000); % 100*10 +% assert(all(fbin==10)); +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(3,4); + +if ~iscolumn(idepth) || isempty(idepth) + errormsg('First argument is invalid. Instrument depth should be a non-empty column vector') +end +if ~iscolumn(bin_dist) + errormsg('Second argument is invalid. ADCP bin distances should be a non-empty column vector') +end + +upward_looking = strcmpi(beam_face_config,'up'); +downward_looking = strcmpi(beam_face_config,'down'); +if ~upward_looking && ~downward_looking + errormsg('Third Argument is invalid. ADCP beam face config should be ''up'' or ''down''') +end +if downward_looking && nargin<4 + errormsg('Not enough arguments. Down-looking ADCPs requires an extra argument representing the local bathymetry.') +end + +if nargin>3 && (isempty(bathy) || numel(bathy)~=1) + errormsg('Fourth argument is invalid. Bathymetry should be a non-empty scalar.') +end + +nz = numel(bin_dist); +bin_distance = reshape(abs(bin_dist), 1, nz); + +%follow convention that depth is positive down. +if upward_looking + visible_water_column_height = idepth; % + zeta +else + if isempty(bathy) + errormsg('No site_depth in for %s. Cannot estimate visible water column height.', sample_data.toolbox_input_file) + end + visible_water_column_height = bathy - idepth; % + zeta +end +wb = bin_distance < visible_water_column_height; + +if nargout > 1 + fbin = sum(wb, 2); +end + +end diff --git a/Util/+IMOS/+adcp/contains_adcp_dimensions.m b/Util/+IMOS/+adcp/contains_adcp_dimensions.m new file mode 100644 index 000000000..e340d7a60 --- /dev/null +++ b/Util/+IMOS/+adcp/contains_adcp_dimensions.m @@ -0,0 +1,48 @@ +function [bool,dname] = contains_adcp_dimensions(sample_data) +% function [bool,dname] = contains_adcp_dimensions(sample_data) +% +% Detect if the sample data contains adcp dimensions. +% +% Inputs: +% +% sample_data - toolbox struct. +% +% Outputs: +% +% bool - True or False for ADCP data. +% dname - the dimension associated with the transducer bins. +% +% Example: +% +% %basic usage +% z.dimensions = IMOS.gen_dimensions('adcp'); +% [b,name] = IMOS.adcp.contains_adcp_dimensions(z); +% assert(b); +% assert(strcmpi(name,'DIST_ALONG_BEAMS')); +% z.dimensions{3}.name = 'HEIGHT_ABOVE_SENSOR'; +% [b,name] = IMOS.adcp.contains_adcp_dimensions(z); +% assert(b); +% %HEIGHT is checked first. +% assert(strcmpi(name,'HEIGHT_ABOVE_SENSOR')); +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(1, 1) +bool = false; +dname = ''; + +avail_dimensions = IMOS.get(sample_data.dimensions, 'name'); + +height = 'HEIGHT_ABOVE_SENSOR'; +beam = 'DIST_ALONG_BEAMS'; + +if inCell(avail_dimensions,height) + dname = height; + bool = true; +elseif inCell(avail_dimensions,beam) + dname = beam; + bool = true; +end + +end diff --git a/Util/+IMOS/+adcp/is_along_beam.m b/Util/+IMOS/+adcp/is_along_beam.m new file mode 100644 index 000000000..2b1dd7bf9 --- /dev/null +++ b/Util/+IMOS/+adcp/is_along_beam.m @@ -0,0 +1,29 @@ +function [bool] = is_along_beam(sample_data) +% function [bool] = is_along_beam(sample_data) +% +% Check if a dataset is defined along beam coordinates. +% +% Inputs: +% +% sample_data [struct] - a toolbox dataset. +% +% Outputs: +% +% bool - True if 'DIST_ALONG_BEAMS' is defined. +% +% Example: +% +% %basic usage +% x.dimensions = IMOS.gen_dimensions('adcp'); +% assert(IMOS.adcp.is_along_beam(x)); +% x.dimensions{3}.name = 'HEIGHT_ABOVE_SENSOR'; +% assert(IMOS.adcp.is_along_beam(x)); +% x.dimensions{2}.name = 'HEIGHT_ABOVE_SENSOR'; +% assert(~IMOS.adcp.is_along_beam(x)); +% +% +% author: hugo.oliveira@utas.edu.au +% +avail_dims = IMOS.get(sample_data.dimensions,'name'); +bool = inCell(avail_dims,'DIST_ALONG_BEAMS'); +end diff --git a/Util/+IMOS/+adcp/number_of_beams.m b/Util/+IMOS/+adcp/number_of_beams.m new file mode 100644 index 000000000..b9aee48fa --- /dev/null +++ b/Util/+IMOS/+adcp/number_of_beams.m @@ -0,0 +1,57 @@ +function [number_of_beams] = number_of_beams(sample_data) +% function [n] = number_of_beams(sample_data) +% +% Detect the number of beams in a toolbox +% struct. +% +% Inputs: +% +% sample_data [struct] - the toolbox struct. +% +% Outputs: +% +% number_of_beams [numeric] - ditto. +% +% Example: +% +% %basic usage +% s.meta.adcp_info.number_of_beams = 4; +% assert(IMOS.adcp.number_of_beams(s)==4); +% x.meta.nBeams = 3; +% assert(IMOS.adcp.number_of_beams(x)==3); +% y.variables = {struct('name','WCUR_2')}; +% assert(IMOS.adcp.number_of_beams(y)==4); +% z.variables = {struct('name','ECUR')}; +% assert(IMOS.adcp.number_of_beams(y)==4); +% q.variables = {struct('name','WCUR')}; +% assert(IMOS.adcp.number_of_beams(q)==3); +% q = struct(); +% assert(IMOS.adcp.number_of_beams(q)==0); +% +% +% author: hugo.oliveira@utas.edu.au +% +try + number_of_beams = sample_data.meta.adcp_info.number_of_beams; +catch + + try + number_of_beams = sample_data.meta.nBeams; + catch + + try + vars = IMOS.get(sample_data.variables, 'name'); + + if inCell(vars, 'WCUR_2') || inCell(vars, 'ECUR') + number_of_beams = 4; + else + number_of_beams = 3; + end + + catch + number_of_beams = 0; + end + + end + +end diff --git a/Util/+IMOS/+meta/beam_face_config.m b/Util/+IMOS/+meta/beam_face_config.m new file mode 100644 index 000000000..99b658085 --- /dev/null +++ b/Util/+IMOS/+meta/beam_face_config.m @@ -0,0 +1,44 @@ +function [beam_face_config] = beam_face_config(sample_data) +% function [beam_face_config] = beam_face_config(sample_data) +% +% Get the beam face config of an ADCP. +% +% Inputs: +% +% sample_data [struct] - the toolbox dataset. +% +% Outputs: +% +% beam_face_config [char] - 'up','down' or empty. +% +% Example: +% +% %basic usage +% +% +% author: hugo.oliveira@utas.edu.au +% +try + beam_face_config = sample_data.meta.adcp_info.beam_face_config; + return +catch + try + bin_dist = IMOS.get_data(sample_data.dimensions,'HEIGHT_ABOVE_SENSOR'); + catch + try + bin_dist = IMOS.get_data(sample_data.dimensions,'DIST_ALONG_BEAMS'); + catch + beam_face_config = ''; + return + end + end + if all(bin_dist>0) + beam_face_config = 'up'; + elseif all(bin_dist<0) + beam_face_config = 'down'; + else + beam_face_config = ''; + end +end + +end diff --git a/Util/+IMOS/+meta/site_bathymetry.m b/Util/+IMOS/+meta/site_bathymetry.m new file mode 100644 index 000000000..74a982fbf --- /dev/null +++ b/Util/+IMOS/+meta/site_bathymetry.m @@ -0,0 +1,68 @@ +function [depth] = site_bathymetry(sample_data, oper) +%function [depth] = site_bathymetry(sample_data, oper) +% +% Disambiguate the site bathymetry from a sample +% +% Inputs: +% +% sample_data [struct] - the toolbox struct +% oper [function handle] - a function handle to operate on +% the different estimates if any. +% Default: @min +% +% Outputs: +% +% depth [numeric] - the bathymetry value +% +% +% Example: +% sample_data.site_nominal_depth = []; +% sample_data.site_depth_at_deployment = []; +% assert(isempty(IMOS.meta.get_site_bathymetry(sample_data))) +% sample_data.site_nominal_depth = 20; +% assert(IMOS.meta.get_site_bathymetry(sample_data)==20) +% sample_data.site_depth_at_deployment = 10; +% assert(IMOS.meta.get_site_bathymetry(sample_data)==10) +% +% author: hugo.oliveira@utas.edu.au +% +% +narginchk(1, 2) + +if ~isstruct(sample_data) + errormsg('Not a struct') +end + +if nargin == 1 + oper = @min; +else + + if ~isfunctionhandle(oper) + errormsg('Second argument not a function handle') + end + +end + +try + b1 = sample_data.site_depth_at_deployment; +catch + b1 = []; +end + +try + b2 = sample_data.site_nominal_depth; +catch + b2 = []; +end + +if isempty(b1) + b1 = b2; +end + +if isempty(b2) + b2 = b1; +end + +depth = oper(b1, b2); + +end diff --git a/Util/+IMOS/+meta/velocity_variables.m b/Util/+IMOS/+meta/velocity_variables.m new file mode 100644 index 000000000..dec5a94f2 --- /dev/null +++ b/Util/+IMOS/+meta/velocity_variables.m @@ -0,0 +1,53 @@ +function [vnames] = velocity_variables(sample_data) +% function [vnames] = velocity_variables(sample_data) +% +% Get the velocity variable names (UCUR,VCUR,and variants) +% from a dataset. +% +% Inputs: +% +% sample_data [struct] - the toolbox dataset. +% +% Outputs: +% +% vnames [cell{str}] - A cell with velocity variable names. +% +% Example: +% +% %basic usage +% dims = IMOS.gen_dimensions('adcp'); +% vars = IMOS.gen_variables(dims,{'VEL1','VEL2','UCUR','VCUR'}); +% x.variables = vars; +% vnames = IMOS.meta.velocity_variables(x); +% assert(isempty(setdiff(vnames,{'UCUR','VCUR','VEL1','VEL2'}))); +% x.variables = IMOS.gen_variables(dims,{'UCUR_MAG','VCUR_MAG','WCUR'}); +% vnames = IMOS.meta.velocity_variables(x); +% assert(isempty(setdiff(vnames,{'UCUR_MAG','VCUR_MAG','WCUR'}))); +% +% % ambiguity is not handled. +% x.variables = IMOS.gen_variables(dims,{'UCUR','UCUR_MAG'}); +% vnames = IMOS.meta.velocity_variables(x); +% assert(isempty(setdiff(vnames,{'UCUR','UCUR_MAG'}))) +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(1,1); +try + varcell = sample_data.variables; +catch + errormsg('No variable fieldname available.') +end +avail_variables = IMOS.get(varcell,'name'); +allowed_variables = {'UCUR_MAG','UCUR','VCUR_MAG','VCUR','WCUR','WCUR_2','ECUR','VEL1','VEL2','VEL3','VEL4'}; +vnames=cell(1,numel(allowed_variables)); +c=0; +for k= 1:numel(avail_variables) + vname = avail_variables{k}; + if inCell(allowed_variables,vname) + c=c+1; + vnames{c} = vname; + end +end +vnames=vnames(1:c); +end diff --git a/Util/+IMOS/+resolve/function_parameters.m b/Util/+IMOS/+resolve/function_parameters.m new file mode 100644 index 000000000..626d2f4cd --- /dev/null +++ b/Util/+IMOS/+resolve/function_parameters.m @@ -0,0 +1,53 @@ +function [opts] = function_parameters() +% function [opts] = function_parameters() +% +% Load the current function parameters +% This typically involve loading a txt file with +% the same name as the running function, and +% converting the arguments. +% +% Inputs: +% +% +% Outputs: +% +% opts - the options loaded as map, converted +% to numbers of a number. +% +% Example: +% %see imosSurfaceDetectionQC.m +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(0, 0) +opts = containers.Map(); + +called_by_function = callerName(); +if isempty(called_by_function) + return +end +fpath = which(called_by_function); +[folder, fname] = fileparts(fpath); +ofile_path = [fullfile(folder, fname) '.txt']; + +missing_file = ~exist(ofile_path, 'file'); +if missing_file + return +end + +delimiter = detectMappingDelimiter(ofile_path); +opts = readMappings(ofile_path,delimiter); +keys = opts.keys; + +for k = 1:numel(keys) + name = keys{k}; + nvalue = str2double(opts(name)); + + if strcmpi(opts(name), 'nan') || ~isnan(nvalue) + opts(name) = nvalue; + end + +end + +end diff --git a/Util/+IMOS/+schemas/imosEchoIntensitySetQC.m b/Util/+IMOS/+schemas/imosEchoIntensitySetQC.m new file mode 100644 index 000000000..fd28f3a71 --- /dev/null +++ b/Util/+IMOS/+schemas/imosEchoIntensitySetQC.m @@ -0,0 +1,50 @@ +function [bool, reason] = imosEchoIntensitySetQC(sample_data) +% function [bool,reason] = imosEchoIntensitySetQC(sample_data) +% +% Check if sample_data is a valid input for SurfaceDetectionByEchoIntensityQC. +% +% Inputs: +% +% sample_data [struct] - A toolbox dataset. +% +% Outputs: +% +% bool - True if dataset is valid. False otherwise. +% reason - The reasons why the dataset is invalid. +% +% Example: +% +% %see test_imosEchoIntensitySetQC. +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(1,1) +reason = {}; + +if ~IMOS.adcp.contains_adcp_dimensions(sample_data) + reason{1} = 'Not an adcp file.'; +end + +avail_variables = IMOS.get(sample_data.variables,'name'); +absic_counter = sum(contains(avail_variables,'ABSIC')); +if absic_counter == 0 + reason{end+1} = 'Missing ABSIC variables.'; +end + +vel_vars = IMOS.meta.velocity_variables(sample_data); +if numel(vel_vars) == 0 + reason{end+1} = 'Missing at leat one velocity variable to flag.'; +end + + +if absic_counter > 4 + reason{end+1} = 'Number of ABSIC variables is invalid'; +end + +if isempty(reason) + bool = true; +else + bool = false; +end + +end diff --git a/Util/+IMOS/+schemas/imosSurfaceDetectionByDepthSetQC.m b/Util/+IMOS/+schemas/imosSurfaceDetectionByDepthSetQC.m new file mode 100644 index 000000000..6a4f62f8e --- /dev/null +++ b/Util/+IMOS/+schemas/imosSurfaceDetectionByDepthSetQC.m @@ -0,0 +1,89 @@ +function [bool, reason] = imosSurfaceDetectionByDepthSetQC(sample_data) +% function [bool,reason] = imosSurfaceDetectionByDepthSetQC(sample_data) +% +% Check if sample_data is a valid input for SurfaceDetectionByDepthSetQC. +% +% Inputs: +% +% sample_data [struct] - A toolbox dataset. +% +% Outputs: +% +% bool - True if dataset is valid. False otherwise. +% reason - The reasons why the dataset is invalid. +% +% Example: +% +% %see test_imosSurfaceDetectionByDepthSetQC.m +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(1,1) +reason = {}; + +if ~IMOS.adcp.contains_adcp_dimensions(sample_data) + reason{1} = 'Not an adcp file.'; +end + +try + 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 && ~along_beams + reason{end + 1} = 'Missing adcp bin dimensions'; + end + + time = IMOS.get_data(sample_data.dimensions, 'TIME'); + + if isempty(time) + reason{end + 1} = 'No TIME dimension.'; + end + + 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'); + else + bin_dist = []; + end + + if isempty(bin_dist) + reason{end + 1} = 'No bin distance dimensions.'; + end + +catch + reason{end + 1} = 'No dimensions in dataset.'; +end + +try + idepth = IMOS.get_data(sample_data.variables, 'DEPTH'); + + if isempty(idepth) + reason{end + 1} = 'no DEPTH variable.'; + end + +catch + reason{end + 1} = 'No variables in dataset.'; +end + +bathy = IMOS.meta.site_bathymetry(sample_data); + +if isempty(bathy) + reason{end + 1} = 'no bathymetry metadata.'; +end + +beam_face_config = IMOS.meta.beam_face_config(sample_data); + +if ~strcmpi(beam_face_config, 'up') && ~strcmpi(beam_face_config, 'down') + reason{end + 1} = 'unknown ADCP beam face config or inconsistent sensor height values.'; +end + +if isempty(reason) + bool = true; +else + bool = false; +end + +end diff --git a/Util/+IMOS/+update/function_parameters.m b/Util/+IMOS/+update/function_parameters.m new file mode 100644 index 000000000..9b68b8ea0 --- /dev/null +++ b/Util/+IMOS/+update/function_parameters.m @@ -0,0 +1,38 @@ +function function_parameters(optname,optvalue) +% function function_parameters(optname,optvalue) +% +% Update the current function parameters +% This opens a txt file with the same name +% as the running function, and +% update the respective option with the respective value. +% +% Inputs: +% +% optname - the name of the option. +% optvalue - the value of the option. +% +% Outputs: +% +% +% Example: +% %see usage in imosEchoIntensityQC. +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(2, 2) +called_by_function = callerName(); +if isempty(called_by_function) + return +end +fpath = which(called_by_function); +[folder, fname] = fileparts(fpath); +ofile_path = [fullfile(folder, fname) '.txt']; + +missing_file = ~exist(ofile_path, 'file'); +if missing_file + errormsg('no parameter file for %s',fname); +end + +updateMappings(ofile_path,optname,optvalue); +end diff --git a/Util/+IMOS/validate_dataset.m b/Util/+IMOS/validate_dataset.m new file mode 100644 index 000000000..70fa156c8 --- /dev/null +++ b/Util/+IMOS/validate_dataset.m @@ -0,0 +1,49 @@ +function [bool, reason] = validate_dataset(sample_data, fname) +% function [bool,reason] = validate_dataset(sample_data,fname) +% +% Validate a toolbox dataset against a function name. +% +% The idea of this function is to execute the associated +% `fname` schema function, which verifies the feasability +% of the sample_data content for use within `fname` function. +% +% Think all the checks you do at a function but wrapped in +% one call. +% +% Inputs: +% +% sample_data [struct] - a toolbox sample data. +% fname [string] = the function name. +% +% Outputs: +% +% bool[logical] - True if dataset is accepted by fname. +% reason [cell{str}] - The reasons why the dataset is invalid for fname. +% +% Example: +% +% %basic usage +% +% +% author: hugo.oliveira@utas.edu.au +% +func = ['IMOS.schemas.' fname]; +try + feval(func) + missing_schema = true; +catch me + missing_schema = ~strcmpi(me.message,'Not enough input arguments.'); +end + +if missing_schema + errormsg('Schema %s not available',fname); +end + +try + [bool, reason] = feval(func, sample_data); +catch + bool = false; + reason{1} = sprintf('%s failed.', func); +end + +end diff --git a/Util/+IMOS/variables_with_dimensions.m b/Util/+IMOS/variables_with_dimensions.m new file mode 100644 index 000000000..894849965 --- /dev/null +++ b/Util/+IMOS/variables_with_dimensions.m @@ -0,0 +1,66 @@ +function [varnames] = variables_with_dimensions(dims, vars, dnamecell) +% function [varnames] = variables_with_dimensions(dims, vars, dnamecell) +% +% Return variable names with certain named dimensions, in order. +% +% Inputs: +% +% dims [cell{struct}] - The toolbox dimensions cell. +% vars [cell{struct}] - The toolbox variables cell. +% dnamecell [cell{str}] - A cell with dimension names. +% +% Outputs: +% +% varnames [cell{str}] - A cell with variable names. +% +% Example: +% +% %basic usage +% +% dims = IMOS.gen_dimensions('timeSeries',3,{'x','y','z'},{@double,@double,@double},{(1:5)',(1:10)',(1:15)'}); +% vars = IMOS.gen_variables(dims,{'a','b','c'},{@double,@double,@double},{zeros(5,10),zeros(1,10),zeros(5,15)}); +% a = IMOS.variables_with_dimensions(dims,vars,{'x','y'}); +% assert(isequal(a,{'a'})) +% b = IMOS.variables_with_dimensions(dims,vars,{'y'}); +% assert(isequal(b,{'b'})) +% c = IMOS.variables_with_dimensions(dims,vars,{'x','z'}); +% assert(isequal(c,{'c'})) +% none = IMOS.variables_with_dimensions(dims,vars,{'x'}); +% assert(isempty(none)) +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(3, 3) + +if ~IMOS.is_toolbox_dimcell(dims) + errormsg('First argument is not a toolbox dimensions cell') +end + +if ~IMOS.is_toolbox_varcell(vars) + errormsg('Second argument is not a toolbox variables cell') +end + +if ~iscellstr(dnamecell) + errormsg('Third argument is not a cell of strings') +end + +allvars = IMOS.get(vars, 'name'); +nvars = numel(allvars); +varnames = cell(1, nvars); +c = 0; + +for k = 1:numel(allvars) + vname = allvars{k}; + vdims = IMOS.get_dimension_names(dims, vars, vname); + contain_dim = isequal(vdims, dnamecell); + + if contain_dim + c = c + 1; + varnames{c} = vname; + end + +end + +varnames = varnames(1:c); +end diff --git a/Util/callerName.m b/Util/callerName.m new file mode 100644 index 000000000..47be6f30c --- /dev/null +++ b/Util/callerName.m @@ -0,0 +1,31 @@ +function [name] = callerName() +% function [name] = callerName() +% +% Return the name of the function that +% called this function. +% +% Inputs: +% +% +% Outputs: +% +% name - the function name. +% +% Example: +% %only works here since this is evaluated. +% name = callerName(); +% assert(strcmpi(callerName,'testDocstring')); +% %outside any function +% %assert(callerName(),'user-interaction'); +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(0,0); +s = dbstack(2); +if isempty(s) + name = ''; +else + name = s(1).name; +end +end diff --git a/Util/detectMappingDelimiter.m b/Util/detectMappingDelimiter.m new file mode 100644 index 000000000..29b8da894 --- /dev/null +++ b/Util/detectMappingDelimiter.m @@ -0,0 +1,60 @@ +function [delimiter] = detectMappingDelimiter(fname) +% function [delimiter] = detectMappingDelimiter(fname) +% +% A simple delimiter detection for toolbox mapping/Property files. +% +% Inputs: +% +% fname - the text file. +% +% Outputs: +% +% delimiter - The delimiter used. +% +% Example: +% +% %basic usage +% comma_file = [toolboxRootPath 'GUI/instrumentAliases.txt']; +% assert(isequal(detectMappingDelimiter(comma_file),',')) +% equal_file = [toolboxRootPath 'AutomaticQC/imosEchoIntensityQC.txt']; +% assert(isequal(detectMappingDelimiter(equal_file),'=')) +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(1,1) + +delimiter = ''; +comma_ok = false; +equal_ok = false; +if ~exist(fname,'file') + errormsg('%s doesn''t exist.',fname) +end + + +try + comma_opts = readMappings(fname,','); + comma_ok = true; +catch +end + +try + equal_opts = readMappings(fname,'='); + equal_ok = true; +catch +end + +if comma_ok && equal_ok + if all(contains(comma_opts.keys,'=')) + delimiter = '='; + elseif all(contains(equal_opts.keys,',')) + delimiter = ','; + end + return +elseif comma_ok + delimiter = ','; +elseif equal_ok + delimiter = '='; +end + +end diff --git a/Util/dispmsg.m b/Util/dispmsg.m new file mode 100644 index 000000000..112a3b254 --- /dev/null +++ b/Util/dispmsg.m @@ -0,0 +1,29 @@ +function dispmsg(msg,varargin) +% function dispmsg(msg,varargin) +% +% A Wrapper to disp, but showing where +% the function is being called. +% +% Inputs: +% +% msg - A message string. +% varargin - sprintf/further arguments (e.g.fields for msg). +% +% +% Example: +% +% %basic usage +% +% +% author: hugo.oliveira@utas.edu.au +% +if nargin < 1 + error('Error at dispmsg(line 21): msg argument is compulsory') +end + +cstack = dbstack(1); +name = cstack(1).name; +line = cstack(1).line; +actualmsg = sprintf([name '(' num2str(line) '): ' msg],varargin{:}); +disp(actualmsg) +end diff --git a/Util/errormsg.m b/Util/errormsg.m index f365bc8ce..6e11525c3 100644 --- a/Util/errormsg.m +++ b/Util/errormsg.m @@ -23,12 +23,14 @@ function errormsg(msg,varargin) cstack = dbstack; if numel(cstack) == 1 - name = cstack(1).name; + n = 1; else - name = cstack(2).name; + n = 2; end +name = cstack(n).name; +line = cstack(n).line; -actualmsg = sprintf(msg,varargin{:}); +actualmsg = sprintf([name '(' num2str(line) '):' msg],varargin{:}); dotsplit = split(name,'.'); if numel(dotsplit) == 1 name_as_id = sprintf('%s:%s',dotsplit{1},dotsplit{1}); diff --git a/Util/genSampleDataDesc.m b/Util/genSampleDataDesc.m index 9728ec9be..d93b0e02d 100644 --- a/Util/genSampleDataDesc.m +++ b/Util/genSampleDataDesc.m @@ -17,24 +17,6 @@ % % Author: Paul McCarthy % - -% -% 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 . -% narginchk(1,2); if ~isstruct(sam), error('sam must be a struct'); end @@ -70,18 +52,36 @@ timeFmt = readProperty('toolbox.timeFormat'); -timeRange = ['from ' datestr(sam.time_coverage_start, timeFmt) ' to ' ... - datestr(sam.time_coverage_end, timeFmt)]; +try + time_coverage_start = sam.time_coverage_start; +catch + time_coverage_start = 0; +end + +try + time_coverage_end = sam.time_coverage_end; +catch + time_coverage_end = 0; +end + +timeRange = ['from ' datestr(time_coverage_start, timeFmt) ' to ' ... + datestr(time_coverage_end, timeFmt)]; + +try + filename = sam.toolbox_input_file; +catch + filename = 'Unknown'; +end [~, fName, fSuffix] = fileparts(sam.toolbox_input_file); fName = [fName fSuffix]; -alias_file = ''; try alias_file = readProperty('toolbox.instrumentAliases'); catch + alias_file = ''; end instrument_entry = [sam.meta.instrument_make ' ' sam.meta.instrument_model]; @@ -93,29 +93,47 @@ end end +try + depth = sam.meta.depth; +catch + depth = NaN; +end + +try + instrument_serial_no = sam.meta.instrument_serial_no; +catch + instrument_serial_no = 'Unknown'; +end + +try + site_id = sam.meta.site_id; +catch + site_id = 'Unknown'; +end + switch detailLevel case 'name-only' - desc = [ instrument_entry ]; + desc = instrument_entry ; case 'short' - desc = [ instrument_entry ' @' num2str(sam.meta.depth) 'm']; + desc = [ instrument_entry ' @' num2str(depth) 'm']; case 'medium' desc = [ instrument_entry ... - ' SN=' sam.meta.instrument_serial_no ... - ' @' num2str(sam.meta.depth) 'm' ... + ' SN=' instrument_serial_no ... + ' @' num2str(depth) 'm' ... ' (' fName ')']; case 'id' - desc = [ '(' fName ')' ' SN=' sam.meta.instrument_serial_no ' @' num2str(sam.meta.depth) 'm']; + desc = [ '(' fName ')' ' SN=' instrument_serial_no ' @' num2str(depth) 'm']; otherwise % full details - desc = [ sam.meta.site_id ... + desc = [ site_id ... ' - ' instrument_entry ... - ' SN=' sam.meta.instrument_serial_no ... - ' @' num2str(sam.meta.depth) 'm' ... + ' SN=' instrument_serial_no ... + ' @' num2str(depth) 'm' ... ' ' timeRange ... ' (' fName ')']; end diff --git a/Util/isinside.m b/Util/isinside.m index db317694b..7a8c73185 100644 --- a/Util/isinside.m +++ b/Util/isinside.m @@ -21,7 +21,18 @@ % assert(isinside({'a','b','c'},{'a','b','c'})) % assert(~isinside({'a'},{'a','b','c'})) % +% %ignore matlab inconsistent comparison against cell of different shapes +% assert(isinside({'a','b','c'},{'a','b'}')) +% % author: hugo.oliveira@utas.edu.au % -bool = isequal(unique(intersect(obj,items)),unique(items)); +a = unique(intersect(obj,items)); +b = unique(items); + +if iscell(obj) || iscell(items) + bool = isequal(a(:),b(:)); +else + bool = isequal(a,b); +end + end diff --git a/Util/readDatasetParameter.m b/Util/readDatasetParameter.m index be51391d7..5aafae833 100644 --- a/Util/readDatasetParameter.m +++ b/Util/readDatasetParameter.m @@ -68,7 +68,7 @@ % for .pqc files. [pPath, oldPFile, ~] = fileparts(rawDataFile); oldPFile = fullfile(pPath, [oldPFile, '.', pType]); -if exist(oldPFile, 'file') +if exist(oldPFile, 'file') && ~strcmpi(oldPFile,pFile) movefile(oldPFile, pFile); end diff --git a/Util/readMappings.m b/Util/readMappings.m index b8ada8bbe..eeb85a63d 100644 --- a/Util/readMappings.m +++ b/Util/readMappings.m @@ -37,13 +37,12 @@ nf = fopen(file, 'r'); raw_read = textscan(nf, '%s', 'Delimiter', delimiter, 'CommentStyle', '%'); raw_read = raw_read{1}; +fclose(nf); try - mappings = containers.Map(raw_read(1:2:end), raw_read(2:2:end)); + mappings = containers.Map(strip(raw_read(1:2:end)), strip(raw_read(2:2:end)), 'UniformValues', false); catch error('Mapping file %s is incomplete.', file); end -fclose(nf); - end diff --git a/Util/updateMappings.m b/Util/updateMappings.m new file mode 100644 index 000000000..bd3b4b6b1 --- /dev/null +++ b/Util/updateMappings.m @@ -0,0 +1,123 @@ +function updateMappings(fname, name, value) +% function updateMappings(fname, name, value) +% +% This updates a mapping/Property in a mapping file, +% preserving the previous state of the file, +% including comments. +% +% A mapping/property file is a text file with the +% the first column of the file is a key, +% the second column is a value. +% +% The function will fail if the property can't be found. +% +% +% Inputs: +% +% fname [char] - the filename path. +% name [char] - an option/mapping name. +% value [char] - an option/mapping value. +% +% Outputs: +% +% +% Example: +% +% file = [toolboxRootPath 'test_updateMappings']; +% nf = fopen(file,'w'); +% fprintf(nf,'%s\n','%1st line is a comment'); +% fprintf(nf,'%s\n','zero = a'); +% fprintf(nf,'%s\n','%3rd line is another comment'); +% fprintf(nf,'%s\n','one = b'); +% fprintf(nf,'%s\n','%5th line is a comment'); +% fprintf(nf,'%s\n','%6th line is also comment'); +% fprintf(nf,'%s\n','spam = good'); +% fprintf(nf,'%s\n','spam_and_spam = better'); +% fprintf(nf,'%s\n','spam_and_eggs = heaven'); + +% fclose(nf); +% %fail at undefined name. +% try; updateMappings(file,'two','2');catch f=true;end +% assert(f) +% +% %pass - but file is untouched. +% updateMappings(file,'zero','a'); +% +% %update pass - file is rewritten, comments are kept. +% updateMappings(file,'spam_and_eggs','and_spam'); +% updateMappings(file,'spam_and_spam','and_eggs'); +% updateMappings(file,'one',1); +% updateMappings(file,'zero',0); +% opts = readMappings(file,'='); +% assert(strcmpi(opts('zero'),'0')); +% assert(strcmpi(opts('one'),'1')); +% assert(strcmpi(opts('spam'),'good')); +% assert(strcmpi(opts('spam_and_eggs'),'and_spam')); +% assert(strcmpi(opts('spam_and_spam'),'and_eggs')); +% nf = fopen(file,'r'); +% text = textscan(nf,'%s','Delimiter','\n'); +% text = strip(text{1}); +% fclose(nf); +% assert(strcmpi(text{3},'%3rd line is another comment')); +% assert(strcmpi(text{6},'%6th line is also comment')); +% +% +% author: hugo.oliveira@utas.edu.au +% +narginchk(3,3) +if ~ischar(fname) + errormsg('First argument not a char') +end +if ~ischar(name) + errormsg('Second argument not a char') +end + +delimiter = detectMappingDelimiter(fname); +nf = fopen(fname,'r'); +raw_text = textscan(nf,'%s','Delimiter',delimiter); +if isempty(raw_text) + errormsg('File %s is empty.',fname); +end +raw_text = raw_text{1}; +text = strip(raw_text); +text_value = strip(num2str(value)); + +[found,cind] = inCell(text,name); +fclose(nf); + +if ~found + errormsg('option %s not in parameter file %s.',name,fname) +end + +missing_value = cind+1 > length(text) || strcmpi(text{cind+1}(1),'%'); +if missing_value + errormsg('option %s is missing at file %s.',name,fname); +else + update_value = cind+1 <= length(text) && ~strcmpi(text{cind+1},text_value); +end + +if update_value + nf = fopen(fname,'w'); + is_name = false; + for k=1:length(text) + is_comment = strcmpi(text{k}(1),'%'); + if is_comment + fprintf(nf,'%s\n',text{k}); + is_name = false; + elseif ~is_name + fprintf(nf,['%s ' delimiter ' '],text{k}); + is_name = true; + else + is_update = k == cind+1; + if is_update + fprintf(nf,'%s\n',text_value); + else + fprintf(nf,'%s\n',text{k}); + end + is_name = false; + end + end + fclose(nf); +end + +end diff --git a/Util/writeDatasetParameter.m b/Util/writeDatasetParameter.m index 4f85d3b17..ecc9dbea9 100644 --- a/Util/writeDatasetParameter.m +++ b/Util/writeDatasetParameter.m @@ -60,12 +60,11 @@ function writeDatasetParameter(rawDataFile, routine, param, value) return; end pFile = [rawDataFile, '.', pType]; - % we need to migrate any remnants of the old file naming convention % for .pqc files. [pPath, oldPFile, ~] = fileparts(rawDataFile); oldPFile = fullfile(pPath, [oldPFile, '.', pType]); -if exist(oldPFile, 'file') +if exist(oldPFile, 'file') && ~strcmpi(oldPFile,pFile) movefile(oldPFile, pFile); end diff --git a/test/AutomaticQC/test_imosEchoIntensitySetQC.m b/test/AutomaticQC/test_imosEchoIntensitySetQC.m new file mode 100644 index 000000000..34920acfa --- /dev/null +++ b/test/AutomaticQC/test_imosEchoIntensitySetQC.m @@ -0,0 +1,265 @@ +classdef test_imosEchoIntensitySetQC < matlab.unittest.TestCase + + properties (TestParameter) + zdim = {'HEIGHT_ABOVE_SENSOR','DIST_ALONG_BEAMS'} + end + + methods (TestMethodTeardown) + function reset_default_options(~) + file = which('imosEchoIntensitySetQC'); + file(end:end+2) = 'txt'; + updateMappings(file,'echo_intensity_threshold',50); + updateMappings(file,'bound_by_depth',0); + updateMappings(file,'bound_by_index',0); + updateMappings(file,'bound_value',99999); + updateMappings(file,'propagate',0); + end + end + + methods (Test) + + function test_simple_detection(~,zdim) + % only data 500:600,11 is bad. + sample_data = create_simple_data(zdim); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:13)==1,'all')) + assert(all(flag(500:600,14)==4,'all')) + assert(all(flag(500:600,15:end)==1,'all')) + end + + function test_simple_detection_upward_propagate(~,zdim) + % only data 500:600,11 is bad. + sample_data = create_simple_data(zdim); + switch_option('propagate',1); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:13)==1,'all')) + assert(all(flag(500:600,14:end)==4,'all')) + end + + function test_simple_detection_downward_propagate(~,zdim) + % only data 500:600,11 is bad. + sample_data = create_simple_data(zdim,-1*(10:10:200)'); + switch_option('propagate',1); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:13)==1,'all')) + assert(all(flag(500:600,14:end)==4,'all')) + end + + function test_bound_by_depth_upward_clean_marked(~,zdim) + depth_constraint = 30; %only shallower than 30m is allowed markings. + site_nominal_depth = 170; % 14th bad bin will be at 30m, so will be clear out + sample_data = create_simple_data(zdim,(10:10:200)',site_nominal_depth); + switch_option('bound_by_depth',1); + switch_option('bound_value',depth_constraint); + switch_option('propagate',0); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag==1,'all')) + end + + function test_bound_by_depth_downward_clean_marked(~,zdim) + depth_constraint = 140; %only deeper than 140m is allowed markings. + site_nominal_depth = 0; % 14th bin will be clear. + sample_data = create_simple_data(zdim,-1*(10:10:200)',site_nominal_depth); + switch_option('bound_by_depth',1); + switch_option('bound_value',depth_constraint); + switch_option('propagate',0); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag==1,'all')) + end + + function test_bound_by_depth_upward_simple(~,zdim) + depth_constraint = 40; %above 40m to be marked. + site_nominal_depth = 140 + depth_constraint; % bindepth of 14th bin + constraint + zrange = 0; + sample_data = create_simple_data(zdim,(10:10:200)',site_nominal_depth,zrange); + echo_ind = IMOS.find(sample_data.variables,'ABSIC1'); + sample_data.variables{echo_ind}.data(500:600,14) = 205; % increase gradient so 15th bin is also marked. + % we now unmark all bins deeper than 40m. and keep everything as detected above 40m. + switch_option('bound_by_depth',1); + switch_option('bound_value',depth_constraint); + switch_option('propagate',0); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:14)==1,'all')) + assert(all(flag(500:600,15)==4,'all')) + assert(all(flag(500:600,16:end)==1,'all')) + end + + function test_bound_by_depth_upward_zeta(~,zdim) + depth_constraint = 40; %above 40m to be marked. + site_nominal_depth = 140 + depth_constraint; % bindepth of 14th bin + constraint + zrange = 2; + sample_data = create_simple_data(zdim,(10:10:200)',site_nominal_depth,zrange); + echo_ind = IMOS.find(sample_data.variables,'ABSIC1'); + sample_data.variables{echo_ind}.data(500:600,14) = 205; % increase gradient so 15th bin is also marked. + % we now unmark some bins deeper than 40m. and keep everything as detected shallower than 40m. + switch_option('bound_by_depth',1); + switch_option('bound_value',depth_constraint); + switch_option('propagate',0); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:13)==1,'all')) + assert(any(flag(500:600,1:14)==4,'all')) + assert(all(flag(500:600,15)==4,'all')) + assert(all(flag(500:600,16:end)==1,'all')) + switch_option('propagate',1) + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:13)==1,'all')) + assert(any(flag(500:600,1:14)==4,'all')) + assert(all(flag(500:600,15:end)==4,'all')) + + end + + function test_bound_by_depth_downward_simple(~,zdim) + depth_constraint = 140; %only deeper than 600m is allowed markings. + site_nominal_depth = 300; % bindepth of 14th bin will be 440m + sample_data = create_simple_data(zdim,-1*(10:10:200)',site_nominal_depth); + echo_ind = IMOS.find(sample_data.variables,'ABSIC1'); + sample_data.variables{echo_ind}.data(500:600,14) = 205; % increase gradient so 15th bin is also marked. + % now unmark bins at 440m but keep all detected towards the bottom. + switch_option('bound_by_depth',1); + switch_option('bound_value',depth_constraint); + switch_option('propagate',0); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:14)==1,'all')) + assert(all(flag(500:600,15)==4,'all')) + assert(all(flag(500:600,16:end)==1,'all')) + end + + function test_bound_by_depth_downward_zeta(~,zdim) + depth_constraint = 140; %only deeper than 600m is allowed markings. + site_nominal_depth = 300; % bindepth of 14th bin will be 440m + zrange = 2; + sample_data = create_simple_data(zdim,-1*(10:10:200)',site_nominal_depth,zrange); + echo_ind = IMOS.find(sample_data.variables,'ABSIC1'); + sample_data.variables{echo_ind}.data(500:600,14) = 205; % increase gradient so 15th bin is also marked. + % now unmark some bins around 40m but keep all beyond 40m to the surface. + switch_option('bound_by_depth',1); + switch_option('bound_value',depth_constraint); + switch_option('propagate',0); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:13)==1,'all')) + assert(any(flag(500:600,1:14)==4,'all')) + assert(all(flag(500:600,15)==4,'all')) + assert(all(flag(500:600,16:end)==1,'all')) + switch_option('propagate',1); + new = imosEchoIntensitySetQC(sample_data); + flag = new.variables{end}.flags; + assert(all(flag(500:600,1:13)==1,'all')) + assert(any(flag(500:600,1:14)==4,'all')) + assert(all(flag(500:600,15:end)==4,'all')) + end + + function test_realdata_bound_index(~) + % this file is clipped to out of water points. + % There are several above the threshold at the first bin interface and two points at bin26. + % We will ignore the first bin interface gradients with a bound index, and match the two bad data. + adcp_file = fullfile(toolboxRootPath, 'data/testfiles/Teledyne/workhorse/v000/beam/16072000.000.reduced'); + sample_data = load_binmapped_sample_data(adcp_file); + time = IMOS.get_data(sample_data.dimensions,'TIME'); + height_above_sensor = IMOS.get_data(sample_data.dimensions, 'HEIGHT_ABOVE_SENSOR'); + nbins = numel(height_above_sensor); + beam_height = IMOS.adcp.beam_height(height_above_sensor); + idepth = transpose(randomBetween(100,200,numel(time))); + sdepth = IMOS.gen_variables(sample_data.dimensions, {'DEPTH'}, {@double}, {idepth}, 'positive', 'down'); + sample_data.variables{end + 1} = sdepth{1}; + sample_data.site_nominal_depth = beam_height - mode(idepth); + + new = imosEchoIntensitySetQC(sample_data); + u_flag = new.variables{IMOS.find(new.variables,'VEL1')}.flags; + assert(all(u_flag(:,1)==1,'all')) + assert(any(u_flag(:,2)==4,'all')) + assert(all(u_flag(:,3:25)==1,'all')) + assert(sum(u_flag(:,26)==4)==2) + assert(all(u_flag(:,27:end)==1,'all')) + + switch_option('bound_by_index',1) + switch_option('bound_value',25);% 1:25 always good. + + new = imosEchoIntensitySetQC(sample_data); + u_flag = new.variables{IMOS.find(new.variables,'VEL1')}.flags; + v_flag = new.variables{IMOS.find(new.variables,'VEL2')}.flags; + w_flag = new.variables{IMOS.find(new.variables,'VEL3')}.flags; + e_flag = new.variables{IMOS.find(new.variables,'VEL4')}.flags; + a1_flag = new.variables{IMOS.find(new.variables,'ABSIC1')}.flags; + a2_flag = new.variables{IMOS.find(new.variables,'ABSIC2')}.flags; + a3_flag = new.variables{IMOS.find(new.variables,'ABSIC3')}.flags; + a4_flag = new.variables{IMOS.find(new.variables,'ABSIC4')}.flags; + assert(isequal(u_flag,v_flag,w_flag,e_flag,a1_flag,a2_flag,a3_flag,a4_flag)); + assert(all(u_flag(:,1:25)==1,'all')) + assert(sum(u_flag(:,26)==4)==2) + assert(all(u_flag(:,27:end)==1,'all')) + + converted = adcpWorkhorseVelocityBeam2EnuPP({sample_data},''); + switch_option('bound_by_index',1) + switch_option('bound_value',26) + new = imosEchoIntensitySetQC(converted{1}); + u_flag = new.variables{IMOS.find(new.variables,'UCUR_MAG')}.flags; + v_flag = new.variables{IMOS.find(new.variables,'VCUR_MAG')}.flags; + w_flag = new.variables{IMOS.find(new.variables,'WCUR')}.flags; + e_flag = new.variables{IMOS.find(new.variables,'ECUR')}.flags; + assert(isequal(u_flag,v_flag,w_flag,e_flag)); + assert(all(u_flag==1,'all')) + + end + end + +end + +function switch_option(optname,optvalue) + file = which('imosEchoIntensitySetQC'); + file(end:end+2) = 'txt'; + updateMappings(file,optname,optvalue); +end + +function [binmapped] = load_binmapped_sample_data(file) + data = workhorseParse({file}, ''); + binmapped = adcpBinMappingPP({data}, ''); + binmapped = binmapped{1}; +end + +function [sample_data] = create_simple_data(zdim_name,bin_dist, site_nominal_depth,zrange) + % Create a simple test adcp dataset + % The instrument depth is at the bottom for upward looking, and at the surface + % for downward looking. + % The depth is contaminated with oscillations of ~[-zrange,+zrange] dbar, + % and echo intensity is randomized between 100-150 counts. + % A detectable echo intensity count is found between 500-600 timestamps + % at the 14th bin with magnitude == 51 counts, which is above + % the default threshold of the test. + time = (1:1000)'; + if nargin<2 + bin_dist = (10:10:200)'; + end + if nargin<3 + site_nominal_depth = NaN; + end + if nargin<4 + zrange = 0; + end + dims = IMOS.gen_dimensions('adcp', 2, {'TIME', zdim_name}, {@double, @double}, {time, bin_dist}); %upward looking == positive depth + zeta = randomBetween(-zrange,zrange,numel(time)); + if all(bin_dist>0) + instrument_depth = transpose(site_nominal_depth + zeta); + else + instrument_depth = transpose(zeta); + end + vars_0 = randomBetween(100,150,1000*20); + vars_0 = reshape(vars_0,1000,20); + vars_0(500:600,14) = 201; % default threshold is 50, so fail at index 14. + vars_0(500:600,15) = 151; %drop off to only actually mark at index 14. + vars_0(500:600,16) = 141; %as above + vars = IMOS.gen_variables(dims, {'DEPTH','UCUR', 'ABSIC1'}, {@double, @double,@double}, {instrument_depth, vars_0, vars_0}); + sample_data.dimensions = dims; + sample_data.variables = vars; + sample_data.site_nominal_depth = site_nominal_depth; + sample_data.toolbox_input_file = 'testdriven'; +end diff --git a/test/AutomaticQC/test_imosSurfaceDetectionByDepthSetQC.m b/test/AutomaticQC/test_imosSurfaceDetectionByDepthSetQC.m new file mode 100644 index 000000000..0f623f717 --- /dev/null +++ b/test/AutomaticQC/test_imosSurfaceDetectionByDepthSetQC.m @@ -0,0 +1,184 @@ +classdef test_imosSurfaceDetectionByDepthSetQC < matlab.unittest.TestCase + + methods (Test) + + function test_simple_binmapped_detection_minimal_input(~) + sample_data = create_simple_data('HEIGHT_ABOVE_SENSOR'); + new = imosSurfaceDetectionByDepthSetQC(sample_data); + flag = new.variables{1}.flags; + first_bin_at_surface = find(sum(flag == 1, 1) == 0, 1, 'first'); + assert(first_bin_at_surface == 10); + end + + function test_simple_raw_detection_minimal_input(~) + sample_data = create_simple_data('DIST_ALONG_BEAMS'); + new = imosSurfaceDetectionByDepthSetQC(sample_data); + flag = new.variables{1}.flags; + first_bin_at_surface = find(sum(flag == 1, 1) == 0, 1, 'first'); + assert(first_bin_at_surface == 10); + end + + function test_raw_file_downward_adcp_bottom_detection(~) + adcp_file = fullfile(toolboxRootPath, 'data/testfiles/Teledyne/workhorse/v000/beam/16072000.000.reduced'); + sample_data = workhorseParse({adcp_file},''); + time = IMOS.get_data(sample_data.dimensions, 'TIME'); + height_above_sensor = IMOS.get_data(sample_data.dimensions, 'DIST_ALONG_BEAMS'); + beam_height = IMOS.adcp.beam_height(height_above_sensor); + + % use a instrument deployment depth between 100~200 + % which is the original file depth range. + min_allowed_depth = 100; + max_allowed_depth = 200; + idepth = transpose(randomBetween(min_allowed_depth, max_allowed_depth, numel(time))); %1d vars are stored as col vectors. + sdepth = IMOS.gen_variables(sample_data.dimensions, {'DEPTH'}, {@double}, {idepth}, 'positive', 'down'); + sample_data.variables{end + 1} = sdepth{1}; + sample_data.site_nominal_depth = beam_height - mode(idepth); + + new = imosSurfaceDetectionByDepthSetQC(sample_data); + + vars_to_check = IMOS.variables_with_dimensions(new.dimensions, new.variables, {'TIME', 'DIST_ALONG_BEAMS'}); + vars_inds = IMOS.find(new.variables, vars_to_check); + + for k = 1:length(vars_inds) + assert(isfield(new.variables{k}, 'flags')) + last_water_bin = find(sum(new.variables{k}.flags == 1, 1), 1, 'last'); % look at good ones + last_waterbin_depth = abs(height_above_sensor(last_water_bin)); + assert(last_waterbin_depth + min_allowed_depth < sample_data.site_nominal_depth) + assert(last_waterbin_depth + max_allowed_depth < sample_data.site_nominal_depth) + first_badbin_depth = abs(height_above_sensor(last_water_bin + 1)); + assert(first_badbin_depth - min_allowed_depth > min_allowed_depth) + end + + end + + function test_binmapped_file_downward_adcp_bottom_detection(~) + adcp_file = fullfile(toolboxRootPath, 'data/testfiles/Teledyne/workhorse/v000/beam/16072000.000.reduced'); + sample_data = load_binmapped_sample_data(adcp_file); + time = IMOS.get_data(sample_data.dimensions, 'TIME'); + height_above_sensor = IMOS.get_data(sample_data.dimensions, 'HEIGHT_ABOVE_SENSOR'); + beam_height = IMOS.adcp.beam_height(height_above_sensor); + + % use a instrument deployment depth between 100~200 + % which is the original file depth range. + min_allowed_depth = 100; + max_allowed_depth = 200; + idepth = transpose(randomBetween(min_allowed_depth, max_allowed_depth, numel(time))); %1d vars are stored as col vectors. + sdepth = IMOS.gen_variables(sample_data.dimensions, {'DEPTH'}, {@double}, {idepth}, 'positive', 'down'); + sample_data.variables{end + 1} = sdepth{1}; + sample_data.site_nominal_depth = beam_height - mode(idepth); + + new = imosSurfaceDetectionByDepthSetQC(sample_data); + + vars_to_check = IMOS.variables_with_dimensions(new.dimensions, new.variables, {'TIME', 'HEIGHT_ABOVE_SENSOR'}); + vars_inds = IMOS.find(new.variables, vars_to_check); + + for k = 1:length(vars_inds) + assert(isfield(new.variables{k}, 'flags')) + last_water_bin = find(sum(new.variables{k}.flags == 1, 1), 1, 'last'); % look at good ones + last_waterbin_depth = abs(height_above_sensor(last_water_bin)); + assert(last_waterbin_depth + min_allowed_depth < sample_data.site_nominal_depth) + assert(last_waterbin_depth + max_allowed_depth < sample_data.site_nominal_depth) + first_badbin_depth = abs(height_above_sensor(last_water_bin + 1)); + assert(first_badbin_depth - min_allowed_depth > min_allowed_depth) + end + + end + + function test_raw_file_upward_adcp_surface_detection(~) + adcp_file = fullfile(toolboxRootPath, 'data/testfiles/Teledyne/workhorse/v000/beam/1759001.000.reduced'); + sample_data = workhorseParse({adcp_file},''); + time = IMOS.get_data(sample_data.dimensions, 'TIME'); + height_above_sensor = IMOS.get_data(sample_data.dimensions, 'DIST_ALONG_BEAMS'); + beam_height = IMOS.adcp.beam_height(height_above_sensor); + + n_bad_at_surface = 3; + first_bad_bin = length(height_above_sensor) - n_bad_at_surface + 1; + bin_cell_height = mode(diff(height_above_sensor)); + + % use beam height as rigid surface with + % some oscillations with cell_height modulated magnitudes. + min_allowed_zeta = -bin_cell_height * n_bad_at_surface; + max_allowed_zeta = +bin_cell_height * n_bad_at_surface; + zeta = transpose(randomBetween(min_allowed_zeta, max_allowed_zeta, numel(time))); %1d vars are stored as col vectors. + idepth = beam_height + zeta; + sdepth = IMOS.gen_variables(sample_data.dimensions, {'DEPTH'}, {@double}, {idepth}, 'positive', 'down'); + sample_data.variables{end + 1} = sdepth{1}; + sample_data.site_nominal_depth = beam_height; + + new = imosSurfaceDetectionByDepthSetQC(sample_data); + + vars_to_check = IMOS.variables_with_dimensions(new.dimensions, new.variables, {'TIME', 'DIST_ALONG_BEAMS'}); + vars_inds = IMOS.find(new.variables, vars_to_check); + + for k = 1:length(vars_inds) + assert(isfield(new.variables{k}, 'flags')) + surface_bins = find(sum(new.variables{k}.flags == 4, 1)); % look at bad ones + assert(numel(surface_bins) == n_bad_at_surface); + first_surface_bin = min(surface_bins); + assert(first_surface_bin == first_bad_bin); + surface_bin_depth_variation = max(idepth) - beam_height; + assert(isequal_tol(n_bad_at_surface * bin_cell_height, surface_bin_depth_variation, 1)) + end + + end + + + function test_binmapped_file_upward_adcp_surface_detection(~) + adcp_file = fullfile(toolboxRootPath, 'data/testfiles/Teledyne/workhorse/v000/beam/1759001.000.reduced'); + sample_data = load_binmapped_sample_data(adcp_file); + time = IMOS.get_data(sample_data.dimensions, 'TIME'); + height_above_sensor = IMOS.get_data(sample_data.dimensions, 'HEIGHT_ABOVE_SENSOR'); + beam_height = IMOS.adcp.beam_height(height_above_sensor); + + n_bad_at_surface = 3; + first_bad_bin = length(height_above_sensor) - n_bad_at_surface + 1; + bin_cell_height = mode(diff(height_above_sensor)); + + % use beam height as rigid surface with + % some oscillations with cell_height modulated magnitudes. + min_allowed_zeta = -bin_cell_height * n_bad_at_surface; + max_allowed_zeta = +bin_cell_height * n_bad_at_surface; + zeta = transpose(randomBetween(min_allowed_zeta, max_allowed_zeta, numel(time))); %1d vars are stored as col vectors. + idepth = beam_height + zeta; + sdepth = IMOS.gen_variables(sample_data.dimensions, {'DEPTH'}, {@double}, {idepth}, 'positive', 'down'); + sample_data.variables{end + 1} = sdepth{1}; + sample_data.site_nominal_depth = beam_height; + + new = imosSurfaceDetectionByDepthSetQC(sample_data); + + vars_to_check = IMOS.variables_with_dimensions(new.dimensions, new.variables, {'TIME', 'HEIGHT_ABOVE_SENSOR'}); + vars_inds = IMOS.find(new.variables, vars_to_check); + + for k = 1:length(vars_inds) + assert(isfield(new.variables{k}, 'flags')) + surface_bins = find(sum(new.variables{k}.flags == 4, 1)); % look at bad ones + assert(numel(surface_bins) == n_bad_at_surface); + first_surface_bin = min(surface_bins); + assert(first_surface_bin == first_bad_bin); + surface_bin_depth_variation = max(idepth) - beam_height; + assert(isequal_tol(n_bad_at_surface * bin_cell_height, surface_bin_depth_variation, 1)) + end + + end + + end + +end + +function [binmapped] = load_binmapped_sample_data(file) + data = workhorseParse({file}, ''); + binmapped = adcpBinMappingPP({data}, ''); + binmapped = binmapped{1}; +end + +function [sample_data] = create_simple_data(zdim_name) + dims = IMOS.gen_dimensions('adcp', 2, {'TIME', zdim_name}, {@double, @double}, {(1:1000)', (10:10:200)'}); %upward looking == positive depth + site_nominal_depth = 91; + zeta = randomBetween(-1, 1, 1000); + instrument_depth = transpose(site_nominal_depth + zeta); + vars_0 = zeros(1000, 20); + vars = IMOS.gen_variables(dims, {'UCUR', 'DEPTH'}, {@double, @double}, {vars_0, instrument_depth}); + sample_data.dimensions = dims; + sample_data.variables = vars; + sample_data.site_nominal_depth = site_nominal_depth; +end