...

Commits (3)
math/boxav.m 0 → 100644
 function qA = boxav(x,y,qx,qy,q) % Synthax : qA = boxav(x,y,qx,qy,q) % % Box averaging on a rectangular but not necessarily regular grid. Grid % coordinates are given by matrices 'x' and 'y'. Averaged variable 'q' is % located at 'qx' and 'qy' coordinates and can be vectors or matrices. % Averaged variable 'qA', the size of 'x' and 'y' is the output. % Initialize output qA = nan(size(x)) ; % Conduct averaging for ii = 1:numel(x(:,1)) % longitude loop for jj = 1:numel(y(1,:)) % latitude loop % calculate distance to neighbor points try dxl = abs( x(ii,jj)-x(ii,jj-1) ); catch; dxl = abs( x(ii,jj)-x(ii,jj+1) ); end try dxr = abs( x(ii,jj)-x(ii,jj+1) ); catch; dxr = abs( x(ii,jj)-x(ii,jj-1) ); end try dyu = abs( y(ii,jj)-y(ii+1,jj) ); catch; dyu = abs( y(ii,jj)-y(ii-1,jj) ); end try dyd = abs( y(ii,jj)-y(ii-1,jj) ); catch; dyd = abs( y(ii,jj)-y(ii+1,jj) ); end % set logical condition to find points of interest I = find(qx <= x(ii,jj) + dxr./2 & ... qx >= x(ii,jj) - dxl./2 & ... qy <= y(ii,jj) + dyu./2 & ... qy >= y(ii,jj) - dyd./2) ; % average if ~isempty(I) qA(ii,jj) = mean(q(I),'omitnan') ; else qA(ii,jj) = nan ; end end end end \ No newline at end of file
math/genopt.m 0 → 100644
 function MOM = genopt(Q,X,dom,NP,NI,pc,pm) % Synthax: MOM = genopt(@Q,@X,dom,NP,NI,pc,pm) % % Optimisation routine following a genetic algorithm. Returns the best parameter % set which solves the problem expressed by the anonymous function 'Q', when % either the criteria expressed by the anonymous function 'X' has been met or % NI iterations have been completed. 'X' should be formulated in terms of the % values returned by 'Q'. The 'dom' variable is a matrix of size [p 2], where % p is the number of parameters passed to the function 'Q'. It's rows are the % domains from which the initial population will be generated for each parameter. % Variable 'NP' is the size of the population and must be an even number. Variables % 'pc' and 'pm' are respectively the probability of mixing genetic material and % probability of mutation. Recommended values are 0.8 and 0.2 . %INIT POPULATION NA = numel(dom(:,1)) ; PP = nan(NP,NA) ; SC = nan(NP,1) ; RK = nan(NP,1) ; for ii = 1:NP for jj = 1:NA PP(ii,jj) = dom(jj,1) + (dom(jj,2)-dom(jj,1))*rand(1) ; end end for kk = 1:NI %EVALUATION for ii = 1:NP ; SC(ii) = Q(PP(ii,:)) ; end %RANKING [HS,RK] = sort(SC) ; PP = PP(RK,:) ; %TEST if X(HS(1)) ; break; end %SELECTION DICE = randi([2 4]) ; MOM = PP(1, :) ; DAD = PP(DICE,:) ; %REPRODUCTION PP(1,:) = MOM ; PP(2,:) = DAD ; for ii = 3:2:NP for jj = 1:NA % CROISEMENT if rand(1) < pc DICE = rand(1) ; PP(ii,jj) = MOM(jj)*DICE + DAD(jj)*(1-DICE) ; PP(ii+1,jj) = DAD(jj)*DICE + MOM(jj)*(1-DICE) ; else PP(ii,jj) = MOM(jj) ; PP(ii+1,jj) = DAD(jj) ; end % MUTATION if rand(1) < pm PP(ii,jj) = dom(jj,1) + (dom(jj,2)-dom(jj,1))*rand(1) ; end end end end end
 ... ... @@ -27,6 +27,7 @@ elseif sz(1) == 1 & sz(2) > 1 | sz(2) == 1 & sz(1) > 1 start = pars ; end % Minimize parameters A = fminsearch(@(a)sum( (y - rhs(a,x,varargin{:})).^2 ),start) ; yfit = rhs(A,x,varargin{:}) ; ... ...
