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
%
%//////////////////////////////////////////////////////////////////////////
% 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 = 2;

% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ...
deploymentNumber = 3;

% SET CORRECTION PARAMETERS STRUCTURE...
correctionParams = [0.1587 0.0214 6.5316 1.5969];

%
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



% add paths for required files...
addpath('MATLAB/util/');
%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, '_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');