Merge pull request #483 from aodn/2.5.36

2.5.36

AutomaticQC/imosDensityInversionSetQC.m

@@ -0,0 +1,179 @@
+function [sample_data, varChecked, paramsLog] = imosDensityInversionSetQC( sample_data, auto )
+%DENSITYINVERSIONSETQC Flags any PSAL + TEMP + CNDC + DENS value that shows an inversion in density 
+% or an increase/decrease in density that is greater than a threshold set in 
+% imosDensityInversionQC.txt.
+%
+% Density inversion test on profiles data from ARGO which finds and flags 
+% good any data which value Vn passes the tests:
+% - density difference when going down <= threshold
+% - density difference when going up >= -threshold
+%
+% Inputs:
+%   sample_data - struct containing the data set.
+%
+%   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 <[email protected]>
+%
+
+%
+% Copyright (C) 2017, Australian Ocean Data Network (AODN) and Integrated 
+% Marine Observing System (IMOS).
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation version 3 of the License.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+
+% You should have received a copy of the GNU General Public License
+% along with this program.
+% If not, see <https://www.gnu.org/licenses/gpl-3.0.en.html>.
+%
+
+narginchk(1, 2);
+if ~isstruct(sample_data),              error('sample_data must be a struct');      end
+
+% auto logical in input to enable running under batch processing
+if nargin<5, auto=false; end
+
+varChecked = {};
+paramsLog  = [];
+
+% this test only applies to profile mode
+mode = readProperty('toolbox.mode');
+if ~strcmpi(mode, 'profile')
+    return;
+end
+
+% this test only applies on PSAL when TEMP and pressure data is available
+tempIdx = getVar(sample_data.variables, 'TEMP');
+psalIdx = getVar(sample_data.variables, 'PSAL');
+
+[presRel, zName, ~] = getPresRelForGSW(sample_data);
+[lat, lon] = getLatLonForGSW(sample_data);
+
+if ~(psalIdx && tempIdx && ~isempty(zName) && ~isempty(lat) && ~isempty(lon)), return; end
+
+temp = sample_data.variables{tempIdx}.data;
+psal = sample_data.variables{psalIdx}.data;
+
+% read threshold value from properties file
+threshold = str2double(readProperty('threshold', fullfile('AutomaticQC', 'imosDensityInversionSetQC.txt')));
+
+% read dataset QC parameters if exist and override previous parameters file
+currentQCtest = mfilename;
+threshold = readDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'threshold', 0.03); % default threshold value of 0.03kg/m3
+
+qcSet = str2double(readProperty('toolbox.qc_set'));
+rawFlag         = imosQCFlag('raw',         qcSet, 'flag');
+goodFlag        = imosQCFlag('good',        qcSet, 'flag');
+probGoodFlag    = imosQCFlag('probablyGood',qcSet, 'flag');
+probBadFlag     = imosQCFlag('probablyBad', qcSet, 'flag');
+badFlag         = imosQCFlag('bad',         qcSet, 'flag');
+
+badFlags = [probBadFlag, badFlag];
+
+paramsLog = ['threshold=' num2str(threshold)];
+
+% matrix case, we unfold the matrix in one vector for profile study
+% purpose
+isMatrix = size(temp, 1)>1 & size(temp, 2)>1;
+if isMatrix
+    len1 = size(temp, 1);
+    len2 = size(temp, 2);
+    len3 = size(temp, 3);
+    temp = temp(:);
+end
+
+lenData = length(temp);
+
+% we don't consider already bad data in the current test
+tempFlags = sample_data.variables{tempIdx}.flags;
+psalFlags = sample_data.variables{psalIdx}.flags;
+
+zIdx = getVar(sample_data.variables, zName);
+if zIdx == 0
+    % we either have DEPTH as a dimension or we have NOMINAL_DEPTH: in both
+    % cases we don't have any flag information, then it is safe to duplicate TEMP flags
+    % for example.
+    presFlags = tempFlags;
+else
+    presFlags = sample_data.variables{zIdx}.flags;
+end
+
+iBadData = ismember(tempFlags, badFlags) | ismember(psalFlags, badFlags) | ismember(presFlags, badFlags);
+
+% we need at least 2 valid measurements to perform this test
+if sum(~iBadData) <= 1, return; end
+
+temp    = temp(~iBadData);
+psal    = psal(~iBadData);
+presRel = presRel(~iBadData);
+
+lenDataTested = length(temp);
+
+flags = ones(lenData, 1, 'int8')*rawFlag;
+flagsTested = ones(lenDataTested, 1, 'int8')*rawFlag;
+
+I = true(lenDataTested, 1);
+I(end) = false;
+
+Ip1 = [false; I(1:end-1)];
+
+presRelRef = (presRel(I) + presRel(Ip1))/2;
+
+SA = gsw_SA_from_SP(psal, presRel, lon, lat);
+CT = gsw_CT_from_t(SA, temp, presRel);
+rhoPotI   = gsw_rho(SA(I), CT(I), presRelRef); % potential density referenced to the mid-point pressure
+rhoPotIp1 = gsw_rho(SA(Ip1), CT(Ip1), presRelRef);
+
+deltaUpBottom = [rhoPotI - rhoPotIp1; NaN];
+deltaBottomUp = [NaN; rhoPotIp1 - rhoPotI];
+
+threshold = ones(lenDataTested, 1) .* threshold;
+
+iDensInv = deltaUpBottom >= threshold | deltaBottomUp <= -threshold;
+
+flagsTested(iDensInv) = badFlag;
+flagsTested(~iDensInv) = goodFlag;
+flags(~iBadData) = flagsTested;
+
+if isMatrix
+    % we fold the vector back into a matrix
+    flags = reshape(flags, [len1, len2, len3]);
+end
+
+sample_data.variables{tempIdx}.flags = flags;
+sample_data.variables{psalIdx}.flags = flags;
+
+varChecked = {sample_data.variables{tempIdx}.name, ...
+    sample_data.variables{psalIdx}.name};
+
+densIdx = getVar(sample_data.variables, 'DENS');
+if densIdx
+    sample_data.variables{densIdx}.flags = flags;
+    varChecked = [varChecked, {sample_data.variables{densIdx}.name}];
+end
+
+cndcIdx = getVar(sample_data.variables, 'CNDC');
+if cndcIdx
+    if any(ismember(sample_data.variables{cndcIdx}.flags, [goodFlag, probGoodFlag]))
+        sample_data.variables{cndcIdx}.flags = flags;
+        varChecked = [varChecked, {sample_data.variables{cndcIdx}.name}];
+    end
+end
+
+% write/update dataset QC parameters
+writeDatasetParameter(sample_data.toolbox_input_file, currentQCtest, 'threshold', threshold);

