Changeset Navigation

← Previous Changeset
Next Changeset →

Changeset 508

Timestamp:

07/03/13 11:53:18

Author:

cbc

Message:

Add Harvey's CTD QC code for Ramses.

Files:

gliderproc/trunk/gliderCTD_Generate_L1_Data.m (modified) (3 diffs)

Legend:

: Unmodified
: Added
: Removed
: Modified
: Copied
: Moved

gliderproc/trunk/gliderCTD_Generate_L1_Data.m

r498	r508
21	21	% Created: 21 May 2012
22	22	% reviewed, modified by HES December 2012, March 2013, called v1
	23	% June 2013, implemented some QC procedures for ramses, HES
	24	% 20130703 - Iterate throuth gliders and deployments. - CBC
23	25	%
24	26	%//////////////////////////////////////////////////////////////////////////
…	…
44	46	clear all;
45	47
46
47		~~%!!!!!! SET PRIOR TO RUNNING CODE !!!!!!!!!!!!!!!!!!!!!!!!~~
48		~~%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!~~
49		%
50		~~% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ...~~
51		~~gliderIndex = 2;~~
52
53		~~% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ...~~
54		~~deploymentNumber = 3;~~
55
56		~~% SET CORRECTION PARAMETERS STRUCTURE...~~
57		~~correctionParams = [0.1587 0.0214 6.5316 1.5969];~~
58
59		%
60		~~%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!~~
61		~~%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!~~
62
63
64
65	48	% add paths for required files...
66	49	addpath('MATLAB/util/');
…	…
68	51	addpath('MATLAB/seawater/');
69	52	%addpath('MATLAB/plots/');
70		%addpath('MATLAB/strfun/');
	53	addpath('MATLAB/strfun/');
71	54	%addpath('MATLAB/opnml/');
72	55	%addpath('MATLAB/opnml/FEM/');
73	56
74
75
76		% glider name string...
77		if (gliderIndex==1)
78		strGliderName = 'Pelagia';
79		else
80		strGliderName = 'Ramses';
81		end
82
83		% populate arrays for the deployment start and end dates...
84		% ex. strStart(2, 3) is start date for Ramses, Deployment 3
85		strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'};
86		strEnd = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'};
87
88		% deployment number string...
89		strDeploymentNumber = num2str(deploymentNumber);
90
91		% deployment start date string...
92		strStartDate = strStart(gliderIndex, deploymentNumber);
93
94		% deployment end date string...
95		strEndDate = strEnd(gliderIndex, deploymentNumber);
96
97		% define the path to the glider ascii files...
98		%datadir = strcat('/Users/haloboy/Documents/MASC/MATLAB/CTD_data_correction/GLIDER_CTD_DATA_LEVEL0/',...
99		datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, '_Deployment', strDeploymentNumber, '/');
100
101		% define default bounds for use in plots...
102		switch gliderIndex
103		case 1 % Pelagia
104		switch deploymentNumber
105		case 1 % Deployment 1
106		tempBounds = [17.0 24.0];
107		salinBounds = [36.0 36.4];
108		densBounds = [1025.0 1026.6];
109		chlorBounds = [0.0 4.0];
110
111		case 2 % Deployment 2
112		tempBounds = [17.0 24.0];
113		salinBounds = [36.0 36.5];
114		densBounds = [1024.5 1026.8];
115		chlorBounds = [0.0 4.0];
116
117		case 3 % Deployment 3
118		tempBounds = [17.0 24.0];
119		salinBounds = [35.9 36.7];
120		densBounds = [1024.4 1026.4];
121		chlorBounds = [0.0 4.0];
	57	% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ...
	58	for gliderIndex=1:2
	59
	60	% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ...
	61	for deploymentNumber=1:3
	62
	63	clearvars -except gliderIndex deploymentNumber;
	64
	65	% glider name string...
	66	if (gliderIndex==1)
	67	strGliderName = 'Pelagia';
	68	else
	69	strGliderName = 'Ramses';