math/rho_t.m 0 → 100644
 function [rho] = rho_t(x,y,degrad) % [rho] = rho_t(x,y,degrad) % % Explicitly calculates the correlation coefficient rho_t detailed in % Fisher and Lee (1983), between the two circular variables 'x' and 'y'. % Variable 'degrad' is a string and should be "deg" if 'x' and 'y' are % in degrees or 'rad' if they are in radians. % % Reference: Fisher, N.I., LEE, A.J.,(1983), A correlation coefficient % for circular data, Biometrika, 70, 2, pp. 327-32 numerator = 0 ; denA = 0 ; denB = 0 ; switch degrad case 'deg' x = x*pi/180 ; y = -y*pi/180 ; case 'rad' % do nothing end for ii = 1:numel(x)-1 numerator = numerator + sum( sin(x(ii)-x(ii+1:end)).*sin(y(ii)-y(ii+1:end)) ) ; denA = denA + sum( sin(x(ii)-x(ii+1:end)).^2 ); denB = denB + sum( sin(y(ii)-y(ii+1:end)).^2 ) ; end rho = numerator./( sqrt(denA).*sqrt(denB) ) ; end
util/bin_data.m 0 → 100644
 function [yb] = bin_data(x,y,xb) % Synthaxe : [yb] = bin_data(x,y,xb) % Cette fonction produit un vecteur 'yb' bined aux bins 'xb' dx = xb(2) - xb(1) ; yb = nan(size(xb)); for i = 1:length(xb) j = find(x <= xb(i) + 0.5*dx & x >= xb(i)-0.5*dx) ; yb(i) = mean(y(j),'omitnan') ; end end
