Skip to content

Commit

Permalink
PWLF Time Drift option added
Browse files Browse the repository at this point in the history
  • Loading branch information
@treynaud
treynaud committed Apr 26, 2024
1 parent 41b9915 commit 560bc36
Show file tree
Hide file tree
Showing 12 changed files with 1,108 additions and 864 deletions.
17 changes: 14 additions & 3 deletions doxy_corr/DOXY_argo_write.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,11 +104,14 @@
% new repertories names to save final netcdf files, takes
% into account PSAT or DOXY correction in the
% directories names
% 21/03/2019 Marine GALLIAN, Altran Ouest,
% v4 21/03/2019 Marine GALLIAN, Altran Ouest,
% new repertories names to save final netcdf files, takes
% into account drift computed on WOA or NCEP in the
% directories names and if a constant correction was
% applied
% v5 26.04.2024 PWLF option added for Time Drift Gain Correction (T. Reynaud)
% kwrite = 1 ==> SLOPE(1) and DRIFT(1) all cases+PWLF linear segment 1
% kwrite = 2 ==> SLOPE(2) and DRIFT(2) PWLF linear segment 2


function [] = DOXY_argo_write(varargin)
Expand Down Expand Up @@ -146,8 +149,16 @@

% Updates scientific calib fields, history fields,
% data_state_indicator, data_mode, date_update.

[monoProf, Dim] = DOXY_update_fields(Work,monoProf,Dim,REF_ARGO,k);

% Modified by T.Reynaud 26.04.2024
diffday=monoProf.juld.data(1)-Work.launchdate;
kwrite=1;
if Work.drift_PWLF & diffday>=Work.PPOX_DRIFTCORR_SEG(2)
kwrite=2;
end
%[monoProf, Dim] = DOXY_update_fields(Work,monoProf,Dim,REF_ARGO,k); % Commented by TR 26.04.2024
[monoProf, Dim] = DOXY_update_fields(Work,monoProf,Dim,REF_ARGO,kwrite);

monoProf = fillValue_2_nan(monoProf);
monoProf.hereDoxyFull.name = 'hereDoxy';
monoProf.hereDoxyFull.dim = {'N_PROF'};
Expand Down
17 changes: 11 additions & 6 deletions doxy_corr/DOXY_corr_prepare.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,24 +101,25 @@
% HISTORY
% $created: 14/01/2015 $author: Emilie Brion, Altran Ouest
% $Revision: version $Date: $author:
% v2 26/01/2017 Emilie Brion, Altran Ouest
% v2 26/01/2017 Emilie Brion, Altran Ouest
% argo - ameliorations 2017 phase 1
% v2.1 09/11/2017 Emilie Brion, Altran Ouest
% v2.1 09/11/2017 Emilie Brion, Altran Ouest
% drift correction is now re-computed using
% polyval or feval, instead of a default case
% ycorrected = raw - a*daydiff => becomming
% ycorrected = raw - fitRegression and
% fitRegression = polyval(driftCoef,daydiff)
% or feval(function_fit,daydiff).
% 26/10/2018 Marine GALLIAN, Altran Ouest :
% v2.2 26/10/2018 Marine GALLIAN, Altran Ouest :
% Add the possibility to apply a
% constant drift from a certain day
% 25/03/2019 Marine GALLIAN, Altran Ouest
% v2.3 25/03/2019 Marine GALLIAN, Altran Ouest
% Fix unit probleme for applying drift :
% Convert PPOX data when drift is computed on
% WOA
% 05/02/2020 Add conversion from mumol/L to mumol/kg when
% v2.4 05/02/2020 Add conversion from mumol/L to mumol/kg when
% drift is computed on PPOX. V. Thierry
% v2.5 26.04.2024 PWLF option added for Time Drift Gain Correction (T. Reynaud)

function [CORR2] = DOXY_corr_prepare(CORR, Work, argoStruct, argoTrajWork)