122	70	end
123		case 2 % Ramses
124		switch deploymentNumber
125		case 1 % Deployment 1
126		tempBounds = [8.0 23.0];
127		salinBounds = [35.0 36.4];
128		densBounds = [1024.5 1027.5];
129		chlorBounds = [0.0 4.0];
130
131		case 2 % Deployment 2
132		tempBounds = [9.0 25.0];
133		salinBounds = [35.2 36.6];
134		densBounds = [1024.0 1027.5];
135		chlorBounds = [0.0 4.0];
136
137		case 3 % Deployment 3
138		tempBounds = [10.0 24.5];
139		salinBounds = [35.3 36.7];
140		densBounds = [1024.4 1027.4];
141		chlorBounds = [0.0 4.0];
	71
	72	disp(['Generating Level 1 CTD data for ', strGliderName, ' Deployment ', num2str(deploymentNumber)]);
	73
	74	% SET CORRECTION PARAMETERS STRUCTURE...
	75	correctionParams = [0.1587 0.0214 6.5316 1.5969];
	76
	77	% populate arrays for the deployment start and end dates...
	78	% ex. strStart(2, 3) is start date for Ramses, Deployment 3
	79	strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'};
	80	strEnd = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'};
	81
	82	% deployment number string...
	83	strDeploymentNumber = num2str(deploymentNumber);
	84
	85	% deployment start date string...
	86	strStartDate = strStart(gliderIndex, deploymentNumber);
	87
	88	% deployment end date string...
	89	strEndDate = strEnd(gliderIndex, deploymentNumber);
	90
	91	% define the path to the glider ascii files...
	92	%datadir = strcat('/Users/haloboy/Documents/MASC/MATLAB/CTD_data_correction/GLIDER_CTD_DATA_LEVEL0/',...
	93	datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, '_Deployment', strDeploymentNumber, '/');
	94
	95	% define default bounds for use in plots...
	96	switch gliderIndex
	97	case 1 % Pelagia
	98	switch deploymentNumber
	99	case 1 % Deployment 1
	100	tempBounds = [17.0 24.0];
	101	salinBounds = [36.0 36.4];
	102	densBounds = [1025.0 1026.6];
	103	chlorBounds = [0.0 4.0];
	104
	105	case 2 % Deployment 2
	106	tempBounds = [17.0 24.0];
	107	salinBounds = [36.0 36.5];
	108	densBounds = [1024.5 1026.8];
	109	chlorBounds = [0.0 4.0];
	110
	111	case 3 % Deployment 3
	112	tempBounds = [17.0 24.0];
	113	salinBounds = [35.9 36.7];
	114	densBounds = [1024.4 1026.4];
	115	chlorBounds = [0.0 4.0];
	116	end
	117	case 2 % Ramses
	118	switch deploymentNumber
	119	case 1 % Deployment 1
	120	tempBounds = [8.0 23.0];
	121	salinBounds = [35.0 36.4];
	122	densBounds = [1024.5 1027.5];
	123	chlorBounds = [0.0 4.0];
	124
	125	case 2 % Deployment 2
	126	tempBounds = [9.0 25.0];
	127	salinBounds = [35.2 36.6];
	128	densBounds = [1024.0 1027.5];
	129	chlorBounds = [0.0 4.0];
	130
	131	case 3 % Deployment 3
	132	tempBounds = [10.0 24.5];
	133	salinBounds = [35.3 36.7];
	134	densBounds = [1024.4 1027.4];
	135	chlorBounds = [0.0 4.0];
	136	end
142	137	end
143		end
144
145
146		%##########################################################################################
147
148
149
150
151
152		%* READ IN EBD DATA ***************************************************
153		% declare variables for storing data...
154		temp=[];
155		cond=[];
156		pres=[];
157		ctd_time=[];
158		%chlor=[];
159		ptime=[];
160		% mtime=[]; scioxy=[]; scibb=[]; scicdom=[]; scichlor=[]; scibbam=[];
161
162		% try to load all *.ebdasc files at once...
163		[files, Dstruct] = wilddir(datadir, '.ebdasc');
164		nfile = size(files, 1);
165
166		for i=1:nfile-1
167		% protect against empty ebd file
168		if(Dstruct(i).bytes>0)
169		data = read_gliderasc2([datadir, files(i,:)]);
170		%data = read_gliderasc3([datadir, files(i,:)]);
171
172		% if the number of values (in data.data) is less than the number of
173		% vars (in data.vars), this means that the data were not completely read
174		% in. To correct this, pad data.data with NaNs until its length
175		% equals that of data.vars...
176		if (length(data.data) < length(data.vars))
177		data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post');
	138
	139
	140	%##########################################################################################
	141
	142
	143
	144
	145
	146	%* READ IN EBD DATA ***************************************************
	147	% declare variables for storing data...
	148	temp=[];
	149	cond=[];
	150	pres=[];
	151	ctd_time=[];
	152	%chlor=[];
	153	ptime=[];
	154	% mtime=[]; scioxy=[]; scibb=[]; scicdom=[]; scichlor=[]; scibbam=[];
	155
	156	% try to load all *.ebdasc files at once...
	157	[files, Dstruct] = wilddir(datadir, '.ebdasc');
	158	nfile = size(files, 1);
	159
	160	for i=1:nfile-1
	161	% protect against empty ebd file
	162	if(Dstruct(i).bytes>0)
	163	data = read_gliderasc2([datadir, files(i,:)]);
	164	%data = read_gliderasc3([datadir, files(i,:)]);
	165
	166	% if the number of values (in data.data) is less than the number of
	167	% vars (in data.vars), this means that the data were not completely read
	168	% in. To correct this, pad data.data with NaNs until its length
	169	% equals that of data.vars...
	170	if (length(data.data) < length(data.vars))
	171	data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post');
	172	end
	173
	174	% populate variables with data...
	175	if(~isempty(data.data))
	176	temp = [temp;data.data(:,strmatch('sci_water_temp', data.vars, 'exact'))]; % temperature
	177	cond = [cond;data.data(:,strmatch('sci_water_cond', data.vars, 'exact'))]; % conductivity
	178	pres = [pres;data.data(:,strmatch('sci_water_pressure', data.vars, 'exact'))]; % pressure (measure of depth) science bay
	179	ctd_time = [ctd_time;data.data(:,strmatch('sci_ctd41cp_timestamp', data.vars, 'exact'))]; % ctd timestamp
	180	% switch gliderIndex
	181	% case 1 % this is Pelagia...
	182	% chlor = [chlor;data.data(:,strmatch('sci_bbfl2s_chlor_scaled', data.vars, 'exact'))]; % chlorophyll
	183	% case 2 % this is Ramses...
	184	% chlor = [chlor;data.data(:,strmatch('sci_flbbcd_chlor_units', data.vars, 'exact'))]; % chlorophyll
	185	% end
	186	ptime = [ptime;data.data(:,strmatch('sci_m_present_time', data.vars, 'exact'))]; % present time
	187	end
	188
	189	data = [];
	190	end
