Index: gliderproc/trunk/gliderCTD_2DPlots.m =================================================================== --- gliderproc/trunk/gliderCTD_2DPlots.m (revision 495) +++ gliderproc/trunk/gliderCTD_2DPlots.m (revision 496) @@ -59,12 +59,4 @@ ylabel('Temperature (°C)', 'fontsize', 18); xlabel('Salinity (psu)', 'fontsize', 18); -% strLegend = []; -% for i=2:lastLeg -% if i==2 -% strLegend = ['Leg ', i-1]; -% else -% strLegend = strcat(strLegend, [', Leg ', i-1]); -% end -% legend(strLegend); title([strGliderName, ' - All legs for ', strDeployment, ' '], 'fontsize', 20, 'fontweight', 'bold'); x = [ax(1):0.1:ax(2)]; Index: gliderproc/trunk/gliderCTD_3DPlots.m =================================================================== --- gliderproc/trunk/gliderCTD_3DPlots.m (revision 495) +++ gliderproc/trunk/gliderCTD_3DPlots.m (revision 496) @@ -128,9 +128,9 @@ hold on; % vertical marker lines at mooring locations - plot3(LB1_X, LB1_Y, LB1_Z, 'k', 'linewidth', [5]); - hold on; - plot3(LB2_X, LB2_Y, LB2_Z, 'k', 'linewidth', [5]); - hold on; - plot3(LB3_X, LB3_Y, LB3_Z, 'k', 'linewidth', [5]); + plot3(LB1_X, LB1_Y, LB1_Z, 'k', 'linewidth', 5); + hold on; + plot3(LB2_X, LB2_Y, LB2_Z, 'k', 'linewidth', 5); + hold on; + plot3(LB3_X, LB3_Y, LB3_Z, 'k', 'linewidth', 5); hold on; % surface marker symbol at mooring locations @@ -191,9 +191,9 @@ hold on; % vertical marker lines at mooring locations - plot3(LB1_X, LB1_Y, LB1_Z, 'k', 'linewidth', [5]); - hold on; - plot3(LB2_X, LB2_Y, LB2_Z, 'k', 'linewidth', [5]); - hold on; - plot3(LB3_X, LB3_Y, LB3_Z, 'k', 'linewidth', [5]); + plot3(LB1_X, LB1_Y, LB1_Z, 'k', 'linewidth', 5); + hold on; + plot3(LB2_X, LB2_Y, LB2_Z, 'k', 'linewidth', 5); + hold on; + plot3(LB3_X, LB3_Y, LB3_Z, 'k', 'linewidth', 5); hold on; % surface marker symbol at mooring locations @@ -236,9 +236,9 @@ hold on; % vertical marker lines at mooring locations - plot3(LB1_X, LB1_Y, LB1_Z, 'k', 'linewidth', [5]); - hold on; - plot3(LB2_X, LB2_Y, LB2_Z, 'k', 'linewidth', [5]); - hold on; - plot3(LB3_X, LB3_Y, LB3_Z, 'k', 'linewidth', [5]); + plot3(LB1_X, LB1_Y, LB1_Z, 'k', 'linewidth', 5); + hold on; + plot3(LB2_X, LB2_Y, LB2_Z, 'k', 'linewidth', 5); + hold on; + plot3(LB3_X, LB3_Y, LB3_Z, 'k', 'linewidth', 5); hold on; % surface marker symbol at mooring locations @@ -281,9 +281,9 @@ hold on; % vertical marker lines at mooring locations - plot3(LB1_X, LB1_Y, LB1_Z, 'k', 'linewidth', [5]); - hold on; - plot3(LB2_X, LB2_Y, LB2_Z, 'k', 'linewidth', [5]); - hold on; - plot3(LB3_X, LB3_Y, LB3_Z, 'k', 'linewidth', [5]); + plot3(LB1_X, LB1_Y, LB1_Z, 'k', 'linewidth', 5); + hold on; + plot3(LB2_X, LB2_Y, LB2_Z, 'k', 'linewidth', 5); + hold on; + plot3(LB3_X, LB3_Y, LB3_Z, 'k', 'linewidth', 5); hold on; % surface marker symbol at mooring locations Index: gliderproc/trunk/gliderCTD_CorrectThermalLag.m =================================================================== --- gliderproc/trunk/gliderCTD_CorrectThermalLag.m (revision 495) +++ (revision ) @@ -1,574 +1,0 @@ -% -% gliderCTD_CorrectThermalLag.m -% -% Purpose: Apply thermal lag correction to glider CTD salinity and density -% and generate mat file -% -% NOTE: At top of file, set glider index, deployment number and desired correction parameters -% -% Requires: -% correctThermalLag.m -% MATLAB folder - contains util, matutil, seawater, plots, strfun, opnml -% -% -% Author: William Stark -% Marine Sciences Department -% UNC-Chapel Hill -% -% Created: 21 May 2012 -% -%////////////////////////////////////////////////////////////////////////// - -clear all; - - -%!!!!!! SET PRIOR TO RUNNING CODE !!!!!!!!!!!!!!!!!!!!!!!! -%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -% -% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ... -gliderIndex = 1; - -% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ... -deploymentNumber = 1; - -% SET CORRECTION PARAMETERS STRUCTURE... -correctionParams = [0.1587 0.0214 6.5316 1.5969]; - -% -%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - -% add paths for required files... -addpath('MATLAB/util/'); -addpath('MATLAB/matutil/'); -addpath('MATLAB/seawater/'); -addpath('MATLAB/plots/'); -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_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 - - -%########################################################################################## - -%*** READ IN EBD DATA ***************************************************** -% declare variables for storing data... -temp=[]; -cond=[]; -pres=[]; -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 - 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 -%************************************************************************** - -%*** 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 -%************************************************************************** - - -% 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); -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 nans from EBD data... -i = find(~isnan(temp)); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(pres)); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(cond)); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(chlor)); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); - -% remove conductivity values less than 1... -i = find(cond>=1); -ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(i); - -% remove pressure values less than 0... -i = find(pres>=0); -ptime = ptime(i); temp = temp(i); pres = pres(i); cond = cond(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... -i = find(~isnan(horizontalVelocity)); -ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i); -depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i); -glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -i = find(~isnan(depth)); -ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i); -depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i); -glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -i = find(~isnan(pitch)); -ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i); -depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i); -glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -i = find(~isnan(avgDepthRate)); -ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i); -depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i); -glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -i = find(~isnan(glideAngle)); -ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i); -depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i); -glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -i = find(~isnan(lat)); -ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i); -depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i); -glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -i = find(~isnan(lon)); -ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i); -depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i); -glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); - -% clip horizontalVelocity data at 0.6... -%i = find(horizontalVelocity~=0 & horizontalVelocity<0.6); -i = find(horizontalVelocity<0.6); -ptime1_dbd = ptime1_dbd(i); horizontalVelocity = horizontalVelocity(i); -depth = depth(i); pitch = pitch(i); avgDepthRate = avgDepthRate(i); -glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); - -% interpolate DBD data to align with EBD data... -horizontalVelocity = interp1(ptime1_dbd, horizontalVelocity, ptime); -depth = interp1(ptime1_dbd, depth, ptime); -pitch = interp1(ptime1_dbd, pitch, ptime); -avgDepthRate = interp1(ptime1_dbd, avgDepthRate, ptime); -glideAngle = interp1(ptime1_dbd, glideAngle, ptime); -lat = interp1(ptime1_dbd, lat, ptime); -lon = interp1(ptime1_dbd, lon, ptime); - -% make sure there are no NaNs in the final set of data... -i = find(~isnan(horizontalVelocity)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(depth)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(pitch)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(avgDepthRate)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(glideAngle)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(lat)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(lon)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(temp)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(cond)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(pres)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(~isnan(chlor)); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); - -% eliminate any zeros in temp, cond, pres, chlor, depth, lat and lon -i = find(temp~=0); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(cond~=0); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(pres~=0); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(chlor~=0); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(depth~=0); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(lat~=0); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); -i = find(lon~=0); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(i); - -% eliminate any negative chlorophyll values -i = find(chlor>0); -horizontalVelocity = horizontalVelocity(i); depth = depth(i); pitch = pitch(i); -avgDepthRate = avgDepthRate(i); glideAngle = glideAngle(i); lat = lat(i); lon = lon(i); -ptime = ptime(i); temp = temp(i); cond = cond(i); pres = pres(i); chlor = chlor(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', ptime, '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... -densCorrected = sw_pden(salinCorrected, temp, pres, 0); - - -% convert ptime into datenum style... -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); - -% 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 mat file name... -strMatFileName = strcat(strGliderName, '_CTD_Deployment', strDeploymentNumber, '.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'); - - - Index: gliderproc/trunk/gliderCTD_GeneratePlots.m =================================================================== --- gliderproc/trunk/gliderCTD_GeneratePlots.m (revision 495) +++ gliderproc/trunk/gliderCTD_GeneratePlots.m (revision 496) @@ -47,5 +47,5 @@ % mat file path string... - strMatFilePath = strcat('GLIDER_CTD_DATA_LEVEL1/', strGliderName, '_CTD_Deployment', strDeploymentNumber, '.mat'); + strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, '_Deployment', strDeploymentNumber, '_DataL1.mat'); @@ -142,5 +142,4 @@ nLegEndPoints(i,2) = I(end); end - end Index: gliderproc/trunk/gliderCTD_Generate_L1_Data_v1.m =================================================================== --- (revision ) +++ gliderproc/trunk/gliderCTD_Generate_L1_Data_v1.m (revision 496) @@ -1,0 +1,502 @@ +% +% gliderCTD_Generate_L1_Data_v1.m +% +% Purpose: Apply thermal lag correction to glider CTD salinity and density +% and generate Level 1 data mat files. +% +% NOTE: At top of file, set glider index, deployment number and desired correction parameters +% +% Requires: +% correctThermalLag.m - Code from authors that implements Garau et al. 2011, Thermal +% lag correction on Slocum CTD glider data, JAOT, 28, +% 1065-1071,doi:10.1175/JTECH-D-10-05030.1 +% +% MATLAB folder - contains util, seawater, +% +% +% Author: William Stark +% Marine Sciences Department +% UNC-Chapel Hill +% +% Created: 21 May 2012 +% reviewed, modified by HES December 2012, March 2013, called v1 +% +%////////////////////////////////////////////////////////////////////////// +% TO DO: (done, 3/12/2013) +% compare sci_m_present_time and sci_ctdxx_timestamp - looks like ptime-ctd_time is +% 0.6 +-0.6 second, never goes negative but can be very big positive +% number. Appears to happen at surfacing, and that the CTD throws several +% measurements with the same time stamp, maybe a flushed buffer? Have +% badflagged the values, identified as times when the two timestamps are very +% different. +% +% check final density calculation - done correctly +% +% resolve best determination of lat and lon - HES: looked at this some; +% lat/lon are most common measure but appear to reflect dead reckoning. +% waypoints are just that, intended targets, so chose gps in the end for +% determining position of each CTD point. Commented out code to process +% other position measures for now. +% +% removed chlorophyll from processing, do separately + +clear all; + + +%!!!!!! SET PRIOR TO RUNNING CODE !!!!!!!!!!!!!!!!!!!!!!!! +%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +% +% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ... +gliderIndex = 1; + +% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ... +deploymentNumber = 1; + +% SET CORRECTION PARAMETERS STRUCTURE... +correctionParams = [0.1587 0.0214 6.5316 1.5969]; + +% +%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + +% add paths for required files... +addpath('MATLAB/util/'); +%addpath('MATLAB/matutil/'); +addpath('MATLAB/seawater/'); +%addpath('MATLAB/plots/'); +%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]; + 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 + + +%########################################################################################## + + + + + +%*** 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 +%************************************************************************** + +%*** 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 +%************************************************************************** + + +% 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, '_DataL1.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'); + + + Index: gliderproc/trunk/gliderOptode_Generate_L1_Data.m =================================================================== --- (revision ) +++ gliderproc/trunk/gliderOptode_Generate_L1_Data.m (revision 496) @@ -1,0 +1,396 @@ +% +% gliderOptode_Generate_L1_Data.m +% +% Purpose: Generate Level 1 data mat files that include optode DO data. +% +% NOTE: At top of file, set glider index and deployment number +% +% Requires: +% MATLAB folder - contains util, matutil, seawater, plots, strfun, opnml +% +% +% Author: William Stark +% Marine Sciences Department +% UNC-Chapel Hill +% +% Created: 21 May 2012 +% +%////////////////////////////////////////////////////////////////////////// + +clear all; + + +%!!!!!! SET PRIOR TO RUNNING CODE !!!!!!!!!!!!!!!!!!!!!!!! +%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +% +% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ... +gliderIndex = 1; + +% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ... +deploymentNumber = 1; + +% +%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + +% add paths for required files... +addpath('MATLAB/util/'); +addpath('MATLAB/matutil/'); +addpath('MATLAB/seawater/'); +addpath('MATLAB/plots/'); +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]; +% 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 + + +%########################################################################################## + + + + + +%*** READ IN EBD DATA ***************************************************** +% declare variables for storing data... +oxy_dphase=[]; +ptime_oxy=[]; + +% 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)) + oxy_dphase = [oxy_dphase;data.data(:,strmatch('sci_oxy3835_wphase_dphase', data.vars, 'exact'))]; % oxygen (dphase) + ptime_oxy = [ptime_oxy;data.data(:,strmatch('sci_m_present_time', data.vars, 'exact'))]; % present time + %should read in other oxy values for comparison with + %re-calculated values + end + + data = []; + end +end +%************************************************************************** + + +% load glider data from existing mat file... +strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, '_Deployment', strDeploymentNumber, '_DataL1.mat'); +load(strMatFilePath); + + +% first, apply the sort() function to ptime_oxy to make sure it increases monotonically... +[Y,I] = sort(ptime_oxy); +ptime_oxy = Y; +oxy_dphase = oxy_dphase(I); + + +% remove nans from EBD data... +i = find(~isnan(oxy_dphase)); +oxy_dphase = oxy_dphase(i); ptime_oxy = ptime_oxy(i); + +% interpolate CTD to times of oxy_dphase +tempi = interp1(ptime,temp,ptime_oxy); +salini = interp1(ptime,salinCorrected,ptime_oxy); +depthi = interp1(ptime,depth,ptime_oxy); + +% first, implement temperature-dependent correction to DO concentration +% that utilizes dphase as the input (manual page 30), coefficients from cal +% sheets (a 5x4 matrix of values) and will be glider (and foil) dependent. + +% for ramses, 2 June 2010 calibration sheet + +C0coef=[4.53793e3 -1.62595e2 3.29574 -2.79285e-2]; +C1coef=[-2.50953e2 8.02322 -1.58398e-1 1.31141e-3]; +C2coef=[5.66417 -1.59647e-1 3.07910e-3 -2.46265e-5]; +C3coef=[-5.99449e-2 1.48326e-3 -2.82110e-5 2.15156e-7]; +C4coef=[2.43614e-4 -5.26759e-6 1.00064e-7 -7.14320e-10]; + +C0 = C0coef(1)+C0coef(2)*tempi+C0coef(3)*tempi.^2+C0coef(4)*tempi.^3; +C1 = C1coef(1)+C1coef(2)*tempi+C1coef(3)*tempi.^2+C1coef(4)*tempi.^3; +C2 = C2coef(1)+C2coef(2)*tempi+C2coef(3)*tempi.^2+C2coef(4)*tempi.^3; +C3 = C3coef(1)+C3coef(2)*tempi+C3coef(3)*tempi.^2+C3coef(4)*tempi.^3; +C4 = C4coef(1)+C4coef(2)*tempi+C4coef(3)*tempi.^2+C4coef(4)*tempi.^3; + +o2 = C0+C1*oxy_dphase+C2*oxy_dphase.^2+C3*oxy_dphase.^3+C4*oxy_dphase.^4; + +% second, implement the salinity correction to DO concentration (page 31) +% ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER +% SCALING IS REQUIRED + +B0=-6.24097e-3; +B1=-6.93498e-3; +B2=-6.90358e-3; +B3=-4.29155e-3; +C0=-3.11680e-7; + +Ts = log((298.15 - tempi) ./ (273.15 + tempi)); + +o2_salcorr = o2*exp(salini*(B0+B1*Ts+B2*Ts.*Ts+B3*Ts.*Ts.*Ts)+C0*salini.^2); + +% third, implement the pressure correction to DO concentration (page 32) + +o2_salprescorr = os_salcorr*(1 + 0.04*depthi/1000); + +% use polynomial to calculate DO satuations using the measured temp and +% sal (manual page 30) + +A0 = 2.00856; +A1 = 3.22400; +A2 = 3.99063; +A3 = 4.80299; +A4 = 0.0978188; +A5 = 1.71069; +B0 = -00624097; +B1 = -00693498; +B2 = -00690358; +B3 = -00429155; +C0 = -000000311680; + +% need interpolated temp, salinity, pressure at times of DO obs + +rslt = A0 + A1*Ts + A2*Ts.^2 + A3*Ts.^3 + A4*Ts.^4 + A5*Ts.^5 +... + salini.*(B0 + B1*Ts + B2*Ts.^2 + B3*Ts.^3) + C0*salini.^2; + +O2sol = exp(rslt); + +o2_sat = (o2_salprescorr * 2.2414) ./ O2sol; + + +% convert ptime into datenum style... +ptime_oxy_datenum = ptime_oxy/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); + + + + +lb = datenum('10-Feb-2012 12:22:00'); +ub = datenum('10-Feb-2012 12:23:00'); +%lb = datenum('22-Mar-2012 12:01:00'); +%ub = datenum('23-Mar-2012 12:02:00'); + +% figure('Position', [500,500,1000,1000]); +% subplot(311); +% ind = find(ptime_oxy_datenum >= lb & ptime_oxy_datenum <= ub); +% plot(ptime_oxy_datenum(ind), oxy_dphase(ind), 'b.'); +% datetick; +% xlabel('time', 'fontsize', 12); +% ylabel('DO (dphase)', 'fontsize', 12); +% subplot(312); +% ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub); +% plot(ptime_datenum(ind2), temp(ind2), 'r.'); +% datetick; +% xlabel('time', 'fontsize', 12); +% ylabel('temp', 'fontsize', 12); +% subplot(313); +% ind = find(ptime_oxy_datenum >= lb & ptime_oxy_datenum <= ub); +% ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub); +% plot(ptime_oxy_datenum(ind), ones(length(ind), 1), 'b.', ptime_datenum(ind2), ones(length(ind2), 1), 'r.'); +% datetick; + + +figure('Position', [500,500,1000,1000]); +subplot(211); +ind = find(ptime_datenum >= lb & ptime_datenum <= ub); +plot(ptime_datenum(ind), oxy_dphase(ind), 'b.', ptime_datenum(ind), oxy_dphase_adjusted(ind), 'r.'); +datetick; +xlabel('time', 'fontsize', 12); +ylabel('DO (dphase)', 'fontsize', 12); +subplot(212); +ind2 = find(ptime_datenum >= lb & ptime_datenum <= ub); +plot(ptime_datenum(ind2), temp(ind2), 'r.'); +datetick; +xlabel('time', 'fontsize', 12); +ylabel('temp', 'fontsize', 12); + + + + + + + + +% interpolate DBD data to align with EBD data... +% horizontalVelocity = interp1(ptime1_dbd, horizontalVelocity, ptime); +% depth = interp1(ptime1_dbd, depth, ptime); +% pitch = interp1(ptime1_dbd, pitch, ptime); +% avgDepthRate = interp1(ptime1_dbd, avgDepthRate, ptime); +% glideAngle = interp1(ptime1_dbd, glideAngle, ptime); +% lat = interp1(ptime1_dbd, lat, ptime); +% lon = interp1(ptime1_dbd, lon, ptime); + + + + + + + + + + + + +% 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, '_WithDO_LEVEL1.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'); + + +