% % 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');