178	191	end
179
180		% populate variables with data...
181		if(~isempty(data.data))
182		temp = [temp;data.data(:,strmatch('sci_water_temp', data.vars, 'exact'))]; % temperature
183		cond = [cond;data.data(:,strmatch('sci_water_cond', data.vars, 'exact'))]; % conductivity
184		pres = [pres;data.data(:,strmatch('sci_water_pressure', data.vars, 'exact'))]; % pressure (measure of depth) science bay
185		ctd_time = [ctd_time;data.data(:,strmatch('sci_ctd41cp_timestamp', data.vars, 'exact'))]; % ctd timestamp
186		% switch gliderIndex
187		% case 1 % this is Pelagia...
188		% chlor = [chlor;data.data(:,strmatch('sci_bbfl2s_chlor_scaled', data.vars, 'exact'))]; % chlorophyll
189		% case 2 % this is Ramses...
190		% chlor = [chlor;data.data(:,strmatch('sci_flbbcd_chlor_units', data.vars, 'exact'))]; % chlorophyll
191		% end
192		ptime = [ptime;data.data(:,strmatch('sci_m_present_time', data.vars, 'exact'))]; % present time
	192	%**************************************************************************
	193
	194	%* READ IN DBD DATA ***************************************************
	195	% declare variables for storing data...
	196	ptime_dbd=[];
	197	horizontalVelocity=[];
	198	depth = [];
	199	pitch=[];
	200	avgDepthRate = [];
	201	angleOfAttack = [];
	202	%lat = [];
	203	%lon = [];
	204	gpsLat = [];
	205	gpsLon = [];
	206	%wptLat = [];
	207	%wptLon = [];
	208
	209	% try to load all *.dbdasc files at once...
	210	[files, Dstruct] = wilddir(datadir, '.dbdasc');
	211	nfile = size(files, 1);
	212
	213	clear data;
	214
	215	for i=1:nfile
	216	% protect against empty dbd file
	217	if(Dstruct(i).bytes>0)
	218	data = read_gliderasc2([datadir, files(i,:)]);
	219	%data = read_gliderasc3([datadir, files(i,:)]);
	220
	221	% if the number of values (in data.data) is less than the number of
	222	% vars (in data.vars), this means that the data were not completely read
	223	% in. To correct this, pad data.data with NaNs until its length
	224	% equals that of data.vars...
	225	if (length(data.data) < length(data.vars))
	226	data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post');
	227	end
	228
	229	% populate variables with data...
	230	if(~isempty(data.data))
	231	ptime_dbd = [ptime_dbd; data.data(:,strmatch('m_present_time', data.vars, 'exact'))]; % present time
	232	horizontalVelocity = [horizontalVelocity; data.data(:,strmatch('m_speed', data.vars, 'exact'))]; % horizontal glider velocity
	233	depth = [depth; data.data(:,strmatch('m_depth', data.vars, 'exact'))]; % depth
	234	pitch = [pitch; data.data(:,strmatch('m_pitch', data.vars, 'exact'))]; % pitch (radians)
	235	avgDepthRate = [avgDepthRate; data.data(:,strmatch('m_avg_depth_rate', data.vars, 'exact'))]; % avg depth rate
	236	angleOfAttack = [angleOfAttack; data.data(:,strmatch('u_angle_of_attack', data.vars, 'exact'))]; % angle of attack (radians)
	237	% wptLat = [wptLat; data.data(:,strmatch('c_wpt_lat', data.vars, 'exact'))]; % Waypoint latitude
	238	% wptLon = [wptLon; data.data(:,strmatch('c_wpt_lon', data.vars, 'exact'))]; % Waypoint longitude
	239	gpsLat = [gpsLat; data.data(:,strmatch('m_gps_lat', data.vars, 'exact'))]; % GPS latitude
	240	gpsLon = [gpsLon; data.data(:,strmatch('m_gps_lon', data.vars, 'exact'))]; % GPS longitude
	241	% lat = [lat; data.data(:,strmatch('m_lat', data.vars, 'exact'))]; % latitude
	242	% lon = [lon; data.data(:,strmatch('m_lon', data.vars, 'exact'))]; % longitude
	243	end
	244
	245	data = [];
	246	end
