Commit cc34cf0a9e9051342e145bab8786f5708dc21189

Authored by Jean-Luc Shaw
1 parent 8f5d4751
Exists in master

Added automated quality control for SeaBird cnv files function seabirdCTD.m

Showing 2 changed files with 254 additions and 0 deletions   Show diff stats
util/bin_data.m 0 → 100644
... ... @@ -0,0 +1,15 @@
  1 +function [yb] = bin_data(x,y,xb)
  2 +
  3 +% Synthaxe : [yb] = bin_data(x,y,xb)
  4 +
  5 +% Cette fonction produit un vecteur 'yb' bined aux bins 'xb'
  6 +
  7 +dx = xb(2) - xb(1) ;
  8 +
  9 +yb = nan(size(xb));
  10 +for i = 1:length(xb)
  11 + j = find(x <= xb(i) + 0.5*dx & x >= xb(i)-0.5*dx) ;
  12 + yb(i) = mean(y(j),'omitnan') ;
  13 +end
  14 +
  15 +end
... ...
util/seabirdCTD.m 0 → 100644
... ... @@ -0,0 +1,239 @@
  1 +function varargout = seabirdCTD(path,varargin)
  2 +% Synthax : varargout = seabirdCTD(path,varargin)
  3 +%
  4 +% Read Seabird electronics CTD .cnv file(s) automatically, select downcast,
  5 +% and bin average the cast profiles. The data is organised into output
  6 +% structure 'S' with time as field 'dnum' and depth as field 'z'. Other
  7 +% fields depend on the data contained in the raw .cnv file. Detected data
  8 +% types are:
  9 +%
  10 +% DATA field name
  11 +%
  12 +% temperature -> S.temp
  13 +% salinity -> S.sal
  14 +% fluorescence -> S.fluo
  15 +% turbidity -> S.turb
  16 +% oxygen -> S.oxy
  17 +% longitude -> S.lon
  18 +% latitude -> S.lat
  19 +% flag -> S.flag
  20 +%
  21 +% Optional parameters should be supplied in parameter/value pairs. Available
  22 +% options are:
  23 +%
  24 +% warmup, [integer] number of minutes allowed for probe warmup
  25 +% when sampling. Defaults to 0.
  26 +% verticalVelocity, [float] minimum vertical velocity kept.
  27 +% defaults to 0.
  28 +%
  29 +% binSize, [float] averaging bin size in (m).
  30 +% defaults to 1 m.
  31 +%
  32 +% showz, [string] plot raw depth against time and cleaned
  33 +% binned depth agains time.
  34 +% Options are: 'y' or 'n'
  35 +% Defaults to 'n'.
  36 +%
  37 +% namest, [cell array] extract information from the file name.
  38 +% The first element of the string array
  39 +% is the delimiter, the remaining are pairs
  40 +% of string and integer
  41 +%
  42 +% gpsgar, [string] path to a Garmin .csv gps track log. Track
  43 +% log is interpolated to the sampling time
  44 +% vector. Probably not general to all Garmin
  45 +% devices. If problematic edit lines 83-92.
  46 +%
  47 +% gsw, [string] calculate absolute salinity, conservative
  48 +% temperature and density using the Gibbs
  49 +% Sea Water matlab package. The package must
  50 +% be on the matlab path for this to work.
  51 +% Added fields to the structure are
  52 +%
  53 +% cons. temperature -> S.ctemp
  54 +% cons. density -> S.pden
  55 +% absolute salinity -> S.asal
  56 +%
  57 +% If the 'path' variable points to a file, only this file is processed and the
  58 +% structure 'S' is returned as a variable. If it points to a folder, all the
  59 +% .cnv files present in this folder are processed and the output suppressed.
  60 +%
  61 +
  62 +
  63 +% MANAGE OPTIONS
  64 +% defaults
  65 +warmup = 0;
  66 +dzthres = 0;
  67 +bsz = 1;
  68 +showz = 'n' ;
  69 +gsw = false;
  70 +
  71 +% user
  72 +if numel(varargin) > 1
  73 + for ii = 1:2:numel(varargin)
  74 + switch varargin{ii}
  75 + case 'warmup'
  76 + warmup = varargin{ii+1}*60 ;
  77 + case 'verticalVelocity'
  78 + dzthres = varargin{ii+1} ;
  79 + case 'binSize'
  80 + bsz = varargin{ii+1} ;
  81 + case 'showz'
  82 + showz = varargin{ii+1} ;
  83 + case 'namestr'
  84 + Nin = numel(varargin{ii+1}) ;
  85 + delim = varargin{ii+1}{1};
  86 + field = cell([(Nin-1)/2 1]) ;
  87 + namef = cell([(Nin-1)/2 1]) ;
  88 + for jj = 1:(Nin-1)/2
  89 + field{jj} = varargin{ii+1}{2*jj} ;
  90 + namef{jj} = varargin{ii+1}{2*jj+1} ;
  91 + end
  92 + case 'gpsgar' % read navigation from a garmin gps .csv
  93 + fid = fopen(varargin{ii+1}) ;
  94 + ln = fgetl(fid);
  95 + while ~startsWith(ln,'ID,trkseg'); ln = fgetl(fid) ; end
  96 + frmt = ['%f%f%f%f%f' repmat('%s',[1 numel(split(ln,','))-5])] ;
  97 + A = textscan(fid,frmt,'delimiter',',');
  98 + lon = A{4};
  99 + lat = A{3};
  100 + tim = datenum(A{6},'yyyy-mm-ddTHH:MM:SS');
  101 + fclose(fid);
  102 + case 'gsw'
  103 + if startsWith(varargin{ii+1},'y'); gsw = true; end
  104 + end
  105 + end
  106 +end
  107 +
  108 +% FIND + ORGANISE FILE CANDIDATES
  109 +if isdir(path)
  110 + dirst = dir([path '/*.cnv']) ;
  111 + fname = {dirst.name}';
  112 + fpath = {dirst.folder}';
  113 +else
  114 + dirst = dir([path]) ;
  115 + fname = {dirst.name}';
  116 + fpath = {dirst.folder}';
  117 +end
  118 +
  119 +for kk = 1:numel(fname)
  120 +
  121 +% DISPLAY FILENAME
  122 +disp(['Processing file : ' fname{kk}])
  123 +
  124 +fid = fopen([fpath{kk} '/' fname{kk}]) ;
  125 +ln = fgetl(fid) ;
  126 +
  127 +sens = {'temperature','salinity','Depth','oxygen','fluorescence','turbidity','flag','latitude','longitude'} ;
  128 +Nsens = 0 ;
  129 +colname = [] ;
  130 +while ~startsWith(ln,'*END*')
  131 + ln = fgetl(fid) ;
  132 + % get time related values
  133 + if startsWith(ln,'# nvalues') ;
  134 + spln = split(ln);
  135 + nvalues = str2num(spln{4}); end
  136 + if startsWith(ln,'# interval') ;
  137 + spln = split(ln) ;
  138 + sr = str2num(spln{end}); end
  139 + if startsWith(ln,'# start_time') ;
  140 + spln = split(ln) ;
  141 + st = datenum([spln{4:7}],'mmmddyyyyHH:MM:SS'); end
  142 +
  143 + % figure out probe channels
  144 + if startsWith(ln,'# name') ;
  145 + Nsens = Nsens + 1 ;
  146 + tmp = regexpi(ln,sens,'match') ;
  147 + I = find(~cellfun('isempty',tmp)) ;
  148 + switch I
  149 + case 1
  150 + colname = [colname; {'temp'}] ;
  151 + case 2
  152 + colname = [colname; {'sal'}] ;
  153 + case 3
  154 + colname = [colname; {'z'}] ; Iz = Nsens ;
  155 + case 4
  156 + colname = [colname; {'oxy'}] ;
  157 + case 5
  158 + colname = [colname; {'fluo'}] ;
  159 + case 6
  160 + colname = [colname; {'turb'}] ;
  161 + case 7
  162 + colname = [colname; {'flag'}] ;
  163 + case 8
  164 + colname = [colname; {'lat'}] ;
  165 + case 9
  166 + colname = [colname; {'lon'}] ;
  167 + end
  168 + end
  169 +end
  170 +
  171 +%% READ DATA
  172 +A = textscan(fid,repmat('%f',[1 Nsens])) ;
  173 +
  174 +%% MAKE TIME VECTOR
  175 +dnum = st+sr*[0:nvalues-1]'/(24*3600) ;
  176 +
  177 +%% GET DOWNCAST
  178 +z = movmean(A{:,Iz},10/sr) ;
  179 +if startsWith(showz,'y'); figure; plot(dnum,z,'b'); hold on; end
  180 +
  181 +dz = gradient(z,dnum*24*3600) ; % compute downward velocity
  182 +z( [1:nvalues]' < warmup/sr | dz<=0.001 ) = nan ; % get rid of start and end values
  183 +[~,Im] = min(z) ; % find min depth
  184 +[~,IM] = max(z) ; % find max depth
  185 +I = find(dz > dzthres & [1:nvalues]' > Im & [1:nvalues]' < IM) ; % select down-going between min z and max z
  186 +if startsWith(showz,'y'); plot(dnum(Im),z(Im),'*'); plot(dnum(IM),z(IM),'o'); end
  187 +
  188 +%% MAKE OUTPUT STRUCTURE
  189 +S = [] ;
  190 +for ii = 1:Nsens
  191 + if ii ~= Iz; S = setfield(S,colname{ii},A{ii}(I)) ; end
  192 +end
  193 +S.dnum = dnum(I) ;
  194 +
  195 +%% BIN AVERAGE
  196 +xb = [0:bsz:max(z(I))+0.5*bsz] ;
  197 +flds = fieldnames(S) ;
  198 +for ii = 1:Nsens
  199 + yb = bin_data(z(I),getfield(S,flds{ii}),xb) ;
  200 + S = setfield(S,flds{ii},yb) ;
  201 +end
  202 +S.z = xb ;
  203 +
  204 +%% EXTRACT INFO FROM FILE NAME
  205 +if exist('delim')
  206 + str = split(fname{kk}(1:end-4),delim);
  207 + for ii = 1:numel(namef)
  208 + S = setfield(S,namef{ii},str{field{ii}}) ;
  209 + end
  210 +end
  211 +
  212 +%% ATTACH GPS COORDINATES
  213 +if exist('lon')
  214 + S.lon = interp1(tim,lon,S.dnum) ;
  215 + S.lat = interp1(tim,lat,S.dnum) ; end
  216 +
  217 +%% USE GSW IF POSSIBLE
  218 +if gsw
  219 + S.asal = gsw_SA_from_SP(S.sal,0,-66.4,49);
  220 + S.ctemp = gsw_CT_from_t(S.asal,S.temp,0);
  221 + S.pden = gsw_sigma0(S.asal,S.ctemp); end
  222 +
  223 +%% SHOW KEPT/DISCARDED DATA
  224 +if startsWith(showz,'y')
  225 + plot(S.dnum,S.z,'.');
  226 + datetick('x','HHMM');
  227 + ylabel('z (m)');
  228 + title(fname{kk},'interpreter','none') ; end
  229 +
  230 +%% OUTPUT
  231 +if nargout<1
  232 + save([fpath{kk} '/' fname{kk}(1:end-4) '.mat'],'S') ;
  233 +elseif nargout==1 & numel(fname)==1
  234 + varargout{1} = S;
  235 +end
  236 +
  237 +end
  238 +
  239 +end
... ...