Newer
Older
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% template_get_adcp_data.m
% -------------------------------
% Author : Jrmie HABASQUE - IRD
% Modification: Pierre ROUSSELOT - IRD
% -------------------------------
% 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'))
% First part --------------------------------------------------------------------------------------------------------------------
% 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');
pierre.rousselot_ird.fr
committed
rawfile='.\FR27_000.000'; % binary file with .000 extension
fpath_output = '.\FR29_bis\';
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
adcp.sn = 15258;
adcp.type = '150 khz Quartermaster'; % Type : Quartermaster, longranger
adcp.direction = 'up'; % upward-looking 'up', downward-looking 'dn'
pierre.rousselot_ird.fr
committed
adcp.instr_depth = 300; % nominal instrument depth% If ADCP was not set up to correct for magnetic deviation internally
% ("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
magnetic_deviation_ini = -9;
magnetic_deviation_end = -8.77;
rot = (magnetic_deviation_ini+magnetic_deviation_end)/2;
raw = read_os3(rawfile,'all');
pierre.rousselot_ird.fr
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));
figure;plot(raw.pressure);
detrend_sdata = detrend(raw.pressure);
trend = raw.pressure - detrend_sdata;
title('Pressure sensor');ylabel('Pressure [dbar]');xlabel('Time index');grid on;
saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Pressure_sensor'],'fig')
figure;plot(raw.temperature);
title('Temperature sensor');ylabel('Temperature [C]');xlabel('Time index');grid on;
% 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;
ea = squeeze(mean(raw.amp(:,:,first:last),2));
figure; colormap jet; pcolor(ea); shading flat; title('Amplitude of the bins'); colorbar;
ylabel('Bins');xlabel('Time index');
saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Amplitude_bins'],'fig')
% Third part --------------------------------------------------------------------------------------------------------------------
% 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
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)];
T = raw.temperature(first:last);
press = raw.pressure(first:last);
mag_dev = mag_dev(first:last);
nbin = raw.config.ncells; % number of bins
bin = 1:nbin;
blen = raw.config.cell; % bin length
blnk = raw.config.blank; % blank distance after transmit
dt = (time(2)-time(1))*24; % Sampling interval in hours
[u(:,ii),v(:,ii)] = uvrot(u2(:,ii), v2(:,ii), -mag_dev(ii)); % Correction of magnetic deviation (over time-P. Rousselot)
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')
dpt = sw_dpth(press,mooring.lat)'; % convert pressure to depth, press needs to have dimension (n x 1)
dpt1 = repmat(dpt,nbin,1);
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_ird.fr
committed
z = dpt1+(binmat-0.5)*blen+blnk;else
error('Bin depth calculation: unknown direction!');
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')
%% 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));
%sbin(i)=30; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% here a manual criterion should be hard-coded if
% adcp_check_surface (below) shows bad velocities close to the
% surface
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;
bins = nmedian(sbin)-4:nmedian(sbin)+2;
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
saveas(figure(4),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Meridional_zonal_velocity'],'fig')
adcp.config = raw.config;
adcp.z_offset = offset;
adcp.ang = ang;
adcp.mag_dev = rot;
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;
save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '.mat'],'adcp','mooring','data','raw','-v7.3');
Z = fliplr(blen/2:blen:max(z(:))+blen);
Zmax = max(Z);
data.Z = Z;
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)]);
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);
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')
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')
save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '_int_filt.mat'],'uifilt','vifilt','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);
pierre.rousselot_ird.fr
committed
data.inttim = inttim;save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '_int_filt_sub.mat'],'adcp','mooring','data','raw');
hf=figure('position', [0, 0, 1400, 1000]);
%u
subplot(2,1,1);
colormap jet
[C,h] = contourf(inttim,Z(bin_start:bin_end),uintfilt(bin_start:bin_end,:),niv_u);
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
title({[mooring.name, ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz']});
pierre.rousselot_ird.fr
committed
[C,h] = contourf(inttim,Z(bin_start:bin_end),vintfilt(bin_start:bin_end,:),niv_v); set(h,'LineColor','none');
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;
title({[mooring.name, ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz']});
graph_name = [fpath_output, mooring.name '_U_V_int_filt_sub'];
pos = get(hf,'Position');
set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
% %% 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');
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);
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);
% rmpath('..\moored_adcp_proc');
clear all; close all;
% -------------------------------------------------------------------------------------------