Commit c4f12544 authored by Dany Dumont's avatar Dany Dumont
Browse files

Ajout de run_process_ADCP.m + nettoyage et documentation de plusieurs fonctions

parent 0441d923
function [adcp,cfg,ens]=rdradcp(name,varargin);
% RDRADCP Read (raw binary) RDI ADCP files,
% ADCP=RDRADCP(NAME) reads the raw binary RDI BB/Workhorse ADCP file NAME and
function [adcp,cfg,ens] = rdradcp(name,varargin)
% RDRADCP - Read (raw binary) RDI ADCP files,
% ADCP = RDRADCP(NAME) reads the raw binary RDI BB/Workhorse ADCP file NAME and
% puts all the relevant configuration and measured data into a data structure
% ADCP (which is self-explanatory). This program is designed for handling data
% recorded by moored instruments (primarily Workhorse-type but can also read
......@@ -83,11 +83,11 @@ century=2000; % ADCP clock does not have century prior to firmware 16.05.
vels='no'; % Default to simple averaging
lv=length(varargin);
if lv>=1 & ~isstr(varargin{1}),
if lv>=1 && ~isstr(varargin{1}),
num_av=varargin{1}; % Block filtering and decimation parameter (# ensembles to block together).
varargin(1)=[];
lv=lv-1;
if lv>=1 & ~isstr(varargin{1}),
if lv>=1 && ~isstr(varargin{1}),
nens=varargin{1};
varargin(1)=[];
lv=lv-1;
......@@ -459,7 +459,7 @@ if num_av<0,
SOURCE=0;
end;
% This reinitializes to whatever length of ens we want to average.
if num_av<0 | isempty(ens),
if num_av<0 || isempty(ens),
FIXOFFSET=0;
n=abs(num_av);
[hdr,pos]=rd_hdr(fd);
......
function S = enx2mat(file,umax,errmax,trans_depth)
% S = enx2mat(file,umax,errmax,trans_depth)
%
% Load RDI raw ADCP data file (.ENX) into a
% ENX2MAT - Load RDI raw ADCP data file (.ENX) into a
% Matlab structure S.
%
% Syntax: S = enx2mat(file,umax,errmax,trans_depth)
%
% Correct ship motion with bottom tracking.
% Remove velocities with errors larger than errmax (m/s)
% and speed larger than umax (m/s).
......@@ -11,7 +11,7 @@ function S = enx2mat(file,umax,errmax,trans_depth)
% with trans_depth (m).
% Cdric Chavanne, updated 2016/11/6.
% ______________________________________________________________________
% load raw data (.ENX file) without ensemble averaging (NUMAV = 1):
S = rdradcp(file,1);
......@@ -31,22 +31,20 @@ else
z = trans_depth + S.config.ranges;
end
% check attitude data:
figure(1)
clf
subplot(211)
% Check attitude data
figure(1); clf
subplot(2,1,1)
plot(t,S.heading)
ylabel('heading')
subplot(212)
subplot(2,1,2)
plot(t,S.pitch,'r',t,S.roll,'b')
ylabel('pitch (r), roll (b)')
% estimate bottom depth from bottom track data:
% Estimate bottom depth from bottom track data:
S.bt_range(find(S.bt_range < z(1))) = nan;
zbot = nanmedian(S.bt_range);
figure(2)
clf
figure(2); clf
plot(t,S.bt_range)
hold on
plot(t,zbot,'k','linewidth',2)
......@@ -55,8 +53,7 @@ plot(t,zbot,'k','linewidth',2)
zgood = zbot * cosd(S.config.beam_angle);
% check bottom track velocities:
figure(3)
clf
figure(3); clf
plot(t,S.bt_vel)
legend('u','v','w','err')
% bottom track vel are in mm/s
......@@ -65,36 +62,37 @@ legend('u','v','w','err')
J = find(abs(S.bt_vel(4,:)) > 1e2);
S.bt_vel(:,J) = nan;
figure(3)
clf
figure(3); clf
plot(t,S.bt_vel)
legend('u','v','w','err')
% check water velocities:
figure(4)
clf
subplot(411)
figure(4); clf
subplot(4,1,1)
pcolor(t-datenum(year,0,0),-z,S.east_vel)
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k')
datetick('x',15)
subplot(412)
subplot(4,1,2)
pcolor(t-datenum(year,0,0),-z,S.north_vel)
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k')
datetick('x',15)
subplot(413)
subplot(4,1,3)
pcolor(t-datenum(year,0,0),-z,S.vert_vel)
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k')
datetick('x',15)
subplot(414)
subplot(4,1,4)
pcolor(t-datenum(year,0,0),-z,S.error_vel)
shading flat
colorbar
......@@ -105,30 +103,32 @@ xlabel(datestr(t(1),1))
% water vel are in m/s and are not corrected for ship motion
% check correlations:
figure(5)
clf
subplot(411)
figure(5); clf
subplot(4,1,1)
pcolor(t-datenum(year,0,0),-z,squeeze(S.corr(:,1,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(412)
subplot(4,1,2)
pcolor(t-datenum(year,0,0),-z,squeeze(S.corr(:,2,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(413)
subplot(4,1,3)
pcolor(t-datenum(year,0,0),-z,squeeze(S.corr(:,3,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(414)
subplot(4,1,4)
pcolor(t-datenum(year,0,0),-z,squeeze(S.corr(:,4,:)))
shading flat
colorbar
......@@ -137,31 +137,33 @@ plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
xlabel(datestr(t(1),1))
% check intensities:
figure(6)
clf
subplot(411)
% Check intensities
figure(6); clf
subplot(4,1,1)
pcolor(t-datenum(year,0,0),-z,squeeze(S.intens(:,1,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(412)
subplot(4,1,2)
pcolor(t-datenum(year,0,0),-z,squeeze(S.intens(:,2,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(413)
subplot(4,1,3)
pcolor(t-datenum(year,0,0),-z,squeeze(S.intens(:,3,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(414)
subplot(4,1,4)
pcolor(t-datenum(year,0,0),-z,squeeze(S.intens(:,4,:)))
shading flat
colorbar
......@@ -170,31 +172,33 @@ plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
xlabel(datestr(t(1),1))
% check status:
figure(7)
clf
subplot(411)
% Check status
figure(7); clf
subplot(4,1,1)
pcolor(t-datenum(year,0,0),-z,squeeze(S.status(:,1,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(412)
subplot(4,1,2)
pcolor(t-datenum(year,0,0),-z,squeeze(S.status(:,2,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(413)
subplot(4,1,3)
pcolor(t-datenum(year,0,0),-z,squeeze(S.status(:,3,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
subplot(414)
subplot(4,1,4)
pcolor(t-datenum(year,0,0),-z,squeeze(S.status(:,4,:)))
shading flat
colorbar
......@@ -218,9 +222,8 @@ for j = 1:length(zgood)
err(I,j) = nan;
end
figure(5)
clf
subplot(411)
figure(5); clf
subplot(4,1,1)
pcolor(t-datenum(year,0,0),-z,u)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
......@@ -229,7 +232,8 @@ shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(zgood) 0])
subplot(412)
subplot(4,1,2)
pcolor(t-datenum(year,0,0),-z,v)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
......@@ -238,7 +242,8 @@ shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(zgood) 0])
subplot(413)
subplot(4,1,3)
pcolor(t-datenum(year,0,0),-z,w)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
......@@ -247,7 +252,8 @@ shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(zgood) 0])
subplot(414)
subplot(4,1,4)
pcolor(t-datenum(year,0,0),-z,err)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
......@@ -258,9 +264,8 @@ datetick('x',15)
xlabel(datestr(t(1),1))
set(gca,'ylim',[-max(zgood) 0])
% remove velocities with large errors:
figure(8)
clf
% Remove velocities with large errors
figure(8); clf
hist(abs(err(:)),128)
I = find(abs(err) > errmax);
......@@ -270,23 +275,23 @@ w(I) = nan;
err(I) = nan;
% check histograms of velocities:
figure(9)
clf
subplot(311)
figure(9); clf
subplot(3,1,1)
hist(u(:),128)
xlabel('u (m/s)')
subplot(312)
subplot(3,1,2)
hist(v(:),128)
xlabel('v (m/s)')
subplot(313)
subplot(3,1,3)
hist(w(:),128)
xlabel('w (m/s)')
% remove large velocities:
speed = sqrt(u.^2+v.^2+w.^2);
figure(10)
clf
figure(10); clf
hist(speed(:),128)
I = find(speed > umax);
......@@ -295,22 +300,22 @@ v(I) = nan;
w(I) = nan;
err(I) = nan;
figure(9)
clf
subplot(311)
figure(9); clf
subplot(3,1,1)
hist(u(:),128)
xlabel('u (m/s)')
subplot(312)
subplot(3,1,2)
hist(v(:),128)
xlabel('v (m/s)')
subplot(313)
subplot(3,1,3)
hist(w(:),128)
xlabel('w (m/s)')
% plot velocities:
figure(5)
clf
subplot(411)
% Plot velocities
figure(5); clf
subplot(4,1,1)
pcolor(t-datenum(year,0,0),-z,u)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
......@@ -321,7 +326,7 @@ datetick('x',15)
set(gca,'ylim',[-max(zgood) 0],'xticklabel',{})
title('u')
subplot(412)
subplot(4,1,2)
pcolor(t-datenum(year,0,0),-z,v)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
......@@ -332,7 +337,7 @@ datetick('x',15)
set(gca,'ylim',[-max(zgood) 0],'xticklabel',{})
title('v')
subplot(413)
subplot(4,1,3)
pcolor(t-datenum(year,0,0),-z,w)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
......@@ -343,7 +348,7 @@ datetick('x',15)
set(gca,'ylim',[-max(zgood) 0],'xticklabel',{})
title('w')
subplot(414)
subplot(4,1,4)
pcolor(t-datenum(year,0,0),-z,err)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
......@@ -355,17 +360,14 @@ xlabel(datestr(t(1),1))
set(gca,'ylim',[-max(zgood) 0])
title('err')
% output variables:
S = [];
S.u = u;
S.v = v;
S.w = w;
S.err = err;
S.z = z;
S.t = t;
S.lon = lon;
S.lat = lat;
% Output variables
S = [];
S.u = u;
S.v = v;
S.w = w;
S.err = err;
S.z = z;
S.t = t;
S.lon = lon;
S.lat = lat;
S.zbot = zbot;
% Load RDI 150-kHz ADCP data for each station
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function load_ADCP(idctd,idadcp)
% LOAD_ADCP - Load RDI 150-kHz ADCP data for each station
%
% Syntax: load_ADCP(idctd,idadcp)
%
% Inputs:
% idctd - String for station number
% idadcp - String for file number corresponding to station
%
% Other m-files required: RDRADCP, ENX2MAT
% Subfunctions: none
% MAT-files required:
% Other files required: ADCP and CTD files
% Cdric Chavanne
% ______________________________________________________________________
rootdir = '/Users/dany/Dropbox/cours/OCE-62014_Experimentale/MS2019';
addpath /Users/dany/Dropbox/cours/OCE-62014_Experimentale/MS2019/mfiles/ -end
addpath /Users/dany/Dropbox/cours/OCE-62014_Experimentale/MS2019/mfiles/ReadADCP/ -end
%addpath /Users/dany/Dropbox/cours/OCE-62014_Experimentale/MS2019/mfiles/RDADCP_Jan18v0/ -end
umax = 2.0; % max water velocity (m/s)
errmax = 0.2; % max error velocity (m/s)
%% Leg 1, day 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
% max water velocity (m/s):
umax = 2;
% max error velocity (m/s):
errmax = 0.2;
% file ID numbers for stations:
idctd = {'F9'};
idadcp = {'001_000000'};
ctdfile=['../data/ctd/mat/MS2019_' idctd '.mat'];
if ~exist(ctdfile)
warning('CTD file not found')
else
load(ctdfile)
CTD = S; clear S
end
for i = 1:length(idctd)
ctdfile=['../data/ctd/proc/MS2019_' idctd{i} '.mat'];
if ~exist(ctdfile)
warning('CTD file not found')
else
load(ctdfile)
CTD = S; clear S
end
adcpfile=['../data/adcp/raw/MS2019-150kHz_20190914T112422_' idadcp{i} '.ENX'];
if ~exist(adcpfile)
warning('ADCP file not found')
adcpfile=['../data/adcp/raw/MS2019-150kHz_' idadcp '.ENX'];
if ~exist(adcpfile)
error('ADCP file not found')
else
ADCP = enx2mat(adcpfile,umax,errmax);
% Select period around CTD cast
I = find(ADCP.t >= CTD.t(1)-5/(24*60) & ADCP.t <= CTD.t(end)+5/(24*60));
if isempty(I)
warning('ADCP time does not match CTD time')
else
ADCP = enx2mat(adcpfile,umax,errmax);
% select period around CTD cast:
I = find(ADCP.t >= CTD.t(1)-5/(24*60) & ADCP.t <= CTD.t(end)+5/(24*60));
if isempty(I)
warning('ADCP time does not match CTD time')
else
% plot velocities:
figure(11)
clf
subplot(411)
pcolor(ADCP.t(I)-ADCP.t(I(1)),-ADCP.z,ADCP.u(:,I))
hold on
plot(ADCP.t(I)-ADCP.t(I(1)),-ADCP.zbot(I),'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(ADCP.zbot(I)) 0],'xticklabel',{})
title('u')
subplot(412)
pcolor(ADCP.t(I)-ADCP.t(I(1)),-ADCP.z,ADCP.v(:,I))
hold on
plot(ADCP.t(I)-ADCP.t(I(1)),-ADCP.zbot(I),'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(ADCP.zbot(I)) 0],'xticklabel',{})
title('v')
subplot(413)
pcolor(ADCP.t(I)-ADCP.t(I(1)),-ADCP.z,ADCP.w(:,I))
hold on
plot(ADCP.t(I)-ADCP.t(I(1)),-ADCP.zbot(I),'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(ADCP.zbot(I)) 0],'xticklabel',{})
title('w')
subplot(414)
pcolor(ADCP.t(I)-ADCP.t(I(1)),-ADCP.z,ADCP.err(:,I))
hold on
plot(ADCP.t(I)-ADCP.t(I(1)),-ADCP.zbot(I),'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(ADCP.zbot(I)) 0],'xticklabel',{})
title('err')
end
% Plot velocities
figure(11); clf
subplot(4,1,1)
pcolor(ADCP.t(I)-ADCP.t(I(1)),-ADCP.z,ADCP.u(:,I))
hold on
plot(ADCP.t(I)-ADCP.t(I(1)),-ADCP.zbot(I),'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(ADCP.zbot(I)) 0],'xticklabel',{})
title('u')
% time-average currents
u = nanmedian(ADCP.u(:,I),2);
v = nanmedian(ADCP.v(:,I),2);
lon = nanmedian(ADCP.lon(I));
lat = nanmedian(ADCP.lat(I));
t = nanmedian(ADCP.t(I));
z = ADCP.z;
nobs = sum(~isnan(ADCP.u(:,I)),2);
subplot(4,1,2)
pcolor(ADCP.t(I)-ADCP.t(I(1)),-ADCP.z,ADCP.v(:,I))
hold on
plot(ADCP.t(I)-ADCP.t(I(1)),-ADCP.zbot(I),'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(ADCP.zbot(I)) 0],'xticklabel',{})
title('v')
figure(12)
clf
plot(nobs,-z)
subplot(4,1,3)
pcolor(ADCP.t(I)-ADCP.t(I(1)),-ADCP.z,ADCP.w(:,I))
hold on
plot(ADCP.t(I)-ADCP.t(I(1)),-ADCP.zbot(I),'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(ADCP.zbot(I)) 0],'xticklabel',{})
title('w')
% remove data with insufficient observations
J = find(nobs < 0.5*max(nobs));
u(J) = nan;
v(J) = nan;
subplot(4,1,4)
pcolor(ADCP.t(I)-ADCP.t(I(1)),-ADCP.z,ADCP.err(:,I))
hold on
plot(ADCP.t(I)-ADCP.t(I(1)),-ADCP.zbot(I),'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(ADCP.zbot(I)) 0],'xticklabel',{})
title('err')
figure(13)
clf
plot(u,-z,'b',v,-z,'r')
set(gca,'ylim',[-max(ADCP.zbot(I)) 0])
% save results
fileout = ['../data/adcp/proc/MS2019-150kHz_20190914T112422_' idctd{i}];
save(fileout,'u','v','lon','lat','z','t')
pause
end
% time-average currents
u = nanmedian(ADCP.u(:,I),2);