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

root/gliderproc/trunk/gliderOptode_Generate_L1_Data.m

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

Regularize corrected variable names.

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         % 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                     case 3 % Deployment 3;
72                            % replaced foil on 12 March 2012;
73                            % 18 August 2010 calibration date;
74                            % batch 1023
75                         C0coef=[4.27019336E+03,-1.32723585E+02,2.15629751E+00,-1.40275831E-02];
76                         C1coef=[-2.29729690E+02,5.74242078E+00,-6.85357898E-02,1.88612346E-04];
77                         C2coef=[5.06401550E+00,-9.62084932E-02,5.22180779E-04,7.70889717E-06];
78                         C3coef=[-5.26332308E-02,7.15467419E-04,3.31185072E-06,-1.86124024E-07];
79                         C4coef=[2.10916841E-04,-1.84087896E-06,-4.28645540E-08,1.11120317E-09];
80                 end
81             case 2 % Ramses all deployments;
82                    % 2 June 2010 calibration date;
83                    % batch 5009
84                 C0coef=[4.53793e3 -1.62595e2 3.29574 -2.79285e-2];
85                 C1coef=[-2.50953e2 8.02322 -1.58398e-1 1.31141e-3];
86                 C2coef=[5.66417 -1.59647e-1 3.07910e-3 -2.46265e-5];
87                 C3coef=[-5.99449e-2 1.48326e-3 -2.82110e-5 2.15156e-7];
88                 C4coef=[2.43614e-4 -5.26759e-6 1.00064e-7 -7.14320e-10];
89         end
90
91         %*** READ IN EBD DATA ****
92         % declare arrays for accumulating data
93         oxyw_oxygen = [];
94         oxyw_saturation = [];
95         oxyw_temp = [];
96         oxyw_dphase = [];
97         % oxyw_bphase = [];
98         % oxyw_rphase = [];
99         % oxyw_bamp = [];
100         % oxyw_bpot = [];
101         % oxyw_ramp = [];
102         % oxyw_rawtemp = [];
103         % oxyw_time = [];
104         % oxyw_installed = [];
105         ptime_ebd = [];
106
107         % try to load all *.ebdasc files at once...
108         [files, Dstruct] = wilddir(datadir, '.ebdasc');
109         nfile = size(files, 1);
110
111         for i=1:nfile-1
112             % protect against empty ebd file
113             if(Dstruct(i).bytes>0)
114                 data = read_gliderasc2([datadir, files(i,:)]);
115
116                 % if the number of values (in data.data) is less than the number
117                 % of vars (in data.vars), this means that the data were not
118                 % completely read in. To correct this, pad data.data with NaNs
119                 % until its length equals that of data.vars...
120                 if (length(data.data) < length(data.vars))
121                     data.data = padarray(data.data, ...
122                         [0 length(data.vars)-length(data.data)], ...
123                         NaN, 'post');
124                 end
125
126                 % concatenate variables with data...
127                 if(~isempty(data.data))
128                     oxyw_oxygen = [oxyw_oxygen; ...
129                         data.data(:,strmatch('sci_oxy3835_wphase_oxygen',...
130                                              data.vars, 'exact'))];
131                     oxyw_saturation = [oxyw_saturation; ...
132                         data.data(:,strmatch('sci_oxy3835_wphase_saturation',...
133                                              data.vars, 'exact'))];
134                     oxyw_temp = [oxyw_temp; ...
135                         data.data(:,strmatch('sci_oxy3835_wphase_temp',...
136                                              data.vars, 'exact'))];
137                     oxyw_dphase = [oxyw_dphase; ...
138                         data.data(:,strmatch('sci_oxy3835_wphase_dphase',...
139                                              data.vars, 'exact'))];
140                     % oxyw_bphase = [oxyw_bphase; ...
141                     %     data.data(:,strmatch('sci_oxy3835_wphase_bphase',...
142                     %                          data.vars, 'exact'))];
143                     % oxyw_rphase = [oxyw_rphase; ...
144                     %     data.data(:,strmatch('sci_oxy3835_wphase_rphase',...
145                     %                          data.vars, 'exact'))];
146                     % oxyw_bamp = [oxyw_bamp; ...
147                     %     data.data(:,strmatch('sci_oxy3835_wphase_bamp',...
148                     %                          data.vars, 'exact'))];
149                     % oxyw_bpot = [oxyw_bpot; ...
150                     %     data.data(:,strmatch('sci_oxy3835_wphase_bpot',...
151                     %                          data.vars, 'exact'))];
152                     % oxyw_ramp = [oxyw_ramp; ...
153                     %     data.data(:,strmatch('sci_oxy3835_wphase_ramp',...
154                     %                          data.vars, 'exact'))];
155                     % oxyw_rawtemp = [oxyw_rawtemp; ...
156                     %     data.data(:,strmatch('sci_oxy3835_wphase_rawtemp',...
157                     %                          data.vars, 'exact'))];
158                     % oxyw_time = [oxyw_time; ...
159                     %     data.data(:,strmatch('sci_oxy3835_wphase_timestamp',...
160                     %                          data.vars, 'exact'))];
161                     % oxyw_installed = [oxyw_installed; ...
162                     %     data.data(:,strmatch('sci_oxy3835_wphase_is_installed',...
163                     %                          data.vars, 'exact'))];
164                     ptime_ebd = [ptime_ebd; ...
165                         data.data(:,strmatch('sci_m_present_time',...
166                                              data.vars, 'exact'))];
167                 end
168
169                 data = [];
170             end 
171         end
172         %*** END READ IN EBD DATA ****
173
174         % remove nans from EBD data...
175         i = find(~isnan(oxyw_dphase));
176         oxyw_dphase = oxyw_dphase(i);
177         ptime_ebd = ptime_ebd(i);
178         oxyw_oxygen = oxyw_oxygen(i);
179         oxyw_saturation = oxyw_saturation(i);
180         oxyw_temp = oxyw_temp(i);
181         % oxyw_bphase = oxyw_bphase(i);
182         % oxyw_rphase = oxyw_rphase(i);
183         % oxyw_bamp = oxyw_bamp(i);
184         % oxyw_bpot = oxyw_bpot(i);
185         % oxyw_ramp = oxyw_ramp(i);
186         % oxyw_rawtemp = oxyw_rawtemp(i);
187         % oxyw_time = oxyw_time(i);
188         % oxyw_installed = oxyw_installed(i);
189
190         % apply the sort() function to ptime_oxy
191         % to make sure it increases monotonically...
192         [Y,I] = sort(ptime_ebd);
193         ptime_ebd = Y;
194         oxyw_oxygen = oxyw_oxygen(I);
195         oxyw_saturation = oxyw_saturation(I);
196         oxyw_temp = oxyw_temp(I);
197         oxyw_dphase = oxyw_dphase(I);
198         % oxyw_bphase = oxyw_bphase(I);
199         % oxyw_rphase = oxyw_rphase(I);
200         % oxyw_bamp = oxyw_bamp(I);
201         % oxyw_bpot = oxyw_bpot(I);
202         % oxyw_ramp = oxyw_ramp(I);
203         % oxyw_rawtemp = oxyw_rawtemp(I);
204         % oxyw_time = oxyw_time(I);
205         % oxyw_installed = oxyw_installed(I);
206
207         % load CTD data from existing mat file...
208         strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, ...
209             '_Deployment', strDeploymentNumber, '_CTD_L1.mat');
210         load(strMatFilePath);
211
212         % interpolate CTD to times of oxyw_dphase
213         tempi = interp1(ptime,temp,ptime_ebd);
214         salini = interp1(ptime,salinCorrected,ptime_ebd);
215         depthi = interp1(ptime,depth,ptime_ebd);
216
217         % first, implement temperature-dependent correction to DO concentration
218         % that utilizes dphase as the input (manual page 30), coefficients from cal
219         % sheets (a 5x4 matrix of values) and will be glider (and foil) dependent.
220         C0 = C0coef(1) + (C0coef(2) .* tempi) + ...
221             (C0coef(3) .* (tempi.^2)) + (C0coef(4) .* (tempi.^3));
222         C1 = C1coef(1) + (C1coef(2) .* tempi) + ...
223             (C1coef(3) .* (tempi.^2)) + (C1coef(4) .* (tempi.^3));
224         C2 = C2coef(1) + (C2coef(2) .* tempi) + ...
225             (C2coef(3) .* (tempi.^2)) + (C2coef(4) .* (tempi.^3));
226         C3 = C3coef(1) + (C3coef(2) .* tempi) + ...
227             (C3coef(3) .* (tempi.^2)) + (C3coef(4) .* (tempi.^3));
228         C4 = C4coef(1) + (C4coef(2) .* tempi) + ...
229             (C4coef(3) .* (tempi.^2)) + (C4coef(4) .* (tempi.^3));
230
231         o2_tcorr = C0 + (C1 .* oxyw_dphase) + (C2 .* (oxyw_dphase.^2)) + ...
232             (C3 .* (oxyw_dphase.^3)) + (C4 .* (oxyw_dphase.^4));
233
234         % second, implement the salinity correction to DO concentration (page 31)
235         % ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER
236         % SCALING IS REQUIRED
237         B0=-6.24097e-3;
238         B1=-6.93498e-3;
239         B2=-6.90358e-3;
240         B3=-4.29155e-3;
241         C0=-3.11680e-7;
242
243         Ts = log((298.15 - tempi) ./ (273.15 + tempi));
244
245         o2_tscorr = o2_tcorr .* ...
246             exp((salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + ...
247                             (B3 .* (Ts.^3))))  + ...
248             (C0 .* (salini.^2)));
249
250         % third, implement the pressure correction to DO concentration (page 32)
251         o2_tspcorr = o2_tscorr .* (1 + (0.04 .* depthi ./ 1000));
252
253         % use polynomial to calculate DO satuations using the measured temp and
254         % sal (manual page 30)
255         A0 = 2.00856;
256         A1 = 3.22400;
257         A2 = 3.99063;
258         A3 = 4.80299;
259         A4 = 9.78188e-1;
260         A5 = 1.71069;
261         B0 = -6.24097e-3;
262         B1 = -6.93498e-3;
263         B2 = -6.90358e-3;
264         B3 = -4.29155e-3;
265         C0 = -3.11680e-7;
266
267         % need interpolated temp, salinity, pressure at times of DO obs
268         rslt = (A0 + (A1 .* Ts) + (A2 .* (Ts.^2)) + (A3 .* (Ts.^3)) + ...
269                (A4 .* (Ts.^4)) + (A5 .* (Ts.^5))) + ...
270             (salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + (B3 .*(Ts.^3)))) + ...
271             (C0 .* (salini.^2));
272
273         o2_sol = exp(rslt);
274
275         o2_sat = (o2_tspcorr .* 2.2414) ./ o2_sol;
276
277         % convert ptime into datenum style...
278         ptime_ebd_datenum = (ptime_ebd/3600/24) + datenum(1970, 1, 1, 0, 0, 0);
279
280         % create configuration struct...
281         units = struct( ...
282             'oxyw_oxygen', '10e-6 mol/dm3', ...
283             'oxyw_saturation', 'percent', ...
284             'oxyw_temp', 'degrees C', ...
285             'oxyw_dphase', 'degrees', ...
286             ... % 'oxyw_bphase', 'degrees', ...
287             ... % 'oxyw_rphase', 'degrees', ...
288             ... % 'oxyw_bamp', 'mA', ...
289             ... % 'oxyw_bpot', 'mV', ...
290             ... % 'oxyw_ramp', 'mA', ...
291             ... % 'oxyw_rawtemp', 'degrees C', ...
292             ... % 'oxyw_time', 'timestamp', ...
293             ... % 'oxyw_installed', 'bool', ...
294             'ptime_ebd', 'seconds since 0000-01-01T00:00', ...
295             'ptime_ebd_datenum', 'days since 1970-01-01T00:00', ...
296             'tempi', 'degrees C', ...
297             'salini', 'psu', ...
298             'depthi', 'meters', ...
299             'o2_tcorr', '10e-6 mol/dm3', ...
300             'o2_tscorr', '10e-6 mol/dm3', ...
301             'o2_tspcorr', '10e-6 mol/dm3', ...
302             'o2_sol', 'cm3/liter at 1031 hPa', ...
303             'o2_sat', 'percent');
304        
305         variable_description = struct( ...
306             'oxyw_oxygen', 'dissolved oxygen', ...
307             'oxyw_saturation', 'dissolved oxygen saturation', ...
308             'oxyw_temp', 'water temperature', ...
309             'oxyw_dphase', 'phase difference', ...
310             ... % 'oxyw_bphase', 'blue phase', ...
311             ... % 'oxyw_rphase', 'red phase', ...
312             ... % 'oxyw_bamp', 'blue current bias', ...
313             ... % 'oxyw_bpot', 'blue voltage bias', ...
314             ... % 'oxyw_ramp', 'red current bias', ...
315             ... % 'oxyw_rawtemp', 'raw water temperature', ...
316             ... % 'oxyw_time', 'optode timestampe', ...
317             ... % 'oxyw_installed', 'bool', ...
318             'ptime_ebd', 'science computer time', ...
319             'ptime_ebd_datenum', 'science computer date', ...
320             'tempi', 'interpolated CTD water temperature', ...
321             'salini', 'interpoloted CTD salinity', ...
322             'depthi', 'interpolated CTD depth', ...
323             'o2_tcorr', 'temperature corrected dissolved oxygen', ...
324             'o2_tscorr', 'temperature and salinity corrected dissolved oxygen', ...
325             'o2_tspcorr', 'temperature, salinity, and pressure corrected dissolved oxygen', ...
326             'o2_sol', 'corrected oxygen solubility', ...
327             'o2_sat', 'corrected oxygen saturation');
328
329         correction_coefficients = struct('C0coef', C0coef, ...
330                                          'C1coef', C1coef, ...
331                                          'C2coef', C2coef, ...
332                                          'C3coef', C3coef, ...
333                                          'C4coef', C4coef);
334
335         config = struct('glider_name', strGliderName,...
336                         'deployment_number', strDeploymentNumber,...
337                         'start_date', strStartDate,...
338                         'end_date', strEndDate,...
339                         'correction_coefficients', correction_coefficients,...
340                         'var_descriptions', variable_description,...
341                         'var_units', units);
342
343         % set Level 1 data mat file name...
344         strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_DO_L1.mat');
345
346         % save flight data to mat file...
347         save(strMatFileName,...
348              'config', ...
349              'oxyw_oxygen', ...
350              'oxyw_saturation', ...
351              'oxyw_temp', ...
352              'oxyw_dphase', ...
353              ... % 'oxyw_bphase', ...
354              ... % 'oxyw_rphase', ...
355              ... % 'oxyw_bamp', ...
356              ... % 'oxyw_bpot', ...
357              ... % 'oxyw_ramp', ...
358              ... % 'oxyw_rawtemp', ...
359              ... % 'oxyw_time', ...
360              ... % 'oxyw_installed', ...
361              'ptime_ebd', ...
362              'ptime_ebd_datenum', ...
363              'tempi', ...
364              'salini', ...
365              'depthi', ...
366              'o2_tcorr', ...
367              'o2_tscorr', ...
368              'o2_tspcorr', ...
369              'o2_sol', ...
370              'o2_sat');
371
372     end
373 end
Note: See TracBrowser for help on using the browser.