Compare View

switch
from
...
to
 
Commits (3)
math/boxav.m 0 → 100644
... ... @@ -0,0 +1,38 @@
  1 +function qA = boxav(x,y,qx,qy,q)
  2 +% Synthax : qA = boxav(x,y,qx,qy,q)
  3 +%
  4 +% Box averaging on a rectangular but not necessarily regular grid. Grid
  5 +% coordinates are given by matrices 'x' and 'y'. Averaged variable 'q' is
  6 +% located at 'qx' and 'qy' coordinates and can be vectors or matrices.
  7 +% Averaged variable 'qA', the size of 'x' and 'y' is the output.
  8 +
  9 +% Initialize output
  10 +qA = nan(size(x)) ;
  11 +
  12 +% Conduct averaging
  13 +for ii = 1:numel(x(:,1)) % longitude loop
  14 + for jj = 1:numel(y(1,:)) % latitude loop
  15 +
  16 + % calculate distance to neighbor points
  17 + try dxl = abs( x(ii,jj)-x(ii,jj-1) ); catch; dxl = abs( x(ii,jj)-x(ii,jj+1) ); end
  18 + try dxr = abs( x(ii,jj)-x(ii,jj+1) ); catch; dxr = abs( x(ii,jj)-x(ii,jj-1) ); end
  19 + try dyu = abs( y(ii,jj)-y(ii+1,jj) ); catch; dyu = abs( y(ii,jj)-y(ii-1,jj) ); end
  20 + try dyd = abs( y(ii,jj)-y(ii-1,jj) ); catch; dyd = abs( y(ii,jj)-y(ii+1,jj) ); end
  21 +
  22 + % set logical condition to find points of interest
  23 + I = find(qx <= x(ii,jj) + dxr./2 & ...
  24 + qx >= x(ii,jj) - dxl./2 & ...
  25 + qy <= y(ii,jj) + dyu./2 & ...
  26 + qy >= y(ii,jj) - dyd./2) ;
  27 +
  28 + % average
  29 + if ~isempty(I)
  30 + qA(ii,jj) = mean(q(I),'omitnan') ;
  31 + else
  32 + qA(ii,jj) = nan ;
  33 + end
  34 +
  35 + end
  36 +end
  37 +
  38 +end
0 39 \ No newline at end of file
... ...
math/genopt.m 0 → 100644
... ... @@ -0,0 +1,66 @@
  1 +function MOM = genopt(Q,X,dom,NP,NI,pc,pm)
  2 +% Synthax: MOM = genopt(@Q,@X,dom,NP,NI,pc,pm)
  3 +%
  4 +% Optimisation routine following a genetic algorithm. Returns the best parameter
  5 +% set which solves the problem expressed by the anonymous function 'Q', when
  6 +% either the criteria expressed by the anonymous function 'X' has been met or
  7 +% NI iterations have been completed. 'X' should be formulated in terms of the
  8 +% values returned by 'Q'. The 'dom' variable is a matrix of size [p 2], where
  9 +% p is the number of parameters passed to the function 'Q'. It's rows are the
  10 +% domains from which the initial population will be generated for each parameter.
  11 +% Variable 'NP' is the size of the population and must be an even number. Variables
  12 +% 'pc' and 'pm' are respectively the probability of mixing genetic material and
  13 +% probability of mutation. Recommended values are 0.8 and 0.2 .
  14 +
  15 +%INIT POPULATION
  16 +NA = numel(dom(:,1)) ;
  17 +PP = nan(NP,NA) ;
  18 +SC = nan(NP,1) ;
  19 +RK = nan(NP,1) ;
  20 +for ii = 1:NP
  21 + for jj = 1:NA
  22 + PP(ii,jj) = dom(jj,1) + (dom(jj,2)-dom(jj,1))*rand(1) ;
  23 + end
  24 +end
  25 +
  26 +for kk = 1:NI
  27 +
  28 + %EVALUATION
  29 + for ii = 1:NP ; SC(ii) = Q(PP(ii,:)) ; end
  30 +
  31 + %RANKING
  32 + [HS,RK] = sort(SC) ;
  33 + PP = PP(RK,:) ;
  34 +
  35 + %TEST
  36 + if X(HS(1)) ; break; end
  37 +
  38 + %SELECTION
  39 + DICE = randi([2 4]) ;
  40 + MOM = PP(1, :) ;
  41 + DAD = PP(DICE,:) ;
  42 +
  43 + %REPRODUCTION
  44 + PP(1,:) = MOM ;
  45 + PP(2,:) = DAD ;
  46 + for ii = 3:2:NP
  47 + for jj = 1:NA
  48 + % CROISEMENT
  49 + if rand(1) < pc
  50 + DICE = rand(1) ;
  51 + PP(ii,jj) = MOM(jj)*DICE + DAD(jj)*(1-DICE) ;
  52 + PP(ii+1,jj) = DAD(jj)*DICE + MOM(jj)*(1-DICE) ;
  53 + else
  54 + PP(ii,jj) = MOM(jj) ;
  55 + PP(ii+1,jj) = DAD(jj) ;
  56 + end
  57 +
  58 + % MUTATION
  59 + if rand(1) < pm
  60 + PP(ii,jj) = dom(jj,1) + (dom(jj,2)-dom(jj,1))*rand(1) ;
  61 + end
  62 + end
  63 + end
  64 +end
  65 +
  66 +end
... ...
math/nlinfit.m
... ... @@ -27,6 +27,7 @@ elseif sz(1) == 1 &amp; sz(2) &gt; 1 | sz(2) == 1 &amp; sz(1) &gt; 1
27 27 start = pars ;
28 28 end
29 29  
  30 +
30 31 % Minimize parameters
31 32 A = fminsearch(@(a)sum( (y - rhs(a,x,varargin{:})).^2 ),start) ;
32 33 yfit = rhs(A,x,varargin{:}) ;
... ...
math/rho_t.m 0 → 100644
... ... @@ -0,0 +1,33 @@
  1 +function [rho] = rho_t(x,y,degrad)
  2 +% [rho] = rho_t(x,y,degrad)
  3 +%
  4 +% Explicitly calculates the correlation coefficient rho_t detailed in
  5 +% Fisher and Lee (1983), between the two circular variables 'x' and 'y'.
  6 +% Variable 'degrad' is a string and should be "deg" if 'x' and 'y' are
  7 +% in degrees or 'rad' if they are in radians.
  8 +%
  9 +% Reference: Fisher, N.I., LEE, A.J.,(1983), A correlation coefficient
  10 +% for circular data, Biometrika, 70, 2, pp. 327-32
  11 +
  12 +
  13 +numerator = 0 ;
  14 +denA = 0 ;
  15 +denB = 0 ;
  16 +
  17 +switch degrad
  18 + case 'deg'
  19 + x = x*pi/180 ;
  20 + y = -y*pi/180 ;
  21 + case 'rad'
  22 + % do nothing
  23 +end
  24 +
  25 +for ii = 1:numel(x)-1
  26 + numerator = numerator + sum( sin(x(ii)-x(ii+1:end)).*sin(y(ii)-y(ii+1:end)) ) ;
  27 + denA = denA + sum( sin(x(ii)-x(ii+1:end)).^2 );
  28 + denB = denB + sum( sin(y(ii)-y(ii+1:end)).^2 ) ;
  29 +end
  30 +
  31 +rho = numerator./( sqrt(denA).*sqrt(denB) ) ;
  32 +
  33 +end
... ...
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.005 ) = 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
... ...