AutomaticQC/imosDensityInversionSetQC.txt

@@ -0,0 +1 @@
+threshold = 0.03

AutomaticQC/imosVerticalSpikeQC.m

@@ -59,6 +59,12 @@
 paramsLog = [];
 flags     = [];
 
+% this test only applies to profile mode
+mode = readProperty('toolbox.mode');
+if ~strcmpi(mode, 'profile')
+    return;
+end
+
 % read all values from imosSpikeQC properties file
 values = readProperty('*', fullfile('AutomaticQC', 'imosVerticalSpikeQC.txt'));
 
@@ -91,6 +97,7 @@
     passFlag = imosQCFlag('good', qcSet, 'flag');
     failFlag = imosQCFlag('probablyBad',  qcSet, 'flag');
     badFlag  = imosQCFlag('bad',  qcSet, 'flag');
+    badFlags = [failFlag, badFlag];
 
     paramsLog = ['threshold=' thresholdExpr{iParam}];
 
@@ -105,7 +112,7 @@
     end
 
     % we don't consider already bad data in the current test
-    iBadData = sample_data.variables{k}.flags == badFlag;
+    iBadData = ismember(sample_data.variables{k}.flags, badFlags);
     dataTested = data(~iBadData);
 
     if isempty(dataTested), return; end

Graph/DepthProfile/highlightDepthProfileGeneric.m

@@ -19,7 +19,8 @@
 %               data. If no data points lie within the highlight region, 
 %               an empty matrix is returned.
 %
