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
........  0 canceled by user
  %       ........ -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
    % or no TSG measurement at sampling time
    % -------------------------------------------------------------
    ind2 = find(~isnan(tsg.EXT_DIF(ind)));
    
    if ~isempty(ind2)
      
      % Compute mean and standard deviation of the TSG/SAMPLE difference
      % that are suggested as default value for bias and error
      % ----------------------------------------------------------------
      if 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 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 length(ind2) == 1
        
        meanDif = tsg.EXT_DIF(ind(ind2));
        stdDif = -1;
        
      end
      
      defaultValue = {num2str(meanDif),num2str(stdDif)};
      
    else
      
      defaultValue = {'0','0'};
      
    end
    
  else
    
    defaultValue = {'0','0'};
    
  end
  
  
  % Enter the bias that will be applied to PARA{1}
  % ----------------------------------------------
  prompt = {['Constant value to be applied to the ' PARA{1} ' time series:'],...
    ['Error value to be applied to the ' PARA{1} ' time series:']};
  answer = inputdlg(prompt,'Bias adjustment',1,defaultValue);
  
  % if the user clic on cancel, answer is empty, do nothing
  % --------------------------------------------------------
  if ~isempty( answer )
    
    % everything OK
    % -------------
    error = 1;
    a = answer(1);
    b = answer(2);
    
    % 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
    
  else
    
    % canceled by user
    % ------------------
    error = 0;
  end
  
else
  
  % DateMax <= DateMin
  % ------------------
  error = -1;
  
end