util/seabirdCTD.m 0 → 100644
 function varargout = seabirdCTD(path,varargin) % Synthax : varargout = seabirdCTD(path,varargin) % % Read Seabird electronics CTD .cnv file(s) automatically, select downcast, % and bin average the cast profiles. The data is organised into output % structure 'S' with time as field 'dnum' and depth as field 'z'. Other % fields depend on the data contained in the raw .cnv file. Detected data % types are: % % DATA field name % % temperature -> S.temp % salinity -> S.sal % fluorescence -> S.fluo % turbidity -> S.turb % oxygen -> S.oxy % longitude -> S.lon % latitude -> S.lat % flag -> S.flag % % Optional parameters should be supplied in parameter/value pairs. Available % options are: % % warmup, [integer] number of minutes allowed for probe warmup % when sampling. Defaults to 0. % verticalVelocity, [float] minimum vertical velocity kept. % defaults to 0. % % binSize, [float] averaging bin size in (m). % defaults to 1 m. % % showz, [string] plot raw depth against time and cleaned % binned depth agains time. % Options are: 'y' or 'n' % Defaults to 'n'. % % namest, [cell array] extract information from the file name. % The first element of the string array % is the delimiter, the remaining are pairs % of string and integer % % gpsgar, [string] path to a Garmin .csv gps track log. Track % log is interpolated to the sampling time % vector. Probably not general to all Garmin % devices. If problematic edit lines 83-92. % % gsw, [string] calculate absolute salinity, conservative % temperature and density using the Gibbs % Sea Water matlab package. The package must % be on the matlab path for this to work. % Added fields to the structure are % % cons. temperature -> S.ctemp % cons. density -> S.pden % absolute salinity -> S.asal % % If the 'path' variable points to a file, only this file is processed and the % structure 'S' is returned as a variable. If it points to a folder, all the % .cnv files present in this folder are processed and the output suppressed. % % MANAGE OPTIONS % defaults warmup = 0; dzthres = 0; bsz = 1; showz = 'n' ; gsw = false; % user if numel(varargin) > 1 for ii = 1:2:numel(varargin) switch varargin{ii} case 'warmup' warmup = varargin{ii+1}*60 ; case 'verticalVelocity' dzthres = varargin{ii+1} ; case 'binSize' bsz = varargin{ii+1} ; case 'showz' showz = varargin{ii+1} ; case 'namestr' Nin = numel(varargin{ii+1}) ; delim = varargin{ii+1}{1}; field = cell([(Nin-1)/2 1]) ; namef = cell([(Nin-1)/2 1]) ; for jj = 1:(Nin-1)/2 field{jj} = varargin{ii+1}{2*jj} ; namef{jj} = varargin{ii+1}{2*jj+1} ; end case 'gpsgar' % read navigation from a garmin gps .csv fid = fopen(varargin{ii+1}) ; ln = fgetl(fid); while ~startsWith(ln,'ID,trkseg'); ln = fgetl(fid) ; end frmt = ['%f%f%f%f%f' repmat('%s',[1 numel(split(ln,','))-5])] ; A = textscan(fid,frmt,'delimiter',','); lon = A{4}; lat = A{3}; tim = datenum(A{6},'yyyy-mm-ddTHH:MM:SS'); fclose(fid); case 'gsw' if startsWith(varargin{ii+1},'y'); gsw = true; end end end end % FIND + ORGANISE FILE CANDIDATES if isdir(path) dirst = dir([path '/*.cnv']) ; fname = {dirst.name}'; fpath = {dirst.folder}'; else dirst = dir([path]) ; fname = {dirst.name}'; fpath = {dirst.folder}'; end for kk = 1:numel(fname) % DISPLAY FILENAME disp(['Processing file : ' fname{kk}]) fid = fopen([fpath{kk} '/' fname{kk}]) ; ln = fgetl(fid) ; sens = {'temperature','salinity','Depth','oxygen','fluorescence','turbidity','flag','latitude','longitude'} ; Nsens = 0 ; colname = [] ; while ~startsWith(ln,'*END*') ln = fgetl(fid) ; % get time related values if startsWith(ln,'# nvalues') ; spln = split(ln); nvalues = str2num(spln{4}); end if startsWith(ln,'# interval') ; spln = split(ln) ; sr = str2num(spln{end}); end if startsWith(ln,'# start_time') ; spln = split(ln) ; st = datenum([spln{4:7}],'mmmddyyyyHH:MM:SS'); end % figure out probe channels if startsWith(ln,'# name') ; Nsens = Nsens + 1 ; tmp = regexpi(ln,sens,'match') ; I = find(~cellfun('isempty',tmp)) ; switch I case 1 colname = [colname; {'temp'}] ; case 2 colname = [colname; {'sal'}] ; case 3 colname = [colname; {'z'}] ; Iz = Nsens ; case 4 colname = [colname; {'oxy'}] ; case 5 colname = [colname; {'fluo'}] ; case 6 colname = [colname; {'turb'}] ; case 7 colname = [colname; {'flag'}] ; case 8 colname = [colname; {'lat'}] ; case 9 colname = [colname; {'lon'}] ; end end end %% READ DATA A = textscan(fid,repmat('%f',[1 Nsens])) ; %% MAKE TIME VECTOR dnum = st+sr*[0:nvalues-1]'/(24*3600) ; %% GET DOWNCAST z = movmean(A{:,Iz},10/sr) ; if startsWith(showz,'y'); figure; plot(dnum,z,'b'); hold on; end dz = gradient(z,dnum*24*3600) ; % compute downward velocity z( [1:nvalues]' < warmup/sr | dz<=0.005 ) = nan ; % get rid of start and end values [~,Im] = min(z) ; % find min depth [~,IM] = max(z) ; % find max depth I = find(dz > dzthres & [1:nvalues]' > Im & [1:nvalues]' < IM) ; % select down-going between min z and max z if startsWith(showz,'y'); plot(dnum(Im),z(Im),'*'); plot(dnum(IM),z(IM),'o'); end %% MAKE OUTPUT STRUCTURE S = [] ; for ii = 1:Nsens if ii ~= Iz; S = setfield(S,colname{ii},A{ii}(I)) ; end end S.dnum = dnum(I) ; %% BIN AVERAGE xb = [0:bsz:max(z(I))+0.5*bsz] ; flds = fieldnames(S) ; for ii = 1:Nsens yb = bin_data(z(I),getfield(S,flds{ii}),xb) ; S = setfield(S,flds{ii},yb) ; end S.z = xb ; %% EXTRACT INFO FROM FILE NAME if exist('delim') str = split(fname{kk}(1:end-4),delim); for ii = 1:numel(namef) S = setfield(S,namef{ii},str{field{ii}}) ; end end %% ATTACH GPS COORDINATES if exist('lon') S.lon = interp1(tim,lon,S.dnum) ; S.lat = interp1(tim,lat,S.dnum) ; end %% USE GSW IF POSSIBLE if gsw S.asal = gsw_SA_from_SP(S.sal,0,-66.4,49); S.ctemp = gsw_CT_from_t(S.asal,S.temp,0); S.pden = gsw_sigma0(S.asal,S.ctemp); end %% SHOW KEPT/DISCARDED DATA if startsWith(showz,'y') plot(S.dnum,S.z,'.'); datetick('x','HHMM'); ylabel('z (m)'); title(fname{kk},'interpreter','none') ; end %% OUTPUT if nargout<1 save([fpath{kk} '/' fname{kk}(1:end-4) '.mat'],'S') ; elseif nargout==1 & numel(fname)==1 varargout{1} = S; end end end