-% Author: Paul McCarthy <[email protected]>
+% Author:       Paul McCarthy <[email protected]>
+% Contributor:  Guillaume Galibert <[email protected]>
 %
 
 %
@@ -39,4 +40,48 @@
 % along with this program.
 % If not, see <https://www.gnu.org/licenses/gpl-3.0.en.html>.
 %
-highlight = highlightTimeSeriesGeneric(region, data, variable, type);
+narginchk(4, 4);
+
+if ~isnumeric(region) || ~isvector(region) || length(region) ~= 4
+  error('region must be a numeric vector of length 4');
+end
+
+if ~ishandle(data),     error('data must be a graphics handle'); end
+if ~isstruct(variable), error('variable must be a struct');      end
+if ~ischar(type),       error('type must be a string');          end
+
+xdata = get(data(1), 'XData'); % data(1) retrieves the current graphic handles only, in case extra sample is selected
+ydata = get(data(1), 'YData');
+
+if iscell(xdata)
+   xdata = cell2mat(xdata)'; 
+   ydata = cell2mat(ydata)';
+end
+
+if strcmp(type, 'alt') % right click
+    % figure out indices of all data points outside the X range but within
+    % the Y range
+    xidx  = (xdata <  region(1)) | (xdata >  region(3));
+    yidx  = (ydata >= region(2)) & (ydata <= region(4));
+else
+    % figure out indices of all data points within the X and Y range
+    xidx  = (xdata >= region(1)) & (xdata <= region(3));
+    yidx  = (ydata >= region(2)) & (ydata <= region(4));
+end
+
+% figure out indices of all the points to be highlighted
+idx = xidx & yidx;
+
+if ~any(idx)
+    % return nothing if no points to plot
+    highlight = [];
+else
+    % create the highlight
+    highlight = line(xdata(idx),ydata(idx), ...
+        'UserData',        idx, ...
+        'Parent',          gca, ...
+        'LineStyle',       'none', ...
+        'Marker',          'o', ...
+        'MarkerEdgeColor', 'white', ...
+        'MarkerFaceColor', 'white');
+end

Graph/TimeSeries/highlightTimeSeriesGeneric.m

@@ -55,20 +55,16 @@
    ydata = cell2mat(ydata)';
 end
 
-% on right click highlight, only highlight 
-% unflagged data points in the region
-if strcmp(type, 'alt')    
-    f = variable.flags;
-    f = f == 0;
-
-    xdata = xdata(f);
-    ydata = ydata(f);
+if strcmp(type, 'alt') % right click
+    % figure out indices of all data points within the X range but outside the Y range
+    xidx  = (xdata >= region(1)) & (xdata <= region(3));
+    yidx  = (ydata <  region(2)) | (ydata >  region(4));
+else
+    % figure out indices of all data points within the X and Y range
+    xidx  = (xdata >= region(1)) & (xdata <= region(3));
+    yidx  = (ydata >= region(2)) & (ydata <= region(4));
 end
 
-% figure out indices of all data points within the range
-xidx  = (xdata >= region(1) & xdata <= region(3));
-yidx  = (ydata >= region(2) & ydata <= region(4));
-
 % figure out indices of all the points to be highlighted
 idx = xidx & yidx;
 

Graph/TimeSeries/highlightTimeSeriesTimeDepth.m