Expand Down Expand Up @@ -197,7 +198,11 @@
O2DriftCoef = 'PPOX';
end
if Work.drift_spec == 0
fitRegression = polyval(Work.([O2DriftCoef '_DRIFTCORR_COEF']),double(dayjul));
if ~Work.drift_PWLF % Added T.Reynaud 18.04.2024
fitRegression = polyval(Work.([O2DriftCoef '_DRIFTCORR_COEF']),double(dayjul));
else
fitRegression=polyval_PWLF(Work.([O2DriftCoef '_DRIFTCORR_COEF']),Work.PPOX_DRIFTCORR_SEG,double(dayjul));
end
if iscolumn(fitRegression)
fitRegression = fitRegression';
end
Expand Down
102 changes: 81 additions & 21 deletions doxy_corr/DOXY_drift.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@
% extremum day.
% 11.04.2022 Bug corrected for Work.launchdate (T.Reynaud)
% 31.10.2023 Online questions versus Dialog Box choice added (T. Reynaud and C. Kermabon)

% 26.04.2024 PWLF option added for Time Drift Gain Correction (T. Reynaud)

function [Work, DRIFT] = DOXY_drift(Work, argoWork, argo, argoTrajWork)
% =========================================================================
Expand Down Expand Up @@ -221,7 +221,8 @@
% Nb of day after deployment
%daydiff = (argo.juld.data - argo.juld.data(1));
% Nb of day since launch date
daydiff = (argo.juld.data - Work.launchdate);
%daydiff = (argo.juld.data - Work.launchdate);% TR 08.04.2024 -commented
daydiff = (argo.juld.data - argo.launchdate);% TR 08.04.2024 - added



Expand Down Expand Up @@ -361,11 +362,20 @@
Work.drift_fitPolynomialDegree = str2num(char(inputval)); %#ok<ST2NM>
end
end
%Added by T.Reynaud 05.04.2024 PWLF Time Drift
if Work.drift_PWLF
num_seg=Work.drift_PWLF_N;
else
num_seg=1;
end

if Work.drift_spec == 0
coeffFit = NaN(Work.drift_fitPolynomialDegree+1,length(Zq));
coeffFit = NaN(num_seg,Work.drift_fitPolynomialDegree+1,length(Zq));
else
coeffFit = NaN(numcoeffs(Work.drift_fittype),length(Zq));
coeffFit = NaN(num_seg,numcoeffs(Work.drift_fittype),length(Zq));
end
%

cpt = 0;
noDriftComputation = false;

Expand All @@ -377,18 +387,33 @@
isok = ~isnan(daydiff) & ~isnan(DRIFT.gainppox(:,z));
if sum(isok) >= nbr_lin_reg && any(~isnan(DRIFT.gainppox(:,z)))
if Work.drift_spec == 0
fitResult = polyfit(double(daydiff(isok)),double(DRIFT.gainppox(isok,z)),Work.drift_fitPolynomialDegree);
coeffFit(:,z) = fitResult;
if ~Work.drift_PWLF % Added T.Reynaud 05.04.2024
fitResult = polyfit(double(daydiff(isok)),double(DRIFT.gainppox(isok,z)),Work.drift_fitPolynomialDegree);
coeffFit(:,:,z) = fitResult;
else
x=double(daydiff(isok));
y=double(DRIFT.gainppox(isok,z));
if ~exist('XI','var')
[XI,YI] = ginput(1);
XI=[0,XI,x(end)];
end
fitResult_PWLF=polyfit_PWLF(x,y,XI,Work.drift_PWLF_N);
coeffFit(:,:,z)=fitResult_PWLF;
end
else
fitResult = fit(double(daydiff(isok)),double(DRIFT.gainppox(isok,z)),Work.drift_fittype);
coeffFit(:,z) = coeffvalues(fitResult);
coeffFit(num_seg,:,z) = coeffvalues(fitResult);
end

