dev_correctBucket.m

clc
clear
close all

% Get TSG and bucket data
% -----------------------
% tsg    = getappdata( hMainFig, 'tsg_data');
% bucket = getappdata( hMainFig, 'bucket_data');

% Lecture fichier TSG
% -------------------
filename = 'F:\work\M_TsgQc\tsg_data\past0601.txt';
[tsg, error] = dev_readTsg( filename);

% Lecture fichier bucket
% -----------------------
filename = 'F:\work\M_TsgQc\tsg_data\past0601.btl';
[bucket, error] = dev_readBucket( filename);

% Salinity in one 1 hour interval, should not depart the average for more
% than SAL_STD_LIM standard deviation
% -----------------------------------------------------------------------
SAL_STD_LIM = 0.1;

% 1 hour interval expressed in MATLAB serial Number
% -------------------------------------------------
INTERVAL_SMOOTHING = datenum(0, 0, 0, 1, 0 , 0);

% dt between 2 tsg measurements
% -----------------------------
TSG_SAMPLING_TIME = datenum(0, 0, 0, 0, 5 , 0);
% ou plus general
TSG_SAMPLING_TIME = tsg.TIME(2) - tsg.TIME(1);

% Correction Windows in days
% --------------------------
TIME_WINDOWS = 10;


% *************************************************************************
%                             First Step
%
% Running average over INTERVAL_SMOOTHING hour.
% Salinity in one INTERVAL_SMOOTHING hour interval, should not depart the
% average for more than SAL_STD_LIM
%
% *************************************************************************

[psal_smooth, nval] = ...
      dev_moveaverage(tsg.TIME, tsg.PSAL, INTERVAL_SMOOTHING, SAL_STD_LIM);

  
% In the future we could use additional samples (i.e. CTD data)
% These data should be merged with the bucket measurement.
% All the samples would be stored in a structure 'sample'
% Function to be built
% ---------------------------------------------------------------------
sample.TIME        = bucket.TIME_WS;
sample.LATITUDE    = bucket.LATITUDE_WS;
sample.LONGITUDE   = bucket.LONGITUDE_WS;
sample.PSAL        = bucket.PSAL_WS;
sample.PSAL_QC     = bucket.PSAL_QC_WS;
sample.PSAL_DIF    = zeros(size(sample.PSAL));
sample.PSAL_SMOOTH = zeros(size(sample.PSAL));
  
% *************************************************************************
%                             Second Step
% 
% Co-location of samples and TSG measurements
% Compute the sample-TSG differences

sample = dev_diffTsgSample(tsg, psal_smooth, sample, TSG_SAMPLING_TIME);


% Exclusion des points aberrants (ecart trop grand)
% -------------------------------------------------
%A = find( sal_diff  < (median(sal_diff) + 3 * std(sal_diff)));
%B = find( sal_diff  > (median(sal_diff) - 3 * std(sal_diff)));
%C = find( sal_diff == (median(sal_diff) + 3 * std(sal_diff)));


% *************************************************************************

%COR_mediane    = [0 0];
%erreur_mediane = [0 0];

%for i = 1:length(sample.PSAL_DIF)
    
%    is = find( sample.TIME > sample.TIME(i) - TIME_WINDOWS/2 & ...
%               sample.TIME < sample.TIME(i) + TIME_WINDOWS/2);
           
%    if ~isempty(is)
%        COR_mediane    = [COR_mediane;...
%                         [ sample.TIME(i) median(sample.PSAL_DIF(is))]];
%        erreur_mediane = [erreur_mediane;...
%                         [std(sample.PSAL_DIF(is))/sqrt(length(is)) length(is)]];
%    end
%end

%COR_mediane    = COR_mediane(2:length(COR_mediane),:);
%erreur_mediane = erreur_mediane(2:length(erreur_mediane),:);

%erreur_mediane( find(erreur_mediane(:,2) < 4), 1) = 1;

% Complete les series de correction PAS TRES BON
%if COR_mediane(1,1) ~= tsg.TIME(1)
%    COR_mediane = [[tsg.TIME(1) COR_mediane(1,2)];COR_mediane];
%end

% Pas trs bon : La correction ne doit pas tre applique + de 10 jours
% avant ou aprs le dernier echantillon
% --------------------------------------------------------------------
%if COR_mediane(length(COR_mediane),1) ~= tsg.TIME(end)
%    COR_mediane = [COR_mediane; ...
%                  [tsg.TIME(end) ...
%                  COR_mediane(length(COR_mediane),2)]];
%end

%COR = interp1(COR_mediane(:,1), COR_mediane(:,2), tsg.TIME);

%tsg.PSAL_ADJ = tsg.PSAL + COR;   %salinite initiale corrige

dateMin = tsg.TIME(5000);
dateMax = tsg.TIME(8000);
tsg = dev_corMethod1(tsg, sample, dateMin, dateMax, TIME_WINDOWS);


% ******************************************************************
%                                   Trac
% ******************************************************************    
figure
plot( tsg.TIME, tsg.PSAL, 'k' );
hold on
plot( tsg.TIME, psal_smooth, 'r' );

figure
hist(nval, 1:13)

figure
plot( tsg.TIME, psal_smooth, 'k' );
hold on
ind = find( sample.PSAL_QC == 1 ); 
plot( sample.TIME(ind), sample.PSAL(ind), '.b' );

ind = find( sample.PSAL_QC == 0 ); 
if ~isempty(ind)
    plot( sample.TIME(ind), sample.PSAL(ind), '.r' );
end

figure
plot( sample.TIME, sample.PSAL_DIF, '.k' );

figure
plot( tsg.TIME, tsg.PSAL, 'k' );
hold on
plot( tsg.TIME, tsg.PSAL_ADJ, 'r' );
plot( sample.TIME, sample.PSAL, '.b', 'Markersize',  6 );