Commit 0441d923 authored by Dany Dumont's avatar Dany Dumont
Browse files

Premier depot

parents
% ADCPDEMO Demo script for reading ADCP data
%
%******************Demo for ADCP reading programs*******************
%
% This package consists of two functions that read binary data generated by
% RDI ADCPs into a matlab data structure. Some averaging and subsampling
% can be performed but after that data analysis and visualization is
% up to you.
%
% Two kinds of RDI files can be read with this software:
%
% 1) "Raw" files (usually downloaded from moored instruments) using
% the BBDF format (compatible with BB and WH ADCPs) for (hopefully
% all) firmware versions <=16.19.
%
% Note that slight modifications of this "raw" format are used
% by WINRIVER and VMDAS.
%
% 2) "P-files" (processed files generated by TRANSECT or WINRIVER after
% processing data acquired in real-time. P-files are preferable
% in this case because navigation data and configuration data is
% merged with the raw data). NB, BB, and WH instrument files can
% be read.
%
%
% This demo script will provide some examples of how to use the code.
%
% WARNING: RDI is constantly improving their products, and so the format
% of their binary files is constantly changing. I have done the
% best I can to track down all firmware and program version
% dependencies and this code has been tested on datasets from
% a large-ish number of researchers and instruments but it
% may still fail on your file. If it does, please let me know
% so I can figure out why!
%
% Author: R. Pawlowicz (rich@eos.ubc.ca, www.eos.ubc.ca/~rich/research.html)
% [some support for parts of this work was provided by RESCAN
% ENVIRONMENTAL SERVICES www.rescan.com]
% Changes - clc instead of escape code wizardry for clearing screen , thanks to F. Bahr
more on
help adcpdemo
fprintf('--> Now we being the demonstration.....\n');
disp('press any key to continue');
pause;
clc;
fprintf('--> Default read (5-ping averages of whole file) \n');
fprintf('--> [This data segment was provided by R. Dewey] \n');
fprintf('\n\necho on\n');
echo on;
adcp=rdradcp('mooredwh-bbdf.000');
echo off;
adcp
fprintf('--> Once read this is what the ADCP data structure looks like!\n');
disp('press any key to continue');
pause;
clc;
adcp.config
fprintf('--> Configuration data is stored in a sub-structure "adcp.config"\n');
fprintf('--> Note that this particular dataset was recorded in\n');
fprintf('--> Beam coordinates by someone who preferred doing the\n');
fprintf('--> velocity mapping themselves for this deployment\n');
disp('press any key to continue');
pause;
clc;
fprintf('echo on\n');
echo on;
clf
subplot(311);
plot(adcp.number,adcp.roll);
ylabel('Roll angle');
title('Deployment of ADCP');
subplot(312);
plot(adcp.number,adcp.pitch);
ylabel('Pitch angle');
subplot(313);
pcolor(adcp.number,adcp.config.ranges,squeeze(mean(adcp.intens,2)));
shading flat;
ylabel('Range (m)');
shg;
echo off;
fprintf('--> A simple check of some of the raw data to see how\n');
fprintf('--> the instrument behaved on deployment \n');
fprintf('--> (the surface can be seen in mean backscatter) \n\n\n');
fprintf('--> Now let''s load some real-time data - we shall read\n');
fprintf('--> the first 100 pings of the file, NOT use bottom or\n');
fprintf('--> navigation as a reference, and perform a de-spiked 10-ping average. \n');
disp('press any key to continue');
pause;
fprintf('\n\necho on\n');
echo on;
adcp=rdpadcp('realtimewh-p.000',10,100,'ref','none','despike','yes');
clf
plot(adcp.north_vel,adcp.config.range);
set(gca,'ydir','reverse');
xlabel('Velocity (m/s)');ylabel('depth (m)');
title('ADCP NORTH VELOCITY','fontsize',16);
shg;
echo off
fprintf('--> ...And we plot the north velocity of all pings.\n');
fprintf('--> ..Simple!\n');
fprintf('--> - Enjoy, RP. \n');
more off
fname = '1001_000.mat';
umin = -1.0;
umax = 1.0;
wmin = -0.15;
wmax = 0.1;
load(fname);
% Set the z axis
z0 = 0.8; % Depth of the transducer. This info can come from the ADCP.depth data
z = ADCP.config.ranges + z0;
% Set the time axis
mtime = ADCP.mtime;
dvec = datevec(mtime(1));
mtime0 = datenum(dvec(1),dvec(2),dvec(3),0,0,0);
% Velocity
U = ADCP.east_vel;
V = ADCP.north_vel;
W = ADCP.vert_vel;
% Backscatter intensity
I1 = ADCP.intens(:,1,:);
I2 = ADCP.intens(:,2,:);
I3 = ADCP.intens(:,3,:);
I4 = ADCP.intens(:,4,:);
I1 = squeeze(I1);
I2 = squeeze(I2);
I3 = squeeze(I3);
I4 = squeeze(I4);
Imean = (I1 + I2 + I3 + I4)/4;
figure(1);
pcolor(mtime-mtime0,z,I1);
shading('interp');
set(gca,'ydir','reverse');
set(gca,'XTickLabel',[]);
ylabel('Depth (m)');
datetick('x',15);
xlabel(['Time on ',datestr(mtime0)]);
colorbar;
figure(2);
pcolor(mtime-mtime0,z,I2);
shading('interp');
set(gca,'ydir','reverse');
set(gca,'XTickLabel',[]);
ylabel('Depth (m)');
datetick('x',15);
xlabel(['Time on ',datestr(mtime0)]);
colorbar;
figure(3);
pcolor(mtime-mtime0,z,I3);
shading('interp');
set(gca,'ydir','reverse');
set(gca,'XTickLabel',[]);
ylabel('Depth (m)');
datetick('x',15);
xlabel(['Time on ',datestr(mtime0)]);
colorbar;
figure(4);
pcolor(mtime-mtime0,z,I4);
shading('interp');
set(gca,'ydir','reverse');
set(gca,'XTickLabel',[]);
ylabel('Depth (m)');
datetick('x',15);
xlabel(['Time on ',datestr(mtime0)]);
colorbar;
figure(5)
pcolor(mtime-mtime0,z,U);
shading('interp');
caxis([umin umax]);
set(gca,'ydir','reverse');
set(gca,'XTickLabel',[]);
ylabel('Depth (m)');
datetick('x',15);
xlabel(['Time on ',datestr(mtime0)]);
colorbar;
figure(6)
pcolor(mtime-mtime0,z,V);
shading('interp');
caxis([umin umax]);
set(gca,'ydir','reverse');
set(gca,'XTickLabel',[]);
ylabel('Depth (m)');
datetick('x',15);
xlabel(['Time on ',datestr(mtime0)]);
colorbar;
figure(7);
pcolor(mtime-mtime0,z,W);
shading('interp');
caxis([wmin wmax]);
set(gca,'ydir','reverse');
ylabel('Depth (m)');
datetick('x',15);
xlabel(['Time on ',datestr(mtime0)]);
colorbar;
This diff is collapsed.
This diff is collapsed.
% calibrate O2 sensor with O2 titrations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
% O2 titrations %%%%%%%%%%%%%%%%%%%%%%%
T.stn = {'ES10','ES10','ES10','ES8','ES8','ES8','ES12','ES12','ES7','ES5'};
T.dep = [2,8,13,44,28,2,3,11,88,93];
T.oxy = [273.527885,276.057051,276.251603,275.667949,272.263303,276.202965,273.090145,271.776924,273.187421,222.020448];
% O2 CTD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S.stn = {'ES10','ES10','ES10','ES8','ES8','ES8','ES12','ES12','ES7','ES5'};
S.dep = [2,8,13,44,28,2,3,11,88,93];
S.oxy1 = [253.874,254.182,253.174,249.667,249.925,253.445,255.923,255.148,247.822,242.779];
S.oxy2 = [253.445,253.949,253.394,249.735,249.882,254.124,255.765,255.095,247.737,242.944];
% plot
OT = T.oxy;
OS = (S.oxy1+S.oxy2)/2;
omin = min([OT,OS]);
omax = max([OT,OS]);
figure(1)
clf
plot(OT,OS,'.','markersize',15)
hold on
plot([omin omax],[omin,omax],'k--')
axis image
xlabel('Oxygen from titration (\mumoles/L)')
ylabel('Oxygen from Rosette (\mumoles/kg)')
print -f1 -dpng ../figs/calib_O2
% calibrate fluorometer with chla samples
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
% chla samples %%%%%%%%%%%%%%%%%%%%%%%
R.stn = {'ES8','ES8','ES8','ES10','ES10','ES10','ES12','ES12'};
R.dep = [44,28,2,13,8,2,11,3];
R.chl = [0.639,0.639,0.601,0.444,0.441,0.726,8.131,7.396];
% fluorometer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S.stn = {'ES8','ES8','ES8','ES10','ES10','ES10','ES12','ES12'};
S.dep = [44,28,2,13,8,2,11,3];
S.chl1 = [5.4721e-02,6.8326e-02,5.1110e-01,6.2398e-01,6.2421e-01,7.3365e-01,3.0981e+00,3.2165e+00];
S.chl2 = [5.6814e-02,5.2927e-02,4.9570e-01,6.1823e-01,6.1957e-01,7.4112e-01,3.0948e+00,3.4058e+00];
% regression
CR = R.chl;
CS = (S.chl1+S.chl2)/2;
P = polyfit(CR,CS,1);
cmin = min([CR,CS]);
cmax = max([CR,CS]);
Cfit = polyval(P,[cmin,cmax]);
% plot
figure(1)
clf
plot(CR,CS,'.','markersize',15)
hold on
plot([cmin cmax],[cmin,cmax],'k--')
plot([cmin cmax],Cfit,'r')
text((cmin+cmax)/2,mean(Cfit),['y = ',num2str(P(1),'%3.2f'),' x + ',num2str(P(2),'%3.2f')])
axis image
xlabel('Chla from samples (\mug/L)')
ylabel('Chla from sensor (\mug/L)')
print -f1 -dpng ../figs/calib_chla
function S = enx2mat(file,umax,errmax,trans_depth)
% S = enx2mat(file,umax,errmax,trans_depth)
%
% Load RDI raw ADCP data file (.ENX) into a
% Matlab structure S.
% Correct ship motion with bottom tracking.
% Remove velocities with errors larger than errmax (m/s)
% and speed larger than umax (m/s).
% If transducer depth was not set in VmDAS, specify it
% with trans_depth (m).
% Cdric Chavanne, updated 2016/11/6.
% load raw data (.ENX file) without ensemble averaging (NUMAV = 1):
S = rdradcp(file,1);
% compute central position of each ping:
lon = (S.nav_slongitude + S.nav_elongitude)/2;
lat = (S.nav_slatitude + S.nav_elatitude)/2;
% convert time to datenum format:
t = datenum(2000,0,0) + S.mtime;
[year,month,day] = datevec(t(1));
% bins depths:
if nargin<4
z = nanmedian(S.depth) + S.config.ranges;
else
z = trans_depth + S.config.ranges;
end
% check attitude data:
figure(1)
clf
subplot(211)
plot(t,S.heading)
ylabel('heading')
subplot(212)
plot(t,S.pitch,'r',t,S.roll,'b')
ylabel('pitch (r), roll (b)')
% 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
plot(t,S.bt_range)
hold on
plot(t,zbot,'k','linewidth',2)
% shallowest depth not contaminated by side lobes:
zgood = zbot * cosd(S.config.beam_angle);
% check bottom track velocities:
figure(3)
clf
plot(t,S.bt_vel)
legend('u','v','w','err')
% bottom track vel are in mm/s
% clean bottom track velocities:
J = find(abs(S.bt_vel(4,:)) > 1e2);
S.bt_vel(:,J) = nan;
figure(3)
clf
plot(t,S.bt_vel)
legend('u','v','w','err')
% check water velocities:
figure(4)
clf
subplot(411)
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)
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)
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)
pcolor(t-datenum(year,0,0),-z,S.error_vel)
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k')
datetick('x',15)
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)
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)
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)
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)
pcolor(t-datenum(year,0,0),-z,squeeze(S.corr(:,4,:)))
shading flat
colorbar
hold on
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)
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)
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)
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)
pcolor(t-datenum(year,0,0),-z,squeeze(S.intens(:,4,:)))
shading flat
colorbar
hold on
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)
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)
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)
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)
pcolor(t-datenum(year,0,0),-z,squeeze(S.status(:,4,:)))
shading flat
colorbar
hold on
plot(t-datenum(year,0,0),-zgood,'k','linewidth',2)
datetick('x',15)
xlabel(datestr(t(1),1))
% correct water vel for ship motion:
u = S.east_vel - repmat(S.bt_vel(1,:)/1000,length(z),1);
v = S.north_vel - repmat(S.bt_vel(2,:)/1000,length(z),1);
w = S.vert_vel - repmat(S.bt_vel(3,:)/1000,length(z),1);
err = S.error_vel;
% remove depth cells contaminated by side lobes:
for j = 1:length(zgood)
I = find(z > zgood(j));
u(I,j) = nan;
v(I,j) = nan;
w(I,j) = nan;
err(I,j) = nan;
end
figure(5)
clf
subplot(411)
pcolor(t-datenum(year,0,0),-z,u)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(zgood) 0])
subplot(412)
pcolor(t-datenum(year,0,0),-z,v)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(zgood) 0])
subplot(413)
pcolor(t-datenum(year,0,0),-z,w)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
hold off
shading flat
colorbar
datetick('x',15)
set(gca,'ylim',[-max(zgood) 0])
subplot(414)
pcolor(t-datenum(year,0,0),-z,err)
hold on
plot(t-datenum(year,0,0),-zgood,'k')
hold off
shading flat
colorbar
datetick('x',15)
xlabel(datestr(t(1),1))
set(gca,'ylim',[-max(zgood) 0])
% remove velocities with large errors:
figure(8)
clf
hist(abs(err(:)),128)
I = find(abs(err) > errmax);
u(I) = nan;
v(I) = nan;
w(I) = nan;
err(I) = nan;
% check histograms of velocities:
figure(9)
clf
subplot(311)
hist(u(:),128)
xlabel('u (m/s)')
subplot(312)
hist(v(:),128)
xlabel('v (m/s)')
subplot(313)
hist(w(:),128)
xlabel('w (m/s)')
% remove large velocities:
speed = sqrt(u.^2+v.^2+w.^2);
figure(10)
clf
hist(speed(:),128)
I = find(speed > umax);
u(I) = nan;
v(I) = nan;
w(I) = nan;
err(I) = nan;
figure(9)
clf
subplot(311)
hist(u(:),128)</