193	247	end
194
195		data = [];
196		end
197		end
198		%**************************************************************************
199
200		%* READ IN DBD DATA ***************************************************
201		% declare variables for storing data...
202		ptime_dbd=[];
203		horizontalVelocity=[];
204		depth = [];
205		pitch=[];
206		avgDepthRate = [];
207		angleOfAttack = [];
208		%lat = [];
209		%lon = [];
210		gpsLat = [];
211		gpsLon = [];
212		%wptLat = [];
213		%wptLon = [];
214
215		% try to load all *.dbdasc files at once...
216		[files, Dstruct] = wilddir(datadir, '.dbdasc');
217		nfile = size(files, 1);
218
219		clear data;
220
221		for i=1:nfile
222		% protect against empty dbd file
223		if(Dstruct(i).bytes>0)
224		data = read_gliderasc2([datadir, files(i,:)]);
225		%data = read_gliderasc3([datadir, files(i,:)]);
226
227		% if the number of values (in data.data) is less than the number of
228		% vars (in data.vars), this means that the data were not completely read
229		% in. To correct this, pad data.data with NaNs until its length
230		% equals that of data.vars...
231		if (length(data.data) < length(data.vars))
232		data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post');
	248	%**************************************************************************
	249
	250
	251	% first, apply the sort() function to make sure that values in the time vectors
	252	% (ptime and ptime_dbd) increase monotonically...
	253	[Y,I] = sort(ptime);
	254	ptime = Y;
	255	temp = temp(I);
	256	cond = cond(I);
	257	pres = pres(I);
	258	ctd_time=ctd_time(I);
	259
	260	[Y,I] = sort(ptime_dbd);
	261	ptime_dbd = Y;
	262	horizontalVelocity = horizontalVelocity(I);
	263	depth = depth(I);
	264	pitch = pitch(I);
	265	avgDepthRate = avgDepthRate(I);
	266	angleOfAttack = angleOfAttack(I);
	267	%wptLat = wptLat(I);
	268	%wptLon = wptLon(I);
	269	gpsLat = gpsLat(I);
	270	gpsLon = gpsLon(I);
	271	%lat = lat(I);
	272	%lon = lon(I);
	273
	274
	275	% remove ctd measurements when ptime and ctd_time are a lot different
	276	iweird=find(ptime-ctd_time > 10);
	277	ptime(iweird)=NaN; temp(iweird)=NaN; pres(iweird)=NaN; cond(iweird)=NaN;
	278	ctd_time(iweird)=NaN;
	279
	280	% remove nans from EBD data...
	281	% HES - need full triplet from CTD
	282	i = find(~isnan(temp) & ~isnan(pres) & ~isnan(cond));
	283	ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i);
	284	ctd_time = ctd_time(i);
	285
	286	% remove conductivity values less than 1, must be at surface or bad
	287	i = find(cond>=1);
	288	ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i);
	289	ctd_time = ctd_time(i);
	290
	291	% remove pressure values less than 0
	292	i = find(pres>=0);
	293	ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i);
	294	ctd_time = ctd_time(i);
	295
	296	% some QC - use diff to ID spikes in temperature and conductivity and
	297	% remove them
	298	% these values chosen from looking at ramses, may need different set for
	299	% pelagia
	300
	301	if(gliderIndex == 2) % apply only to ramses
	302	ib=find(abs(diff(temp))>1.);
	303	ib2=find(abs(diff(cond))>0.15);
	304	ibb=union(ib,ib2);
	305	temp(ibb+1)=NaN;
	306	cond(ibb+1)=NaN;
	307	i=find(~isnan(temp));
	308	ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i);
	309	ctd_time = ctd_time(i);