% control plot : check the coefficient and the linear regression at
% a local depth
if Work.makePlot
if Work.drift_spec == 0
DRIFT.poliv_c = polyval(coeffFit(:,z), double(daydiff));
if ~Work.drift_PWLF % Added T.Reynaud 05.04.2024
DRIFT.poliv_c = polyval(squeeze(coeffFit(num_seg,:,z)), double(daydiff));
else
DRIFT.poliv_c = polyval_PWLF(squeeze(coeffFit(:,:,z)),XI,double(daydiff));
end
else
DRIFT.poliv_c = fitResult(double(daydiff));
end
Expand Down Expand Up @@ -416,8 +441,17 @@
% else
if sum(isok) >= nbr_lin_reg %|| any(~cellfun(@isnan,DRIFT.gainppox,'UniformOutput',false))
if ~Work.drift_spec
fitResult = polyfit(daydiff(isok),data(isok),Work.drift_fitPolynomialDegree);
coeffFit(:,z) = fitResult;
if ~Work.drift_PWLF % Added T.Reynaud 04.04.2024
fitResult = polyfit(daydiff(isok),data(isok),Work.drift_fitPolynomialDegree);
coeffFit(:,:,z) = fitResult;
else
x=daydiff(isok);
y=data(isok);
[XI,YI] = ginput(1);
XI=[0,XI,x(end)];
fitResult_PWLF=polyfit_PWLF(x,y,XI,Work.drift_PWLF_N);
coeffFit(:,:,z)=fitResult_PWLF;
end
else
% fit doesn't accept row vectors.
if ~iscolumn(daydiff)
Expand All @@ -427,13 +461,17 @@
data = data';
end
fitResult = fit(daydiff(isok),data(isok),Work.drift_fittype);
coeffFit(:,z) = coeffvalues(fitResult);
coeffFit(:,:,z) = coeffvalues(fitResult);
end
% control plot : check the coefficient and the linear regression at
% a local depth
if Work.makePlot
if Work.drift_spec == 0
DRIFT.poliv_c = polyval(coeffFit(:,z), double(daydiff));
if ~Work.drift_PWLF % Added T.Reynaud 04.04.2024
DRIFT.poliv_c = polyval(coeffFit(:,z), double(daydiff));%Pas de isok??
else
DRIFT.poliv_c = polyval_PWLF(squeeze(coeffFit(:,:,z)),XI,double(daydiff));
end
else
DRIFT.poliv_c = fitResult(double(daydiff));
end
Expand All @@ -457,8 +495,9 @@
% -------------------------------------------------------------------------
% Slope averaging
% -------------------------------------------------------------------------
coeffFit_mean = nanmean(coeffFit,2);
DRIFT.coeffFit_mean=coeffFit_mean;
% Added T.Reynaud 05.04.2024
coeffFit_mean = nanmean(coeffFit,3);
DRIFT.coeffFit_mean=coeffFit_mean;% For all Time Drift

if Work.drift_fitPolynomialDegree == 2
first_derivative = polyder(DRIFT.coeffFit_mean);%compute first derivative
Expand All @@ -474,7 +513,11 @@
Work.savePlot = 1;
end
if Work.drift_spec == 0
DRIFT.fitRegression = polyval(coeffFit_mean,double(daydiff));
if ~Work.drift_PWLF % Added T.Reynaud 05.04.2024
DRIFT.fitRegression = polyval(coeffFit_mean,double(daydiff));%For all Time Drift
else
DRIFT.fitRegression = polyval_PWLF(coeffFit_mean,XI,double(daydiff));%For INAIR ONLY
end
else
listCoef = [];
for i=1:length(coeffFit_mean)
Expand All @@ -495,8 +538,14 @@
daydiff = (argo.juld.data - argo.launchdate);

