corTsgBias.m

function [error] = corTsgBias(hMainFig, PARA, dateMin, dateMax)
% Correct the TSG time series with constant value, a bias.
% 
% Input
% hMainFig ..... Handle to the main GUI
% PARA ..........Cell array
%                   PARA{1} contains the characters (SSP, SSJT, SSTP)
%                   PARA{2} contains either the cahracters (SSPS, SSJT, SSTP)
%                           or (SSPS_CAL, SSJT_CAL, SSTP_CAL)
% dateMin ...... the correction is applied between dateMin and date Max
% dateMax ...... the correction is applied between dateMin and date Max
%
% Output
% Error ........  1 everything OK
%       ........ -1 dateMax <= date Min
%
% $Id$

% Get application data
% --------------------
tsg    = getappdata( hMainFig, 'tsg_data');
SAMPLE = tsg.plot.sample;

% -------------------------------------------------------------------------
% Get from the checkbox the QC code on which the correction will be applied
% -------------------------------------------------------------------------

% get list of keys from hashtable tsg.qc.hash, defined inside
% tsg_initialisation.m
% -----------------------------------------------------------
qc_list = keys(tsg.qc.hash);

% TODO: define size of keptCode
% -----------------------------
%keptCode = zeros(numel(qc_list), 1);

% iterate (loop) on each key store inside hastable
% ------------------------------------------------
keptCode = [];

nKeptCode = 0;
for key = qc_list
    
    % get handle of checkbox
    % ----------------------
    hCb = findobj(hMainFig, 'tag', ['TAG_CHECK_CORRECTION_' char(key)]);
    
    if get( hCb, 'value' )
        nKeptCode = nKeptCode + 1;
        keptCode(nKeptCode) = tsg.qc.hash.(key).code;
    end
end

% Get PROBABLY_GOOD, PROBABLY_BAD and VALUE_CHANGED codes
% -------------------------------------------------------
PROBABLY_GOOD = tsg.qc.hash.PROBABLY_GOOD.code;
PROBABLY_BAD  = tsg.qc.hash.PROBABLY_BAD.code;
VALUE_CHANGED = tsg.qc.hash.VALUE_CHANGED.code;

% Intialisation
% 01/09/2009 : intialisation to NaN for real and 0 for byte (QC)
% BE CAREFUL:
% netcdf toolbox failed with assertion when we write NaN to ncbyte variable
% -------------------------------------------------------------------------
if isempty( tsg.([PARA{1} '_ADJUSTED_ERROR']) )
    tsg.([PARA{1} '_ADJUSTED'])       = NaN*ones(size(tsg.(PARA{1})));
    tsg.([PARA{1} '_ADJUSTED_QC'])    = zeros(size(tsg.([PARA{1} '_QC'])));
    tsg.([PARA{1} '_ADJUSTED_ERROR']) = NaN*ones(size(tsg.(PARA{1})));
end

if dateMax > dateMin
    
    
    % Find samples within TIME_WINDOWS with Good, probably Good, QC
    % -------------------------------------------------------------
    ind = find( tsg.DAYD_EXT    >= dateMin &  tsg.DAYD_EXT    <= dateMax &...
        tsg.([SAMPLE '_EXT_QC']) <= PROBABLY_GOOD);
    
    if ~isempty(ind)
        
        % detect NaN in sample.SSPS_DIF due to bad QC code for tsg.SSPS
        % -------------------------------------------------------------
        ind2 = find(~isnan(tsg.EXT_DIF(ind)));
        
        % Compute mean and standard deviation of the TSG/SAMPLE difference
        % that are suggested as default value for bias and error
        % ----------------------------------------------------------------
        if ~isempty(ind2) && length(ind2) > 2
            
            meanDif = mean(tsg.EXT_DIF(ind(ind2)));
            stdDif = std(tsg.EXT_DIF(ind(ind2)));
            
            % Case with 2 samples only: the suggested bias is the
            % mean TSG/SAMPLE difference and the suggested error is the standard
            % deviation of the TSG/SAMPLE difference, the latter negative to warn that
            % the correction is not done with a significant number of samples
            % --------------------------------------------------------------------------
        elseif ~isempty(ind2) && length(ind2) == 2
            
            meanDif = mean( tsg.EXT_DIF(ind(ind2)) );
            stdDif = -std( tsg.EXT_DIF(ind(ind2)) );
            
            % Case with 1 sample only: the suggested bias is the
            % TSG/SAMPLE difference and the suggested error is -1
            % --------------------------------------------------------------------------
        elseif ~isempty(ind2) && length(ind2) == 1
            
            meanDif = tsg.EXT_DIF(ind(ind2));
            stdDif = -1;
            
        end
        
        defaultValueBias = {num2str(meanDif)};
        defaultValueError = {num2str(stdDif)};
        
    else
        
        defaultValueBias = {'0'};
        defaultValueError = {'0'};
        
    end
    
    
    % Enter the bias that will be applied to PARA{1}
    % ----------------------------------------------
    prompt = ['Constant value to be applied to the ' PARA{1} ' time series'];
    a = inputdlg(prompt,'Bias Correction',1,defaultValueBias);
    prompt = ['Error value to be applied to the ' PARA{1} ' time series'];
    b = inputdlg(prompt,'Bias Error',1,defaultValueError);
    
    % everything OK
    % -------------
    error = 1;
    if ~isempty( a )
        
        % If necessary replace a comma by a point
        % ---------------------------------------
        bias = regexprep(a, ',', '.');
        biasError = regexprep(b, ',', '.');
        
        % If bias not a numeric, str2doublereturn a NaN
        % ------------------------------------
        bias = str2double( bias );
        biasError = str2double( biasError );
        
        if isnumeric( bias ) && ~isnan( bias)
            
            
            %       if bias ~= 0
            
            % The correction is applied to the TSG between dateMin and dateMax
            % only to measurements with keptCode Quality Codes
            % ------------------------------------------------------------------------
            dtTsgQCkept=[];
            for icode = 1 : length( keptCode )
                dtTsg = find( tsg.DAYD    >= dateMin  & tsg.DAYD <= dateMax &...
                    tsg.([PARA{1} '_QC']) == keptCode( icode ));
                
                if ~isempty( dtTsg )
                    
                    dtTsgQCkept=[dtTsgQCkept; dtTsg];
                    
                    % Compute the corrected value : orignal value + correction
                    % --------------------------------------------------------
                    tsg.([PARA{1} '_ADJUSTED'])(dtTsg) = tsg.(PARA{2})(dtTsg) + bias;
                    
                    % Attribute an error
                    % ------------------
                    tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsg) = biasError;
                    
                    % Transfer the QC
                    % ---------------
                    tsg.([PARA{1} '_ADJUSTED_QC'])(dtTsg) = tsg.([PARA{1} '_QC'])(dtTsg);
                end
                
            end
            
            % Update tsg application data
            % ---------------------------
            setappdata( hMainFig, 'tsg_data', tsg);
            
            % everything OK
            % -------------
            error = 1;
        end
    end
    
else
    
    % DateMax <= DateMin
    % ------------------
    error = -1;

end