233	310	end
234	311
235		% populate variables with data...
236		if(~isempty(data.data))
237		ptime_dbd = [ptime_dbd; data.data(:,strmatch('m_present_time', data.vars, 'exact'))]; % present time
238		horizontalVelocity = [horizontalVelocity; data.data(:,strmatch('m_speed', data.vars, 'exact'))]; % horizontal glider velocity
239		depth = [depth; data.data(:,strmatch('m_depth', data.vars, 'exact'))]; % depth
240		pitch = [pitch; data.data(:,strmatch('m_pitch', data.vars, 'exact'))]; % pitch (radians)
241		avgDepthRate = [avgDepthRate; data.data(:,strmatch('m_avg_depth_rate', data.vars, 'exact'))]; % avg depth rate
242		angleOfAttack = [angleOfAttack; data.data(:,strmatch('u_angle_of_attack', data.vars, 'exact'))]; % angle of attack (radians)
243		% wptLat = [wptLat; data.data(:,strmatch('c_wpt_lat', data.vars, 'exact'))]; % Waypoint latitude
244		% wptLon = [wptLon; data.data(:,strmatch('c_wpt_lon', data.vars, 'exact'))]; % Waypoint longitude
245		gpsLat = [gpsLat; data.data(:,strmatch('m_gps_lat', data.vars, 'exact'))]; % GPS latitude
246		gpsLon = [gpsLon; data.data(:,strmatch('m_gps_lon', data.vars, 'exact'))]; % GPS longitude
247		% lat = [lat; data.data(:,strmatch('m_lat', data.vars, 'exact'))]; % latitude
248		% lon = [lon; data.data(:,strmatch('m_lon', data.vars, 'exact'))]; % longitude
	312	% convert pitch and angle of attack from radians to degrees...
	313	pitch = pitch*180/pi;
	314	angleOfAttack = angleOfAttack*180/pi;
	315
	316	% compute actual glide angle = pitch + angle of attack...
	317	glideAngle = pitch + angleOfAttack;
	318
	319	% make copy of dbd time stamp vector for use in salinity/density correction...
	320	ptime1_dbd = ptime_dbd;
	321
	322	% remove nans from DBD data...HES - re-wrote this to interpolate each
	323	% variable to ebd timebase using all existing values. Includes threshold
	324	% on horizontal velocity of greater than zero (really important, fair number
	325	% of values <0, not sure where these happen but think at surface) and >0.6
	326	% m/s (less certain of this, could be looked at further)
	327	i = find(~isnan(horizontalVelocity));
	328	% use hv, interpolated horizontal speed BEFORE thresholding,
	329	% for removing poor salinity after processing - still includes
	330	% crazy speeds
	331	hv = interp1(ptime1_dbd(i), horizontalVelocity(i), ctd_time);
	332
	333	i = find(~isnan(horizontalVelocity)&(horizontalVelocity>0.1 & horizontalVelocity<0.6));
	334	horizontalVelocity = interp1(ptime1_dbd(i), horizontalVelocity(i), ctd_time);
	335
	336	i = find(~isnan(depth));
	337	depth = interp1(ptime1_dbd(i), depth(i), ctd_time);
	338
	339	i = find(~isnan(pitch));
	340	pitch = interp1(ptime1_dbd(i), pitch(i), ctd_time);
	341
	342	i = find(~isnan(avgDepthRate));
	343	avgDepthRate = interp1(ptime1_dbd(i), avgDepthRate(i), ctd_time);
	344
	345	i = find(~isnan(glideAngle));
	346	glideAngle = interp1(ptime1_dbd(i), glideAngle(i), ctd_time);
	347
	348	% make sure there are no NaNs in the final set of data...HES: important for
	349	% horizontalVelocity - a start-up problem - just zaps opening points??
	350	i = find(~isnan(horizontalVelocity));
	351	horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
	352	avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i);
	353	ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i);
	354	ctd_time = ctd_time(i);
	355
	356	% scale up the pressure...
	357	pres = pres*10;
	358
	359	% calculate salinity (without correction)...
	360	salin = sw_salt(10*cond/sw_c3515, temp, pres);
	361
	362	% calculate density (without correction)...
	363	dens = sw_pden(salin, temp, pres, 0);
	364
	365	% calculate glider velocity using horizontal velocity (m_speed) and average depth rate (m_avg_depth_rate)...
	366	gliderVelocity = sqrt(horizontalVelocity.^2 + avgDepthRate.^2);
	367
	368	% pass the correction parameters into the correctThermalLags function, which
	369	% returns the corrected profile structure (with corrected temp and cond added)...
	370	%
	371	% profileStructure: ptime: Present time instant at which this row was collected
	372	% depth: Depth (pressure in decibars) measured by the CTD
	373	% temp: Temperature measured by the CTD
	374	% cond: Conductivity measured by the CTD
	375	% pitch: Pitch angle of the glider (optional)
	376	%
	377	%
	378
	379	% CORRECTION SCHEME: Pass in a glider velocity vector, so that correctThermalLag() will return
	380	% profile data that has been corrected using flow speed equal to this glider velocity.
	381	% Correction parameters are calculated using passed-in glider velocity.
	382	profileStructure = struct('ptime', ctd_time, 'depth', pres, 'temp', temp, 'cond', cond, 'pitch', glideAngle);
	383	[correctedProfileData, varargout] = correctThermalLag(profileStructure, correctionParams, gliderVelocity);
	384
	385	% get the corrected temperature...
	386	tempCorrected = correctedProfileData.tempInCell;
	387
	388	% calculate the corrected salinity using the corrected temperature...
	389	salinCorrected = sw_salt(10*cond/sw_c3515, tempCorrected, pres);
	390
	391	% implement some clean-up here? Will eliminate a lot of points. Issue is with salinities when glider is at
	392	% top or bottom of profiles. Use original velocity measure (hv) and pitch
	393	% to identify points for exclusion
	394	% these values set by look at ramses; may need alternate set for pelagia
	395
	396	if(gliderIndex == 2)
	397	iv = find(hv<0.1 \| hv > 0.7);
	398	ip = find (abs(pitch) < 10.);
	399	ib = union(iv,ip);
	400	salinCorrected(ib) = NaN;
	401
	402	% the step above likely removes many points, need to make sure that dataset
	403	% is consistent, so use salinity to ID valid times going forward
	404
	405	i = find(~isnan(salinCorrected));
	406	ptime=ptime(i); temp=temp(i); tempCorrected=tempCorrected(i); salin=salin(i);
	407	salinCorrected=salinCorrected(i); pres=pres(i); dens=dens(i);
