Skip to content
Snippets Groups Projects
template_get_adcp_data.m 13.3 KiB
Newer Older
habasque's avatar
habasque committed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% template_get_adcp_data.m
% -------------------------------
% Author : Jrmie HABASQUE - IRD
% Modification: Pierre ROUSSELOT - IRD
habasque's avatar
habasque committed
% -------------------------------
% INPUTS:
% - binary raw file with .000 extension
% OUTPUTS:
% - U and V fields interpolated on a regulard grid, filtered and subsampled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('C:\Users\proussel\Documents\outils\TOOLS\climato\nansuite')
addpath(genpath('C:\Users\proussel\Documents\outils\ADCP\ADCP_mooring_data_processing\tools'))
Pierre Rousselot's avatar
Pierre Rousselot committed

% First part --------------------------------------------------------------------------------------------------------------------
habasque's avatar
habasque committed
close all
clear all
 
Pierre Rousselot's avatar
Pierre Rousselot committed
% Path
addpath('.\moored_adcp_proc'); % ou par exemple ('C:\Users\IRD_US_IMAGO\TRAITEMENTS\ADCP_MOUILLAGE\01_DATA_PROCESSING\moored_adcp_proc');
addpath('.\backscatter'); % (Optionnel) / ou par exemple ('C:\Users\IRD_US_IMAGO\TRAITEMENTS\ADCP_MOUILLAGE\01_DATA_PROCESSING\backscatter');

habasque's avatar
habasque committed
% Location rawfile
fpath = '';
rawfile='.\FR27_000.000'; % binary file with .000 extension 
% Directory for outputs
fpath_output = '.\FR29_bis\';
Pierre Rousselot's avatar
Pierre Rousselot committed
% Cruise/mooring info
cruise.name  = 'PIRATA-FR29';
mooring.name = '0N10W'; % 0N10W par exemple
mooring.lat  = 00+00.15/60; %latitude en degrs dcimaux
mooring.lon  = -09-53.15/60; %longitude en degrs dcimaux
clock_drift  = 208/3600; % convert into hrs
Pierre Rousselot's avatar
Pierre Rousselot committed
% ADCP info
adcp.sn          = 15258;
adcp.type        = '150 khz Quartermaster'; % Type : Quartermaster, longranger
adcp.direction   = 'up';        % upward-looking 'up', downward-looking 'dn'
adcp.instr_depth = 300;       % nominal instrument depth% If ADCP was not set up to correct for magnetic deviation internally
habasque's avatar
habasque committed
% ("EA0" code in configuration file), use http://www.ngdc.noaa.gov/geomag-web/#declination
% Magnetic deviation: Mean of deviations at time of deployment and time of recovery 

Pierre Rousselot's avatar
Pierre Rousselot committed
% Magnetic deviation values
magnetic_deviation_ini = -9;
magnetic_deviation_end = -8.77;
rot                    = (magnetic_deviation_ini+magnetic_deviation_end)/2;  
Pierre Rousselot's avatar
Pierre Rousselot committed
% Read rawfile
fprintf('Read %s\n', rawfile);
raw            = read_os3(rawfile,'all');
Pierre Rousselot's avatar
Pierre Rousselot committed
% Correct clock drift
time0          = julian(raw.juliandate);
clockd         = linspace(0, clock_drift, length(time0));
raw.juliandate = raw.juliandate - clockd / 24;                 % P. Rousselot - change clock drift /24

% Magnetic devitation over time %P . Rousselot
mag_dev         = linspace(magnetic_deviation_ini, magnetic_deviation_end, length(time0)); 
Pierre Rousselot's avatar
Pierre Rousselot committed

figure;plot(raw.pressure);
detrend_sdata = detrend(raw.pressure);
trend         = raw.pressure - detrend_sdata;
Pierre Rousselot's avatar
Pierre Rousselot committed
hold on
plot(trend, 'r--')
hold off
title('Pressure sensor');ylabel('Pressure [dbar]');xlabel('Time index');grid on; 
Pierre Rousselot's avatar
Pierre Rousselot committed
saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Pressure_sensor'],'fig')

Pierre Rousselot's avatar
Pierre Rousselot committed
figure;plot(raw.temperature);
title('Temperature sensor');ylabel('Temperature [C]');xlabel('Time index');grid on; 
Pierre Rousselot's avatar
Pierre Rousselot committed
% Second part --------------------------------------------------------------------------------------------------------------------

