1 |
% |
---|
2 |
% gliderOptode_Generate_L1_Data.m |
---|
3 |
% |
---|
4 |
% Purpose: Generate Level 1 data mat files that include optode DO data. |
---|
5 |
% |
---|
6 |
% NOTE: At top of file, set glider index and deployment number |
---|
7 |
% |
---|
8 |
% Requires: |
---|
9 |
% MATLAB folder - contains util, matutil, seawater, plots, strfun, opnml |
---|
10 |
% |
---|
11 |
% |
---|
12 |
% Author: William Stark |
---|
13 |
% Marine Sciences Department |
---|
14 |
% UNC-Chapel Hill |
---|
15 |
% |
---|
16 |
% Created: 21 May 2012 |
---|
17 |
% |
---|
18 |
%////////////////////////////////////////////////////////////////////////// |
---|
19 |
|
---|
20 |
clear all; |
---|
21 |
|
---|
22 |
|
---|
23 |
%!!!!!! SET PRIOR TO RUNNING CODE !!!!!!!!!!!!!!!!!!!!!!!! |
---|
24 |
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
25 |
% |
---|
26 |
% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ... |
---|
27 |
gliderIndex = 1; |
---|
28 |
|
---|
29 |
% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ... |
---|
30 |
deploymentNumber = 1; |
---|
31 |
|
---|
32 |
% |
---|
33 |
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
34 |
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
35 |
|
---|
36 |
|
---|
37 |
|
---|
38 |
% add paths for required files... |
---|
39 |
addpath('MATLAB/util/'); |
---|
40 |
addpath('MATLAB/matutil/'); |
---|
41 |
addpath('MATLAB/seawater/'); |
---|
42 |
addpath('MATLAB/plots/'); |
---|
43 |
addpath('MATLAB/strfun/'); |
---|
44 |
addpath('MATLAB/opnml/'); |
---|
45 |
addpath('MATLAB/opnml/FEM/'); |
---|
46 |
|
---|
47 |
|
---|
48 |
|
---|
49 |
% glider name string... |
---|
50 |
if (gliderIndex==1) |
---|
51 |
strGliderName = 'Pelagia'; |
---|
52 |
else |
---|
53 |
strGliderName = 'Ramses'; |
---|
54 |
end |
---|
55 |
|
---|
56 |
% populate arrays for the deployment start and end dates... |
---|
57 |
% ex. strStart(2, 3) is start date for Ramses, Deployment 3 |
---|
58 |
strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'}; |
---|
59 |
strEnd = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'}; |
---|
60 |
|
---|
61 |
% deployment number string... |
---|
62 |
strDeploymentNumber = num2str(deploymentNumber); |
---|
63 |
|
---|
64 |
% deployment start date string... |
---|
65 |
strStartDate = strStart(gliderIndex, deploymentNumber); |
---|
66 |
|
---|
67 |
% deployment end date string... |
---|
68 |
strEndDate = strEnd(gliderIndex, deploymentNumber); |
---|
69 |
|
---|
70 |
% define the path to the glider ascii files... |
---|
71 |
%datadir = strcat('/Users/haloboy/Documents/MASC/MATLAB/CTD_data_correction/GLIDER_CTD_DATA_LEVEL0/',... |
---|
72 |
datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, '_Deployment', strDeploymentNumber, '/'); |
---|
73 |
|
---|
74 |
|
---|
75 |
|
---|
76 |
% define default bounds for use in plots... |
---|
77 |
% switch gliderIndex |
---|
78 |
% case 1 % Pelagia |
---|
79 |
% switch deploymentNumber |
---|
80 |
% case 1 % Deployment 1 |
---|
81 |
% tempBounds = [17.0 24.0]; |
---|
82 |
% salinBounds = [36.0 36.4]; |
---|
83 |
% densBounds = [1025.0 1026.6]; |
---|
84 |
% chlorBounds = [0.0 4.0]; |
---|
85 |
% |
---|
86 |
% case 2 % Deployment 2 |
---|
87 |
% tempBounds = [17.0 24.0]; |
---|
88 |
% salinBounds = [36.0 36.5]; |
---|
89 |
% densBounds = [1024.5 1026.8]; |
---|
90 |
% chlorBounds = [0.0 4.0]; |
---|
91 |
% |
---|
92 |
% case 3 % Deployment 3 |
---|
93 |
% tempBounds = [17.0 24.0]; |
---|
94 |
% salinBounds = [35.9 36.7]; |
---|
95 |
% densBounds = [1024.4 1026.4]; |
---|
96 |
% chlorBounds = [0.0 4.0]; |
---|
97 |
% end |
---|
98 |
% case 2 % Ramses |
---|
99 |
% switch deploymentNumber |
---|
100 |
% case 1 % Deployment 1 |
---|
101 |
% tempBounds = [8.0 23.0]; |
---|
102 |
% salinBounds = [35.0 36.4]; |
---|
103 |
% densBounds = [1024.5 1027.5]; |
---|
104 |
% chlorBounds = [0.0 4.0]; |
---|
105 |
% |
---|
106 |
% case 2 % Deployment 2 |
---|
107 |
% tempBounds = [9.0 25.0]; |
---|
108 |
% salinBounds = [35.2 36.6]; |
---|
109 |
% densBounds = [1024.0 1027.5]; |
---|
110 |
% chlorBounds = [0.0 4.0]; |
---|
111 |
% |
---|
112 |
% case 3 % Deployment 3 |
---|
113 |
% tempBounds = [10.0 24.5]; |
---|
114 |
% salinBounds = [35.3 36.7]; |
---|
115 |
% densBounds = [1024.4 1027.4]; |
---|
116 |
% chlorBounds = [0.0 4.0]; |
---|
117 |
% end |
---|
118 |
% end |
---|
119 |
|
---|
120 |
|
---|
121 |
%########################################################################################## |
---|
122 |
|
---|
123 |
|
---|
124 |
|
---|
125 |
|
---|
126 |
|
---|
127 |
%*** READ IN EBD DATA ***************************************************** |
---|
128 |
% declare variables for storing data... |
---|
129 |
oxy_dphase=[]; |
---|
130 |
ptime_oxy=[]; |
---|
131 |
|
---|
132 |
% try to load all *.ebdasc files at once... |
---|
133 |
[files, Dstruct] = wilddir(datadir, '.ebdasc'); |
---|
134 |
nfile = size(files, 1); |
---|
135 |
|
---|
136 |
for i=1:nfile-1 |
---|
137 |
% protect against empty ebd file |
---|
138 |
if(Dstruct(i).bytes>0) |
---|
139 |
data = read_gliderasc2([datadir, files(i,:)]); |
---|
140 |
%data = read_gliderasc3([datadir, files(i,:)]); |
---|
141 |
|
---|
142 |
% if the number of values (in data.data) is less than the number of |
---|
143 |
% vars (in data.vars), this means that the data were not completely read |
---|
144 |
% in. To correct this, pad data.data with NaNs until its length |
---|
145 |
% equals that of data.vars... |
---|
146 |
if (length(data.data) < length(data.vars)) |
---|
147 |
data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post'); |
---|
148 |
end |
---|
149 |
|
---|
150 |
% populate variables with data... |
---|
151 |
if(~isempty(data.data)) |
---|
152 |
oxy_dphase = [oxy_dphase;data.data(:,strmatch('sci_oxy3835_wphase_dphase', data.vars, 'exact'))]; % oxygen (dphase) |
---|
153 |
ptime_oxy = [ptime_oxy;data.data(:,strmatch('sci_m_present_time', data.vars, 'exact'))]; % present time |
---|
154 |
%should read in other oxy values for comparison with |
---|
155 |
%re-calculated values |
---|
156 |
end |
---|
157 |
|
---|
158 |
data = []; |
---|
159 |
end |
---|
160 |
end |
---|
161 |
%************************************************************************** |
---|
162 |
|
---|
163 |
|
---|
164 |
% load glider data from existing mat file... |
---|
165 |
strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, '_Deployment', strDeploymentNumber, '_DataL1.mat'); |
---|
166 |
load(strMatFilePath); |
---|
167 |
|
---|
168 |
|
---|
169 |
% first, apply the sort() function to ptime_oxy to make sure it increases monotonically... |
---|
170 |
[Y,I] = sort(ptime_oxy); |
---|
171 |
ptime_oxy = Y; |
---|
172 |
oxy_dphase = oxy_dphase(I); |
---|
173 |
|
---|
174 |
|
---|
175 |
% remove nans from EBD data... |
---|
176 |
i = find(~isnan(oxy_dphase)); |
---|
177 |
oxy_dphase = oxy_dphase(i); ptime_oxy = ptime_oxy(i); |
---|
178 |
|
---|
179 |
% interpolate CTD to times of oxy_dphase |
---|
180 |
tempi = interp1(ptime,temp,ptime_oxy); |
---|
181 |
salini = interp1(ptime,salinCorrected,ptime_oxy); |
---|
182 |
depthi = interp1(ptime,depth,ptime_oxy); |
---|
183 |
|
---|
184 |
% first, implement temperature-dependent correction to DO concentration |
---|
185 |
% that utilizes dphase as the input (manual page 30), coefficients from cal |
---|
186 |
% sheets (a 5x4 matrix of values) and will be glider (and foil) dependent. |
---|
187 |
|
---|
188 |
% for ramses, 2 June 2010 calibration sheet |
---|
189 |
|
---|
190 |
C0coef=[4.53793e3 -1.62595e2 3.29574 -2.79285e-2]; |
---|
191 |
C1coef=[-2.50953e2 8.02322 -1.58398e-1 1.31141e-3]; |
---|
192 |
C2coef=[5.66417 -1.59647e-1 3.07910e-3 -2.46265e-5]; |
---|
193 |
C3coef=[-5.99449e-2 1.48326e-3 -2.82110e-5 2.15156e-7]; |
---|
194 |
C4coef=[2.43614e-4 -5.26759e-6 1.00064e-7 -7.14320e-10]; |
---|
195 |
|
---|
196 |
C0 = C0coef(1)+C0coef(2)*tempi+C0coef(3)*tempi.^2+C0coef(4)*tempi.^3; |
---|
197 |
C1 = C1coef(1)+C1coef(2)*tempi+C1coef(3)*tempi.^2+C1coef(4)*tempi.^3; |
---|
198 |
C2 = C2coef(1)+C2coef(2)*tempi+C2coef(3)*tempi.^2+C2coef(4)*tempi.^3; |
---|
199 |
C3 = C3coef(1)+C3coef(2)*tempi+C3coef(3)*tempi.^2+C3coef(4)*tempi.^3; |
---|
200 |
C4 = C4coef(1)+C4coef(2)*tempi+C4coef(3)*tempi.^2+C4coef(4)*tempi.^3; |
---|
201 |
|
---|
202 |
o2 = C0+C1.*oxy_dphase+C2.*oxy_dphase.^2+C3.*oxy_dphase.^3+C4.*oxy_dphase.^4; |
---|
203 |
|
---|
204 |
% second, implement the salinity correction to DO concentration (page 31) |
---|
205 |
% ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER |
---|
206 |
% SCALING IS REQUIRED |
---|
207 |
|
---|
208 |
B0=-6.24097e-3; |
---|
209 |
B1=-6.93498e-3; |
---|
210 |
B2=-6.90358e-3; |
---|
211 |
B3=-4.29155e-3; |
---|
212 |
C0=-3.11680e-7; |
---|
213 |
|
---|
214 |
Ts = log((298.15 - tempi) ./ (273.15 + tempi)); |
---|
215 |
|
---|
216 |
o2_salcorr = o2*exp(salini*(B0+B1*Ts+B2*Ts.*Ts+B3*Ts.*Ts.*Ts)+C0*salini.^2); |
---|
217 |
|
---|
218 |
% third, implement the pressure correction to DO concentration (page 32) |
---|
219 |
|
---|
220 |
o2_salprescorr = os_salcorr*(1 + 0.04*depthi/1000); |
---|
221 |
|
---|
222 |
% use polynomial to calculate DO satuations using the measured temp and |
---|
223 |
% sal (manual page 30) |
---|
224 |
|
---|
225 |
A0 = 2.00856; |
---|
226 |
A1 = 3.22400; |
---|
227 |
A2 = 3.99063; |
---|
228 |
A3 = 4.80299; |
---|
229 |
A4 = 0.0978188; |
---|
230 |
A5 = 1.71069; |
---|
231 |
B0 = -00624097; |
---|
232 |
B1 = -00693498; |
---|
233 |
B2 = -00690358; |
---|
234 |
B3 = -00429155; |
---|
235 |
C0 = -000000311680; |
---|
236 |
|
---|
237 |
% need interpolated temp, salinity, pressure at times of DO obs |
---|
238 |
|
---|
239 |
rslt = A0 + A1*Ts + A2*Ts.^2 + A3*Ts.^3 + A4*Ts.^4 + A5*Ts.^5 +... |
---|
240 |
salini.*(B0 + B1*Ts + B2*Ts.^2 + B3*Ts.^3) + C0*salini.^2; |
---|
241 |
|
---|
242 |
O2sol = exp(rslt); |
---|
243 |
|
---|
244 |
o2_sat = (o2_salprescorr * 2.2414) ./ O2sol; |
---|
245 |
|
---|
246 |
|
---|
247 |
% convert ptime into datenum style... |
---|
248 |
ptime_oxy_datenum = ptime_oxy/3600/24+datenum(1970, 1, 1, 0, 0, 0); |
---|
249 |
% ptimehrly = fix(ptime_datenum*24)/24; |
---|
250 |
% ptimedaily = fix(ptime_datenum); |
---|
251 |
% ptimedaily2 = unique(ptimedaily); |
---|
252 |
% ptimedaily2 = ptimedaily2(1:2:end); |
---|
253 |
|
---|
254 |
|
---|
255 |
|
---|
256 |
|
---|
257 |
lb = datenum('10-Feb-2012 12:22:00'); |
---|
258 |
ub = datenum('10-Feb-2012 12:23:00'); |
---|
259 |
%lb = datenum('22-Mar-2012 12:01:00'); |
---|
260 |
%ub = datenum('23-Mar-2012 12:02:00'); |
---|
261 |
|
---|
262 |
% figure('Position', [500,500,1000,1000]); |
---|
263 |
% subplot(311); |
---|
264 |
% ind = find(ptime_oxy_datenum >= lb & ptime_oxy_datenum <= ub); |
---|
265 |
% plot(ptime_oxy_datenum(ind), oxy_dphase(ind), 'b.'); |
---|
266 |
% datetick; |
---|
267 |
% xlabel('time', 'fontsize', 12); |
---|
268 |
% ylabel('DO (dphase)', 'fontsize', 12); |
---|
269 |
% subplot(312); |
---|
270 |
% ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub); |
---|
271 |
% plot(ptime_datenum(ind2), temp(ind2), 'r.'); |
---|
272 |
% datetick; |
---|
273 |
% xlabel('time', 'fontsize', 12); |
---|
274 |
% ylabel('temp', 'fontsize', 12); |
---|
275 |
% subplot(313); |
---|
276 |
% ind = find(ptime_oxy_datenum >= lb & ptime_oxy_datenum <= ub); |
---|
277 |
% ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub); |
---|
278 |
% plot(ptime_oxy_datenum(ind), ones(length(ind), 1), 'b.', ptime_datenum(ind2), ones(length(ind2), 1), 'r.'); |
---|
279 |
% datetick; |
---|
280 |
|
---|
281 |
|
---|
282 |
figure('Position', [500,500,1000,1000]); |
---|
283 |
subplot(211); |
---|
284 |
ind = find(ptime_datenum >= lb & ptime_datenum <= ub); |
---|
285 |
plot(ptime_datenum(ind), oxy_dphase(ind), 'b.', ptime_datenum(ind), oxy_dphase_adjusted(ind), 'r.'); |
---|
286 |
datetick; |
---|
287 |
xlabel('time', 'fontsize', 12); |
---|
288 |
ylabel('DO (dphase)', 'fontsize', 12); |
---|
289 |
subplot(212); |
---|
290 |
ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub); |
---|
291 |
plot(ptime_datenum(ind2), temp(ind2), 'r.'); |
---|
292 |
datetick; |
---|
293 |
xlabel('time', 'fontsize', 12); |
---|
294 |
ylabel('temp', 'fontsize', 12); |
---|
295 |
|
---|
296 |
|
---|
297 |
|
---|
298 |
|
---|
299 |
|
---|
300 |
|
---|
301 |
|
---|
302 |
|
---|
303 |
% interpolate DBD data to align with EBD data... |
---|
304 |
% horizontalVelocity = interp1(ptime1_dbd, horizontalVelocity, ptime); |
---|
305 |
% depth = interp1(ptime1_dbd, depth, ptime); |
---|
306 |
% pitch = interp1(ptime1_dbd, pitch, ptime); |
---|
307 |
% avgDepthRate = interp1(ptime1_dbd, avgDepthRate, ptime); |
---|
308 |
% glideAngle = interp1(ptime1_dbd, glideAngle, ptime); |
---|
309 |
% lat = interp1(ptime1_dbd, lat, ptime); |
---|
310 |
% lon = interp1(ptime1_dbd, lon, ptime); |
---|
311 |
|
---|
312 |
|
---|
313 |
|
---|
314 |
|
---|
315 |
|
---|
316 |
|
---|
317 |
|
---|
318 |
|
---|
319 |
|
---|
320 |
|
---|
321 |
|
---|
322 |
|
---|
323 |
% create configuration struct... |
---|
324 |
% units = struct('chlor', 'micrograms liter-1',... |
---|
325 |
% 'dens', 'kg m-3',... |
---|
326 |
% 'densCorrected', 'kg m-3',... |
---|
327 |
% 'depth', 'm',... |
---|
328 |
% 'gpsLat', 'decimal degrees',... |
---|
329 |
% 'gpsLon', 'decimal degrees',... |
---|
330 |
% 'pres', 'decibars',... |
---|
331 |
% 'ptime', 'seconds since 0000-01-01T00:00',... |
---|
332 |
% 'ptime_datenum', 'seconds since 0000-01-01T00:00',... |
---|
333 |
% 'salin', 'psu',... |
---|
334 |
% 'salinCorrected', 'psu',... |
---|
335 |
% 'temp', 'deg C'); |
---|
336 |
|
---|
337 |
% variable_description = struct('chlor', 'chlorophyll measured by glider',... |
---|
338 |
% 'chlorBounds', 'default chlorophyll bounds for plots',... |
---|
339 |
% 'dens', 'density measured by glider',... |
---|
340 |
% 'densBounds', 'default density bounds for plots',... |
---|
341 |
% 'densCorrected', 'density corrected for thermal lag',... |
---|
342 |
% 'depth', 'depth calculated as function of pressure and position latitude',... |
---|
343 |
% 'gpsLat', 'position latitude measured by glider GPS',... |
---|
344 |
% 'gpsLon', 'position longitude measured by glider GPS',... |
---|
345 |
% 'pres', 'pressure measured by glider',... |
---|
346 |
% 'ptime', 'time vector reported by glider',... |
---|
347 |
% 'ptime_datenum', 'Serial Date Number string',... |
---|
348 |
% 'salin', 'salinity measured by glider',... |
---|
349 |
% 'salinBounds', 'default salinity bounds for plots',... |
---|
350 |
% 'salinCorrected', 'salinity corrected for thermal lag',... |
---|
351 |
% 'temp', 'temperature measured by glider',... |
---|
352 |
% 'tempBounds', 'default temperature bounds for plots'); |
---|
353 |
|
---|
354 |
% correction_parameters = struct('alpha_offset', correctionParams(1),... |
---|
355 |
% 'alpha_slope', correctionParams(2),... |
---|
356 |
% 'tau_offset', correctionParams(3),... |
---|
357 |
% 'tau_slope', correctionParams(4)); |
---|
358 |
|
---|
359 |
% config = struct('glider_name', strGliderName,... |
---|
360 |
% 'deployment_number', strDeploymentNumber,... |
---|
361 |
% 'start_date', strStartDate,... |
---|
362 |
% 'end_date', strEndDate,... |
---|
363 |
% 'thermal_lag_correction_parameters', correction_parameters,... |
---|
364 |
% 'var_descriptions', variable_description,... |
---|
365 |
% 'var_units', units); |
---|
366 |
|
---|
367 |
|
---|
368 |
|
---|
369 |
|
---|
370 |
% set Level 1 data mat file name... |
---|
371 |
%strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_WithDO_LEVEL1.mat'); |
---|
372 |
|
---|
373 |
|
---|
374 |
|
---|
375 |
% save glider/deployment data to mat file... |
---|
376 |
% save(strMatFileName,... |
---|
377 |
% 'config',... |
---|
378 |
% 'chlor',... |
---|
379 |
% 'chlorBounds',... |
---|
380 |
% 'dens',... |
---|
381 |
% 'densBounds',... |
---|
382 |
% 'densCorrected',... |
---|
383 |
% 'depth',... |
---|
384 |
% 'gpsLat',... |
---|
385 |
% 'gpsLon',... |
---|
386 |
% 'pres',... |
---|
387 |
% 'ptime',... |
---|
388 |
% 'ptime_datenum',... |
---|
389 |
% 'salin',... |
---|
390 |
% 'salinBounds',... |
---|
391 |
% 'salinCorrected',... |
---|
392 |
% 'temp',... |
---|
393 |
% 'tempBounds'); |
---|
394 |
|
---|
395 |
|
---|
396 |
|
---|