Merge pull request #744 from aodn/hotfix_parser_qc_development

Hotfix parser qc development

.gitignore

@@ -1,4 +1,6 @@
 *.ppp
+*.mqc
+*.pqc
 .mypy_cache/*
 Java/bin
 *.pyc
@@ -8,3 +10,10 @@ data
 dist
 data/*
 dist/*
+Geomag/sample_coords.txt
+Geomag/sample_out_IGRF13.txt
+.python-version
+.standalone_canonical_version
+.build.py@*
+*.nc
+tmp/test_*

AutomaticQC/SpikeClassifiers/imosSpikeClassifierBurstHampel.m

@@ -10,7 +10,7 @@
 % signal - A Burst 1-d signal.
 % use_burst_window  - a boolean to consider the half_window_width applied at burst scale. If false, will consider the burst series as continuous.
 % half_window_width - The half_window_width burst range (the full window size will be 1+2*half_window_width)
-% madfactor -  A multipling scale for the MAD.
+% madfactor - a scale for the mean absolute deviation which values above will be marked as a spike.
 % repeated_only - a boolean to mark spikes in burst only if they are detected more than one time (only for half_window_width>0).
 % lower_mad_limit - a lower threshold for the MAD values, which values below will be ignored.
 %
@@ -23,13 +23,11 @@
 %
 % % simple spikes
 % x = randn(1,100)*1e-2;
-% spikes = [3,7,33,92,99]
+% spikes = [3,7,33,92,99];
 % x(spikes) = 1000;
 % bduration = 6;
 % v = 1:bduration:length(x)+bduration;
-% for k=1:length(v)-1;
-% bind{k} = [v(k),min(v(k+1)-1,length(x))];
-% end
+% for k=1:length(v)-1; bind{k} = [v(k),min(v(k+1)-1,length(x))]; end
 % [dspikes] = imosSpikeClassifierBurstHampel(bind,x);
 % assert(isequal(dspikes,spikes));
 % % equal to Hampel
@@ -38,12 +36,18 @@
 % assert(isequal(dspikes,spikes));
 % assert(isequal(dspikes,dspikes2));
 %
-% % detecting entire burst as spike
+% % by ignoring the burst nature, we need to be aware
+% % that entire bursts as spike will be missing.
+%
 % x = randn(1,100)*1e-2;
-% fullburst_spiked = [3,7,13,14,15,16,17,18,33,92,99]
+% fullburst_spiked = [3,7,13,14,15,16,17,18,33,92,99];
 % x(fullburst_spiked) = 1000;
-% [dspikes] = imosSpikeClassifierBurstHampel(bind,x);
+% half_window_width = 1; % 3 point window
+% madfactor = mad([1000,0,0]); %pick only big spikes. 
+% [dspikes] = imosSpikeClassifierBurstHampel(bind,x,false,1,madfactor);
 % assert(isequal(dspikes,[3,7,33,92,99]));
+%
+% %considering the burst nature of sampling fix this
 % [dspikes2] = imosSpikeClassifierBurstHampel(bind,x,true,2);
 % assert(isequal(dspikes2,fullburst_spiked))
 %

AutomaticQC/SpikeClassifiers/imosSpikeClassifierBurstRunningStats.m

@@ -25,8 +25,7 @@
 % z = randn(1,10);
 % z(10) = 10000;
 % [spikes] = imosSpikeClassifierRunningStats(z,@mean,@std,2);
-% assert(spikes(10)==1)
-% assert(spikes(1:9)==0)
+% assert(spikes==10)
 %
 % author: [email protected]
 %

AutomaticQC/SpikeClassifiers/imosSpikeClassifierNonBurstSavGolOTSU.m

@@ -26,7 +26,7 @@
 % signal = 10*sin(omega*t);
 % signal([20]) = max(signal).^3;
 % [spikes] = imosSpikeClassifierNonBurstSavGolOTSU(signal,5,2,100,5);
-% assert(spikes,20)
+% assert(spikes==20)
 %
 % author: Unknown
 % author: Ken Ridgway

AutomaticQC/SpikeClassifiers/imosSpikeClassifierOTSU.m

@@ -21,24 +21,26 @@
 %
 % signal = [0,1,0];
 % [spikes] = imosSpikeClassifierOTSU(signal,3,1,false);
-% assert(spikes==2)
+% %limitation of otsu method - signal difference influence markups
+% assert(isequal(spikes,[1,2]))
 %
 % %non-centralized
-% signal = [0,0,0,1,0,1,0,0,0,0]
+% signal = [0,0,0,1,0,1,0,0,0,0];
 % [spikes] = imosSpikeClassifierOTSU(signal,3,1,false);
-% assert(spikes==[3,4,5,6])
-% centralized
+% assert(isequal(spikes,[3,4,5,6]))
+%
+% %centralized
 % [spikes] = imosSpikeClassifierOTSU(signal,3,1,true);
 % assert(~isequal(spikes,[4,6])) % fail
-% assert(spikes==[4]) % fail
+% assert(spikes==5) % fail
 %
 % %harmonic
 % omega = 2*pi/86400;
 % t = 0:3600:86400*5;
-% signal = 10*sin(omega*t)
+% signal = 10*sin(omega*t);
 % signal([10,20]) = max(signal)^3;
 % [spikes] = imosSpikeClassifierOTSU(signal,100,1,true);
-% assert(isequal(spikes,[10,20])
+% assert(isequal(spikes,[10,20]))
 %
 % %
 % omega = 2*pi/86400;

AutomaticQC/SpikeClassifiers/imosSpikeClassifierRunningStats.m

@@ -24,8 +24,7 @@
 % z = randn(1,10);
 % z(10) = 10000;
 % [spikes] = imosSpikeClassifierRunningStats(z,@mean,@std,2);
-% assert(spikes(10)==1)
-% assert(spikes(1:9)==0)
+% assert(spikes==10)
 %
 % author: [email protected]
 %

AutomaticQC/SpikeClassifiers/imosSpikeClassifiersLimits.m

@@ -13,11 +13,11 @@
 %
 % Example:
 %
-% len = 100
-% [limits] = imosTimeSeriesSpikeQCLimits(len)
+% len = 100;
+% [limits] = imosSpikeClassifiersLimits(len);
 % assert(limits.hampel_half_window_width.min==1)
 % assert(limits.hampel_half_window_width.max==len/2)
-% assert(limits.hampel_madfactor.max = 20)
+% assert(limits.hampel_madfactor.max==20)
 %
 % author: [email protected]
 %

AutomaticQC/SpikeClassifiers/otsu_threshold.m

@@ -18,9 +18,30 @@
 %
 % Example:
 %
+% % find different samplings
+% dsignal = diff([0:10:100, 101:110]);
+% [t] = otsu_threshold(dsignal,2);
+% assert(t>1);
+% assert(isequal(find(dsignal>t),1:10))
+% [t] = otsu_threshold(dsignal,100);
+% assert(t>1);
+%
+% % spike simple case
+% x = randn(1,100)*10;
+% spikes = [3,30:40,90];
+% x(spikes) = 100;
+% dsignal = diff(x);
+% t = otsu_threshold(dsignal,2);
+% assert(t<100);
+% assert(abs(x(3)-x(2))>t)
+% assert(abs(x(30)-x(29))>t)
+% assert(abs(x(90)-x(89))>t)
+% 
+% % simple threshold only detects the first
+% % occurence, and with a shifted index given
+% % the difference operation
+% assert(isequal(find(dsignal>t),[2,29,89]))
 %
-% [threshold] = compute_otsu_threshold(signal,nbins)
-% assert()
 %
 % author: [email protected]
 %

AutomaticQC/imosEchoIntensitySetQC.m

@@ -15,7 +15,7 @@
 % 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.
 %
 %
@@ -61,17 +61,15 @@
 
 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);
+nbeams = numel(absic_vars);
+flag_vars = IMOS.adcp.echo_intensity_variables(sample_data);
 
 switch nbeams
     case 3
@@ -163,7 +161,6 @@
 
 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);
@@ -176,7 +173,5 @@
 further_bad_bin = find(nbadbins, 1, 'last');
 
 paramsLog = ['echo_intensity_threshold=' threshold, 'propagate=', propagate, 'near_bad_bin=' near_bad_bin, 'further_bad_bin=' further_bad_bin];
-
 writeDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'echo_intensity_threshold', threshold);
-
 end

AutomaticQC/imosSurfaceDetectionByDepthSetQC.m

@@ -7,7 +7,7 @@
 % observable water column height. Any bin beyond this height
 % is marked as bad.
 %
-% This test only works for datasets with available TIME dimension, 
+% 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.
 %
@@ -51,11 +51,11 @@
 flags = ones(size(wbins), 'int8') * badFlag;
 flags(wbins) = goodFlag;
 
-if along_beams
+if along_height
+    dims_tz = {'TIME','HEIGHT_ABOVE_SENSOR'};
+else
     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);

Parser/+Workhorse/import_mappings.m

@@ -22,78 +22,116 @@
 % Examples:
 %
 % sensors = struct('Temperature',true,'Depth',true,'Conductivity',false,'SpeedOfSound',true,'Pitch',true,'Roll',true,'Heading',true);
-%
 % imap = Workhorse.import_mappings(sensors,4,'','beam');
 % fields = fieldnames(imap);
 % assert(inCell(fields,'VEL1'))
 % assert(~inCell(fields,'VCUR'))
 % assert(inCell(fields,'HEADING'))
 % assert(inCell(fields,'ABSIC4'))
+% assert(isequal(imap.('VEL1'),{'velocity','velocity1'}))
+% assert(isequal(imap.('VEL2'),{'velocity','velocity2'}))
+% assert(isequal(imap.('VEL3'),{'velocity','velocity3'}))
+% assert(isequal(imap.('VEL4'),{'velocity','velocity4'}))
 %
 % imap = Workhorse.import_mappings(sensors,4,'','earth');
 % fields = fieldnames(imap);
 % assert(inCell(fields,'VCUR'))
 % assert(~inCell(fields,'VEL1'))
 % assert(inCell(fields,'HEADING'))
 % assert(inCell(fields,'ABSIC4'))
+% assert(isequal(imap.('UCUR'),{'velocity','velocity1'}))
+% assert(isequal(imap.('VCUR'),{'velocity','velocity2'}))
+% assert(isequal(imap.('WCUR'),{'velocity','velocity3'}))
+% assert(isequal(imap.('ECUR'),{'velocity','velocity4'}))
 %
 % %MAG declination renaming.
 %
 % imap = Workhorse.import_mappings(sensors,4,'_MAG','earth');
 % fields = fieldnames(imap);
 %
 % assert(~inCell(fields,'VEL1_MAG'))
-% assert(inCell(fields,'VCUR_MAG'))
-% assert(isequal(imap.('VCUR_MAG'),{'velocity','velocity1'}))
 % assert(inCell(fields,'UCUR_MAG'))
-% assert(isequal(imap.('UCUR_MAG'),{'velocity','velocity2'}))
-%
-% assert(inCell(fields,'ECUR'))
+% assert(isequal(imap.('UCUR_MAG'),{'velocity','velocity1'}))
+% assert(inCell(fields,'VCUR_MAG'))
+% assert(isequal(imap.('VCUR_MAG'),{'velocity','velocity2'}))
 % assert(inCell(fields,'WCUR'))
-%
-% assert(inCell(fields,'CMAG1'))
-% assert(isequal(imap.('CMAG1'),{'corrMag','field1'}))
-%
-% assert(inCell(fields,'PERG4'))
-% assert(isequal(imap.('PERG4'),{'percentGood','field4'}))
+% assert(isequal(imap.('WCUR'),{'velocity','velocity3'}))
+% assert(inCell(fields,'ECUR'))
+% assert(isequal(imap.('ECUR'),{'velocity','velocity4'}))
 %
 % assert(inCell(fields,'ABSIC4'))
 % assert(isequal(imap.('ABSIC4'),{'echoIntensity','field4'}))
+% assert(inCell(fields,'CMAG3'))
+% assert(isequal(imap.('CMAG3'),{'corrMag','field3'}))
+% assert(inCell(fields,'PERG2'))
+% assert(isequal(imap.('PERG2'),{'percentGood','field2'}))
 %
 % assert(inCell(fields,'HEADING_MAG'))
 % assert(isequal(imap.('HEADING_MAG'),{'variableLeader','heading'}))
 % assert(inCell(fields,'TX_VOLT'))
 % assert(isequal(imap.('TX_VOLT'),{'variableLeader','adcChannel1'}))
 %
+% %test beam_vars
+% [~,vel_vars,beam_vars,ts_vars] = Workhorse.import_mappings(sensors,3,'_MAG','earth');
+% assert(inCell(vel_vars,'UCUR_MAG'))
+% assert(inCell(vel_vars,'WCUR'))
+% assert(inCell(beam_vars,'ABSIC1'))
+% assert(inCell(beam_vars,'CMAG2'))
+% assert(inCell(beam_vars,'PERG3'))
+% assert(inCell(ts_vars,'HEADING_MAG'))
+% assert(~inCell(ts_vars,'PSAL'))
+%
+%
 % author: [email protected]
 %
 
 narginchk(4, 4);
+
+imap = struct();
+
 switch frame_of_reference
     case 'earth'
-        vel_vars = {['VCUR' name_extension], ['UCUR' name_extension], 'WCUR'};
-        if num_beams>3
+        imap.(['UCUR' name_extension]) = {'velocity', 'velocity1'};
+        imap.(['VCUR' name_extension]) = {'velocity', 'velocity2'};
+        imap.('WCUR') = {'velocity', 'velocity3'};
+        vel_vars = {['UCUR' name_extension], ['VCUR' name_extension], 'WCUR'};
+        if num_beams == 4 
+            imap.('ECUR') = {'velocity', 'velocity4'};
             vel_vars{end+1} = 'ECUR';
         end
     case 'beam'
+        imap.('VEL1') = {'velocity', 'velocity1'};
+        imap.('VEL2') = {'velocity', 'velocity2'};
+        imap.('VEL3') = {'velocity', 'velocity3'};
         vel_vars = {'VEL1','VEL2','VEL3'};
-        if num_beams>3
+        if num_beams == 4 
+            imap.('VEL4') = {'velocity', 'velocity4'};
             vel_vars{end+1} = 'VEL4';
         end
+    otherwise
+        errormsg('Frame of reference %s not implemented.',frame_of_reference)
 end
 
-beam_vars = {'ABSIC1', 'ABSIC2', 'ABSIC3', 'ABSIC4', 'CMAG1', 'CMAG2', 'CMAG3', 'CMAG4', 'PERG1', 'PERG2', 'PERG3', 'PERG4'};
-ts_vars = {};
-
-imap = struct();
-
-for k = 1:num_beams
-    imap.(vel_vars{k}) = {'velocity', ['velocity' num2str(k)]};
+%follow the order of previous definitions for compatibility
+beam_vars = cell(1,num_beams*3); %ABSIC,CMAG,PERG
+c=0;
+for k=1:num_beams
+    c=c+1;
     imap.(['ABSIC' num2str(k)]) = {'echoIntensity', ['field' num2str(k)]}; %backscatter
+    beam_vars{c} = ['ABSIC' num2str(k)];
+end
+for k=1:num_beams
+    c=c+1;
     imap.(['CMAG' num2str(k)]) = {'corrMag', ['field' num2str(k)]}; %correlation
+    beam_vars{c} = ['CMAG' num2str(k)];
+end
+for k=1:num_beams
+    c=c+1;
     imap.(['PERG' num2str(k)]) = {'percentGood', ['field' num2str(k)]}; %percentGood
+    beam_vars{c} = ['PERG' num2str(k)];
 end
 
+ts_vars = {};
 if sensors.Temperature
     ts_vars = [ts_vars, 'TEMP'];
     imap.('TEMP') = {'variableLeader', 'temperature'};

Parser/GenericParser/InstrumentParsers/StaroddiParser.m

@@ -26,7 +26,7 @@
 %
 % Example:
 %
-% idata = StaroddiParser({'0-4T3769.DAT'},'timeSeries');
+% %see test/Parser/testStaroddi.m
 %
 % author: [email protected]
 %