if Work.drift_spec == 0
DRIFT.fitRegression = polyval(coeffFit_mean,double(daydiff));
DRIFT.daydiff=daydiff;% Added by T. Reynaud 16/07/2021
if ~Work.drift_PWLF % Added T.Reynaud 05.04.2024
DRIFT.fitRegression = polyval(coeffFit_mean,double(daydiff));
DRIFT.daydiff=daydiff;% Added by T. Reynaud 16/07/2021
else
DRIFT.fitRegression = polyval_PWLF(coeffFit_mean,XI,double(daydiff)); % Added T.Reynaud 08.04.2024
DRIFT.daydiff=daydiff;
DRIFT.daydiff_Seg=XI;
end
else
listCoef = [];
for i=1:length(coeffFit_mean)
Expand All @@ -512,8 +561,14 @@
end
else
if Work.drift_spec == 0
DRIFT.fitRegression = polyval(coeffFit_mean,double(daydiff));
DRIFT.daydiff=daydiff;% Added by T. Reynaud 16/07/2021
if ~Work.drift_PWLF % Added T.Reynaud 17.04.2024
DRIFT.fitRegression = polyval(coeffFit_mean,double(daydiff));
DRIFT.daydiff=daydiff;
else
DRIFT.fitRegression = polyval_PWLF(coeffFit_mean,XI,double(daydiff));
DRIFT.daydiff=daydiff;
DRIFT.daydiff_Seg=XI;
end
else
listCoef = [];
for i=1:length(coeffFit_mean)
Expand Down Expand Up @@ -543,8 +598,13 @@
end
%Recompute drift
if Work.drift_spec == 0
DRIFT.fitRegression_airwater = polyval(coeffFit_mean,double(daydiff_airwater));
DRIFT.daydiff_airwater=daydiff_airwater;%Added by T. Reynaud 16/07/2021
if ~Work.drift_PWLF % Added T.Reynaud 08.04.2024
DRIFT.fitRegression_airwater = polyval(coeffFit_mean,double(daydiff_airwater));
DRIFT.daydiff_airwater=daydiff_airwater;%Added by T. Reynaud 16/07/2021
else
DRIFT.fitRegression_airwater = polyval_PWLF(coeffFit_mean,XI,double(daydiff_airwater)); % Added T.Reynaud 08.04.2024
DRIFT.daydiff_airwater=daydiff_airwater;% Added T.Reynaud 08.04.2024
end
else
listCoef_airwater = [];
for i=1:length(coeffFit_mean)
Expand Down
17 changes: 16 additions & 1 deletion doxy_corr/DOXY_drift_apply.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,8 @@
% from Work
% v3.5 09.07.2021 For NCEP time drift the constant time drift
% index recalculated
% v3.6 26.04.2024 PWLF option added for Time Drift Gain Correction (T. Reynaud)


function [Work, DRIFT] = DOXY_drift_apply(Work, argoWork, argo, argoTrajWork, DRIFT)

Expand Down Expand Up @@ -195,6 +197,15 @@

Work.PPOX_DRIFTCORR_COEF = DRIFT.coeffFit_mean;
Work.PPOX_DRIFTCORR = DRIFT.ppox_corrected;
Work.PPOX_fitRegression=DRIFT.fitRegression;
Work.PPOX_daydiff=DRIFT.daydiff;

if Work.drift_PWLF % Added T.Reynaud 16.04.2024
Work.PPOX_DRIFTCORR_SEG=DRIFT.daydiff_Seg;
num_seg=Work.drift_PWLF_N;
else
num_seg=1;
end

if strcmp(Work.whichCorr,'INAIR')

Expand All @@ -211,7 +222,11 @@
end
%Recompute drift
if Work.drift_spec == 0
DRIFT.fitRegression_airwater = polyval(DRIFT.coeffFit_mean,double(daydiff_airwater));
if ~Work.drift_PWLF % Added T.Reynaud 16.04.2024
DRIFT.fitRegression_airwater = polyval(DRIFT.coeffFit_mean,double(daydiff_airwater));
else
DRIFT.fitRegression_airwater = polyval_PWLF(DRIFT.coeffFit_mean,DRIFT.daydiff_Seg,double(daydiff_airwater));
end
else
listCoef_airwater = [];
for i=1:length(DRIFT.coeffFit_mean)
Expand Down
Loading

0 comments on commit 560bc36

Please sign in to comment.