-
jacques.grelet_ird.fr authoredjacques.grelet_ird.fr authored
corTsgLinear.m 8.32 KiB
function [error] = corTsgLinear(hMainFig, PARA, dateMin, dateMax)
% Correct the TSG salinity time series with the Water sample.
% Use a linear fit to the water sample/tsg difference
%
% 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
%
% 2009/04/15 G. Alory
% If the correction is done with 2 samples only, the error is
% the Std.Dev. of the TSG-sample differences but negative
% The error minimum cannot be lower than the sensor accuracy
%
% 2009/03/23 Y. Gouriou
% 1) Compute linear fit with a minimum of 3 points.
% 2) If only 2 points: no linear fit, but compute the average difference
% with the 2 points.
%
% $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 linear fit of the TSG/SAMPLE difference
% -----------------------------------------------
if ~isempty(ind2) && length(ind2) > 2
% if ~isempty(tsg.EXT_DIF(ind(ind2))) Is this line useful?
% Linear fit applied to the difference tsg-sample
% -----------------------------------------------
X = tsg.DAYD_EXT(ind(ind2));
Y = tsg.EXT_DIF(ind(ind2));
[p, S, mu] = polyfit( X, Y, 1);
% The correction is applied to the TSG between dateMin and dateMax using
% a linear interpolation only on 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 correction + the error
% ----------------------------------
[tsg.([PARA{1} '_ADJUSTED'])(dtTsg),...
tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsg)] =...
polyval( p, tsg.DAYD(dtTsg), S, mu);
% Compute the corrected value : orignal value + correction
% --------------------------------------------------------
tsg.([PARA{1} '_ADJUSTED'])(dtTsg) =...
tsg.(PARA{2})(dtTsg) + tsg.([PARA{1} '_ADJUSTED'])(dtTsg);
% Transfer the QC
% ---------------
tsg.([PARA{1} '_ADJUSTED_QC'])(dtTsg) = tsg.([PARA{1} '_QC'])(dtTsg);
end
end
% Case with 2 samples only: the TSG time series is shifted by the
% mean TSG/SAMPLE difference and the error is computed as the standard
% deviation of the TSG/SAMPLE difference,but 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)) );
% The correction is applied to the TSG between dateMin and dateMax using
% a linear interpolation only on 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];
tsg.([PARA{1} '_ADJUSTED'])(dtTsg) = tsg.(PARA{2})(dtTsg) + meanDif;
tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsg) = -stdDif;
tsg.([PARA{1} '_ADJUSTED_QC'])(dtTsg) = tsg.([PARA{1} '_QC'])(dtTsg);
end
end
% Case with 1 sample only: the TSG time series is shifted by the
% TSG/SAMPLE difference and the error is put to -1
% --------------------------------------------------------------------------
elseif ~isempty(ind2) && length(ind2) == 1
% The correction is applied to the TSG between dateMin and dateMax
% only on 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];
tsg.([PARA{1} '_ADJUSTED'])(dtTsg) = tsg.(PARA{2})(dtTsg) + tsg.EXT_DIF(ind(ind2));
tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsg) = -1;
tsg.([PARA{1} '_ADJUSTED_QC'])(dtTsg) = tsg.([PARA{1} '_QC'])(dtTsg);
end
end
end
% The error minimum cannot be lower than the sensor accuracy
% ----------------------------------------------------------
if ~isempty(ind2) && length(ind2) >= 2
accuracy= tsg_accuracy(hMainFig, PARA{1},dtTsgQCkept );
error_too_small=find(abs(tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept)) < accuracy);
if ~isempty(error_too_small)
if length(ind2) > 2
tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept(error_too_small))=accuracy(error_too_small);
else
tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept(error_too_small))=-accuracy(error_too_small);
end
end
end
end
% Update tsg application data
% ---------------------------
setappdata( hMainFig, 'tsg_data', tsg);
% everything OK
% -------------
error = 1;
else
% DateMax <= DateMin
% ------------------
error = -1;
end