root/gliderproc/trunk/gliderCTD_Generate_L1_Data.m

%
%  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
%   June 2013, implemented some QC procedures for ramses, HES
%   20130703 - Iterate throuth gliders and deployments. - CBC
%   20130725 - Implemented some QC procedures for pelagia. - HES
%
%//////////////////////////////////////////////////////////////////////////
% 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;

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

% 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
        
        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


        %##########################################################################################





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

        [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);
        elseif(gliderIndex == 1)
            ib=find(abs(diff(temp))>1.5);
            ib2=find(abs(diff(cond))>0.1);
            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

        % 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); hv = hv(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 (now set based on deployment 1)

        if(gliderIndex == 2)
            iv = find(hv<0.1 | hv > 0.7);
            ip = find (abs(pitch) < 10.);
            ib = union(iv,ip);
            salinCorrected(ib) = NaN;
        elseif(gliderIndex == 1)
            iv = find(hv<0.1);
            ip = find(pitch>5 & pitch < 15);
            ib = union(iv,ip);
            salinCorrected(ib) = NaN;
        end
            
            % 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);
        

        % 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