From f49a895211c0c784d44ff49dbb9596aa3343fa76 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Wed, 10 Mar 2021 14:22:08 +1100 Subject: [PATCH 01/22] cleanup: remove geomag tmp files --- Geomag/sample_coords.txt | 6 ------ Geomag/sample_out_IGRF12.txt | 7 ------- 2 files changed, 13 deletions(-) delete mode 100644 Geomag/sample_coords.txt delete mode 100644 Geomag/sample_out_IGRF12.txt 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 From b2f87084d7adf9d24ee2d397facd48325b2a5f82 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:33:27 +1100 Subject: [PATCH 02/22] fix(Util): only move files if their name differs Avoid errors when trying to move files with the same source/destination names. --- Util/readDatasetParameter.m | 2 +- Util/writeDatasetParameter.m | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) 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/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 From 10de06b75d69a240facdc21c425e3247c4497fde Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:34:01 +1100 Subject: [PATCH 03/22] fix(Util): isinside support to diff. cells lens --- Util/isinside.m | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) 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 From d78fc69197488b11826463576e1bd027b6d1996f Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:36:37 +1100 Subject: [PATCH 04/22] fix(errormsg): put line error in single line --- Util/errormsg.m | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) 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}); From 0ff821e3f7e50738cf808371ff48cb82ec4d8889 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:38:36 +1100 Subject: [PATCH 05/22] fix(readMappings): close file id before error --- Util/readMappings.m | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Util/readMappings.m b/Util/readMappings.m index a1ce6d2aa..eeb85a63d 100644 --- a/Util/readMappings.m +++ b/Util/readMappings.m @@ -37,9 +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 + +end From 4f11f093a27d2d1749cc103f69491d53ea520f1f Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:39:31 +1100 Subject: [PATCH 06/22] fix(genSampleDataDesc): remove license body --- Util/genSampleDataDesc.m | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/Util/genSampleDataDesc.m b/Util/genSampleDataDesc.m index 9728ec9be..e8a378361 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 From 5e28104f6224915dd0e28fc1d55fcbfdd3cc94e4 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:40:15 +1100 Subject: [PATCH 07/22] feat(genSampleDataDesc): support missing fields The updates are to support quick testing and to avoid the UI being broken by a certain datasets with some missing fields. --- Util/genSampleDataDesc.m | 58 ++++++++++++++++++++++++++++++++-------- 1 file changed, 47 insertions(+), 11 deletions(-) diff --git a/Util/genSampleDataDesc.m b/Util/genSampleDataDesc.m index e8a378361..d93b0e02d 100644 --- a/Util/genSampleDataDesc.m +++ b/Util/genSampleDataDesc.m @@ -52,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]; @@ -75,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 From aeb6523619aa465e4a63ffc8875b97e5dfd3eb4d Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:44:18 +1100 Subject: [PATCH 08/22] fix(adcpBinMappingPP): allow downward adcp This fix a bug for downlooking adcp introduced when we changed the interpolation function and cleanup a variable not used. A constraint for the interpolation function is that the values along all dimensions should be increasing. This fix just shift the value within the function closure so we can perform interpolation of downlooking bins. --- Preprocessing/adcpBinMappingPP.m | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) 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); From 41de5f8528bd190c4d48e6b49181d368c444b919 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:49:35 +1100 Subject: [PATCH 09/22] fix(readWorkhorseEnsembles): turn on memory saving In a recent commit, the memory saving option to use memory mapped files was commented out. This turn the feature on for good. --- Parser/readWorkhorseEnsembles.m | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Parser/readWorkhorseEnsembles.m b/Parser/readWorkhorseEnsembles.m index 153461080..cdb1e4479 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' From 7a9a1fd47029eb887e278834c5c72fc67cc50704 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 15:10:48 +1100 Subject: [PATCH 10/22] fix(readWorkhorseEnsembles): remove typo --- Parser/readWorkhorseEnsembles.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Parser/readWorkhorseEnsembles.m b/Parser/readWorkhorseEnsembles.m index cdb1e4479..99fa8d745 100644 --- a/Parser/readWorkhorseEnsembles.m +++ b/Parser/readWorkhorseEnsembles.m @@ -363,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. From dd47547ded807a69b388df5a7ce4416034c9481f Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 15:11:08 +1100 Subject: [PATCH 11/22] fix(readWorkhorseEnsembles): compatibility revert This fix a compatibility problem introduced recently by using the de2bi function. This function may not be avaiable to everyone. Hence, the old bit wrangling was restore for all variables affected. --- Parser/readWorkhorseEnsembles.m | 36 ++++++++++++++++++++++++++++----- 1 file changed, 31 insertions(+), 5 deletions(-) diff --git a/Parser/readWorkhorseEnsembles.m b/Parser/readWorkhorseEnsembles.m index 99fa8d745..54891aa53 100644 --- a/Parser/readWorkhorseEnsembles.m +++ b/Parser/readWorkhorseEnsembles.m @@ -383,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)); @@ -410,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); From 25433303c4924c6868636a1c877fff29e35a7791 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 15:12:53 +1100 Subject: [PATCH 12/22] fix:(adcpWorkhorseVelocityBeam2EnuPP): comments --- Preprocessing/adcpWorkhorseVelocityBeam2EnuPP.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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. From 28baf98da76a2542fa563167b9df283f1459d80a Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 15:55:50 +1100 Subject: [PATCH 13/22] feat(Util): allow update config files This commit introduces functionality to allow easy update to Property/Mapping files while keeping their comment structure. `Util/detectMappingsDelimiter` detects the delimiter of a Parameter/Mapping file. `Util/callerName.m` detects the function calling it from the current workspace. `Util/updateMappings` perform the update for a given configuration name and value. `Util/+IMOS/+resolve/function_parameters` automatically resolve/load the parameters of the current running function. This function substitutes `readProperty.m`. `Util/+IMOS/+update/function_parameters` as above, but update the file with parameters as arguments. This function, in particular, is not used yet but may be a good substitute to `writeProperty.m`. For usage examples, see imosEchoIntensitySetQC, imosSurfaceDetectionByDepthSetQC, and associated tests. --- Util/+IMOS/+resolve/function_parameters.m | 53 ++++++++++ Util/+IMOS/+update/function_parameters.m | 38 +++++++ Util/callerName.m | 31 ++++++ Util/detectMappingDelimiter.m | 60 +++++++++++ Util/updateMappings.m | 123 ++++++++++++++++++++++ 5 files changed, 305 insertions(+) create mode 100644 Util/+IMOS/+resolve/function_parameters.m create mode 100644 Util/+IMOS/+update/function_parameters.m create mode 100644 Util/callerName.m create mode 100644 Util/detectMappingDelimiter.m create mode 100644 Util/updateMappings.m 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/+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/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/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 From bad34c6a907ae60b509869ce314cc35143221fcd Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 16:08:04 +1100 Subject: [PATCH 14/22] feat(Util): disp that shows where the call came --- Util/dispmsg.m | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 Util/dispmsg.m 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 From cccffc881fee6084cc2e7350263a73f3c4fdfefc Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 16:09:22 +1100 Subject: [PATCH 15/22] feat(IMOS): detect variables with given dimensions `IMOS.variables_with_dimensions` implements a filterby/groupby action on variables based on dimensions names, in order. --- Util/+IMOS/variables_with_dimensions.m | 66 ++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 Util/+IMOS/variables_with_dimensions.m 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 From d994c7f09d961842601ed6006b43d3ea6b6cf93a Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 16:17:41 +1100 Subject: [PATCH 16/22] feat(IMOS): adcp utilities This provide some basic adcp query functions. `IMOS.adcp.beam_height` - computes the total beam height/length of adcp. `IMOS.adcp.bin_in_water` - computes the points of the adcp profile that are within the water, based on the beam_face config, beam_height, instrument depth, and local bathymetry. `IMOS.adcp.contains_adcp_dimensions` - check if a dataset contains typical adcp dimensions. `IMOS.adcp.is_along_beam` - check if the data is collected in beam coordinates. `IMOS.adcp.number_of_beams` - return the number of beams based on variable content. --- Util/+IMOS/+adcp/beam_height.m | 32 ++++++ Util/+IMOS/+adcp/bin_in_water.m | 115 ++++++++++++++++++++ Util/+IMOS/+adcp/contains_adcp_dimensions.m | 48 ++++++++ Util/+IMOS/+adcp/is_along_beam.m | 29 +++++ Util/+IMOS/+adcp/number_of_beams.m | 57 ++++++++++ 5 files changed, 281 insertions(+) create mode 100644 Util/+IMOS/+adcp/beam_height.m create mode 100644 Util/+IMOS/+adcp/bin_in_water.m create mode 100644 Util/+IMOS/+adcp/contains_adcp_dimensions.m create mode 100644 Util/+IMOS/+adcp/is_along_beam.m create mode 100644 Util/+IMOS/+adcp/number_of_beams.m 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 From fc362f05d8b6c096e385e50ba57181aa530b147d Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 16:22:42 +1100 Subject: [PATCH 17/22] feat(IMOS): metadata extraction queries Includes some repeatable patterns of query/meta data access. `IMOS.meta.beam_face_config` - discover the beam face config of an ADCP. `IMOS.meta.site_bathymetry` - discover the local bathymetry. `IMOS.meta.velocity_variables` - discover all velocity variables in a dataset. --- Util/+IMOS/+meta/beam_face_config.m | 44 +++++++++++++++++ Util/+IMOS/+meta/site_bathymetry.m | 68 +++++++++++++++++++++++++++ Util/+IMOS/+meta/velocity_variables.m | 53 +++++++++++++++++++++ 3 files changed, 165 insertions(+) create mode 100644 Util/+IMOS/+meta/beam_face_config.m create mode 100644 Util/+IMOS/+meta/site_bathymetry.m create mode 100644 Util/+IMOS/+meta/velocity_variables.m 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 From 97d1146dcb27808ed374c2008b7ad8933d9f3753 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 16:24:28 +1100 Subject: [PATCH 18/22] feat(IMOS): schema evaluation for PP/QC This implements a simple function based schema validation for PP/QC routines. The idea is to externalise some long verification calls and inspections of individual routines and turn the actual code that perform the PP/QC cleaner. The schemas must be named exactly as functions within the `+schema` folder. `IMOS.validate_dataset` perform a schema validation of a function. For practical examples, see imosEchoIntensitySetQC & imosSurfaceDetectionByDepthSetQC. --- Util/+IMOS/+schemas/imosEchoIntensitySetQC.m | 50 +++++++++++ .../imosSurfaceDetectionByDepthSetQC.m | 89 +++++++++++++++++++ Util/+IMOS/validate_dataset.m | 49 ++++++++++ 3 files changed, 188 insertions(+) create mode 100644 Util/+IMOS/+schemas/imosEchoIntensitySetQC.m create mode 100644 Util/+IMOS/+schemas/imosSurfaceDetectionByDepthSetQC.m create mode 100644 Util/+IMOS/validate_dataset.m 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/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 From 9bafc2a4f4f9dca8ce49c5a2270d606c61f7f89c Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 16:28:11 +1100 Subject: [PATCH 19/22] feat(QC): a customisable echoIntensity QC test. This re-implements the EchoIntensityQC test with the objective of deprecating EchoIntensityVelocitySetQC. The major advantages are a cleaner code and customisable behaviour: 1. Propagate - to propagate a echo gradient further down the beam. 2. Bound by index - to limit markings by bin index. 3. Bound by depth - to limit markings by depth value. See also tests and documentation. --- AutomaticQC/imosEchoIntensitySetQC.m | 184 ++++++++++++ AutomaticQC/imosEchoIntensitySetQC.txt | 10 + .../AutomaticQC/test_imosEchoIntensitySetQC.m | 265 ++++++++++++++++++ 3 files changed, 459 insertions(+) create mode 100644 AutomaticQC/imosEchoIntensitySetQC.m create mode 100644 AutomaticQC/imosEchoIntensitySetQC.txt create mode 100644 test/AutomaticQC/test_imosEchoIntensitySetQC.m 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/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 From 02658c7c2e628998397dc1649ee23b439aa89083 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 16:32:08 +1100 Subject: [PATCH 20/22] feat(QC): adcp surface detection by Depth The `imosSurfaceDetectionByDepthSetQC` function implements a simple surface detection based on the instrument depth and the beam height of the adcp. Supports both up-looking and down-looking adcps. --- .../imosSurfaceDetectionByDepthSetQC.m | 71 +++++++ .../test_imosSurfaceDetectionByDepthSetQC.m | 184 ++++++++++++++++++ 2 files changed, 255 insertions(+) create mode 100644 AutomaticQC/imosSurfaceDetectionByDepthSetQC.m create mode 100644 test/AutomaticQC/test_imosSurfaceDetectionByDepthSetQC.m 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/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 From 21bb7381fb58a3d9551f31e1b389fc63f74bee86 Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 14:23:16 +1100 Subject: [PATCH 21/22] fix: rm surfacedetectionsetQC old code This was actually a very old code imported when we merged the 675 branch. Current QC functions like `EchoIntensityVelocity`, the reimplemented one `EchoIntensity`, or the new `SurfaceDetectionByDepth` may achieve equivalent results. --- AutomaticQC/imosSurfaceDetectionSetQC.m | 270 ------------------------ 1 file changed, 270 deletions(-) delete mode 100644 AutomaticQC/imosSurfaceDetectionSetQC.m 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 From 7c26add135243fa2ba339878f6dd07c5c19998bc Mon Sep 17 00:00:00 2001 From: Hugo Oliveira Date: Thu, 1 Apr 2021 17:08:22 +1100 Subject: [PATCH 22/22] fix(imosEchoRangeSetQC): enu detection & clean This clean up the function from the license block and fix enu detection in the test given the changes done at the workhorseParse function. --- AutomaticQC/imosEchoRangeSetQC.m | 33 +++++++++++++------------------- 1 file changed, 13 insertions(+), 20 deletions(-) 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)