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

root/gliderproc/trunk/gliderCTD_CorrectThermalLag.m

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

Initial import of Stark code.

Line 
1 %
2 %  gliderCTD_CorrectThermalLag.m
3 %
4 %  Purpose: Apply thermal lag correction to glider CTD salinity and density
5 %           and generate mat file
6 %
7 %  NOTE: At top of file, set glider index, deployment number and desired correction parameters
8 %
9 %  Requires:
10 %  correctThermalLag.m
11 %  MATLAB folder - contains util, matutil, seawater, plots, strfun, opnml
12 %
13 %
14 %  Author:  William Stark
15 %           Marine Sciences Department
16 %           UNC-Chapel Hill
17 %
18 %  Created: 21 May 2012
19 %
20 %//////////////////////////////////////////////////////////////////////////
21
22 clear all;
23
24
25 %!!!!!!  SET PRIOR TO RUNNING CODE  !!!!!!!!!!!!!!!!!!!!!!!!
26 %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
27 %
28 % SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ...
29 gliderIndex = 1;
30
31 % SET THE DEPLOYMENT NUMBER (1, 2 or 3) ...
32 deploymentNumber = 1;
33
34 % SET CORRECTION PARAMETERS STRUCTURE...
35 correctionParams = [0.1587 0.0214 6.5316 1.5969];
36
37 %
38 %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
39 %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
40
41
42
43 % add paths for required files...
44 addpath('MATLAB/util/');
45 addpath('MATLAB/matutil/');
46 addpath('MATLAB/seawater/');
47 addpath('MATLAB/plots/');
48 addpath('MATLAB/strfun/');
49 addpath('MATLAB/opnml/');
50 addpath('MATLAB/opnml/FEM/');
51
52
53
54 % glider name string...
55 if (gliderIndex==1)
56     strGliderName = 'Pelagia';
57 else
58     strGliderName = 'Ramses';
59 end
60
61 % populate arrays for the deployment start and end dates...
62 % ex. strStart(2, 3) is start date for Ramses, Deployment 3
63 strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'};
64 strEnd   = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'};
65
66 % deployment number string...
67 strDeploymentNumber = num2str(deploymentNumber);
68
69 % deployment start date string...
70 strStartDate = strStart(gliderIndex, deploymentNumber);
71
72 % deployment end date string...
73 strEndDate = strEnd(gliderIndex, deploymentNumber);
74
75 % define the path to the glider ascii files...
76 datadir = strcat('/Users/haloboy/Documents/MASC/MATLAB/CTD_data_correction/GLIDER_DATA_LEVEL0/',...
77                  strGliderName, '_Deployment', strDeploymentNumber, '/');
78
79 % define default bounds for use in plots...
80 switch gliderIndex
81     case 1  % Pelagia
82         switch deploymentNumber
83             case 1  % Deployment 1
84                 tempBounds =  [17.0 24.0];
85                 salinBounds = [36.0 36.4];
86                 densBounds =  [1025.0 1026.6];
87                 chlorBounds = [0.0 4.0];
88                
89             case 2  % Deployment 2
90                 tempBounds =  [17.0 24.0];
91                 salinBounds = [36.0 36.5];
92                 densBounds =  [1024.5 1026.8];
93                 chlorBounds = [0.0 4.0];
94                
95             case 3  % Deployment 3
96                 tempBounds =  [17.0 24.0];
97                 salinBounds = [35.9 36.7];
98                 densBounds =  [1024.4 1026.4];
99                 chlorBounds = [0.0 4.0];
100         end
101     case 2  % Ramses
102         switch deploymentNumber
103             case 1  % Deployment 1
104                 tempBounds =  [8.0 23.0];
105                 salinBounds = [35.0 36.4];
106                 densBounds =  [1024.5 1027.5];
107                 chlorBounds = [0.0 4.0];
108            
109             case 2  % Deployment 2
110                 tempBounds =  [9.0 25.0];
111                 salinBounds = [35.2 36.6];
112                 densBounds =  [1024.0 1027.5];
113                 chlorBounds = [0.0 4.0];
114            
115             case 3  % Deployment 3
116                 tempBounds =  [10.0 24.5];
117                 salinBounds = [35.3 36.7];
118                 densBounds =  [1024.4 1027.4];
119                 chlorBounds = [0.0 4.0];
120         end
121 end
122
123
124 %##########################################################################################
125    
126 %*** READ IN EBD DATA *****************************************************
127 % declare variables for storing data...
128 temp=[];
129 cond=[];
130 pres=[];
131 chlor=[];
132 ptime=[];
133 % mtime=[]; scioxy=[]; scibb=[]; scicdom=[]; scichlor=[]; scibbam=[];
134
135 % try to load all *.ebdasc files at once...
136 [files, Dstruct] = wilddir(datadir, '.ebdasc');
137 nfile = size(files, 1);
138
139 for i=1:nfile-1
140     % protect against empty ebd file
141     if(Dstruct(i).bytes>0)
142         data = read_gliderasc2([datadir, files(i,:)]);
143         %data = read_gliderasc3([datadir, files(i,:)]);
144
145         % if the number of values (in data.data) is less than the number of
146         % vars (in data.vars), this means that the data were not completely read
147         % in.  To correct this, pad data.data with NaNs until its length
148         % equals that of data.vars...
149         if (length(data.data) < length(data.vars))
150             data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post');
151         end
152
153         % populate variables with data...
154         if(~isempty(data.data))
155             temp = [temp;data.data(:,strmatch('sci_water_temp', data.vars, 'exact'))];             % temperature
156             cond = [cond;data.data(:,strmatch('sci_water_cond', data.vars, 'exact'))];             % conductivity
157             pres = [pres;data.data(:,strmatch('sci_water_pressure', data.vars, 'exact'))];         % pressure (measure of depth) science bay
158             switch gliderIndex
159                 case 1  % this is Pelagia...
160                     chlor = [chlor;data.data(:,strmatch('sci_bbfl2s_chlor_scaled', data.vars, 'exact'))];  % chlorophyll
161                 case 2  % this is Ramses...
162                     chlor = [chlor;data.data(:,strmatch('sci_flbbcd_chlor_units', data.vars, 'exact'))];  % chlorophyll
163             end
164             ptime = [ptime;data.data(:,strmatch('sci_m_present_time', data.vars, 'exact'))];       % present time
165         end
166
167         data = [];
168     end 
169 end
170 %**************************************************************************
171
172 %*** READ IN DBD DATA *****************************************************
173 % declare variables for storing data...
174 ptime_dbd=[];
175 horizontalVelocity=[];
176 depth = [];
177 pitch=[];
178 avgDepthRate = [];
179 angleOfAttack = [];
180 lat = [];
181 lon = [];
182 gpsLat = [];
183 gpsLon = [];
184 wptLat = [];
185 wptLon = [];
186
187 % try to load all *.dbdasc files at once...
188 [files, Dstruct] = wilddir(datadir, '.dbdasc');
189 nfile = size(files, 1);
190
191 clear data;
192
193 for i=1:nfile
194     % protect against empty dbd file
195     if(Dstruct(i).bytes>0)
196         data = read_gliderasc2([datadir, files(i,:)]);
197         %data = read_gliderasc3([datadir, files(i,:)]);
198
199         % if the number of values (in data.data) is less than the number of
200         % vars (in data.vars), this means that the data were not completely read
201         % in.  To correct this, pad data.data with NaNs until its length
202         % equals that of data.vars...
203         if (length(data.data) < length(data.vars))
204             data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post');
205         end
206
207         % populate variables with data...
208         if(~isempty(data.data))
209             ptime_dbd = [ptime_dbd; data.data(:,strmatch('m_present_time', data.vars, 'exact'))];               % present time
210             horizontalVelocity = [horizontalVelocity; data.data(:,strmatch('m_speed', data.vars, 'exact'))];    % horizontal glider velocity     
211             depth = [depth; data.data(:,strmatch('m_depth', data.vars, 'exact'))];                              % depth     
212             pitch = [pitch; data.data(:,strmatch('m_pitch', data.vars, 'exact'))];                              % pitch (radians)
213             avgDepthRate = [avgDepthRate; data.data(:,strmatch('m_avg_depth_rate', data.vars, 'exact'))];       % avg depth rate
214             angleOfAttack = [angleOfAttack; data.data(:,strmatch('u_angle_of_attack', data.vars, 'exact'))];    % angle of attack (radians)
215             wptLat = [wptLat; data.data(:,strmatch('c_wpt_lat', data.vars, 'exact'))];                          % Waypoint latitude
216             wptLon = [wptLon; data.data(:,strmatch('c_wpt_lon', data.vars, 'exact'))];                          % Waypoint longitude
217             gpsLat = [gpsLat; data.data(:,strmatch('m_gps_lat', data.vars, 'exact'))];                          % GPS latitude
218             gpsLon = [gpsLon; data.data(:,strmatch('m_gps_lon', data.vars, 'exact'))];                          % GPS longitude
219             lat = [lat; data.data(:,strmatch('m_lat', data.vars, 'exact'))];                                    % latitude
220             lon = [lon; data.data(:,strmatch('m_lon', data.vars, 'exact'))];                                    % longitude
221         end
222
223         data = [];
224     end
225 end
226 %**************************************************************************
227
228
229 % first, apply the sort() function to make sure that values in the time vectors
230 % (ptime and ptime_dbd) increase monotonically...
231 [Y,I] = sort(ptime);
232 ptime = Y;
233 temp = temp(I);
234 cond = cond(I);
235 pres = pres(I);
236 chlor = chlor(I);
237
238 [Y,I] = sort(ptime_dbd);
239 ptime_dbd = Y;
240 horizontalVelocity = horizontalVelocity(I);
241 depth = depth(I);
242 pitch = pitch(I);
243 avgDepthRate = avgDepthRate(I);
244 angleOfAttack = angleOfAttack(I);
245 wptLat = wptLat(I);
246 wptLon = wptLon(I);
247 gpsLat = gpsLat(I);
248 gpsLon = gpsLon(I);
249 lat = lat(I);
250 lon = lon(I);
251
252
253 % remove nans from EBD data...
254 i = find(~isnan(temp));
255 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
256 i = find(~isnan(pres));
257 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
258 i = find(~isnan(cond));
259 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
260 i = find(~isnan(chlor));
261 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
262
263 % remove conductivity values less than 1...
264 i = find(cond>=1);
265 ptime = ptime(i);  temp = temp(i);  pres = pres(i);  cond = cond(i);
266
267 % remove pressure values less than 0...
268 i = find(pres>=0);
269 ptime = ptime(i);  temp = temp(i);  pres = pres(i);  cond = cond(i);
270
271 % convert pitch and angle of attack from radians to degrees...
272 pitch = pitch*180/pi;
273 angleOfAttack = angleOfAttack*180/pi;
274
275 % compute actual glide angle = pitch + angle of attack...
276 glideAngle = pitch + angleOfAttack;
277
278 % make copy of dbd time stamp vector for use in salinity/density correction...
279 ptime1_dbd = ptime_dbd;
280
281
282 % remove nans from DBD data...
283 i = find(~isnan(horizontalVelocity));
284 ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i);
285 depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i);
286 glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
287 i = find(~isnan(depth));
288 ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i);
289 depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i);
290 glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
291 i = find(~isnan(pitch));
292 ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i);
293 depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i);
294 glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
295 i = find(~isnan(avgDepthRate));
296 ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i);
297 depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i);
298 glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
299 i = find(~isnan(glideAngle));
300 ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i);
301 depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i);
302 glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
303 i = find(~isnan(lat));
304 ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i);
305 depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i);
306 glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
307 i = find(~isnan(lon));
308 ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i);
309 depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i);
310 glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
311
312 % clip horizontalVelocity data at 0.6...
313 %i = find(horizontalVelocity~=0 & horizontalVelocity<0.6);
314 i = find(horizontalVelocity<0.6);
315 ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i);
316 depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i);
317 glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
318
319 % interpolate DBD data to align with EBD data...
320 horizontalVelocity = interp1(ptime1_dbd, horizontalVelocity, ptime);
321 depth = interp1(ptime1_dbd, depth, ptime);
322 pitch = interp1(ptime1_dbd, pitch, ptime);
323 avgDepthRate = interp1(ptime1_dbd, avgDepthRate, ptime);
324 glideAngle = interp1(ptime1_dbd, glideAngle, ptime);
325 lat = interp1(ptime1_dbd, lat, ptime);
326 lon = interp1(ptime1_dbd, lon, ptime);
327
328 % make sure there are no NaNs in the final set of data...
329 i = find(~isnan(horizontalVelocity));
330 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
331 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
332 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
333 i = find(~isnan(depth));
334 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
335 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
336 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
337 i = find(~isnan(pitch));
338 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
339 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
340 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
341 i = find(~isnan(avgDepthRate));
342 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
343 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
344 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
345 i = find(~isnan(glideAngle));
346 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
347 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
348 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
349 i = find(~isnan(lat));
350 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
351 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
352 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
353 i = find(~isnan(lon));
354 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
355 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
356 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
357 i = find(~isnan(temp));
358 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
359 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
360 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
361 i = find(~isnan(cond));
362 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
363 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
364 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
365 i = find(~isnan(pres));
366 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
367 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
368 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
369 i = find(~isnan(chlor));
370 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
371 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
372 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
373
374 % eliminate any zeros in temp, cond, pres, chlor, depth, lat and lon
375 i = find(temp~=0);
376 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
377 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
378 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
379 i = find(cond~=0);
380 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
381 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
382 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
383 i = find(pres~=0);
384 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
385 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
386 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
387 i = find(chlor~=0);
388 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
389 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
390 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
391 i = find(depth~=0);
392 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
393 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
394 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
395 i = find(lat~=0);
396 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
397 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
398 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
399 i = find(lon~=0);
400 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
401 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
402 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
403
404 % eliminate any negative chlorophyll values
405 i = find(chlor>0);
406 horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
407 avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i);
408 ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i);
409
410
411
412 % scale up the pressure...
413 pres = pres*10;
414
415 % calculate salinity (without correction)...
416 salin = sw_salt(10*cond/sw_c3515, temp, pres);
417
418 % calculate density (without correction)...
419 dens = sw_pden(salin, temp, pres, 0);
420
421 % calculate glider velocity using horizontal velocity (m_speed) and average depth rate (m_avg_depth_rate)...
422 gliderVelocity = sqrt(horizontalVelocity.^2 + avgDepthRate.^2);
423
424
425
426 % pass the correction parameters into the correctThermalLags function, which
427 % returns the corrected profile structure (with corrected temp and cond added)...
428 %
429 % profileStructure:  ptime: Present time instant at which this row was collected
430 %                    depth: Depth (pressure in decibars) measured by the CTD
431 %                    temp: Temperature measured by the CTD
432 %                    cond: Conductivity measured by the CTD
433 %                    pitch: Pitch angle of the glider (optional)
434 %
435 %
436
437 % CORRECTION SCHEME:  Pass in a glider velocity vector, so that correctThermalLag() will return
438 % profile data that has been corrected using flow speed equal to this glider velocity.
439 % Correction parameters are calculated using passed-in glider velocity.
440 profileStructure = struct('ptime', ptime, 'depth', pres, 'temp', temp, 'cond', cond, 'pitch', glideAngle);
441 [correctedProfileData, varargout] = correctThermalLag(profileStructure, correctionParams, gliderVelocity);
442
443 % get the corrected temperature...
444 tempCorrected = correctedProfileData.tempInCell;
445
446 % calculate the corrected salinity using the corrected temperature...
447 salinCorrected = sw_salt(10*cond/sw_c3515, tempCorrected, pres);
448
449 % calculate density...
450 densCorrected = sw_pden(salinCorrected, temp, pres, 0);
451
452
453 % convert ptime into datenum style...
454 ptime_datenum = ptime/3600/24+datenum(1970, 1, 1, 0, 0, 0);
455 ptime_datenum_dbd = ptime_dbd/3600/24+datenum(1970, 1, 1, 0, 0, 0);
456 ptimehrly = fix(ptime_datenum*24)/24;
457 ptimedaily = fix(ptime_datenum);
458 ptimedaily2 = unique(ptimedaily);
459 ptimedaily2 = ptimedaily2(1:2:end);
460
461
462 % make copies of dbd time stamp vector for use in lat/lon interpolation...
463 ptime_dbd_gps = ptime_dbd;
464 ptime_dbd_wpt = ptime_dbd;
465
466 % convert lats and lons to digital degrees...
467 gpsLat = ddmm2decdeg(gpsLat);
468 gpsLon = ddmm2decdeg(gpsLon);
469 wptLat = ddmm2decdeg(wptLat);
470 wptLon = ddmm2decdeg(wptLon);
471
472 % eliminate outliers in gpsLat, gpsLon...
473 i = find(abs(gpsLat) <= 90.0);
474 gpsLat = gpsLat(i);  gpsLon = gpsLon(i);  ptime_dbd_gps = ptime_dbd_gps(i);
475 i = find(abs(gpsLon) <= 180.0);
476 gpsLat = gpsLat(i);  gpsLon = gpsLon(i);  ptime_dbd_gps = ptime_dbd_gps(i);
477
478 % eliminate outliers in wptLat, wptLon...
479 i = find(abs(wptLat) <= 90.0);
480 wptLat = wptLat(i);  wptLon = wptLon(i);  ptime_dbd_wpt = ptime_dbd_wpt(i);
481 i = find(abs(wptLon) <= 180.0);
482 wptLat = wptLat(i);  wptLon = wptLon(i);  ptime_dbd_wpt = ptime_dbd_wpt(i);
483
484 % eliminate entries where wptLat==0
485 i = find(wptLat~=0);
486 wptLat = wptLat(i);  wptLon = wptLon(i);  ptime_dbd_wpt = ptime_dbd_wpt(i);
487
488 % eliminate nans before interpolating...
489 i = find(~isnan(gpsLat));
490 gpsLat = gpsLat(i);  gpsLon = gpsLon(i);  ptime_dbd_gps = ptime_dbd_gps(i);
491 i = find(~isnan(wptLat));
492 wptLat = wptLat(i);  wptLon = wptLon(i);  ptime_dbd_wpt = ptime_dbd_wpt(i);
493
494 % interpolate DBD lat/lon data to align with EBD data...
495 gpsLat = interp1(ptime_dbd_gps, gpsLat, ptime);
496 gpsLon = interp1(ptime_dbd_gps, gpsLon, ptime);
497 wptLat = interp1(ptime_dbd_wpt, wptLat, ptime);
498 wptLon = interp1(ptime_dbd_wpt, wptLon, ptime);
499
500 % use sw_dpth() to calculate depth from pres...
501 depth = sw_dpth(pres, gpsLat);
502
503
504 % create configuration struct...
505 units = struct('chlor', 'micrograms liter-1',...
506                'dens', 'kg m-3',...
507                'densCorrected', 'kg m-3',...
508                'depth', 'm',...
509                'gpsLat', 'decimal degrees',...
510                'gpsLon', 'decimal degrees',...
511                'pres', 'decibars',...
512                'ptime', 'seconds since 0000-01-01T00:00',...
513                'ptime_datenum', 'seconds since 0000-01-01T00:00',...
514                'salin', 'psu',...
515                'salinCorrected', 'psu',...
516                'temp', 'deg C');
517
518 variable_description = struct('chlor', 'chlorophyll measured by glider',...
519                               'chlorBounds', 'default chlorophyll bounds for plots',...
520                               'dens', 'density measured by glider',...
521                               'densBounds', 'default density bounds for plots',...
522                               'densCorrected', 'density corrected for thermal lag',...
523                               'depth', 'depth calculated as function of pressure and position latitude',...
524                               'gpsLat', 'position latitude measured by glider GPS',...
525                               'gpsLon', 'position longitude measured by glider GPS',...
526                               'pres', 'pressure measured by glider',...
527                               'ptime', 'time vector reported by glider',...
528                               'ptime_datenum', 'Serial Date Number string',...
529                               'salin', 'salinity measured by glider',...
530                               'salinBounds', 'default salinity bounds for plots',...
531                               'salinCorrected', 'salinity corrected for thermal lag',...
532                               'temp', 'temperature measured by glider',...
533                               'tempBounds', 'default temperature bounds for plots');
534
535 correction_parameters = struct('alpha_offset', correctionParams(1),...
536                                'alpha_slope', correctionParams(2),...
537                                'tau_offset', correctionParams(3),...
538                                'tau_slope', correctionParams(4));
539            
540 config = struct('glider_name', strGliderName,...
541                 'deployment_number', strDeploymentNumber,...
542                 'start_date', strStartDate,...
543                 'end_date', strEndDate,...
544                 'thermal_lag_correction_parameters', correction_parameters,...
545                 'var_descriptions', variable_description,...
546                 'var_units', units);
547
548 % set mat file name...
549 strMatFileName = strcat(strGliderName, '_CTD_Deployment', strDeploymentNumber, '.mat');
550
551
552
553 % save glider/deployment data to mat file...
554 save(strMatFileName,...
555      'config',...
556      'chlor',...
557      'chlorBounds',...
558      'dens',...
559      'densBounds',...
560      'densCorrected',...
561      'depth',...
562      'gpsLat',...
563      'gpsLon',...
564      'pres',...
565      'ptime',...
566      'ptime_datenum',...
567      'salin',...
568      'salinBounds',...
569      'salinCorrected',...
570      'temp',...
571      'tempBounds');
572
573
574
Note: See TracBrowser for help on using the browser.