Index: gliderproc/trunk/gliderCTD_Generate_L1_Data.m =================================================================== --- gliderproc/trunk/gliderCTD_Generate_L1_Data.m (revision 498) +++ gliderproc/trunk/gliderCTD_Generate_L1_Data.m (revision 508) @@ -21,4 +21,6 @@ % Created: 21 May 2012 % reviewed, modified by HES December 2012, March 2013, called v1 +% June 2013, implemented some QC procedures for ramses, HES +% 20130703 - Iterate throuth gliders and deployments. - CBC % %////////////////////////////////////////////////////////////////////////// @@ -44,23 +46,4 @@ clear all; - -%!!!!!! SET PRIOR TO RUNNING CODE !!!!!!!!!!!!!!!!!!!!!!!! -%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -% -% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ... -gliderIndex = 2; - -% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ... -deploymentNumber = 3; - -% SET CORRECTION PARAMETERS STRUCTURE... -correctionParams = [0.1587 0.0214 6.5316 1.5969]; - -% -%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - % add paths for required files... addpath('MATLAB/util/'); @@ -68,436 +51,477 @@ addpath('MATLAB/seawater/'); %addpath('MATLAB/plots/'); -%addpath('MATLAB/strfun/'); +addpath('MATLAB/strfun/'); %addpath('MATLAB/opnml/'); %addpath('MATLAB/opnml/FEM/'); - - -% glider name string... -if (gliderIndex==1) - strGliderName = 'Pelagia'; -else - strGliderName = 'Ramses'; -end - -% populate arrays for the deployment start and end dates... -% ex. strStart(2, 3) is start date for Ramses, Deployment 3 -strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'}; -strEnd = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'}; - -% deployment number string... -strDeploymentNumber = num2str(deploymentNumber); - -% deployment start date string... -strStartDate = strStart(gliderIndex, deploymentNumber); - -% deployment end date string... -strEndDate = strEnd(gliderIndex, deploymentNumber); - -% define the path to the glider ascii files... -%datadir = strcat('/Users/haloboy/Documents/MASC/MATLAB/CTD_data_correction/GLIDER_CTD_DATA_LEVEL0/',... -datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, '_Deployment', strDeploymentNumber, '/'); - -% define default bounds for use in plots... -switch gliderIndex - case 1 % Pelagia - switch deploymentNumber - case 1 % Deployment 1 - tempBounds = [17.0 24.0]; - salinBounds = [36.0 36.4]; - densBounds = [1025.0 1026.6]; - chlorBounds = [0.0 4.0]; - - case 2 % Deployment 2 - tempBounds = [17.0 24.0]; - salinBounds = [36.0 36.5]; - densBounds = [1024.5 1026.8]; - chlorBounds = [0.0 4.0]; - - case 3 % Deployment 3 - tempBounds = [17.0 24.0]; - salinBounds = [35.9 36.7]; - densBounds = [1024.4 1026.4]; - chlorBounds = [0.0 4.0]; +% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ... +for gliderIndex=1:2 + + % SET THE DEPLOYMENT NUMBER (1, 2 or 3) ... + for deploymentNumber=1:3 + + clearvars -except gliderIndex deploymentNumber; + + % glider name string... + if (gliderIndex==1) + strGliderName = 'Pelagia'; + else + strGliderName = 'Ramses'; end - case 2 % Ramses - switch deploymentNumber - case 1 % Deployment 1 - tempBounds = [8.0 23.0]; - salinBounds = [35.0 36.4]; - densBounds = [1024.5 1027.5]; - chlorBounds = [0.0 4.0]; - - case 2 % Deployment 2 - tempBounds = [9.0 25.0]; - salinBounds = [35.2 36.6]; - densBounds = [1024.0 1027.5]; - chlorBounds = [0.0 4.0]; - - case 3 % Deployment 3 - tempBounds = [10.0 24.5]; - salinBounds = [35.3 36.7]; - densBounds = [1024.4 1027.4]; - chlorBounds = [0.0 4.0]; + + disp(['Generating Level 1 CTD data for ', strGliderName, ' Deployment ', num2str(deploymentNumber)]); + + % SET CORRECTION PARAMETERS STRUCTURE... + correctionParams = [0.1587 0.0214 6.5316 1.5969]; + + % populate arrays for the deployment start and end dates... + % ex. strStart(2, 3) is start date for Ramses, Deployment 3 + strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'}; + strEnd = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'}; + + % deployment number string... + strDeploymentNumber = num2str(deploymentNumber); + + % deployment start date string... + strStartDate = strStart(gliderIndex, deploymentNumber); + + % deployment end date string... + strEndDate = strEnd(gliderIndex, deploymentNumber); + + % define the path to the glider ascii files... + %datadir = strcat('/Users/haloboy/Documents/MASC/MATLAB/CTD_data_correction/GLIDER_CTD_DATA_LEVEL0/',... + datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, '_Deployment', strDeploymentNumber, '/'); + + % define default bounds for use in plots... + switch gliderIndex + case 1 % Pelagia + switch deploymentNumber + case 1 % Deployment 1 + tempBounds = [17.0 24.0]; + salinBounds = [36.0 36.4]; + densBounds = [1025.0 1026.6]; + chlorBounds = [0.0 4.0]; + + case 2 % Deployment 2 + tempBounds = [17.0 24.0]; + salinBounds = [36.0 36.5]; + densBounds = [1024.5 1026.8]; + chlorBounds = [0.0 4.0]; + + case 3 % Deployment 3 + tempBounds = [17.0 24.0]; + salinBounds = [35.9 36.7]; + densBounds = [1024.4 1026.4]; + chlorBounds = [0.0 4.0]; + end + case 2 % Ramses + switch deploymentNumber + case 1 % Deployment 1 + tempBounds = [8.0 23.0]; + salinBounds = [35.0 36.4]; + densBounds = [1024.5 1027.5]; + chlorBounds = [0.0 4.0]; + + case 2 % Deployment 2 + tempBounds = [9.0 25.0]; + salinBounds = [35.2 36.6]; + densBounds = [1024.0 1027.5]; + chlorBounds = [0.0 4.0]; + + case 3 % Deployment 3 + tempBounds = [10.0 24.5]; + salinBounds = [35.3 36.7]; + densBounds = [1024.4 1027.4]; + chlorBounds = [0.0 4.0]; + end end -end - - -%########################################################################################## - - - - - -%*** READ IN EBD DATA ***************************************************** -% declare variables for storing data... -temp=[]; -cond=[]; -pres=[]; -ctd_time=[]; -%chlor=[]; -ptime=[]; -% mtime=[]; scioxy=[]; scibb=[]; scicdom=[]; scichlor=[]; scibbam=[]; - -% try to load all *.ebdasc files at once... -[files, Dstruct] = wilddir(datadir, '.ebdasc'); -nfile = size(files, 1); - -for i=1:nfile-1 - % protect against empty ebd file - if(Dstruct(i).bytes>0) - data = read_gliderasc2([datadir, files(i,:)]); - %data = read_gliderasc3([datadir, files(i,:)]); - - % if the number of values (in data.data) is less than the number of - % vars (in data.vars), this means that the data were not completely read - % in. To correct this, pad data.data with NaNs until its length - % equals that of data.vars... - if (length(data.data) < length(data.vars)) - data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post'); + + + %########################################################################################## + + + + + + %*** READ IN EBD DATA ***************************************************** + % declare variables for storing data... + temp=[]; + cond=[]; + pres=[]; + ctd_time=[]; + %chlor=[]; + ptime=[]; + % mtime=[]; scioxy=[]; scibb=[]; scicdom=[]; scichlor=[]; scibbam=[]; + + % try to load all *.ebdasc files at once... + [files, Dstruct] = wilddir(datadir, '.ebdasc'); + nfile = size(files, 1); + + for i=1:nfile-1 + % protect against empty ebd file + if(Dstruct(i).bytes>0) + data = read_gliderasc2([datadir, files(i,:)]); + %data = read_gliderasc3([datadir, files(i,:)]); + + % if the number of values (in data.data) is less than the number of + % vars (in data.vars), this means that the data were not completely read + % in. To correct this, pad data.data with NaNs until its length + % equals that of data.vars... + if (length(data.data) < length(data.vars)) + data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post'); + end + + % populate variables with data... + if(~isempty(data.data)) + temp = [temp;data.data(:,strmatch('sci_water_temp', data.vars, 'exact'))]; % temperature + cond = [cond;data.data(:,strmatch('sci_water_cond', data.vars, 'exact'))]; % conductivity + pres = [pres;data.data(:,strmatch('sci_water_pressure', data.vars, 'exact'))]; % pressure (measure of depth) science bay + ctd_time = [ctd_time;data.data(:,strmatch('sci_ctd41cp_timestamp', data.vars, 'exact'))]; % ctd timestamp + % switch gliderIndex + % case 1 % this is Pelagia... + % chlor = [chlor;data.data(:,strmatch('sci_bbfl2s_chlor_scaled', data.vars, 'exact'))]; % chlorophyll + % case 2 % this is Ramses... + % chlor = [chlor;data.data(:,strmatch('sci_flbbcd_chlor_units', data.vars, 'exact'))]; % chlorophyll + % end + ptime = [ptime;data.data(:,strmatch('sci_m_present_time', data.vars, 'exact'))]; % present time + end + + data = []; + end end - - % populate variables with data... - if(~isempty(data.data)) - temp = [temp;data.data(:,strmatch('sci_water_temp', data.vars, 'exact'))]; % temperature - cond = [cond;data.data(:,strmatch('sci_water_cond', data.vars, 'exact'))]; % conductivity - pres = [pres;data.data(:,strmatch('sci_water_pressure', data.vars, 'exact'))]; % pressure (measure of depth) science bay - ctd_time = [ctd_time;data.data(:,strmatch('sci_ctd41cp_timestamp', data.vars, 'exact'))]; % ctd timestamp - % switch gliderIndex - % case 1 % this is Pelagia... - % chlor = [chlor;data.data(:,strmatch('sci_bbfl2s_chlor_scaled', data.vars, 'exact'))]; % chlorophyll - % case 2 % this is Ramses... - % chlor = [chlor;data.data(:,strmatch('sci_flbbcd_chlor_units', data.vars, 'exact'))]; % chlorophyll - % end - ptime = [ptime;data.data(:,strmatch('sci_m_present_time', data.vars, 'exact'))]; % present time + %************************************************************************** + + %*** READ IN DBD DATA ***************************************************** + % declare variables for storing data... + ptime_dbd=[]; + horizontalVelocity=[]; + depth = []; + pitch=[]; + avgDepthRate = []; + angleOfAttack = []; + %lat = []; + %lon = []; + gpsLat = []; + gpsLon = []; + %wptLat = []; + %wptLon = []; + + % try to load all *.dbdasc files at once... + [files, Dstruct] = wilddir(datadir, '.dbdasc'); + nfile = size(files, 1); + + clear data; + + for i=1:nfile + % protect against empty dbd file + if(Dstruct(i).bytes>0) + data = read_gliderasc2([datadir, files(i,:)]); + %data = read_gliderasc3([datadir, files(i,:)]); + + % if the number of values (in data.data) is less than the number of + % vars (in data.vars), this means that the data were not completely read + % in. To correct this, pad data.data with NaNs until its length + % equals that of data.vars... + if (length(data.data) < length(data.vars)) + data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post'); + end + + % populate variables with data... + if(~isempty(data.data)) + ptime_dbd = [ptime_dbd; data.data(:,strmatch('m_present_time', data.vars, 'exact'))]; % present time + horizontalVelocity = [horizontalVelocity; data.data(:,strmatch('m_speed', data.vars, 'exact'))]; % horizontal glider velocity + depth = [depth; data.data(:,strmatch('m_depth', data.vars, 'exact'))]; % depth + pitch = [pitch; data.data(:,strmatch('m_pitch', data.vars, 'exact'))]; % pitch (radians) + avgDepthRate = [avgDepthRate; data.data(:,strmatch('m_avg_depth_rate', data.vars, 'exact'))]; % avg depth rate + angleOfAttack = [angleOfAttack; data.data(:,strmatch('u_angle_of_attack', data.vars, 'exact'))]; % angle of attack (radians) + % wptLat = [wptLat; data.data(:,strmatch('c_wpt_lat', data.vars, 'exact'))]; % Waypoint latitude + % wptLon = [wptLon; data.data(:,strmatch('c_wpt_lon', data.vars, 'exact'))]; % Waypoint longitude + gpsLat = [gpsLat; data.data(:,strmatch('m_gps_lat', data.vars, 'exact'))]; % GPS latitude + gpsLon = [gpsLon; data.data(:,strmatch('m_gps_lon', data.vars, 'exact'))]; % GPS longitude + % lat = [lat; data.data(:,strmatch('m_lat', data.vars, 'exact'))]; % latitude + % lon = [lon; data.data(:,strmatch('m_lon', data.vars, 'exact'))]; % longitude + end + + data = []; + end end - - data = []; - end -end -%************************************************************************** - -%*** READ IN DBD DATA ***************************************************** -% declare variables for storing data... -ptime_dbd=[]; -horizontalVelocity=[]; -depth = []; -pitch=[]; -avgDepthRate = []; -angleOfAttack = []; -%lat = []; -%lon = []; -gpsLat = []; -gpsLon = []; -%wptLat = []; -%wptLon = []; - -% try to load all *.dbdasc files at once... -[files, Dstruct] = wilddir(datadir, '.dbdasc'); -nfile = size(files, 1); - -clear data; - -for i=1:nfile - % protect against empty dbd file - if(Dstruct(i).bytes>0) - data = read_gliderasc2([datadir, files(i,:)]); - %data = read_gliderasc3([datadir, files(i,:)]); - - % if the number of values (in data.data) is less than the number of - % vars (in data.vars), this means that the data were not completely read - % in. To correct this, pad data.data with NaNs until its length - % equals that of data.vars... - if (length(data.data) < length(data.vars)) - data.data = padarray(data.data, [0 length(data.vars)-length(data.data)], NaN, 'post'); + %************************************************************************** + + + % first, apply the sort() function to make sure that values in the time vectors + % (ptime and ptime_dbd) increase monotonically... + [Y,I] = sort(ptime); + ptime = Y; + temp = temp(I); + cond = cond(I); + pres = pres(I); + ctd_time=ctd_time(I); + + [Y,I] = sort(ptime_dbd); + ptime_dbd = Y; + horizontalVelocity = horizontalVelocity(I); + depth = depth(I); + pitch = pitch(I); + avgDepthRate = avgDepthRate(I); + angleOfAttack = angleOfAttack(I); + %wptLat = wptLat(I); + %wptLon = wptLon(I); + gpsLat = gpsLat(I); + gpsLon = gpsLon(I); + %lat = lat(I); + %lon = lon(I); + + + % remove ctd measurements when ptime and ctd_time are a lot different + iweird=find(ptime-ctd_time > 10); + ptime(iweird)=NaN; temp(iweird)=NaN; pres(iweird)=NaN; cond(iweird)=NaN; + ctd_time(iweird)=NaN; + + % remove nans from EBD data... + % HES - need full triplet from CTD + i = find(~isnan(temp) & ~isnan(pres) & ~isnan(cond)); + ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); + ctd_time = ctd_time(i); + + % remove conductivity values less than 1, must be at surface or bad + i = find(cond>=1); + ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i); + ctd_time = ctd_time(i); + + % remove pressure values less than 0 + i = find(pres>=0); + ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i); + ctd_time = ctd_time(i); + + % some QC - use diff to ID spikes in temperature and conductivity and + % remove them + % these values chosen from looking at ramses, may need different set for + % pelagia + + if(gliderIndex == 2) % apply only to ramses + ib=find(abs(diff(temp))>1.); + ib2=find(abs(diff(cond))>0.15); + ibb=union(ib,ib2); + temp(ibb+1)=NaN; + cond(ibb+1)=NaN; + i=find(~isnan(temp)); + ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i); + ctd_time = ctd_time(i); end - % populate variables with data... - if(~isempty(data.data)) - ptime_dbd = [ptime_dbd; data.data(:,strmatch('m_present_time', data.vars, 'exact'))]; % present time - horizontalVelocity = [horizontalVelocity; data.data(:,strmatch('m_speed', data.vars, 'exact'))]; % horizontal glider velocity - depth = [depth; data.data(:,strmatch('m_depth', data.vars, 'exact'))]; % depth - pitch = [pitch; data.data(:,strmatch('m_pitch', data.vars, 'exact'))]; % pitch (radians) - avgDepthRate = [avgDepthRate; data.data(:,strmatch('m_avg_depth_rate', data.vars, 'exact'))]; % avg depth rate - angleOfAttack = [angleOfAttack; data.data(:,strmatch('u_angle_of_attack', data.vars, 'exact'))]; % angle of attack (radians) - % wptLat = [wptLat; data.data(:,strmatch('c_wpt_lat', data.vars, 'exact'))]; % Waypoint latitude - % wptLon = [wptLon; data.data(:,strmatch('c_wpt_lon', data.vars, 'exact'))]; % Waypoint longitude - gpsLat = [gpsLat; data.data(:,strmatch('m_gps_lat', data.vars, 'exact'))]; % GPS latitude - gpsLon = [gpsLon; data.data(:,strmatch('m_gps_lon', data.vars, 'exact'))]; % GPS longitude - % lat = [lat; data.data(:,strmatch('m_lat', data.vars, 'exact'))]; % latitude - % lon = [lon; data.data(:,strmatch('m_lon', data.vars, 'exact'))]; % longitude + % convert pitch and angle of attack from radians to degrees... + pitch = pitch*180/pi; + angleOfAttack = angleOfAttack*180/pi; + + % compute actual glide angle = pitch + angle of attack... + glideAngle = pitch + angleOfAttack; + + % make copy of dbd time stamp vector for use in salinity/density correction... + ptime1_dbd = ptime_dbd; + + % remove nans from DBD data...HES - re-wrote this to interpolate each + % variable to ebd timebase using all existing values. Includes threshold + % on horizontal velocity of greater than zero (really important, fair number + % of values <0, not sure where these happen but think at surface) and >0.6 + % m/s (less certain of this, could be looked at further) + i = find(~isnan(horizontalVelocity)); + % use hv, interpolated horizontal speed BEFORE thresholding, + % for removing poor salinity after processing - still includes + % crazy speeds + hv = interp1(ptime1_dbd(i), horizontalVelocity(i), ctd_time); + + i = find(~isnan(horizontalVelocity)&(horizontalVelocity>0.1 & horizontalVelocity<0.6)); + horizontalVelocity = interp1(ptime1_dbd(i), horizontalVelocity(i), ctd_time); + + i = find(~isnan(depth)); + depth = interp1(ptime1_dbd(i), depth(i), ctd_time); + + i = find(~isnan(pitch)); + pitch = interp1(ptime1_dbd(i), pitch(i), ctd_time); + + i = find(~isnan(avgDepthRate)); + avgDepthRate = interp1(ptime1_dbd(i), avgDepthRate(i), ctd_time); + + i = find(~isnan(glideAngle)); + glideAngle = interp1(ptime1_dbd(i), glideAngle(i), ctd_time); + + % make sure there are no NaNs in the final set of data...HES: important for + % horizontalVelocity - a start-up problem - just zaps opening points?? + i = find(~isnan(horizontalVelocity)); + horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); + avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); + ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); + ctd_time = ctd_time(i); + + % scale up the pressure... + pres = pres*10; + + % calculate salinity (without correction)... + salin = sw_salt(10*cond/sw_c3515, temp, pres); + + % calculate density (without correction)... + dens = sw_pden(salin, temp, pres, 0); + + % calculate glider velocity using horizontal velocity (m_speed) and average depth rate (m_avg_depth_rate)... + gliderVelocity = sqrt(horizontalVelocity.^2 + avgDepthRate.^2); + + % pass the correction parameters into the correctThermalLags function, which + % returns the corrected profile structure (with corrected temp and cond added)... + % + % profileStructure: ptime: Present time instant at which this row was collected + % depth: Depth (pressure in decibars) measured by the CTD + % temp: Temperature measured by the CTD + % cond: Conductivity measured by the CTD + % pitch: Pitch angle of the glider (optional) + % + % + + % CORRECTION SCHEME: Pass in a glider velocity vector, so that correctThermalLag() will return + % profile data that has been corrected using flow speed equal to this glider velocity. + % Correction parameters are calculated using passed-in glider velocity. + profileStructure = struct('ptime', ctd_time, 'depth', pres, 'temp', temp, 'cond', cond, 'pitch', glideAngle); + [correctedProfileData, varargout] = correctThermalLag(profileStructure, correctionParams, gliderVelocity); + + % get the corrected temperature... + tempCorrected = correctedProfileData.tempInCell; + + % calculate the corrected salinity using the corrected temperature... + salinCorrected = sw_salt(10*cond/sw_c3515, tempCorrected, pres); + + % implement some clean-up here? Will eliminate a lot of points. Issue is with salinities when glider is at + % top or bottom of profiles. Use original velocity measure (hv) and pitch + % to identify points for exclusion + % these values set by look at ramses; may need alternate set for pelagia + + if(gliderIndex == 2) + iv = find(hv<0.1 | hv > 0.7); + ip = find (abs(pitch) < 10.); + ib = union(iv,ip); + salinCorrected(ib) = NaN; + + % the step above likely removes many points, need to make sure that dataset + % is consistent, so use salinity to ID valid times going forward + + i = find(~isnan(salinCorrected)); + ptime=ptime(i); temp=temp(i); tempCorrected=tempCorrected(i); salin=salin(i); + salinCorrected=salinCorrected(i); pres=pres(i); dens=dens(i); end - data = []; + % calculate density...should this use temp corrected?? HES - no, now have + % best estimate of salinity to combine with originally measured temp + densCorrected = sw_pden(salinCorrected, temp, pres, 0); + + % convert ptime into datenum style...for plotting I think, commenting out + ptime_datenum = ptime/3600/24+datenum(1970, 1, 1, 0, 0, 0); + %ptime_datenum_dbd = ptime_dbd/3600/24+datenum(1970, 1, 1, 0, 0, 0); + %ptimehrly = fix(ptime_datenum*24)/24; + %ptimedaily = fix(ptime_datenum); + %ptimedaily2 = unique(ptimedaily); + %ptimedaily2 = ptimedaily2(1:2:end); + + + % make copies of dbd time stamp vector for use in lat/lon interpolation... + ptime_dbd_gps = ptime_dbd; + %ptime_dbd_wpt = ptime_dbd; + + % convert lats and lons to digital degrees... + gpsLat = ddmm2decdeg(gpsLat); + gpsLon = ddmm2decdeg(gpsLon); + %wptLat = ddmm2decdeg(wptLat); + %wptLon = ddmm2decdeg(wptLon); + %lat = ddmm2decdeg(lat); + %lon = ddmm2decdeg(lon); + + % eliminate outliers in gpsLat, gpsLon... + i = find(abs(gpsLat) <= 90.0); + gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i); + i = find(abs(gpsLon) <= 180.0); + gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i); + + % eliminate nans before interpolating... + i = find(~isnan(gpsLat)); + gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i); + %i = find(~isnan(wptLat)); + %wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i); + + % interpolate DBD lat/lon data to align with EBD data... + gpsLat = interp1(ptime_dbd_gps, gpsLat, ptime); + gpsLon = interp1(ptime_dbd_gps, gpsLon, ptime); + %wptLat = interp1(ptime_dbd_wpt, wptLat, ptime); + %wptLon = interp1(ptime_dbd_wpt, wptLon, ptime); + + % use sw_dpth() to calculate depth from pres... + depth = sw_dpth(pres, gpsLat); + + + % create configuration struct... + units = struct(... %'chlor', 'micrograms liter-1',... + 'dens', 'kg m-3',... + 'densCorrected', 'kg m-3',... + 'depth', 'm',... + 'gpsLat', 'decimal degrees',... + 'gpsLon', 'decimal degrees',... + 'pres', 'decibars',... + 'ptime', 'seconds since 0000-01-01T00:00',... + 'ptime_datenum', 'seconds since 0000-01-01T00:00',... + 'salin', 'psu',... + 'salinCorrected', 'psu',... + 'temp', 'deg C'); + + variable_description = struct(... %'chlor', 'chlorophyll measured by glider',... + ... %'chlorBounds', 'default chlorophyll bounds for plots',... + 'dens', 'density measured by glider',... + 'densBounds', 'default density bounds for plots',... + 'densCorrected', 'density corrected for thermal lag',... + 'depth', 'depth calculated as function of pressure and position latitude',... + 'gpsLat', 'position latitude measured by glider GPS',... + 'gpsLon', 'position longitude measured by glider GPS',... + 'pres', 'pressure measured by glider',... + 'ptime', 'time vector reported by glider',... + 'ptime_datenum', 'Serial Date Number string',... + 'salin', 'salinity measured by glider',... + 'salinBounds', 'default salinity bounds for plots',... + 'salinCorrected', 'salinity corrected for thermal lag',... + 'temp', 'temperature measured by glider',... + 'tempBounds', 'default temperature bounds for plots'); + + correction_parameters = struct('alpha_offset', correctionParams(1),... + 'alpha_slope', correctionParams(2),... + 'tau_offset', correctionParams(3),... + 'tau_slope', correctionParams(4)); + + config = struct('glider_name', strGliderName,... + 'deployment_number', strDeploymentNumber,... + 'start_date', strStartDate,... + 'end_date', strEndDate,... + 'thermal_lag_correction_parameters', correction_parameters,... + 'var_descriptions', variable_description,... + 'var_units', units); + + % set Level 1 data mat file name... + strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_CTD_L1.mat'); + + + + % save glider/deployment data to mat file... + save(strMatFileName,... + 'config',... + ...% 'chlor',... + ...% 'chlorBounds',... + 'dens',... + 'densBounds',... + 'densCorrected',... + 'depth',... + 'gpsLat',... + 'gpsLon',... + 'pres',... + 'ptime',... + 'ptime_datenum',... + 'salin',... + 'salinBounds',... + 'salinCorrected',... + 'temp',... + 'tempBounds'); end end -%************************************************************************** - - -% first, apply the sort() function to make sure that values in the time vectors -% (ptime and ptime_dbd) increase monotonically... -[Y,I] = sort(ptime); -ptime = Y; -temp = temp(I); -cond = cond(I); -pres = pres(I); -ctd_time=ctd_time(I); -%chlor = chlor(I); - -[Y,I] = sort(ptime_dbd); -ptime_dbd = Y; -horizontalVelocity = horizontalVelocity(I); -depth = depth(I); -pitch = pitch(I); -avgDepthRate = avgDepthRate(I); -angleOfAttack = angleOfAttack(I); -%wptLat = wptLat(I); -%wptLon = wptLon(I); -gpsLat = gpsLat(I); -gpsLon = gpsLon(I); -%lat = lat(I); -%lon = lon(I); - - -% remove ctd measurements when ptime and ctd_time are a lot different -iweird=find(ptime-ctd_time > 10); -ptime(iweird)=NaN; temp(iweird)=NaN; pres(iweird)=NaN; cond(iweird)=NaN; -ctd_time(iweird)=NaN; - -% remove nans from EBD data... -% HES - need full triplet from CTD -i = find(~isnan(temp) & ~isnan(pres) & ~isnan(cond)); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); -ctd_time = ctd_time(i); - -% remove conductivity values less than 1, must be at surface or bad -i = find(cond>=1); -ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i); -ctd_time = ctd_time(i); - -% remove pressure values less than 0 -i = find(pres>=0); -ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i); -ctd_time = ctd_time(i); - -% convert pitch and angle of attack from radians to degrees... -pitch = pitch*180/pi; -angleOfAttack = angleOfAttack*180/pi; - -% compute actual glide angle = pitch + angle of attack... -glideAngle = pitch + angleOfAttack; - -% make copy of dbd time stamp vector for use in salinity/density correction... -ptime1_dbd = ptime_dbd; - -% remove nans from DBD data...HES - re-wrote this to interpolate each -% variable to ebd timebase using all existing values. Includes threshold -% on horizontal velocity of greater than zero (really important, fair number -% of values <0, not sure where these happen but think at surface) and >0.6 -% m/s (less certain of this, could be looked at further) -i = find(~isnan(horizontalVelocity)&(horizontalVelocity>0 & horizontalVelocity<0.6)); -horizontalVelocity = interp1(ptime1_dbd(i), horizontalVelocity(i), ctd_time); - -i = find(~isnan(depth)); -depth = interp1(ptime1_dbd(i), depth(i), ctd_time); - -i = find(~isnan(pitch)); -pitch = interp1(ptime1_dbd(i), pitch(i), ctd_time); - -i = find(~isnan(avgDepthRate)); -avgDepthRate = interp1(ptime1_dbd(i), avgDepthRate(i), ctd_time); - -i = find(~isnan(glideAngle)); -glideAngle = interp1(ptime1_dbd(i), glideAngle(i), ctd_time); - -% make sure there are no NaNs in the final set of data...HES: important for -% horizontalVelocity - a start-up problem - just zaps opening points?? -i = find(~isnan(horizontalVelocity)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); -ctd_time = ctd_time(i); - -% scale up the pressure... -pres = pres*10; - -% calculate salinity (without correction)... -salin = sw_salt(10*cond/sw_c3515, temp, pres); - -% calculate density (without correction)... -dens = sw_pden(salin, temp, pres, 0); - -% calculate glider velocity using horizontal velocity (m_speed) and average depth rate (m_avg_depth_rate)... -gliderVelocity = sqrt(horizontalVelocity.^2 + avgDepthRate.^2); - -% pass the correction parameters into the correctThermalLags function, which -% returns the corrected profile structure (with corrected temp and cond added)... -% -% profileStructure: ptime: Present time instant at which this row was collected -% depth: Depth (pressure in decibars) measured by the CTD -% temp: Temperature measured by the CTD -% cond: Conductivity measured by the CTD -% pitch: Pitch angle of the glider (optional) -% -% - -% CORRECTION SCHEME: Pass in a glider velocity vector, so that correctThermalLag() will return -% profile data that has been corrected using flow speed equal to this glider velocity. -% Correction parameters are calculated using passed-in glider velocity. -profileStructure = struct('ptime', ctd_time, 'depth', pres, 'temp', temp, 'cond', cond, 'pitch', glideAngle); -[correctedProfileData, varargout] = correctThermalLag(profileStructure, correctionParams, gliderVelocity); - -% get the corrected temperature... -tempCorrected = correctedProfileData.tempInCell; - -% calculate the corrected salinity using the corrected temperature... -salinCorrected = sw_salt(10*cond/sw_c3515, tempCorrected, pres); - -% calculate density...should this use temp corrected?? HES - no, now have -% best estimate of salinity to combine with originally measured temp -densCorrected = sw_pden(salinCorrected, temp, pres, 0); - -% convert ptime into datenum style...for plotting I think, commenting out -ptime_datenum = ptime/3600/24+datenum(1970, 1, 1, 0, 0, 0); -%ptime_datenum_dbd = ptime_dbd/3600/24+datenum(1970, 1, 1, 0, 0, 0); -%ptimehrly = fix(ptime_datenum*24)/24; -%ptimedaily = fix(ptime_datenum); -%ptimedaily2 = unique(ptimedaily); -%ptimedaily2 = ptimedaily2(1:2:end); - - -% make copies of dbd time stamp vector for use in lat/lon interpolation... -ptime_dbd_gps = ptime_dbd; -%ptime_dbd_wpt = ptime_dbd; - -% convert lats and lons to digital degrees... -gpsLat = ddmm2decdeg(gpsLat); -gpsLon = ddmm2decdeg(gpsLon); -%wptLat = ddmm2decdeg(wptLat); -%wptLon = ddmm2decdeg(wptLon); -%lat = ddmm2decdeg(lat); -%lon = ddmm2decdeg(lon); - -% eliminate outliers in gpsLat, gpsLon... -i = find(abs(gpsLat) <= 90.0); -gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i); -i = find(abs(gpsLon) <= 180.0); -gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i); - -% eliminate outliers in wptLat, wptLon... -%i = find(abs(wptLat) <= 90.0); -%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i); -%i = find(abs(wptLon) <= 180.0); -%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i); - -% eliminate entries where wptLat==0 -%i = find(wptLat~=0); -%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i); - -% eliminate nans before interpolating... -i = find(~isnan(gpsLat)); -gpsLat = gpsLat(i); gpsLon = gpsLon(i); ptime_dbd_gps = ptime_dbd_gps(i); -%i = find(~isnan(wptLat)); -%wptLat = wptLat(i); wptLon = wptLon(i); ptime_dbd_wpt = ptime_dbd_wpt(i); - -% interpolate DBD lat/lon data to align with EBD data... -gpsLat = interp1(ptime_dbd_gps, gpsLat, ptime); -gpsLon = interp1(ptime_dbd_gps, gpsLon, ptime); -%wptLat = interp1(ptime_dbd_wpt, wptLat, ptime); -%wptLon = interp1(ptime_dbd_wpt, wptLon, ptime); - -% use sw_dpth() to calculate depth from pres... -depth = sw_dpth(pres, gpsLat); - - -% create configuration struct... -units = struct(... %'chlor', 'micrograms liter-1',... - 'dens', 'kg m-3',... - 'densCorrected', 'kg m-3',... - 'depth', 'm',... - 'gpsLat', 'decimal degrees',... - 'gpsLon', 'decimal degrees',... - 'pres', 'decibars',... - 'ptime', 'seconds since 0000-01-01T00:00',... - 'ptime_datenum', 'seconds since 0000-01-01T00:00',... - 'salin', 'psu',... - 'salinCorrected', 'psu',... - 'temp', 'deg C'); - -variable_description = struct(... %'chlor', 'chlorophyll measured by glider',... - ... %'chlorBounds', 'default chlorophyll bounds for plots',... - 'dens', 'density measured by glider',... - 'densBounds', 'default density bounds for plots',... - 'densCorrected', 'density corrected for thermal lag',... - 'depth', 'depth calculated as function of pressure and position latitude',... - 'gpsLat', 'position latitude measured by glider GPS',... - 'gpsLon', 'position longitude measured by glider GPS',... - 'pres', 'pressure measured by glider',... - 'ptime', 'time vector reported by glider',... - 'ptime_datenum', 'Serial Date Number string',... - 'salin', 'salinity measured by glider',... - 'salinBounds', 'default salinity bounds for plots',... - 'salinCorrected', 'salinity corrected for thermal lag',... - 'temp', 'temperature measured by glider',... - 'tempBounds', 'default temperature bounds for plots'); - -correction_parameters = struct('alpha_offset', correctionParams(1),... - 'alpha_slope', correctionParams(2),... - 'tau_offset', correctionParams(3),... - 'tau_slope', correctionParams(4)); - -config = struct('glider_name', strGliderName,... - 'deployment_number', strDeploymentNumber,... - 'start_date', strStartDate,... - 'end_date', strEndDate,... - 'thermal_lag_correction_parameters', correction_parameters,... - 'var_descriptions', variable_description,... - 'var_units', units); - -% set Level 1 data mat file name... -strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_CTD_L1.mat'); - - - -% save glider/deployment data to mat file... -save(strMatFileName,... - 'config',... - ...% 'chlor',... - ...% 'chlorBounds',... - 'dens',... - 'densBounds',... - 'densCorrected',... - 'depth',... - 'gpsLat',... - 'gpsLon',... - 'pres',... - 'ptime',... - 'ptime_datenum',... - 'salin',... - 'salinBounds',... - 'salinCorrected',... - 'temp',... - 'tempBounds'); - - - +