% Determine first and last indiced when instrument was at depth (you can do this by plotting 'raw.pressure' for example           
first = 10; 
last  = 17487; 
Pierre Rousselot's avatar
Pierre Rousselot committed

% amplitude of the bins / Correction ADCP's depth
ea    = squeeze(mean(raw.amp(:,:,first:last),2));   
Pierre Rousselot's avatar
Pierre Rousselot committed
figure; colormap jet; pcolor(ea); shading flat; title('Amplitude of the bins'); colorbar;
Pierre Rousselot's avatar
Pierre Rousselot committed
ylabel('Bins');xlabel('Time index'); 
saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Amplitude_bins'],'fig')

% Third part --------------------------------------------------------------------------------------------------------------------
habasque's avatar
habasque committed

% If upward looking: range of surface bins used for instrument depth correction below!
sbins      = 32:39;%30:35; % here a range of bins is given which cover the surface reflection
habasque's avatar
habasque committed
% Exclude data with percent good below prct_good
habasque's avatar
habasque committed

%% Read data
freq       = raw.config.sysconfig.frequency;
u2         = squeeze(raw.vel(:,1,first:last));
v2         = squeeze(raw.vel(:,2,first:last));
w          = squeeze(raw.vel(:,3,first:last));
vel_err    = squeeze(raw.vel(:,4,first:last));  % the difference in vertical velocity between the two pairs of transducers 
time       = raw.juliandate(first:last);
ang        = [raw.pitch(first:last) raw.roll(first:last) raw.heading(first:last)]; 
habasque's avatar
habasque committed
soundspeed = raw.soundspeed(first:last);
T          = raw.temperature(first:last);
press      = raw.pressure(first:last);
mag_dev    = mag_dev(first:last);
habasque's avatar
habasque committed

nbin       = raw.config.ncells;  % number of bins
bin        = 1:nbin;
blen       = raw.config.cell;    % bin length
blnk       = raw.config.blank;   % blank distance after transmit
habasque's avatar
habasque committed

dt         = (time(2)-time(1))*24;   % Sampling interval in hours
habasque's avatar
habasque committed

Pierre Rousselot's avatar
Pierre Rousselot committed
for ii = 1 : length(mag_dev)
    [u(:,ii),v(:,ii)] = uvrot(u2(:,ii), v2(:,ii), -mag_dev(ii));   % Correction of magnetic deviation (over time-P. Rousselot)
Pierre Rousselot's avatar
Pierre Rousselot committed
end
habasque's avatar
habasque committed

pg            = squeeze(raw.pg(:,4,first:last));  % percent good

igap          = find(pg<prct_good);           % Exclude data with percent good below prct_good
u(igap)       = NaN;
v(igap)       = NaN;
w(igap)       = NaN;
vel_err(igap) = NaN;


%% Control attitude
figure
subplot(3,1,1)
plot(ang(:,1))
ylabel('Pitch []')
subplot(3,1,2)
plot(ang(:,2))
ylabel('Roll []')
subplot(3,1,3)
plot(ang(:,3))
ylabel('Heading []')
axis([0 18000 0 360])
xlabel('Time Index')
habasque's avatar
habasque committed

%% Calculate depth of each bin:
dpt    = sw_dpth(press,mooring.lat)';  % convert pressure to depth, press needs to have dimension (n x 1)
dpt1   = repmat(dpt,nbin,1);
habasque's avatar
habasque committed
binmat = repmat((1:nbin)',1,length(dpt1));

% If ADCP is upward-looking a depth correction can be inferred from the
% surface reflection, which is done in adcp_surface_fit
if strcmp(adcp.direction,'up')  
    [z,dpt1,offset,xnull] = adcp_surface_fit(dpt,ea,sbins,blen,blnk,nbin);
Pierre Rousselot's avatar
Pierre Rousselot committed
elseif strcmp(adcp.direction,'dn')
    error('Bin depth calculation: unknown direction!');
habasque's avatar
habasque committed
end

Pierre Rousselot's avatar
Pierre Rousselot committed
saveas(figure(1),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Hist_diff_orig-depth_recon-depth'],'fig')
saveas(figure(2),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Offset_depth'],'fig')
saveas(figure(3),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Amplitude_bins_2'],'fig')

habasque's avatar
habasque committed
%% Remove bad data if ADCP is looking upward
u1=u; v1=v; w1=w; vel_err1=vel_err; ea1=ea;

if strcmp(adcp.direction,'up')
    for i = 1:length(time)
        sz_dpt(i) = adcp_shadowzone(dpt(i),raw.config.sysconfig.angle); % depending on the instrument depth and the beam angle the shadow zone, i.e. the depth below the surface which is contaminated by the surface reflection is determined
        iz(i)     = find(z(:,i)>sz_dpt(i),1,'last');
        sbin(i)   = bin(iz(i));
habasque's avatar
habasque committed
        
Pierre Rousselot's avatar
Pierre Rousselot committed
        %sbin(i)=30; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
habasque's avatar
habasque committed
        % here a manual criterion should be hard-coded if
        % adcp_check_surface (below) shows bad velocities close to the
        % surface

habasque's avatar
habasque committed
    end

    for i = 1:length(time) 
        u1(sbin(i)+1:end,i)       = NaN;
        v1(sbin(i)+1:end,i)       = NaN;
        w1(sbin(i)+1:end,i)       = NaN;
        vel_err1(sbin(i)+1:end,i) = NaN;
        ea1(sbin(i)+1:end,i)      = NaN;
habasque's avatar
habasque committed
    end

    if 1
        bins = nmedian(sbin)-4:nmedian(sbin)+2;
habasque's avatar
habasque committed
        adcp_check_surface(bins,u,u1,v,v1,time,bin,z); 
        % here the closest bins below the surface are plotted that are supposed to have good velocities, if there are still bad velocities a manual criterion needs to be found
    end
end
Pierre Rousselot's avatar
Pierre Rousselot committed
saveas(figure(4),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Meridional_zonal_velocity'],'fig')
habasque's avatar
habasque committed

%% SAVE DATA
% More meta information
%adcp.comment='';
adcp.config   = raw.config;
adcp.z_offset = offset;
adcp.ang      = ang;
adcp.mag_dev  = rot;
habasque's avatar
habasque committed

% Data structure
data.u      = u1;
data.v      = v1;
data.w      = w1;
data.e      = vel_err1;
data.ea     = ea1;
data.pg     = pg;
data.time   = time;
data.z_bins = z;
data.depth  = dpt;
data.temp   = T;
data.sspd   = soundspeed;
data.lat    = mooring.lat;
data.lon    = mooring.lon;
habasque's avatar
habasque committed

% Save
save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '.mat'],'adcp','mooring','data','raw','-v7.3');
habasque's avatar
habasque committed

%% Interpolate data on a regular vertical grid
Z      = fliplr(blen/2:blen:max(z(:))+blen);
Zmax   = max(Z);
data.Z = Z;
habasque's avatar
habasque committed

u_interp = NaN(length(time),length(Z));
v_interp = NaN(length(time),length(Z));
for i=1:length(time)
    % indice correspondant sur la grille finale Z
    ind = round((Zmax-z(1,i))/blen)+1;
    % filling the grid
    npts = min([length(Z)+ind+1 length(bin)]);
habasque's avatar
habasque committed
    u_interp(i,ind:ind+npts-1) = u1(1:npts,i);
    v_interp(i,ind:ind+npts-1) = v1(1:npts,i);
end
  
%% Horizontal interpolation, filtering and subsampling
[uintfilt,vintfilt,uifilt,vifilt,inttim,utid_baro,vtid_baro] = adcp_filt_sub(data,u_interp',v_interp',1:length(Z),40);
Pierre Rousselot's avatar
Pierre Rousselot committed
saveas(figure(5),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','data_raw_filt_subsampled_1'],'fig')
saveas(figure(6),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','data_raw_filt_subsampled_2'],'fig')
habasque's avatar
habasque committed
saveas(figure(5),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','data_raw_filt_subsampled_1'],'fig')
saveas(figure(6),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','data_raw_filt_subsampled_2'],'fig')
habasque's avatar
habasque committed

save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '_int_filt.mat'],'uifilt','vifilt','data');

habasque's avatar
habasque committed
% Save interpolated data
bin_start     = 4; % bin indice where good interpolated data for the whole dataset start
bin_end       = length(Z);
data.uintfilt = uintfilt(bin_start:bin_end,:);
data.vintfilt = vintfilt(bin_start:bin_end,:);
data.Z        = Z(bin_start:bin_end);
data.inttim   = inttim;save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '_int_filt_sub.mat'],'adcp','mooring','data','raw');
habasque's avatar
habasque committed

%% Figure
Pierre Rousselot's avatar
Pierre Rousselot committed
niv_u = (-1:0.05:1);
niv_v = (-1:0.05:1);
habasque's avatar
habasque committed

hf=figure('position', [0, 0, 1400, 1000]);
%u
subplot(2,1,1);
Pierre Rousselot's avatar
Pierre Rousselot committed
colormap jet
[C,h] = contourf(inttim,Z(bin_start:bin_end),uintfilt(bin_start:bin_end,:),niv_u); 
habasque's avatar
habasque committed
set(h,'LineColor','none');
caxis(niv_u([1 end]));
h=colorbar;
ylabel(h,'U [m s^-^1]');
set(gca,'ydir', 'reverse');
ylabel('Depth (m)');
ylim([0,adcp.instr_depth]);
%change figure label in HH:MM
Pierre Rousselot's avatar
Pierre Rousselot committed
title({[mooring.name, ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz']});
habasque's avatar
habasque committed

%v
subplot(2,1,2);
[C,h] = contourf(inttim,Z(bin_start:bin_end),vintfilt(bin_start:bin_end,:),niv_v); set(h,'LineColor','none');
habasque's avatar
habasque committed
caxis(niv_v([1 end]));
habasque's avatar
habasque committed
ylabel(h,'V [m s^-^1]');
set(gca,'ydir', 'reverse');
ylabel('Depth (m)');
ylim([0,adcp.instr_depth]);
%change figure label in HH:MM
gregtick;
Pierre Rousselot's avatar
Pierre Rousselot committed
title({[mooring.name, ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz']});
habasque's avatar
habasque committed

graph_name = [fpath_output, mooring.name '_U_V_int_filt_sub'];
habasque's avatar
habasque committed
set(hf,'Units','Inches');
pos        = get(hf,'Position');
set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
habasque's avatar
habasque committed
print(hf,graph_name,'-dpdf','-r300');

% %% Plot tide
% hf=figure('position', [0, 0, 1400, 1000]);
% niv_tide = (-5:0.2:5);
% utid_baro = utid_baro * 100;
% vtid_baro = vtid_baro * 100;
% %u
% subplot(2,1,1);
% colormap jet
% [C,h] = contourf(inttim,Z(bin_start:bin_end),utid_baro(bin_start:bin_end,:),niv_tide); 
% set(h,'LineColor','none');
% caxis(niv_tide([1 end]));
% h=colorbar;
% ylabel(h,'U [cm s^-^1]');
% set(gca,'ydir', 'reverse');
% ylabel('Depth (m)');
% ylim([0,adcp.instr_depth]);
% %change figure label in HH:MM
% gregtick;
% title({[mooring.name, ' - ZONAL TIDE VELOCITY - RDI ',num2str(freq),' kHz']});
% 
% %v
% subplot(2,1,2);
% [C,h] = contourf(inttim,Z(bin_start:bin_end),vtid_baro(bin_start:bin_end,:),niv_tide); 
% set(h,'LineColor','none');
% caxis(niv_tide([1 end]));
% h     = colorbar;
% ylabel(h,'V [cm s^-^1]');
% set(gca,'ydir', 'reverse');
% ylabel('Depth (m)');
% ylim([0,adcp.instr_depth]);
% %change figure label in HH:MM
% gregtick;
% title({[mooring.name, ' - MERIDIONAL TIDE VELOCITY - RDI ',num2str(freq),' kHz']});

%%  Write netcdf file           
[yr_start , ~, ~] = gregorian(inttim(1));
[yr_end,  ~, ~]   = gregorian(inttim(length(inttim)));
ncid     = netcdf.create([fpath_output,'ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'_1d.nc'],'NC_WRITE');
 
%create dimension
dimidt   = netcdf.defDim(ncid,'time',length(inttim));
dimidz   = netcdf.defDim(ncid,'depth',length(Z));
%Define IDs for the dimension variables (pressure,time,latitude,...)
time_ID  = netcdf.defVar(ncid,'time','double',dimidt);
depth_ID = netcdf.defVar(ncid,'depth','double',dimidz);
%Define the main variable ()
u_ID     = netcdf.defVar(ncid,'u','double',[dimidt dimidz]);
v_ID     = netcdf.defVar(ncid,'v','double',[dimidt dimidz]);
%We are done defining the NetCdf
netcdf.endDef(ncid);
%Then store the dimension variables in
netcdf.putVar(ncid,time_ID,inttim);
netcdf.putVar(ncid,depth_ID,Z);  
%Then store my main variable
netcdf.putVar(ncid,u_ID,uintfilt');
netcdf.putVar(ncid,v_ID,vintfilt');
%We're done, close the netcdf
netcdf.close(ncid);

Pierre Rousselot's avatar
Pierre Rousselot committed
% rmpath('..\moored_adcp_proc');
clear all; close all;

% -------------------------------------------------------------------------------------------