NCCOOS Trac Projects: Top | Web | Platforms | Processing | Viz | Sprints | Sandbox | (Wind)

root/gliderproc/trunk/gliderOptode_Generate_L1_Data.m

Revision 514 (checked in by cbc, 10 years ago)

Fix optode extraction code to ignore EBD files with no optode data.

Line 
1 %
2 %  gliderOptode_Generate_L1_Data.m
3 %
4 %  Purpose: Generate Level 1 data mat files that include optode DO data.
5 %
6 %  Requires:
7 %  MATLAB folder - contains util
8 %
9 %  Authors:  Adapted from William Stark by Harvey Seim and Chris Calloway
10 %            Marine Sciences Department
11 %            UNC-Chapel Hill
12 %
13 %  Created: March 2013
14 %
15 %///////////////////////////////////////////////////////////////////////////////
16
17 clear all;
18
19 % add paths for required files...
20 addpath('MATLAB/util/');
21
22 % populate arrays for the deployment start and end dates...
23 % ex. strStart(2, 3) is start date for Ramses, Deployment 3
24 strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; ...
25             '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'};
26 strEnd   = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; ...
27             '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'};
28
29 % SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ...
30 for gliderIndex=1:2
31    
32     % SET THE DEPLOYMENT NUMBER (1, 2 or 3) ...
33     for deploymentNumber=1:3
34        
35         clearvars -except gliderIndex deploymentNumber strStart strEnd;
36
37         % glider name string...
38         if (gliderIndex==1)
39             strGliderName = 'Pelagia';
40         else
41             strGliderName = 'Ramses';
42         end
43
44         % deployment number string...
45         strDeploymentNumber = num2str(deploymentNumber);
46
47         % deployment start date string...
48         strStartDate = strStart(gliderIndex, deploymentNumber);
49
50         % deployment end date string...
51         strEndDate = strEnd(gliderIndex, deploymentNumber);
52
53         % define the path to the glider ascii files...
54         datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, ...
55                          '_Deployment', strDeploymentNumber, '/');
56
57         % define optode foil calibration coefficients
58         % and oxygen bounds values by glider and deployment
59         switch gliderIndex
60             case 1 % Pelagia
61                 switch deploymentNumber
62                     case {1, 2} % Deployment 1 or 2;
63                                 % before foil replacement;
64                                 % 21 June 2004 calibration date;
65                                 % batch 2204
66                         C0coef=[3.12899E+03,-1.05764E+02,2.06640E+00,-1.68983E-02];
67                         C1coef=[-1.67671E+02,4.80582E+00,-8.73323E-02,6.61507E-04];
68                         C2coef=[3.73685E+00,-8.78300E-02,1.47560E-03,-9.96701E-06];
69                         C3coef=[-3.96052E-02,7.46930E-04,-1.17804E-05,6.677619E-08];
70                         C4coef=[1.61999E-04,-2.37870E-06,3.63223E-08,-1.62194E-10];
71                         lo_oxy = 50;
72                         hi_oxy = 2000;
73                         lo_sat = 50;
74                         hi_sat = 2000;
75                     case 3 % Deployment 3;
76                            % replaced foil on 12 March 2012;
77                            % 18 August 2010 calibration date;
78                            % batch 1023
79                         C0coef=[4.27019336E+03,-1.32723585E+02,2.15629751E+00,-1.40275831E-02];
80                         C1coef=[-2.29729690E+02,5.74242078E+00,-6.85357898E-02,1.88612346E-04];
81                         C2coef=[5.06401550E+00,-9.62084932E-02,5.22180779E-04,7.70889717E-06];
82                         C3coef=[-5.26332308E-02,7.15467419E-04,3.31185072E-06,-1.86124024E-07];
83                         C4coef=[2.10916841E-04,-1.84087896E-06,-4.28645540E-08,1.11120317E-09];
84                         lo_oxy = 100;
85                         hi_oxy = 300;
86                         lo_sat = 20;
87                         hi_sat = 140;
88                 end
89             case 2 % Ramses all deployments;
90                    % 2 June 2010 calibration date;
91                    % batch 5009
92                 C0coef=[4.53793e3 -1.62595e2 3.29574 -2.79285e-2];
93                 C1coef=[-2.50953e2 8.02322 -1.58398e-1 1.31141e-3];
94                 C2coef=[5.66417 -1.59647e-1 3.07910e-3 -2.46265e-5];
95                 C3coef=[-5.99449e-2 1.48326e-3 -2.82110e-5 2.15156e-7];
96                 C4coef=[2.43614e-4 -5.26759e-6 1.00064e-7 -7.14320e-10];
97                 lo_oxy = 10;
98                 hi_oxy = 250;
99                 lo_sat = 20;
100                 hi_sat = 100;
101         end
102
103         disp(['Processing DO for ', strGliderName, ' Deployment ', strDeploymentNumber]);
104         disp(['lo_oxy filter = ', num2str(lo_oxy)]);
105         disp(['hi_oxy filter = ', num2str(hi_oxy)]);
106         disp(['lo_sat filter = ', num2str(lo_sat)]);
107         disp(['hi_sat filter = ', num2str(hi_sat)]);
108        
109         %*** READ IN EBD DATA ****
110         % declare arrays for accumulating data
111         oxyw_oxygen = [];
112         oxyw_saturation = [];
113         oxyw_temp = [];
114         oxyw_dphase = [];
115         % oxyw_bphase = [];
116         % oxyw_rphase = [];
117         % oxyw_bamp = [];
118         % oxyw_bpot = [];
119         % oxyw_ramp = [];
120         % oxyw_rawtemp = [];
121         % oxyw_time = [];
122         % oxyw_installed = [];
123         ptime_ebd = [];
124
125         % try to load all *.ebdasc files at once...
126         [files, Dstruct] = wilddir(datadir, '.ebdasc');
127         nfile = size(files, 1);
128
129         for i=1:nfile-1
130             % protect against empty ebd file
131             if(Dstruct(i).bytes>0)
132                 data = read_gliderasc2([datadir, files(i,:)]);
133
134                 % if the number of values (in data.data) is less than the number
135                 % of vars (in data.vars), this means that the data were not
136                 % completely read in. To correct this, pad data.data with NaNs
137                 % until its length equals that of data.vars...
138                 if (length(data.data) < length(data.vars))
139                     data.data = padarray(data.data, ...
140                         [0 length(data.vars)-length(data.data)], ...
141                         NaN, 'post');
142                 end
143
144                 % concatenate variables with data...
145                 if(~isempty(data.data))
146                     if (~isempty(strmatch('sci_oxy3835_wphase_oxygen',data.vars,'exact')))
147                         oxyw_oxygen = [oxyw_oxygen; ...
148                             data.data(:,strmatch('sci_oxy3835_wphase_oxygen',...
149                                                  data.vars, 'exact'))];
150                         oxyw_saturation = [oxyw_saturation; ...
151                             data.data(:,strmatch('sci_oxy3835_wphase_saturation',...
152                                                  data.vars, 'exact'))];
153                         oxyw_temp = [oxyw_temp; ...
154                             data.data(:,strmatch('sci_oxy3835_wphase_temp',...
155                                                  data.vars, 'exact'))];
156                         oxyw_dphase = [oxyw_dphase; ...
157                             data.data(:,strmatch('sci_oxy3835_wphase_dphase',...
158                                                  data.vars, 'exact'))];
159                         % oxyw_bphase = [oxyw_bphase; ...
160                         %     data.data(:,strmatch('sci_oxy3835_wphase_bphase',...
161                         %                          data.vars, 'exact'))];
162                         % oxyw_rphase = [oxyw_rphase; ...
163                         %     data.data(:,strmatch('sci_oxy3835_wphase_rphase',...
164                         %                          data.vars, 'exact'))];
165                         % oxyw_bamp = [oxyw_bamp; ...
166                         %     data.data(:,strmatch('sci_oxy3835_wphase_bamp',...
167                         %                          data.vars, 'exact'))];
168                         % oxyw_bpot = [oxyw_bpot; ...
169                         %     data.data(:,strmatch('sci_oxy3835_wphase_bpot',...
170                         %                          data.vars, 'exact'))];
171                         % oxyw_ramp = [oxyw_ramp; ...
172                         %     data.data(:,strmatch('sci_oxy3835_wphase_ramp',...
173                         %                          data.vars, 'exact'))];
174                         % oxyw_rawtemp = [oxyw_rawtemp; ...
175                         %     data.data(:,strmatch('sci_oxy3835_wphase_rawtemp',...
176                         %                          data.vars, 'exact'))];
177                         % oxyw_time = [oxyw_time; ...
178                         %     data.data(:,strmatch('sci_oxy3835_wphase_timestamp',...
179                         %                          data.vars, 'exact'))];
180                         % oxyw_installed = [oxyw_installed; ...
181                         %     data.data(:,strmatch('sci_oxy3835_wphase_is_installed',...
182                         %                          data.vars, 'exact'))];
183                         ptime_ebd = [ptime_ebd; ...
184                             data.data(:,strmatch('sci_m_present_time',...
185                                                  data.vars, 'exact'))];
186                     end
187                 end
188
189                 data = [];
190             end 
191         end
192         %*** END READ IN EBD DATA ****
193
194         % remove nans from EBD data...
195         i = find(~isnan(oxyw_dphase));
196         oxyw_dphase = oxyw_dphase(i);
197         ptime_ebd = ptime_ebd(i);
198         oxyw_oxygen = oxyw_oxygen(i);
199         oxyw_saturation = oxyw_saturation(i);
200         oxyw_temp = oxyw_temp(i);
201         % oxyw_bphase = oxyw_bphase(i);
202         % oxyw_rphase = oxyw_rphase(i);
203         % oxyw_bamp = oxyw_bamp(i);
204         % oxyw_bpot = oxyw_bpot(i);
205         % oxyw_ramp = oxyw_ramp(i);
206         % oxyw_rawtemp = oxyw_rawtemp(i);
207         % oxyw_time = oxyw_time(i);
208         % oxyw_installed = oxyw_installed(i);
209
210         % remove unreasonably low oxygen values from EBD data...
211         i = find(gt(oxyw_oxygen, lo_oxy));
212         oxyw_oxygen = oxyw_oxygen(i);
213         ptime_ebd = ptime_ebd(i);
214         oxyw_saturation = oxyw_saturation(i);
215         oxyw_dphase = oxyw_dphase(i);
216         oxyw_temp = oxyw_temp(i);
217         % oxyw_bphase = oxyw_bphase(i);
218         % oxyw_rphase = oxyw_rphase(i);
219         % oxyw_bamp = oxyw_bamp(i);
220         % oxyw_bpot = oxyw_bpot(i);
221         % oxyw_ramp = oxyw_ramp(i);
222         % oxyw_rawtemp = oxyw_rawtemp(i);
223         % oxyw_time = oxyw_time(i);
224         % oxyw_installed = oxyw_installed(i);
225
226         % remove unreasonably high oxygen values from EBD data...
227         i = find(lt(oxyw_oxygen, hi_oxy));
228         oxyw_oxygen = oxyw_oxygen(i);
229         ptime_ebd = ptime_ebd(i);
230         oxyw_saturation = oxyw_saturation(i);
231         oxyw_dphase = oxyw_dphase(i);
232         oxyw_temp = oxyw_temp(i);
233         % oxyw_bphase = oxyw_bphase(i);
234         % oxyw_rphase = oxyw_rphase(i);
235         % oxyw_bamp = oxyw_bamp(i);
236         % oxyw_bpot = oxyw_bpot(i);
237         % oxyw_ramp = oxyw_ramp(i);
238         % oxyw_rawtemp = oxyw_rawtemp(i);
239         % oxyw_time = oxyw_time(i);
240         % oxyw_installed = oxyw_installed(i);
241
242         % remove unreasonably low saturation values from EBD data...
243         i = find(gt(oxyw_saturation, lo_sat));
244         oxyw_saturation = oxyw_saturation(i);
245         ptime_ebd = ptime_ebd(i);
246         oxyw_oxygen = oxyw_oxygen(i);
247         oxyw_dphase = oxyw_dphase(i);
248         oxyw_temp = oxyw_temp(i);
249         % oxyw_bphase = oxyw_bphase(i);
250         % oxyw_rphase = oxyw_rphase(i);
251         % oxyw_bamp = oxyw_bamp(i);
252         % oxyw_bpot = oxyw_bpot(i);
253         % oxyw_ramp = oxyw_ramp(i);
254         % oxyw_rawtemp = oxyw_rawtemp(i);
255         % oxyw_time = oxyw_time(i);
256         % oxyw_installed = oxyw_installed(i);
257
258         % remove unreasonably high saturation values from EBD data...
259         i = find(lt(oxyw_saturation, hi_sat));
260         oxyw_saturation = oxyw_saturation(i);
261         ptime_ebd = ptime_ebd(i);
262         oxyw_oxygen = oxyw_oxygen(i);
263         oxyw_dphase = oxyw_dphase(i);
264         oxyw_temp = oxyw_temp(i);
265         % oxyw_bphase = oxyw_bphase(i);
266         % oxyw_rphase = oxyw_rphase(i);
267         % oxyw_bamp = oxyw_bamp(i);
268         % oxyw_bpot = oxyw_bpot(i);
269         % oxyw_ramp = oxyw_ramp(i);
270         % oxyw_rawtemp = oxyw_rawtemp(i);
271         % oxyw_time = oxyw_time(i);
272         % oxyw_installed = oxyw_installed(i);
273
274         % apply the sort() function to ptime_oxy
275         % to make sure it increases monotonically...
276         [Y,I] = sort(ptime_ebd);
277         ptime_ebd = Y;
278         oxyw_oxygen = oxyw_oxygen(I);
279         oxyw_saturation = oxyw_saturation(I);
280         oxyw_temp = oxyw_temp(I);
281         oxyw_dphase = oxyw_dphase(I);
282         % oxyw_bphase = oxyw_bphase(I);
283         % oxyw_rphase = oxyw_rphase(I);
284         % oxyw_bamp = oxyw_bamp(I);
285         % oxyw_bpot = oxyw_bpot(I);
286         % oxyw_ramp = oxyw_ramp(I);
287         % oxyw_rawtemp = oxyw_rawtemp(I);
288         % oxyw_time = oxyw_time(I);
289         % oxyw_installed = oxyw_installed(I);
290
291         % load CTD data from existing mat file...
292         strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, ...
293             '_Deployment', strDeploymentNumber, '_CTD_L1.mat');
294         load(strMatFilePath);
295
296         % interpolate CTD to times of oxyw_dphase
297         tempi = interp1(ptime,temp,ptime_ebd);
298         salini = interp1(ptime,salinCorrected,ptime_ebd);
299         depthi = interp1(ptime,depth,ptime_ebd);
300
301         % first, implement temperature-dependent correction to DO concentration
302         % that utilizes dphase as the input (manual page 30), coefficients from cal
303         % sheets (a 5x4 matrix of values) and will be glider (and foil) dependent.
304         C0 = C0coef(1) + (C0coef(2) .* tempi) + ...
305             (C0coef(3) .* (tempi.^2)) + (C0coef(4) .* (tempi.^3));
306         C1 = C1coef(1) + (C1coef(2) .* tempi) + ...
307             (C1coef(3) .* (tempi.^2)) + (C1coef(4) .* (tempi.^3));
308         C2 = C2coef(1) + (C2coef(2) .* tempi) + ...
309             (C2coef(3) .* (tempi.^2)) + (C2coef(4) .* (tempi.^3));
310         C3 = C3coef(1) + (C3coef(2) .* tempi) + ...
311             (C3coef(3) .* (tempi.^2)) + (C3coef(4) .* (tempi.^3));
312         C4 = C4coef(1) + (C4coef(2) .* tempi) + ...
313             (C4coef(3) .* (tempi.^2)) + (C4coef(4) .* (tempi.^3));
314
315         o2_tcorr = C0 + (C1 .* oxyw_dphase) + (C2 .* (oxyw_dphase.^2)) + ...
316             (C3 .* (oxyw_dphase.^3)) + (C4 .* (oxyw_dphase.^4));
317
318         % second, implement the salinity correction to DO concentration (page 31)
319         % ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER
320         % SCALING IS REQUIRED
321         B0=-6.24097e-3;
322         B1=-6.93498e-3;
323         B2=-6.90358e-3;
324         B3=-4.29155e-3;
325         C0=-3.11680e-7;
326
327         Ts = log((298.15 - tempi) ./ (273.15 + tempi));
328
329         o2_tscorr = o2_tcorr .* ...
330             exp((salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + ...
331                             (B3 .* (Ts.^3))))  + ...
332             (C0 .* (salini.^2)));
333
334         % third, implement the pressure correction to DO concentration (page 32)
335         o2_tspcorr = o2_tscorr .* (1 + (0.04 .* depthi ./ 1000));
336
337         % use polynomial to calculate DO satuations using the measured temp and
338         % sal (manual page 30)
339         A0 = 2.00856;
340         A1 = 3.22400;
341         A2 = 3.99063;
342         A3 = 4.80299;
343         A4 = 9.78188e-1;
344         A5 = 1.71069;
345         B0 = -6.24097e-3;
346         B1 = -6.93498e-3;
347         B2 = -6.90358e-3;
348         B3 = -4.29155e-3;
349         C0 = -3.11680e-7;
350
351         % need interpolated temp, salinity, pressure at times of DO obs
352         rslt = (A0 + (A1 .* Ts) + (A2 .* (Ts.^2)) + (A3 .* (Ts.^3)) + ...
353                (A4 .* (Ts.^4)) + (A5 .* (Ts.^5))) + ...
354             (salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + (B3 .*(Ts.^3)))) + ...
355             (C0 .* (salini.^2));
356
357         o2_sol = exp(rslt);
358
359         o2_sat = (o2_tspcorr .* 2.2414) ./ o2_sol;
360
361         % convert ptime into datenum style...
362         ptime_ebd_datenum = (ptime_ebd/3600/24) + datenum(1970, 1, 1, 0, 0, 0);
363
364         % create configuration struct...
365         units = struct( ...
366             'oxyw_oxygen', '10e-6 mol/dm3', ...
367             'oxyw_saturation', 'percent', ...
368             'oxyw_temp', 'degrees C', ...
369             'oxyw_dphase', 'degrees', ...
370             ... % 'oxyw_bphase', 'degrees', ...
371             ... % 'oxyw_rphase', 'degrees', ...
372             ... % 'oxyw_bamp', 'mA', ...
373             ... % 'oxyw_bpot', 'mV', ...
374             ... % 'oxyw_ramp', 'mA', ...
375             ... % 'oxyw_rawtemp', 'degrees C', ...
376             ... % 'oxyw_time', 'timestamp', ...
377             ... % 'oxyw_installed', 'bool', ...
378             'ptime_ebd', 'seconds since 0000-01-01T00:00', ...
379             'ptime_ebd_datenum', 'days since 1970-01-01T00:00', ...
380             'tempi', 'degrees C', ...
381             'salini', 'psu', ...
382             'depthi', 'meters', ...
383             'o2_tcorr', '10e-6 mol/dm3', ...
384             'o2_tscorr', '10e-6 mol/dm3', ...
385             'o2_tspcorr', '10e-6 mol/dm3', ...
386             'o2_sol', 'cm3/liter at 1031 hPa', ...
387             'o2_sat', 'percent');
388        
389         variable_description = struct( ...
390             'oxyw_oxygen', 'dissolved oxygen', ...
391             'oxyw_saturation', 'dissolved oxygen saturation', ...
392             'oxyw_temp', 'water temperature', ...
393             'oxyw_dphase', 'phase difference', ...
394             ... % 'oxyw_bphase', 'blue phase', ...
395             ... % 'oxyw_rphase', 'red phase', ...
396             ... % 'oxyw_bamp', 'blue current bias', ...
397             ... % 'oxyw_bpot', 'blue voltage bias', ...
398             ... % 'oxyw_ramp', 'red current bias', ...
399             ... % 'oxyw_rawtemp', 'raw water temperature', ...
400             ... % 'oxyw_time', 'optode timestampe', ...
401             ... % 'oxyw_installed', 'bool', ...
402             'ptime_ebd', 'science computer time', ...
403             'ptime_ebd_datenum', 'science computer date', ...
404             'tempi', 'interpolated CTD water temperature', ...
405             'salini', 'interpoloted CTD salinity', ...
406             'depthi', 'interpolated CTD depth', ...
407             'o2_tcorr', 'temperature corrected dissolved oxygen', ...
408             'o2_tscorr', 'temperature and salinity corrected dissolved oxygen', ...
409             'o2_tspcorr', 'temperature, salinity, and pressure corrected dissolved oxygen', ...
410             'o2_sol', 'corrected oxygen solubility', ...
411             'o2_sat', 'corrected oxygen saturation');
412
413         correction_coefficients = struct('C0coef', C0coef, ...
414                                          'C1coef', C1coef, ...
415                                          'C2coef', C2coef, ...
416                                          'C3coef', C3coef, ...
417                                          'C4coef', C4coef);
418
419         config = struct('glider_name', strGliderName,...
420                         'deployment_number', strDeploymentNumber,...
421                         'start_date', strStartDate,...
422                         'end_date', strEndDate,...
423                         'correction_coefficients', correction_coefficients,...
424                         'var_descriptions', variable_description,...
425                         'var_units', units);
426
427         % set Level 1 data mat file name...
428         strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_DO_L1.mat');
429
430         % save flight data to mat file...
431         save(strMatFileName,...
432              'config', ...
433              'oxyw_oxygen', ...
434              'oxyw_saturation', ...
435              'oxyw_temp', ...
436              'oxyw_dphase', ...
437              ... % 'oxyw_bphase', ...
438              ... % 'oxyw_rphase', ...
439              ... % 'oxyw_bamp', ...
440              ... % 'oxyw_bpot', ...
441              ... % 'oxyw_ramp', ...
442              ... % 'oxyw_rawtemp', ...
443              ... % 'oxyw_time', ...
444              ... % 'oxyw_installed', ...
445              'ptime_ebd', ...
446              'ptime_ebd_datenum', ...
447              'tempi', ...
448              'salini', ...
449              'depthi', ...
450              'o2_tcorr', ...
451              'o2_tscorr', ...
452              'o2_tspcorr', ...
453              'o2_sol', ...
454              'o2_sat');
455
456     end
457 end
Note: See TracBrowser for help on using the browser.