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

root/gliderproc/trunk/gliderOptode_Generate_L1_Data.m

Revision 507 (checked in by cbc, 11 years ago)

Added hi/lo filters to input 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                     oxyw_oxygen = [oxyw_oxygen; ...
147                         data.data(:,strmatch('sci_oxy3835_wphase_oxygen',...
148                                              data.vars, 'exact'))];
149                     oxyw_saturation = [oxyw_saturation; ...
150                         data.data(:,strmatch('sci_oxy3835_wphase_saturation',...
151                                              data.vars, 'exact'))];
152                     oxyw_temp = [oxyw_temp; ...
153                         data.data(:,strmatch('sci_oxy3835_wphase_temp',...
154                                              data.vars, 'exact'))];
155                     oxyw_dphase = [oxyw_dphase; ...
156                         data.data(:,strmatch('sci_oxy3835_wphase_dphase',...
157                                              data.vars, 'exact'))];
158                     % oxyw_bphase = [oxyw_bphase; ...
159                     %     data.data(:,strmatch('sci_oxy3835_wphase_bphase',...
160                     %                          data.vars, 'exact'))];
161                     % oxyw_rphase = [oxyw_rphase; ...
162                     %     data.data(:,strmatch('sci_oxy3835_wphase_rphase',...
163                     %                          data.vars, 'exact'))];
164                     % oxyw_bamp = [oxyw_bamp; ...
165                     %     data.data(:,strmatch('sci_oxy3835_wphase_bamp',...
166                     %                          data.vars, 'exact'))];
167                     % oxyw_bpot = [oxyw_bpot; ...
168                     %     data.data(:,strmatch('sci_oxy3835_wphase_bpot',...
169                     %                          data.vars, 'exact'))];
170                     % oxyw_ramp = [oxyw_ramp; ...
171                     %     data.data(:,strmatch('sci_oxy3835_wphase_ramp',...
172                     %                          data.vars, 'exact'))];
173                     % oxyw_rawtemp = [oxyw_rawtemp; ...
174                     %     data.data(:,strmatch('sci_oxy3835_wphase_rawtemp',...
175                     %                          data.vars, 'exact'))];
176                     % oxyw_time = [oxyw_time; ...
177                     %     data.data(:,strmatch('sci_oxy3835_wphase_timestamp',...
178                     %                          data.vars, 'exact'))];
179                     % oxyw_installed = [oxyw_installed; ...
180                     %     data.data(:,strmatch('sci_oxy3835_wphase_is_installed',...
181                     %                          data.vars, 'exact'))];
182                     ptime_ebd = [ptime_ebd; ...
183                         data.data(:,strmatch('sci_m_present_time',...
184                                              data.vars, 'exact'))];
185                 end
186
187                 data = [];
188             end 
189         end
190         %*** END READ IN EBD DATA ****
191
192         % remove nans from EBD data...
193         i = find(~isnan(oxyw_dphase));
194         oxyw_dphase = oxyw_dphase(i);
195         ptime_ebd = ptime_ebd(i);
196         oxyw_oxygen = oxyw_oxygen(i);
197         oxyw_saturation = oxyw_saturation(i);
198         oxyw_temp = oxyw_temp(i);
199         % oxyw_bphase = oxyw_bphase(i);
200         % oxyw_rphase = oxyw_rphase(i);
201         % oxyw_bamp = oxyw_bamp(i);
202         % oxyw_bpot = oxyw_bpot(i);
203         % oxyw_ramp = oxyw_ramp(i);
204         % oxyw_rawtemp = oxyw_rawtemp(i);
205         % oxyw_time = oxyw_time(i);
206         % oxyw_installed = oxyw_installed(i);
207
208         % remove unreasonably low oxygen values from EBD data...
209         i = find(gt(oxyw_oxygen, lo_oxy));
210         oxyw_oxygen = oxyw_oxygen(i);
211         ptime_ebd = ptime_ebd(i);
212         oxyw_saturation = oxyw_saturation(i);
213         oxyw_dphase = oxyw_dphase(i);
214         oxyw_temp = oxyw_temp(i);
215         % oxyw_bphase = oxyw_bphase(i);
216         % oxyw_rphase = oxyw_rphase(i);
217         % oxyw_bamp = oxyw_bamp(i);
218         % oxyw_bpot = oxyw_bpot(i);
219         % oxyw_ramp = oxyw_ramp(i);
220         % oxyw_rawtemp = oxyw_rawtemp(i);
221         % oxyw_time = oxyw_time(i);
222         % oxyw_installed = oxyw_installed(i);
223
224         % remove unreasonably high oxygen values from EBD data...
225         i = find(lt(oxyw_oxygen, hi_oxy));
226         oxyw_oxygen = oxyw_oxygen(i);
227         ptime_ebd = ptime_ebd(i);
228         oxyw_saturation = oxyw_saturation(i);
229         oxyw_dphase = oxyw_dphase(i);
230         oxyw_temp = oxyw_temp(i);
231         % oxyw_bphase = oxyw_bphase(i);
232         % oxyw_rphase = oxyw_rphase(i);
233         % oxyw_bamp = oxyw_bamp(i);
234         % oxyw_bpot = oxyw_bpot(i);
235         % oxyw_ramp = oxyw_ramp(i);
236         % oxyw_rawtemp = oxyw_rawtemp(i);
237         % oxyw_time = oxyw_time(i);
238         % oxyw_installed = oxyw_installed(i);
239
240         % remove unreasonably low saturation values from EBD data...
241         i = find(gt(oxyw_saturation, lo_sat));
242         oxyw_saturation = oxyw_saturation(i);
243         ptime_ebd = ptime_ebd(i);
244         oxyw_oxygen = oxyw_oxygen(i);
245         oxyw_dphase = oxyw_dphase(i);
246         oxyw_temp = oxyw_temp(i);
247         % oxyw_bphase = oxyw_bphase(i);
248         % oxyw_rphase = oxyw_rphase(i);
249         % oxyw_bamp = oxyw_bamp(i);
250         % oxyw_bpot = oxyw_bpot(i);
251         % oxyw_ramp = oxyw_ramp(i);
252         % oxyw_rawtemp = oxyw_rawtemp(i);
253         % oxyw_time = oxyw_time(i);
254         % oxyw_installed = oxyw_installed(i);
255
256         % remove unreasonably high saturation values from EBD data...
257         i = find(lt(oxyw_saturation, hi_sat));
258         oxyw_saturation = oxyw_saturation(i);
259         ptime_ebd = ptime_ebd(i);
260         oxyw_oxygen = oxyw_oxygen(i);
261         oxyw_dphase = oxyw_dphase(i);
262         oxyw_temp = oxyw_temp(i);
263         % oxyw_bphase = oxyw_bphase(i);
264         % oxyw_rphase = oxyw_rphase(i);
265         % oxyw_bamp = oxyw_bamp(i);
266         % oxyw_bpot = oxyw_bpot(i);
267         % oxyw_ramp = oxyw_ramp(i);
268         % oxyw_rawtemp = oxyw_rawtemp(i);
269         % oxyw_time = oxyw_time(i);
270         % oxyw_installed = oxyw_installed(i);
271
272         % apply the sort() function to ptime_oxy
273         % to make sure it increases monotonically...
274         [Y,I] = sort(ptime_ebd);
275         ptime_ebd = Y;
276         oxyw_oxygen = oxyw_oxygen(I);
277         oxyw_saturation = oxyw_saturation(I);
278         oxyw_temp = oxyw_temp(I);
279         oxyw_dphase = oxyw_dphase(I);
280         % oxyw_bphase = oxyw_bphase(I);
281         % oxyw_rphase = oxyw_rphase(I);
282         % oxyw_bamp = oxyw_bamp(I);
283         % oxyw_bpot = oxyw_bpot(I);
284         % oxyw_ramp = oxyw_ramp(I);
285         % oxyw_rawtemp = oxyw_rawtemp(I);
286         % oxyw_time = oxyw_time(I);
287         % oxyw_installed = oxyw_installed(I);
288
289         % load CTD data from existing mat file...
290         strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, ...
291             '_Deployment', strDeploymentNumber, '_CTD_L1.mat');
292         load(strMatFilePath);
293
294         % interpolate CTD to times of oxyw_dphase
295         tempi = interp1(ptime,temp,ptime_ebd);
296         salini = interp1(ptime,salinCorrected,ptime_ebd);
297         depthi = interp1(ptime,depth,ptime_ebd);
298
299         % first, implement temperature-dependent correction to DO concentration
300         % that utilizes dphase as the input (manual page 30), coefficients from cal
301         % sheets (a 5x4 matrix of values) and will be glider (and foil) dependent.
302         C0 = C0coef(1) + (C0coef(2) .* tempi) + ...
303             (C0coef(3) .* (tempi.^2)) + (C0coef(4) .* (tempi.^3));
304         C1 = C1coef(1) + (C1coef(2) .* tempi) + ...
305             (C1coef(3) .* (tempi.^2)) + (C1coef(4) .* (tempi.^3));
306         C2 = C2coef(1) + (C2coef(2) .* tempi) + ...
307             (C2coef(3) .* (tempi.^2)) + (C2coef(4) .* (tempi.^3));
308         C3 = C3coef(1) + (C3coef(2) .* tempi) + ...
309             (C3coef(3) .* (tempi.^2)) + (C3coef(4) .* (tempi.^3));
310         C4 = C4coef(1) + (C4coef(2) .* tempi) + ...
311             (C4coef(3) .* (tempi.^2)) + (C4coef(4) .* (tempi.^3));
312
313         o2_tcorr = C0 + (C1 .* oxyw_dphase) + (C2 .* (oxyw_dphase.^2)) + ...
314             (C3 .* (oxyw_dphase.^3)) + (C4 .* (oxyw_dphase.^4));
315
316         % second, implement the salinity correction to DO concentration (page 31)
317         % ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER
318         % SCALING IS REQUIRED
319         B0=-6.24097e-3;
320         B1=-6.93498e-3;
321         B2=-6.90358e-3;
322         B3=-4.29155e-3;
323         C0=-3.11680e-7;
324
325         Ts = log((298.15 - tempi) ./ (273.15 + tempi));
326
327         o2_tscorr = o2_tcorr .* ...
328             exp((salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + ...
329                             (B3 .* (Ts.^3))))  + ...
330             (C0 .* (salini.^2)));
331
332         % third, implement the pressure correction to DO concentration (page 32)
333         o2_tspcorr = o2_tscorr .* (1 + (0.04 .* depthi ./ 1000));
334
335         % use polynomial to calculate DO satuations using the measured temp and
336         % sal (manual page 30)
337         A0 = 2.00856;
338         A1 = 3.22400;
339         A2 = 3.99063;
340         A3 = 4.80299;
341         A4 = 9.78188e-1;
342         A5 = 1.71069;
343         B0 = -6.24097e-3;
344         B1 = -6.93498e-3;
345         B2 = -6.90358e-3;
346         B3 = -4.29155e-3;
347         C0 = -3.11680e-7;
348
349         % need interpolated temp, salinity, pressure at times of DO obs
350         rslt = (A0 + (A1 .* Ts) + (A2 .* (Ts.^2)) + (A3 .* (Ts.^3)) + ...
351                (A4 .* (Ts.^4)) + (A5 .* (Ts.^5))) + ...
352             (salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + (B3 .*(Ts.^3)))) + ...
353             (C0 .* (salini.^2));
354
355         o2_sol = exp(rslt);
356
357         o2_sat = (o2_tspcorr .* 2.2414) ./ o2_sol;
358
359         % convert ptime into datenum style...
360         ptime_ebd_datenum = (ptime_ebd/3600/24) + datenum(1970, 1, 1, 0, 0, 0);
361
362         % create configuration struct...
363         units = struct( ...
364             'oxyw_oxygen', '10e-6 mol/dm3', ...
365             'oxyw_saturation', 'percent', ...
366             'oxyw_temp', 'degrees C', ...
367             'oxyw_dphase', 'degrees', ...
368             ... % 'oxyw_bphase', 'degrees', ...
369             ... % 'oxyw_rphase', 'degrees', ...
370             ... % 'oxyw_bamp', 'mA', ...
371             ... % 'oxyw_bpot', 'mV', ...
372             ... % 'oxyw_ramp', 'mA', ...
373             ... % 'oxyw_rawtemp', 'degrees C', ...
374             ... % 'oxyw_time', 'timestamp', ...
375             ... % 'oxyw_installed', 'bool', ...
376             'ptime_ebd', 'seconds since 0000-01-01T00:00', ...
377             'ptime_ebd_datenum', 'days since 1970-01-01T00:00', ...
378             'tempi', 'degrees C', ...
379             'salini', 'psu', ...
380             'depthi', 'meters', ...
381             'o2_tcorr', '10e-6 mol/dm3', ...
382             'o2_tscorr', '10e-6 mol/dm3', ...
383             'o2_tspcorr', '10e-6 mol/dm3', ...
384             'o2_sol', 'cm3/liter at 1031 hPa', ...
385             'o2_sat', 'percent');
386        
387         variable_description = struct( ...
388             'oxyw_oxygen', 'dissolved oxygen', ...
389             'oxyw_saturation', 'dissolved oxygen saturation', ...
390             'oxyw_temp', 'water temperature', ...
391             'oxyw_dphase', 'phase difference', ...
392             ... % 'oxyw_bphase', 'blue phase', ...
393             ... % 'oxyw_rphase', 'red phase', ...
394             ... % 'oxyw_bamp', 'blue current bias', ...
395             ... % 'oxyw_bpot', 'blue voltage bias', ...
396             ... % 'oxyw_ramp', 'red current bias', ...
397             ... % 'oxyw_rawtemp', 'raw water temperature', ...
398             ... % 'oxyw_time', 'optode timestampe', ...
399             ... % 'oxyw_installed', 'bool', ...
400             'ptime_ebd', 'science computer time', ...
401             'ptime_ebd_datenum', 'science computer date', ...
402             'tempi', 'interpolated CTD water temperature', ...
403             'salini', 'interpoloted CTD salinity', ...
404             'depthi', 'interpolated CTD depth', ...
405             'o2_tcorr', 'temperature corrected dissolved oxygen', ...
406             'o2_tscorr', 'temperature and salinity corrected dissolved oxygen', ...
407             'o2_tspcorr', 'temperature, salinity, and pressure corrected dissolved oxygen', ...
408             'o2_sol', 'corrected oxygen solubility', ...
409             'o2_sat', 'corrected oxygen saturation');
410
411         correction_coefficients = struct('C0coef', C0coef, ...
412                                          'C1coef', C1coef, ...
413                                          'C2coef', C2coef, ...
414                                          'C3coef', C3coef, ...
415                                          'C4coef', C4coef);
416
417         config = struct('glider_name', strGliderName,...
418                         'deployment_number', strDeploymentNumber,...
419                         'start_date', strStartDate,...
420                         'end_date', strEndDate,...
421                         'correction_coefficients', correction_coefficients,...
422                         'var_descriptions', variable_description,...
423                         'var_units', units);
424
425         % set Level 1 data mat file name...
426         strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_DO_L1.mat');
427
428         % save flight data to mat file...
429         save(strMatFileName,...
430              'config', ...
431              'oxyw_oxygen', ...
432              'oxyw_saturation', ...
433              'oxyw_temp', ...
434              'oxyw_dphase', ...
435              ... % 'oxyw_bphase', ...
436              ... % 'oxyw_rphase', ...
437              ... % 'oxyw_bamp', ...
438              ... % 'oxyw_bpot', ...
439              ... % 'oxyw_ramp', ...
440              ... % 'oxyw_rawtemp', ...
441              ... % 'oxyw_time', ...
442              ... % 'oxyw_installed', ...
443              'ptime_ebd', ...
444              'ptime_ebd_datenum', ...
445              'tempi', ...
446              'salini', ...
447              'depthi', ...
448              'o2_tcorr', ...
449              'o2_tscorr', ...
450              'o2_tspcorr', ...
451              'o2_sol', ...
452              'o2_sat');
453
454     end
455 end
Note: See TracBrowser for help on using the browser.