@@ -61,21 +61,18 @@
 Ydata = repmat(ydata', [size(xdata), 1]);
 clear xdata ydata;
 
-% figure out indices of all data points within the range
-X = ((Xdata >= region(1)) & (Xdata <= region(3)));
-Y = ((Ydata >= region(2)) & (Ydata <= region(4)));
+if strcmp(type, 'alt') % right click
+    % figure out indices of all data points within the X range but outside the Y range
+    X = (Xdata >= region(1)) & (Xdata <= region(3));
+    Y = (Ydata <  region(2)) | (Ydata >  region(4));
+else
+    % figure out indices of all data points within the X and Y range
+    X = (Xdata >= region(1)) & (Xdata <= region(3));
+    Y = (Ydata >= region(2)) & (Ydata <= region(4));
+end
 idx = X & Y;
 clear X Y;
 
-% on right click, only highlight unflagged points
-if strcmp(type, 'alt')
-
-  % see if any points in the region haven't been flagged
-  idx = (variable.flags == 0) & idx;
-  if ~any(any(variable.flags)), return; end
-
-end
-
 if ~any(any(idx))
     highlight = [];
 else

Graph/TimeSeries/highlightTimeSeriesTimeFrequencyDirection.m

@@ -50,21 +50,18 @@
 Xdata = get(data, 'XData');
 Ydata = get(data, 'YData');
 
-% figure out indices of all data points within the range
-X = ((Xdata >= region(1)) & (Xdata <= region(3)));
-Y = ((Ydata >= region(2)) & (Ydata <= region(4)));
+if strcmp(type, 'alt') % right click
+    % figure out indices of all data points outside the X and Y range
+    X = ((Xdata < region(1)) | (Xdata > region(3)));
+    Y = ((Ydata < region(2)) | (Ydata > region(4)));
+else
+    % figure out indices of all data points within the X and Y range
+    X = (Xdata >= region(1)) & (Xdata <= region(3));
+    Y = (Ydata >= region(2)) & (Ydata <= region(4));
+end
 idx = X & Y;
 clear X Y;
 
-% on right click, only highlight unflagged points
-if strcmp(type, 'alt')
-
-  % see if any points in the region haven't been flagged
-  idx = (variable.flags == 0) & idx;
-  if ~any(any(variable.flags)), return; end
-
-end
-
 if ~any(any(idx))
     highlight = [];
 else

Graph/TimeSeries/parameters.txt

@@ -35,6 +35,7 @@ NORTEK_QC, TimeDepth
 VEL1,      TimeDepth
 VEL2,      TimeDepth
 VEL3,      TimeDepth
+VEL4,      TimeDepth
 VDEN,      TimeFrequency
 VDEV,      TimeFrequency
 VDEP,      TimeFrequency

Graph/TimeSeries/setTimeSerieColorbarContextMenu.m

@@ -51,7 +51,7 @@
 end
 
 switch upper(var.name(1:4))
-    case {'UCUR', 'VCUR', 'WCUR', 'ECUR', 'VEL1', 'VEL2', 'VEL3'} % 0 centred parameters
+    case {'UCUR', 'VCUR', 'WCUR', 'ECUR', 'VEL1', 'VEL2', 'VEL3', 'VEL4'} % 0 centred parameters
         colormap(r_b);
         cbCLimRange('', '', 'auto, 0 centred', var.data);
 

Graph/diagramMooring1DVarAgainstOther.m

@@ -141,7 +141,7 @@ function diagramMooring1DVarAgainstOther(sample_data, varName, yAxisVarName, isQ
     if iVar && iYAxisVar && ...
             size(sample_data{iSort(i)}.variables{iVar}.data, 2) == 1 && ...
             size(sample_data{iSort(i)}.variables{iYAxisVar}.data, 2) == 1 && ... % we're only plotting 1D variables with depth variable but no current
-            all(~strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3'}, 4))
+            all(~strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3', 'VEL4'}, 4))
         iGood = true(size(sample_data{iSort(i)}.variables{iVar}.data));
         if isQC
             %get time, depth and var QC information

Graph/diagramMooring2DVarAgainstOther.m

@@ -184,7 +184,7 @@ function diagramMooring2DVarAgainstOther(sample_data, varName, yAxisVarName, isQ
         end
 
     elseif iVar && iYAxisVar && ...
-            any(strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3'}, 4)) && ...
+            any(strncmpi(sample_data{iSort(i)}.variables{iVar}.name, {'UCUR', 'VCUR', 'WCUR', 'CDIR', 'CSPD', 'VEL1', 'VEL2', 'VEL3', 'VEL4'}, 4)) && ...
             size(sample_data{iSort(i)}.variables{iVar}.data, 2) == 1 % we're plotting current metre 1D variables with DEPTH variable.
         iGood = true(size(sample_data{iSort(i)}.variables{iVar}.data));
         if isQC