Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
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] = tsg_read( filename);
% Lecture fichier bucket
% -----------------------
filename = 'F:\work\M_TsgQc\tsg_data\past0601.btl';
[bucket, error] = tsg_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] = ...
tsg_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 = tsg_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)));
% *************************************************************************
% Calcul des valeurs medianes de correction: achaque point de
% comparaison, on attribue la mediane des corrections dans une fenetre de
% 10 jours. Puis on interpole entre chaque mediane pour avoir la correction
% a appliquer en chaque point de donnee TSG
% Fenetre de temps = 10 jours pour le calcul des valeurs medianes de
% correction
% ------------------------------------------------------------------
COR_mediane = [0 0];
erreur_mediane = [0 0];
for i = 1:length(sample)
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
% ******************************************************************
% 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 );