Commit 31d166dc authored by Jean-Luc Shaw's avatar Jean-Luc Shaw
Browse files

Initial commit.

parents
function [lond,latd] = dm2deci(londeg,lonmin,lonsec,EW,latdeg,latmin,latsec,NS)
%
% This function converts coordinates from deg/min/sec to decimal
%
% [lond,latd] = dm2deci(londeg,lonmin,lonsec,EW,latdeg,latmin,latsec,NS)
%
% Enter numbers for longitude degrees, minutes and seconds, and a string for cardinal direction ('E' or 'W')
% followed by
% numbers for latitude degrees, minutes and seconds, and a string for cardinal direction ('N' or 'S')
%
% Written by Jean-Luc Shaw
%
if(strmatch(NS,'N')) ; lat_sign = 1. ; end
if(strmatch(NS,'S')) ; lat_sign = -1. ; end
if(strmatch(EW,'E')) ; lon_sign = 1. ; end
if(strmatch(EW,'W')) ; lon_sign = -1. ; end
lond = lon_sign*(londeg+lonmin./60+lonsec./3600) ;
latd = lat_sign*(latdeg+latmin./60+latsec./3600) ;
end
function [LNG,LAT] = globeGrid(lims,sz)
% Synthax : [LNG,LAT] = globeGrid(lims,sz)
%
% Takes as input longitude and latitude limits 'lims', in format :
%
% lims = [lon_min lon_max lat_min lat_max]
%
% and a grid size 'sz' in kilometers, and generates a grid uniform in
% kilometer space, with the center of the x dimension aligned with the
% center of the area specified by 'lims'. This is mostly designed for use
% with small grids (~50-100km^2) where this representation looks almost
% regular in coordinate space and kind of "looks nice"... at least to me it
% does.
%
%
% EARTH RADIUS
R = 6371 ;
% convert to km
lng = [lims(1); lims(1); lims(2); lims(2)] ;
lat = [lims(3); lims(4); lims(3); lims(4)] ;
avl = mean(lng) ;
lng = lng - avl ; % move to prime meridian
% SINUSOIDAL PROJECTION
y = ( lat*R*pi / 180 ) ;
x = ( cosd(lat).*lng*R*pi / 180 ) ;
avy = mean(y) ;
avx = mean(x) ;
% MAKE GRID VECTORS AROUND CENTER IN KM SPACE
xg = [flipud( [avx:-sz:min(x)]' ) ; [avx+sz:sz:max(x)]' ] ;
yg = [flipud( [avy:-sz:min(y)]' ) ; [avy+sz:sz:max(y)]' ] ;
% GENERATE UNIFORM GRID IN KM SPACE
[X,Y] = meshgrid(xg,yg) ;
% CONVERT BACK TO GEOGRAPHICAL COORDINATES
LAT = 180*Y./(R*pi) ;
LNG = 180*X./(R*pi*cosd(LAT)) + avl ; % moved back to average meridian
end
\ No newline at end of file
function A = ll2area(lng,lat,unit)
% Synthax : A = ll2area(lng,lat,unit)
%
% Takes as input the geographical coordinate vectors for longitude 'lng'
% and latitude 'lat' in column format and returns the area inside the
% non-intersecting polygon they form as ouput 'A' in units 'unit'.
% recognized unit options are 'm' for square meters and 'km' for square
% kilometers.
%
% BASED ON:
%
% https://www.periscopedata.com/blog/polygon-area-from-latitude-and-longitude-using-sql
%
% and
%
% http://geomalgorithms.com/a01-_area.html
% PARAMETERS
switch unit
case 'km'
R = 6371 ; % earth radius (km)
case 'm'
R = 6371000 ; % earth radius (m)
end
% CONVERT TO XY
y = ( lat*R*pi / 180 ) ;
x = ( cosd(lat).*lng*R*pi / 180 ) ;
avx = mean(x) ;
avy = mean(y) ;
normx = x - avx ;
normy = y - avy ;
% GET ANGLE FOR EVERY VERTEX (0-360) and convert
ang = atan2d(normy,normx); % + (-1*normx./abs(normx))*90 + 180 ;
ang( ang < 0 ) = ang( ang < 0 ) + 360 ;
ang = ang - 90;
ang = 360 - ang ;
ang(ang < 0) = ang(ang < 0) + 360 ;
ang = mod(ang,360) ;
% SORT TO AVOID SELF INTERSECTION
[~,IA] = sort(ang,'descend') ;
x = x(IA) ;
y = y(IA) ;
% APPEND FOR WRAPPING
x1 = x(1); xn = x(end) ;
y1 = y(1); yn = y(end) ;
x = [xn; x; x1] ;
y = [yn; y; y1] ;
% COMPUTE THE SUM
A = 0 ;
for ii = 2:numel(x)-1
A = A + 0.5*( x(ii).*( y(ii+1) - y(ii-1) ) ) ;
end
end
function [u,v,spd,head] = llt2spd(lon,lat,time)
% Synthax : [u,v,spd,head] = llt2spd(lon,lat,time)
%
% Returns the speed 'spd' and it's eastward/northward components 'u' and
% 'v' in (m/s), and the heading 'head' in degrees with 0 as east increasing
% towards the north (90 degrees).
%
% Values are calculated from a latitude/longitude time series with
% coordinates in decimal degrees and time in days.
%
% Values are estimated between longitude/latitude coordinates and then
% interpolated linearly to the input 'time' vector. Begining and end
% values are extrapolated using the "pchip" method.
%% make everything column vectors
if ~iscolumn(lon) ; lon = lon' ; end
if ~iscolumn(lat) ; lat = lat' ; end
if ~iscolumn(time) ; time = time' ; end
%% make differential vectors
dlon = diff(lon) ;
dlat = diff(lat) ;
dt = diff(time) ;
dl = m_lldist(lon,lat)*1000 ; % (m)
%% calculate norm of speed
spd = dl./(dt*24*60*60) ;
%% calculate heading
head = range_pm180_2_360(atan2d(dlat,dlon)) ;
%% decompose into eastward/northward
u = spd.*cosd(head) ;
v = spd.*sind(head) ;
%% make the time vector
t = time(1:end-1)+0.5*dt(1);
%% interpolate to the input time grid
u = interp1(t,u,time,'linear','extrap') ;
v = interp1(t,v,time,'linear','extrap') ;
spd = interp1(t,spd,time) ;
spd(1) = sqrt(u(1).^2+v(1).^2);
spd(end) = sqrt(u(end).^2+v(end).^2);
head = interp1(t,head,time) ;
end
function output = range_360_2_pm180(input)
I = find(input > 180) ;
input(I) = input(I) - 360 ;
output = input ;
end
function output = range_pm180_2_360(input)
I = find(input < 0) ;
input(I) = input(I) + 360 ;
output = input ;
end
function output = theta2heading(input)
% Synthax : output = theta2heading(input)
%
% Takes as input angle vector in cartesian coordinates running counter-
% clockwise from 0 (east) to 360 degrees and transforms it into a heading
% vector running clockwise from 0 to the north to 360 degrees.
%
input = input - 90;
input = 360 - input ;
input(input < 0) = input(input < 0) + 360 ;
output = mod(input,360) ;
end
function [bins,bs,m,M]=SEplot(dist)
% Synthax : [bins,bs,m,M]=SEplot(dist)
%
% Produces a histogram plot of distribution 'dist' with a color code
% identifying which portions of the distribution correspond to +-1 standard
% error, +-2 standard error, +-3 standard error, and beyond. Limits are
% included to the categories and standard error is calculated with the
% built in matlab std function
sigma = std(dist,'omitnan') ;
mdist = mean(dist,'omitnan') ;
M = max(dist) ;
m = min(dist) ;
bs = sigma./20 ;
bins = sort([mdist:-bs:m mdist+bs:bs:M]);
figure ;
i = find(dist <= mdist + 1*sigma & dist >= mdist - 1.*sigma) ;
histogram(dist(i),bins,'facecolor','b') ;
hold on
i = find(dist <= mdist + 2*sigma & dist > mdist + 1.*sigma | ...
dist < mdist - 1*sigma & dist >= mdist - 2.*sigma) ;
histogram(dist(i),bins,'facecolor','g') ;
i = find(dist <= mdist + 3*sigma & dist > mdist + 2.*sigma | ...
dist < mdist - 2*sigma & dist >= mdist - 3.*sigma) ;
histogram(dist(i),bins,'facecolor','r') ;
i = find(dist <= mdist - 3*sigma | dist >= mdist + 3.*sigma) ;
histogram(dist(i),bins,'facecolor','c') ;
end
function [dist,I] = SEtrunc(dist,thres)
% Synthax : [dist,I] = SEtrunc(dist,thres)
%
% Replaces with nan the values in dist which are outside the range
%
% +-'thres'*std('dist')
%
% where std is matlab built in function. Limits of the range are excluded.
% also returns the index values 'I' of the values nan-ed.
sigma = std(dist,'omitnan') ;
mdist = mean(dist,'omitnan') ;
I = find(dist < mdist - thres.*sigma | dist > mdist + thres.*sigma) ;
dist(I) = nan ;
end
function [vals,vecs,stds] = jlpca(X)
% Synthax : [vals,vecs,stds] = jlpca(X)
%
% Performs principal component analysis on input matrix 'X' . Input matrix
% must be in a form where columns are different features of the data set
% and rows are repetitions of the same measurement.
%
% Also returns the standard deviation of every feature projected in the space
% of the eigenvectors found by the PCA.
%
[m,n] = size(X) ;
%% De-trend
F = repmat(nan,size(X)) ;
for ii = 1:n
F(:,ii) = X(:,ii) - mean(X(:,ii),'omitnan') ;
end
%% Build covariance matrix
R = F'*F ;
%% Find eigenvalues
[C,L] = eig(R) ;
%% sort and normalise
vals = sort(diag(L),'descend') ;
% eigenvectors in growing importance
vecs = fliplr(C) ;
%% Project in the space of the eigenvectors
Fprime = inv(vecs)*F' ;
%% Get the standard error in each "eigendimension"
stds = std(Fprime,[],2) ;
end
function output = nanmean(input,varargin)
%
% Matlab is super cheap for not including this
% in the student license. Seriously.
%
I = find(isfinite(input)) ;
output = mean(input(I),varargin{:}) ;
end
function out = nanmedian(in)
I = find(isfinite(in)) ;
out = median(in(I)) ;
end
function [yfit,A] = nlinfit(rhs,x,y,pars,varargin)
%
% Synthax : [yfit,A] = nlinfit(rhs,x,y,pars,cst1,cst2,...,cstn)
%
% An idea found online and remodeled to make it more modular about how to
% fminsearch.m to perform non-linear, multiple parameter function fitting.
% Variables 'x' and 'y' are the input data, i.e., y = y(x) and must be
% vectors. Variable 'rhs' is a function handle refering to te mathematical
% form of the fitted function. For example, to fit a 1/6 power law, define:
%
% rhs = @(a,x) a(1)*x.^(1/6) + a(2) ;
%
% where 'a' is a vector formed of the parameters to fit and 'x' is the
% input independent variable. If variable 'pars' is entered as an integer
% it is interpreted as the number of parameters to fit, and sets the first
% guess for all parameters to zero. If 'pars' is entered as a vector, then
% it is interpreted as the first guess for all parameters.
%
% The function returns the fitted values at 'x' as well as the optimized
% parameter vector 'A'.
% Is pars first guess for parameters or number of parameters
sz = size(pars) ;
if sz(1) == 1 & sz(2) == 1
start = zeros([pars 1]) ;
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{:}) ;
end
function output = smooth2D(input,hs,vs)
% Synthax : output = smooth2D(input,hs,vs)
%
% 2D moving average. The window size is 'hs' in direction j and 'vs' in
% direction i. Both 'hs' and 'vs' MUST be odd integers. Each element of
% 'output' is calculated from the nanmean of a box of size 'hs' by 'vs'
% centered on the element of corresponding indices in 'input'.
output = repmat(nan,size(input)) ;
hpad = (hs-1)/2 ;
vpad = (vs-1)/2 ;
if vpad < 1 ; vpad = 0 ; end;
if hpad < 1 ; hpad = 0 ; end;
if iseven(vpad) ; vpad = vpad + 1 ; end
if iseven(hpad) ; hpad = hpad + 1 ; end
[II,JJ] = size(input);
for ii = 1:II
for jj = 1:JJ
if ii <= vpad & jj > hpad & jj <= JJ-hpad % TOP
box = input(1:ii+vpad,jj-hpad:jj+hpad) ;
elseif ii > II-vpad & jj > hpad & jj <= JJ-hpad % BOTTOM
box = input(ii-vpad:II,jj-hpad:jj+hpad) ;
elseif jj <= hpad & ii > vpad & ii <= II-vpad % LEFT
box = input(ii-vpad:ii+vpad,1:jj+hpad) ;
elseif jj > JJ-hpad & ii > vpad & ii <= II-vpad % RIGHT
box = input(ii-vpad:ii+vpad,jj-hpad:JJ) ;
elseif ii <= vpad & jj <= hpad % TOP LEFT
box = input(1:ii+vpad,1:jj+hpad) ;
elseif ii > II-vpad & jj <= hpad % BOTTOM LEFT
box = input(ii-vpad:II,1:jj+hpad) ;
elseif ii <= vpad & jj > JJ-hpad % TOP RIGHT
box = input(1:ii+vpad,jj-hpad:JJ) ;
elseif ii > II-vpad & jj > JJ-hpad % BOTTOM RIGHT
box = input(ii-vpad:II,jj-hpad:JJ) ;
else % INSIDE
box = input(ii-vpad:ii+vpad,jj-hpad:jj+hpad) ;
end
output(ii,jj) = mean(box(:),'omitnan') ;
end
end
end
function output = unifrnd(b_min,b_max,d_len,m,n)
% Synthax : output = unifrnd(b_min,b_max,d_len,m,n)
%
% Imitates the unifrnd function from the statistics toolbox. Generates a
% 'm' by 'n' matrix of values chosen from a uniform distribution of minimum
% 'b_min' and maximum 'b_max' . The parameter 'd_len' controls the number
% of elements in the uniform distribution the 'output' values are chosen
% from.
%
dist = linspace(b_min,b_max,d_len)' ;
I = randi([1 d_len],[m n]) ;
output = dist(I) ;
end
\ No newline at end of file
function theta = vvangle(v1,v2)
% Synthax : theta = vvangle(v1,v2)
%
% Get the absolute angle between two line vectors, or between the vectors
% described by the lines v1 and v2.
%
%% make sure the sizes are compatible
if ~( isequal(size(v1),size(v2)) || isvector(v1) && isvector(v2) && numel(v1) == numel(v2) )
disp(' Error : v1 and v2 have different dimensions!') ;
return ;
end
%% get the vector norms
nv1 = sum(abs(v1).^2,2).^(1/2) ;
nv2 = sum(abs(v2).^2,2).^(1/2) ;
%% calculate the angle
theta = acosd( ( dot(v1,v2,2) ) ./ (nv1.*nv2) ) ;
end
\ No newline at end of file
function ADCP_cleanup(ADCP_file_path,ADCP_file_name,BATH_file,CTD_file,GPS_file,strtsamp,stopsamp,theta,corrtag,CTDtag,WD,zerolon,zerolat)
disp('========== 1) ADCP_cleanup ==========') ;
tic
%% PARAMETERS
paperwidth = 20 ;
paperheight = 6 ;
%% OUTPUT FILE TAGS
CC = ['CC' strrep(num2str(theta),'.','') 'd'] ;
%% LOAD data
disp('Loading data ...')
load([WD ADCP_file_path ADCP_file_name]) ;
bathystruc = load(BATH_file) ;
bathyfield = fieldnames(bathystruc) ;
BATHY = getfield(bathystruc,bathyfield{:}) ;
% dataset decimation (for testing)
loaded = 1:1:numel(Sens.dnum) ;
%% CREATE GRAPHICS AND DATA FOLDERS
if 7 == exist([WD 'figs/ADCP_cleanup'])
disp('Graphics folder found.') ;
else
disp('Creating graphics folder : figs/ADCP_cleanup')
mkdir([WD 'figs/ADCP_cleanup']);
end
if 7 == exist([WD 'data'])
disp('data folder found.') ;
else
disp('Creating data folder : data')
mkdir([WD 'data']);
end
%% temps ADCP
I = find(Sens.dnum(loaded) < strtsamp | Sens.dnum(loaded) > stopsamp) ;
loaded(I) = [] ;
time = Sens.dnum(loaded) ;
%% matrices des vitesses
u_raw = Wt.vel(loaded,:,1)' ; % m/s
v_raw = Wt.vel(loaded,:,2)' ; % m/s
w_raw = Wt.vel(loaded,:,3)' ; % m/s
disp('Done')
%% vitesse du bateau
disp('Making boat speed matrices ...')
fid = fopen([WD GPS_file]) ;
frmt = ['%s%s%f%f%s%s' repmat('%s',[1 23])] ;
data = textscan(fid,frmt,'headerlines',49,'delimiter',',') ;
fclose(fid) ;
[boat_tim,IA,~] = unique(datenum(cell2mat(data{6}),'yyyy-mm-ddTHH:MM:SS'));
boat_lon = data{4}(IA) ;
boat_lat = data{3}(IA) ;
% calculer la vitesse du bateau
[boat_u,boat_v,boat_spd,boat_head] = llt2spd(boat_lon,boat_lat,boat_tim) ;
boat_acc = sqrt(gradient(boat_u,boat_tim*(24*3600)).^2 + gradient(boat_v,boat_tim*(24*3600)).^2) ;
% interpoler sur le temps de l ADCP
boat_lon = interp1(boat_tim,boat_lon,time) ;
boat_lat = interp1(boat_tim,boat_lat,time) ;
boat_u = interp1(boat_tim,boat_u,time) ;
boat_v = interp1(boat_tim,boat_v,time) ;
boat_head = interp1(boat_tim,boat_head,time) ;
boat_acc = interp1(boat_tim,boat_acc,time) ;
boat_spd = interp1(boat_tim,boat_spd,time) ;
% creer les matrices de vitesse de bateau de la taille de vitesses raw
boat_u_mat = repmat(boat_u',[numel(u_raw(:,1)) 1]) ;
boat_v_mat = repmat(boat_v',[numel(v_raw(:,1)) 1]) ;
%% positions du centre des bins
% une profondeur du poisson de 1 m est supposee.
z = 1 + Info.cell1 + 0.5*Info.cell ...
+ [0 : Info.cell : 93*Info.cell] ; % ecrit ainsi par souci de clarte
%% raw data depths
depth = BATHY(boat_lon,boat_lat) ;
clear BATHY ;
%% PRINT THE BOAT SPEED MATRICES
disp('Printing boat speed matrices ... ')
figure('visible','off');
pcolor(time,z,boat_u_mat)
hold on
plot(time,depth*cosd(25),'k')
plot(time,depth,'k--')
datetick('x','HH:MM')
xlabel('Time')
ylabel('Depth (m)')
shading flat ;
axis ij;
ylim([0 80]) ; xlim([min(time) max(time)]) ;
caxis([-1 1])
jlcolorbar('u') ;
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperSize',[paperwidth paperheight]);
set(gcf,'PaperPosition',[0 0 paperwidth paperheight]);
print('-dpng', '-r300',[WD 'figs/ADCP_cleanup/u_boat_speed_' CC '_' corrtag '.png']);
close;
figure('visible','off');
pcolor(time,z,boat_v_mat)
hold on
plot(time,depth*cosd(25),'k')
plot(time,depth,'k--')
datetick('x','HH:MM')
xlabel('Time')
ylabel('Depth (m)')
shading flat ;
axis ij;
ylim([0 80]) ; xlim([min(time) max(time)]) ;
caxis([-1 1])
jlcolorbar('v') ;
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperSize',[paperwidth paperheight]);
set(gcf,'PaperPosition',[0 0 paperwidth paperheight]);
print('-dpng', '-r300',[WD 'figs/ADCP_cleanup/v_boat_speed_' CC '_' corrtag '.png']);
close;
disp('Done.') ;
%% PRINT THE RAW DATA
disp('Printing the raw data ... ')
figure('visible','off');
pcolor(time,z,u_raw)
hold on
plot(time,depth*cosd(25),'k')
plot(time,depth,'k--')
datetick('x','HH:MM')
xlabel('Time')
ylabel('Depth (m)')
shading flat ;
axis ij;
ylim([0 80]) ; xlim([min(time) max(time)]) ;
caxis([-1 1])
jlcolorbar('u') ;
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperSize',[paperwidth paperheight]);
set(gcf,'PaperPosition',[0 0 paperwidth paperheight]);
print('-dpng', '-r300',[WD 'figs/ADCP_cleanup/u_raw_' CC '_' corrtag '.png']);
close;
figure('visible','off');