249	408	end
250	409
251		data = [];
	410	% calculate density...should this use temp corrected?? HES - no, now have
	411	% best estimate of salinity to combine with originally measured temp
	412	densCorrected = sw_pden(salinCorrected, temp, pres, 0);
	413
	414	% convert ptime into datenum style...for plotting I think, commenting out
	415	ptime_datenum = ptime/3600/24+datenum(1970, 1, 1, 0, 0, 0);
	416	%ptime_datenum_dbd = ptime_dbd/3600/24+datenum(1970, 1, 1, 0, 0, 0);
	417	%ptimehrly = fix(ptime_datenum*24)/24;
	418	%ptimedaily = fix(ptime_datenum);
	419	%ptimedaily2 = unique(ptimedaily);
	420	%ptimedaily2 = ptimedaily2(1:2:end);
	421
	422
	423	% make copies of dbd time stamp vector for use in lat/lon interpolation...
	424	ptime_dbd_gps = ptime_dbd;
	425	%ptime_dbd_wpt = ptime_dbd;
	426
	427	% convert lats and lons to digital degrees...
	428	gpsLat = ddmm2decdeg(gpsLat);
	429	gpsLon = ddmm2decdeg(gpsLon);
	430	%wptLat = ddmm2decdeg(wptLat);
	431	%wptLon = ddmm2decdeg(wptLon);
	432	%lat = ddmm2decdeg(lat);
	433	%lon = ddmm2decdeg(lon);
	434
	435	% eliminate outliers in gpsLat, gpsLon...
	436	i = find(abs(gpsLat) <= 90.0);
	437	gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i);
	438	i = find(abs(gpsLon) <= 180.0);
	439	gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i);
	440
	441	% eliminate nans before interpolating...
	442	i = find(~isnan(gpsLat));
	443	gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i);
	444	%i = find(~isnan(wptLat));
	445	%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i);
	446
	447	% interpolate DBD lat/lon data to align with EBD data...
	448	gpsLat = interp1(ptime_dbd_gps, gpsLat, ptime);
	449	gpsLon = interp1(ptime_dbd_gps, gpsLon, ptime);
	450	%wptLat = interp1(ptime_dbd_wpt, wptLat, ptime);
	451	%wptLon = interp1(ptime_dbd_wpt, wptLon, ptime);
	452
	453	% use sw_dpth() to calculate depth from pres...
	454	depth = sw_dpth(pres, gpsLat);
	455
	456
	457	% create configuration struct...
	458	units = struct(... %'chlor', 'micrograms liter-1',...
	459	'dens', 'kg m-3',...
	460	'densCorrected', 'kg m-3',...
	461	'depth', 'm',...
	462	'gpsLat', 'decimal degrees',...
	463	'gpsLon', 'decimal degrees',...
	464	'pres', 'decibars',...
	465	'ptime', 'seconds since 0000-01-01T00:00',...
	466	'ptime_datenum', 'seconds since 0000-01-01T00:00',...
	467	'salin', 'psu',...
	468	'salinCorrected', 'psu',...
	469	'temp', 'deg C');
	470
	471	variable_description = struct(... %'chlor', 'chlorophyll measured by glider',...
	472	... %'chlorBounds', 'default chlorophyll bounds for plots',...
	473	'dens', 'density measured by glider',...
	474	'densBounds', 'default density bounds for plots',...
	475	'densCorrected', 'density corrected for thermal lag',...
	476	'depth', 'depth calculated as function of pressure and position latitude',...
	477	'gpsLat', 'position latitude measured by glider GPS',...
	478	'gpsLon', 'position longitude measured by glider GPS',...
	479	'pres', 'pressure measured by glider',...
	480	'ptime', 'time vector reported by glider',...
	481	'ptime_datenum', 'Serial Date Number string',...
	482	'salin', 'salinity measured by glider',...
	483	'salinBounds', 'default salinity bounds for plots',...
	484	'salinCorrected', 'salinity corrected for thermal lag',...
	485	'temp', 'temperature measured by glider',...
	486	'tempBounds', 'default temperature bounds for plots');
	487
	488	correction_parameters = struct('alpha_offset', correctionParams(1),...
	489	'alpha_slope', correctionParams(2),...
	490	'tau_offset', correctionParams(3),...
	491	'tau_slope', correctionParams(4));
	492
	493	config = struct('glider_name', strGliderName,...
	494	'deployment_number', strDeploymentNumber,...
	495	'start_date', strStartDate,...
	496	'end_date', strEndDate,...
	497	'thermal_lag_correction_parameters', correction_parameters,...
	498	'var_descriptions', variable_description,...
	499	'var_units', units);
	500
	501	% set Level 1 data mat file name...
	502	strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_CTD_L1.mat');
	503
	504
	505
	506	% save glider/deployment data to mat file...
	507	save(strMatFileName,...
	508	'config',...
	509	...% 'chlor',...
	510	...% 'chlorBounds',...
	511	'dens',...
	512	'densBounds',...
	513	'densCorrected',...
	514	'depth',...
	515	'gpsLat',...
	516	'gpsLon',...
	517	'pres',...
	518	'ptime',...
	519	'ptime_datenum',...
	520	'salin',...
	521	'salinBounds',...
	522	'salinCorrected',...
	523	'temp',...
	524	'tempBounds');
252	525	end
253	526	end
254		%**************************************************************************
255
256
257		% first, apply the sort() function to make sure that values in the time vectors
258		% (ptime and ptime_dbd) increase monotonically...
259		[Y,I] = sort(ptime);
260		ptime = Y;
261		temp = temp(I);
262		cond = cond(I);
263		pres = pres(I);
264		ctd_time=ctd_time(I);
265		%chlor = chlor(I);
266
267		[Y,I] = sort(ptime_dbd);
268		ptime_dbd = Y;
269		horizontalVelocity = horizontalVelocity(I);
270		depth = depth(I);
271		pitch = pitch(I);
272		avgDepthRate = avgDepthRate(I);
273		angleOfAttack = angleOfAttack(I);
274		%wptLat = wptLat(I);
275		%wptLon = wptLon(I);
276		gpsLat = gpsLat(I);
277		gpsLon = gpsLon(I);
278		%lat = lat(I);
279		%lon = lon(I);
280
281
282		% remove ctd measurements when ptime and ctd_time are a lot different
283		iweird=find(ptime-ctd_time > 10);
284		ptime(iweird)=NaN; temp(iweird)=NaN; pres(iweird)=NaN; cond(iweird)=NaN;
285		ctd_time(iweird)=NaN;
286
287		% remove nans from EBD data...
288		% HES - need full triplet from CTD
289		i = find(~isnan(temp) & ~isnan(pres) & ~isnan(cond));
290		ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i);
291		ctd_time = ctd_time(i);
292
293		% remove conductivity values less than 1, must be at surface or bad
294		i = find(cond>=1);
295		ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i);
296		ctd_time = ctd_time(i);
297
298		% remove pressure values less than 0
299		i = find(pres>=0);
300		ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i);
301		ctd_time = ctd_time(i);
302
303		% convert pitch and angle of attack from radians to degrees...
304		pitch = pitch*180/pi;
305		angleOfAttack = angleOfAttack*180/pi;
306
307		% compute actual glide angle = pitch + angle of attack...
308		glideAngle = pitch + angleOfAttack;
309
310		% make copy of dbd time stamp vector for use in salinity/density correction...
311		ptime1_dbd = ptime_dbd;
312
313		% remove nans from DBD data...HES - re-wrote this to interpolate each
314		% variable to ebd timebase using all existing values. Includes threshold
315		% on horizontal velocity of greater than zero (really important, fair number
316		% of values <0, not sure where these happen but think at surface) and >0.6
317		% m/s (less certain of this, could be looked at further)
318		i = find(~isnan(horizontalVelocity)&(horizontalVelocity>0 & horizontalVelocity<0.6));
319		horizontalVelocity = interp1(ptime1_dbd(i), horizontalVelocity(i), ctd_time);
320
321		i = find(~isnan(depth));
322		depth = interp1(ptime1_dbd(i), depth(i), ctd_time);
323
324		i = find(~isnan(pitch));
325		pitch = interp1(ptime1_dbd(i), pitch(i), ctd_time);
326
327		i = find(~isnan(avgDepthRate));
328		avgDepthRate = interp1(ptime1_dbd(i), avgDepthRate(i), ctd_time);
329
330		i = find(~isnan(glideAngle));
331		glideAngle = interp1(ptime1_dbd(i), glideAngle(i), ctd_time);
332
333		% make sure there are no NaNs in the final set of data...HES: important for
334		% horizontalVelocity - a start-up problem - just zaps opening points??
335		i = find(~isnan(horizontalVelocity));
336		horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i);
337		avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i);
338		ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i);
339		ctd_time = ctd_time(i);
340
341		% scale up the pressure...
342		pres = pres*10;
343
344		% calculate salinity (without correction)...
345		salin = sw_salt(10*cond/sw_c3515, temp, pres);
346
347		% calculate density (without correction)...
348		dens = sw_pden(salin, temp, pres, 0);
349
350		% calculate glider velocity using horizontal velocity (m_speed) and average depth rate (m_avg_depth_rate)...
351		gliderVelocity = sqrt(horizontalVelocity.^2 + avgDepthRate.^2);
352
353		% pass the correction parameters into the correctThermalLags function, which
354		% returns the corrected profile structure (with corrected temp and cond added)...
355		%
356		% profileStructure: ptime: Present time instant at which this row was collected
357		% depth: Depth (pressure in decibars) measured by the CTD
358		% temp: Temperature measured by the CTD
359		% cond: Conductivity measured by the CTD
360		% pitch: Pitch angle of the glider (optional)
361		%
362		%
363
364		% CORRECTION SCHEME: Pass in a glider velocity vector, so that correctThermalLag() will return
365		% profile data that has been corrected using flow speed equal to this glider velocity.
366		% Correction parameters are calculated using passed-in glider velocity.
367		profileStructure = struct('ptime', ctd_time, 'depth', pres, 'temp', temp, 'cond', cond, 'pitch', glideAngle);
368		[correctedProfileData, varargout] = correctThermalLag(profileStructure, correctionParams, gliderVelocity);
369
370		% get the corrected temperature...
371		tempCorrected = correctedProfileData.tempInCell;
372
373		% calculate the corrected salinity using the corrected temperature...
374		salinCorrected = sw_salt(10*cond/sw_c3515, tempCorrected, pres);
375
376		% calculate density...should this use temp corrected?? HES - no, now have
377		% best estimate of salinity to combine with originally measured temp
378		densCorrected = sw_pden(salinCorrected, temp, pres, 0);
379
380		% convert ptime into datenum style...for plotting I think, commenting out
381		ptime_datenum = ptime/3600/24+datenum(1970, 1, 1, 0, 0, 0);
382		%ptime_datenum_dbd = ptime_dbd/3600/24+datenum(1970, 1, 1, 0, 0, 0);
383		%ptimehrly = fix(ptime_datenum*24)/24;
384		%ptimedaily = fix(ptime_datenum);
385		%ptimedaily2 = unique(ptimedaily);
386		%ptimedaily2 = ptimedaily2(1:2:end);
387
388
389		% make copies of dbd time stamp vector for use in lat/lon interpolation...
390		ptime_dbd_gps = ptime_dbd;
391		%ptime_dbd_wpt = ptime_dbd;
392
393		% convert lats and lons to digital degrees...
394		gpsLat = ddmm2decdeg(gpsLat);
395		gpsLon = ddmm2decdeg(gpsLon);
396		%wptLat = ddmm2decdeg(wptLat);
397		%wptLon = ddmm2decdeg(wptLon);
398		%lat = ddmm2decdeg(lat);
399		%lon = ddmm2decdeg(lon);
400
401		% eliminate outliers in gpsLat, gpsLon...
402		i = find(abs(gpsLat) <= 90.0);
403		gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i);
404		i = find(abs(gpsLon) <= 180.0);
405		gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i);
406
407		% eliminate outliers in wptLat, wptLon...
408		%i = find(abs(wptLat) <= 90.0);
409		%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i);
410		%i = find(abs(wptLon) <= 180.0);
411		%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i);
412
413		% eliminate entries where wptLat==0
414		%i = find(wptLat~=0);
415		%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i);
416
417		% eliminate nans before interpolating...
418		i = find(~isnan(gpsLat));
419		gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i);
420		%i = find(~isnan(wptLat));
421		%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i);
422
423		% interpolate DBD lat/lon data to align with EBD data...
424		gpsLat = interp1(ptime_dbd_gps, gpsLat, ptime);
425		gpsLon = interp1(ptime_dbd_gps, gpsLon, ptime);
426		%wptLat = interp1(ptime_dbd_wpt, wptLat, ptime);
427		%wptLon = interp1(ptime_dbd_wpt, wptLon, ptime);
428
429		% use sw_dpth() to calculate depth from pres...
430		depth = sw_dpth(pres, gpsLat);
431
432
433		% create configuration struct...
434		units = struct(... %'chlor', 'micrograms liter-1',...
435		'dens', 'kg m-3',...
436		'densCorrected', 'kg m-3',...
437		'depth', 'm',...
438		'gpsLat', 'decimal degrees',...
439		'gpsLon', 'decimal degrees',...
440		'pres', 'decibars',...
441		'ptime', 'seconds since 0000-01-01T00:00',...
442		'ptime_datenum', 'seconds since 0000-01-01T00:00',...
443		'salin', 'psu',...
444		'salinCorrected', 'psu',...
445		'temp', 'deg C');
446
447		variable_description = struct(... %'chlor', 'chlorophyll measured by glider',...
448		... %'chlorBounds', 'default chlorophyll bounds for plots',...
449		'dens', 'density measured by glider',...
450		'densBounds', 'default density bounds for plots',...
451		'densCorrected', 'density corrected for thermal lag',...
452		'depth', 'depth calculated as function of pressure and position latitude',...
453		'gpsLat', 'position latitude measured by glider GPS',...
454		'gpsLon', 'position longitude measured by glider GPS',...
455		'pres', 'pressure measured by glider',...
456		'ptime', 'time vector reported by glider',...
457		'ptime_datenum', 'Serial Date Number string',...
458		'salin', 'salinity measured by glider',...
459		'salinBounds', 'default salinity bounds for plots',...
460		'salinCorrected', 'salinity corrected for thermal lag',...
461		'temp', 'temperature measured by glider',...
462		'tempBounds', 'default temperature bounds for plots');
463
464		correction_parameters = struct('alpha_offset', correctionParams(1),...
465		'alpha_slope', correctionParams(2),...
466		'tau_offset', correctionParams(3),...
467		'tau_slope', correctionParams(4));
468
469		config = struct('glider_name', strGliderName,...
470		'deployment_number', strDeploymentNumber,...
471		'start_date', strStartDate,...
472		'end_date', strEndDate,...
473		'thermal_lag_correction_parameters', correction_parameters,...
474		'var_descriptions', variable_description,...
475		'var_units', units);
476
477		% set Level 1 data mat file name...
478		strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_CTD_L1.mat');
479
480
481
482		% save glider/deployment data to mat file...
483		save(strMatFileName,...
484		'config',...
485		...% 'chlor',...
486		...% 'chlorBounds',...
487		'dens',...
488		'densBounds',...
489		'densCorrected',...
490		'depth',...
491		'gpsLat',...
492		'gpsLon',...
493		'pres',...
494		'ptime',...
495		'ptime_datenum',...
496		'salin',...
497		'salinBounds',...
498		'salinCorrected',...
499		'temp',...
500		'tempBounds');
501
502
503
	527

Changeset Navigation

Changeset 508

Legend:

gliderproc/trunk/gliderCTD_Generate_L1_Data.m

Download in other formats: