diff --git a/.gitignore b/.gitignore index a0130cfd4953dca0a06d9ece5cee3621001aac4d..d4562b664f03a07cc518cfbcf02d3051150dc820 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,7 @@ *.docx FR*/ IDMX*/ +AMAZOMIX/ old_data/ doc/*.pdf *.jnl diff --git a/moored_adcp_proc/adcp_surface_fit.m b/moored_adcp_proc/adcp_surface_fit.m index 5b640fc179204ae89c336add365270d09adfc5c6..b045a238aabcfe237129c889c0ad8f43f7472c74 100644 --- a/moored_adcp_proc/adcp_surface_fit.m +++ b/moored_adcp_proc/adcp_surface_fit.m @@ -75,7 +75,7 @@ function [zbins,zadcp1,offset,x_null]=adcp_surface_fit(zadcp,ea,surface_bins,ble plot(-((x_null-1)*blen+blen/2+fbind),'r'); legend('Original','Reconstructed from surface reflection'); - if abs(offset)>15 + if any(abs(offset)>15) reply = input('Do you want to overwrite offset? 1/0 [0]:'); if isempty(reply) reply = 0; @@ -83,6 +83,7 @@ function [zbins,zadcp1,offset,x_null]=adcp_surface_fit(zadcp,ea,surface_bins,ble if reply==1 offset=input('->Enter new offset:'); + offset=offset*ones(1,length(zadcp)); else disp(['->Cleaned median filter offset is applied']); end diff --git a/multi_adcp_mooring/merge_data_finalfile.m b/multi_adcp_mooring/merge_data_finalfile.m index 4646dabd28165611f43d4db585b6dfb74e48c1a3..716ed761bb058eb31576687de901dd34e80b9d0f 100644 --- a/multi_adcp_mooring/merge_data_finalfile.m +++ b/multi_adcp_mooring/merge_data_finalfile.m @@ -5,11 +5,11 @@ clear all; close all; %% Variables -FileToMerge1 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/ADCP_0W0N_2016_2020_1d.nc'; %.nc file to merge with -FileToMerge2 = '/media/irdcampagnes/PIRATA/PIRATA-DATA/MOORING-PIRATA/0-0/ADCP_0W0N_2020_2021_1d.nc'; %.nc file to merge +FileToMerge1 = '/media/irdcampagnes/PIRATA-DATA/MOORING-PIRATA/10W/merged_data/ADCP_10W0N_2001_2019.nc'; %.nc file to merge with +FileToMerge2 = '/media/irdcampagnes/PIRATA-DATA/MOORING-PIRATA/10W/ADCP_10W0N_2021_2023_1d.nc'; %.nc file to merge step_subsampling = 1; % 1=daily plot_data = 1; -mooring.name = '0W0N'; +mooring.name = '10W0N'; mooring.lat = 0; mooring.lon = 0; freq = '150 khz Quartermaster'; @@ -35,14 +35,14 @@ ncfile2.v = ncread(FileToMerge2,'VCUR'); %% Create depth matrix %depth = max(vertcat(ncfile1.depth,ncfile2.depth)):ncfile1.depth(2)-ncfile1.depth(1):min(vertcat(ncfile1.depth,ncfile2.depth)); %depth = fliplr(depth); -max(vertcat(ncfile1.depth,ncfile2.depth)) -min(vertcat(ncfile1.depth,ncfile2.depth)) +max(vertcat(ncfile1.depth,ncfile2.depth)); +min(vertcat(ncfile1.depth,ncfile2.depth)); depth = 0:5:320; %% Create time matrix time = vertcat(ncfile1.time,ncfile2.time); [YY,MM,DD,hh,mm,ss] = datevec(time+datenum(1950,01,01)); -time = julian(YY,MM,DD,hh,mm,ss); +%time = julian(YY,MM,DD,hh,mm,ss); %% Create u/v matrix for ii = 1:length(ncfile1.time) @@ -57,10 +57,11 @@ end u = vertcat(u1,u2); v = vertcat(v1,v2); +time = julian(YY,MM,DD,hh,mm,ss); time = time*ones(1,length(depth)); % create a continuous series of daily data, ranging from min(d) to max(d) -final_time = ceil(min(min(time))):0.25:floor(max(max(time))); +final_time = ceil(min(min(time))):1:floor(max(max(time))); ADCP_final.u = NaN(length(final_time),length(depth)); ADCP_final.v = NaN(length(final_time),length(depth)); diff --git a/template_get_adcp_data.m b/template_get_adcp_data.m index 6e76b9e1e4f17714219fd3063df09d0b1e088b7c..225742c9e830d74423ba6081da3ae869eb4a11f0 100644 --- a/template_get_adcp_data.m +++ b/template_get_adcp_data.m @@ -16,26 +16,26 @@ addpath('/home/proussel/Documents/OUTILS/TOOLS/nansuite'); % NaNSuitePath %% META information: % Location rawfile -rawfile = 'FR31/FR29_000.000'; % binary file with .000 extension -fpath_output = './FR31/10Wtest/'; % Output directory +rawfile = 'AMAZOMIX/AMAZ_002.000'; % binary file with .000 extension +fpath_output = './AMAZOMIX/M2/'; % Output directory % Cruise/mooring info -cruise.name = 'PIRATA'; % cruise name -mooring.name = '10W0N'; % '0N10W' -mooring.lat = '00°00.017'; % latitude [°'] -mooring.lon = '-09°53.909'; % longitude [°'] -clock_drift = 418; % [seconds] +cruise.name = 'AMAZOMIX'; % cruise name +mooring.name = 'M2'; % '0N10W' +mooring.lat = '03°56.971'; % latitude [°'] +mooring.lon = '-45°07.518'; % longitude [°'] +clock_drift = 0; % [seconds] % ADCP info -adcp.sn = 24629; % ADCP serial number -adcp.type = '150 khz Quartermaster'; % Type : Quartermaster, longranger +adcp.sn = 'ADCP709'; % ADCP serial number +adcp.type = '75 khz Long Ranger'; % Type : Quartermaster, longranger adcp.direction = 'up'; % upward-looking 'up', downward-looking 'dn' -adcp.instr_depth = 300; % nominal instrument depth +adcp.instr_depth = 400; % nominal instrument depth instr = 1; % this is just for name convention and sorting of all mooring instruments offset_pres_sens = 1; % offset in m between pressure sensore and ADCP % NCFILE info -d_fillvalue = -9999; +d_fillvalue = -9999; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -60,419 +60,454 @@ disp('****') raw_file = [fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '_raw.mat']; if exist(raw_file) fprintf('Read %s\n', raw_file); + prepross = 1; load(raw_file) else fprintf('Read %s\n', rawfile); + prepross = 0; raw = read_os3(rawfile,'all'); save(raw_file,'raw','-v7.3'); end -%% Correct clock drift -time0 = julian(raw.juliandate); -clockd = linspace(0, clock_drift, length(time0)); -raw.juliandate = raw.juliandate - clockd / 24; -disp('****') -disp('Correct clock drift') - -%% Plot pressure and temperature sensor -figure; -subplot(2,1,1) -plot(raw.pressure); -detrend_sdata = detrend(raw.pressure); -trend = raw.pressure - detrend_sdata; -hold on -plot(trend, 'r--') -hold off -title('Pressure sensor'); -ylabel('Depth [m]'); -xlabel('Time index'); -grid on; - -subplot(2,1,2) -plot(raw.temperature); -title('Temperature sensor'); -ylabel('Temperature [cond1C]'); -xlabel('Time index'); -grid on; -saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Pressure_Temp_sensor'],'fig') - -%% Determine first and last indiced when instrument was at depth (you can do this by plotting 'raw.pressure' for example -disp('****') -first = input('Determine first indice when instrument was at depth (with pres/temp plot): '); -disp('****') -last = input('Determine last indice when instrument was at depth (with pres/temp plot): '); - -%% Correct velocity data with external T/S sensor - -%% Extract 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)); -corr = squeeze(mean(raw.cor(:,4,first:last),2)); -ea = squeeze(mean(raw.amp(:,:,first:last),2)); -pg = squeeze(raw.pg(:,4,first:last)); -time = raw.juliandate(first:last); -ang = [raw.pitch(first:last) raw.roll(first:last) raw.heading(first:last)]; -soundspeed = raw.soundspeed(first:last); -temp = raw.temperature(first:last); -press = raw.pressure(first:last); -if press < 0 - press = -press; -end -nbin = raw.config.ncells; % number of bins -bin = 1:nbin; % bin number -blen = raw.config.cell; % bin length -tlen = raw.config.pulse; % pulse length -lag = raw.config.lag; % transmit lag distance -blnk = raw.config.blank; % blank distance after transmit -fbind = raw.config.bin1distance;% middle of first bin distance -dt = (time(2)-time(1))*24; % Sampling interval in hours -EA0 = round(mean(ea(nbin,:))); - -%% Calculate Magnetic deviation values -[a,~] = gregorian(raw.juliandate(1)); -magnetic_deviation_ini = magdev(mooring.lat,mooring.lon,0,a+(raw.juliandate(1)-julian(a,1,1,0,0,0))/365.25); -[a,~] = gregorian(raw.juliandate(end)); -magnetic_deviation_end = magdev(mooring.lat,mooring.lon,0,a+(raw.juliandate(end)-julian(a,1,1,0,0,0))/365.25); -rot = (magnetic_deviation_ini+magnetic_deviation_end)/2; -mag_dev = linspace(magnetic_deviation_ini, magnetic_deviation_end, length(time0)); -mag_dev = mag_dev(first:last); - -%% Correction of magnetic deviation -disp('****') -disp('Correct magnetic deviation') -mag_dev = -mag_dev * pi/180; -for ii = 1 : length(mag_dev) - M(1,1) = cos(mag_dev(ii)); - M(1,2) = -sin(mag_dev(ii)); - M(2,1) = sin(mag_dev(ii)); - M(2,2) = cos(mag_dev(ii)); - vvel(1,:) = u2(:,ii); - vvel(2,:) = v2(:,ii); - vvvel = M * vvel; - u(:,ii) = vvvel(1,:); - v(:,ii) = vvvel(2,:); +if prepross == 1 + cor_data = input('Do you want to apply a new correction to data? 1/0 [1]'); + if isempty(cor_data) + cor_data = 0; + end +else + cor_data = 1; end +if cor_data == 1 - -%% Correct percent good: Exclude data with percent good below prct_good -figure; -colormap jet; -pcolor(pg); -shading flat; -title('Percent good of the bins'); colorbar; -ylabel('Bins'); -xlabel('Time index'); -saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Perceent_good'],'fig') -disp('****') -prct_good = input('Determine percent good threshold (generally 20): '); -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; -ea(igap) = NaN; -corr(igap) = NaN; - -%% Correction data with high attitude (pitch/roll) -figure; -subplot(2,1,1) -plot(abs(ang(:,1))); -hold on -plot(10*ones(1,length(ang(:,1))),'--r'); -hold off -title('Attitude sensor'); -ylabel('Pitch [cond1]'); -xlabel('Time index'); -axis([0 length(ang(:,1)) 0 20]) -grid on; - -subplot(2,1,2) -plot(abs(ang(:,2))); -hold on -plot(10*ones(1,length(ang(:,1))),'--r'); -hold off -ylabel('Roll [cond1]'); -xlabel('Time index'); -axis([0 length(ang(:,1)) 0 20]) -grid on; -saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Attitude'],'fig') -disp('****') -disp('Delete high attitude ADCP data (>=10°)'); - -high_pitch = find(abs(ang(:,1))>=10); -high_roll = find(abs(ang(:,2))>=10); -if ~isempty(high_pitch) || ~isempty(high_roll) - high_att = input('Do you want to delete them? 1/0 [1]'); - if isempty(high_att) - high_att = 1; + %% Correct clock drift + time0 = julian(raw.juliandate); + clockd = linspace(0, clock_drift, length(time0)); + raw.juliandate = raw.juliandate - clockd / 24; + disp('****') + disp('Correct clock drift') + + %% Plot pressure and temperature sensor + figure; + subplot(2,1,1) + plot(raw.pressure); + detrend_sdata = detrend(raw.pressure); + trend = raw.pressure - detrend_sdata; + hold on + plot(trend, 'r--') + hold off + title('Pressure sensor'); + ylabel('Depth [m]'); + xlabel('Time index'); + grid on; + + subplot(2,1,2) + plot(raw.temperature); + title('Temperature sensor'); + ylabel('Temperature [cond1C]'); + xlabel('Time index'); + grid on; + saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Pressure_Temp_sensor'],'fig') + + %% Determine first and last indiced when instrument was at depth (you can do this by plotting 'raw.pressure' for example + disp('****') + first = input('Determine first indice when instrument was at depth (with pres/temp plot): '); + disp('****') + last = input('Determine last indice when instrument was at depth (with pres/temp plot): '); + + %% Correct velocity data with external T/S sensor + + %% Extract 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)); + corr = squeeze(mean(raw.cor(:,4,first:last),2)); + ea = squeeze(mean(raw.amp(:,:,first:last),2)); + pg = squeeze(raw.pg(:,4,first:last)); + time = raw.juliandate(first:last); + ang = [raw.pitch(first:last) raw.roll(first:last) raw.heading(first:last)]; + soundspeed = raw.soundspeed(first:last); + temp = raw.temperature(first:last); + press = raw.pressure(first:last); + if press < 0 + press = -press; end - if high_att == 1 - high_att = union(high_pitch,high_roll); - offset(high_att) = NaN; - ang(high_att,1) = NaN; - ang(high_att,2) = NaN; - ang(high_att,3) = NaN; - mag_dev(high_att) = NaN; - u(high_att) = NaN; - v(high_att) = NaN; - w(high_att) = NaN; - vel_err(high_att) = NaN; - ea(high_att) = NaN; - corr(high_att) = NaN; - pg(high_att) = NaN; - time(high_att) = NaN; - z_bins(high_att) = NaN; - depth(high_att) = NaN; - temp(high_att) = NaN; - soundspeed(high_att) = NaN; + nbin = raw.config.ncells; % number of bins + bin = 1:nbin; % bin number + blen = raw.config.cell; % bin length + tlen = raw.config.pulse; % pulse length + lag = raw.config.lag; % transmit lag distance + blnk = raw.config.blank; % blank distance after transmit + fbind = raw.config.bin1distance;% middle of first bin distance + dt = (time(2)-time(1))*24; % Sampling interval in hours + EA0 = round(mean(ea(nbin,:))); + + %% Calculate Magnetic deviation values + [a,~] = gregorian(raw.juliandate(1)); + magnetic_deviation_ini = magdev(mooring.lat,mooring.lon,0,a+(raw.juliandate(1)-julian(a,1,1,0,0,0))/365.25); + [a,~] = gregorian(raw.juliandate(end)); + magnetic_deviation_end = magdev(mooring.lat,mooring.lon,0,a+(raw.juliandate(end)-julian(a,1,1,0,0,0))/365.25); + rot = (magnetic_deviation_ini+magnetic_deviation_end)/2; + mag_dev = linspace(magnetic_deviation_ini, magnetic_deviation_end, length(time0)); + mag_dev = mag_dev(first:last); + + %% Correction of magnetic deviation + disp('****') + disp('Correct magnetic deviation') + mag_dev = -mag_dev * pi/180; + for ii = 1 : length(mag_dev) + M(1,1) = cos(mag_dev(ii)); + M(1,2) = -sin(mag_dev(ii)); + M(2,1) = sin(mag_dev(ii)); + M(2,2) = cos(mag_dev(ii)); + vvel(1,:) = u2(:,ii); + vvel(2,:) = v2(:,ii); + vvvel = M * vvel; + u(:,ii) = vvvel(1,:); + v(:,ii) = vvvel(2,:); end -end -%% amplitude of the bins / Correction ADCP's depth -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') - -%% Add external pressure sensor choice -disp('****') -pressure_data = input('Do you have external pressure data? 1/0 [0]'); -if isempty(pressure_data) - pressure_data = 0; -end -if pressure_data == 1 - [filename, path, ~] = uigetfile('*.cnv', 'Select external pressure file'); - press_file = fullfile(path, filename); - press_file = read_sbe(press_file); - real_depth = press_file.depSM(first:last); - dpt = sw_dpth(press,mooring.lat)'; % convert pressure to depth, press needs to have dimension (n x 1) - plot(real_depth, 'b') - grid on + + + %% Correct percent good: Exclude data with percent good below prct_good + figure; + colormap jet; + pcolor(pg); + shading flat; + title('Percent good of the bins'); colorbar; + ylabel('Bins'); + xlabel('Time index'); + saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Perceent_good'],'fig') + disp('****') + prct_good = input('Determine percent good threshold (generally 20): '); + 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; + ea(igap) = NaN; + corr(igap) = NaN; + + %% Correction data with high attitude (pitch/roll) + figure; + subplot(2,1,1) + plot(abs(ang(:,1))); hold on - plot(dpt, 'r') - legend('External sensor depth','ADCP sensor depth'); - xlabel('Time Index') - ylabel('Depth [m]') - exsens_data = input('Do you want to use external sensor data? 1/0 [0]'); - if isempty(exsens_data) - exsens_data = 1; - end - if exsens_data - dpt = real_depth-offset_pres_sens; - dpt1 = repmat(dpt,nbin,1); - binmat = repmat((1:nbin)',1,length(dpt1)); - z(1,:) = dpt1(1,:)-fbind; - for ii = 2:length(binmat(:,1)) - z(ii,:) = z(1,:)-(binmat(ii,:)-1.5)*blen; + plot(10*ones(1,length(ang(:,1))),'--r'); + hold off + title('Attitude sensor'); + ylabel('Pitch [cond1]'); + xlabel('Time index'); + axis([0 length(ang(:,1)) 0 20]) + grid on; + + subplot(2,1,2) + plot(abs(ang(:,2))); + hold on + plot(10*ones(1,length(ang(:,1))),'--r'); + hold off + ylabel('Roll [cond1]'); + xlabel('Time index'); + axis([0 length(ang(:,1)) 0 20]) + grid on; + saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Attitude'],'fig') + disp('****') + disp('Delete high attitude ADCP data (>=10°)'); + + high_pitch = find(abs(ang(:,1))>=10); + high_roll = find(abs(ang(:,2))>=10); + if ~isempty(high_pitch) || ~isempty(high_roll) + % high_att = input('Do you want to delete them? 1/0 [1]'); + high_att = 0; + disp('WARNING you have high attitude ADCP data (>=10°) ---- check your data and be carefull'); + disp('Press a key to continue'); + pause + if isempty(high_att) + high_att = 1; + end + if high_att == 1 + high_att = union(high_pitch,high_roll); + offset(high_att) = NaN; + ang(high_att,1) = NaN; + ang(high_att,2) = NaN; + ang(high_att,3) = NaN; + mag_dev(high_att) = NaN; + u(high_att) = NaN; + v(high_att) = NaN; + w(high_att) = NaN; + vel_err(high_att) = NaN; + ea(high_att) = NaN; + corr(high_att) = NaN; + pg(high_att) = NaN; + time(high_att) = NaN; + z_bins(high_att) = NaN; + depth(high_att) = NaN; + temp(high_att) = NaN; + soundspeed(high_att) = NaN; end end -end -%% If upward looking: determine range of surface bins used for instrument depth correction below! -surface_bins = 0; -if pressure_data == 0 + %% amplitude of the bins / Correction ADCP's depth + 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') + + %% Add external pressure sensor choice disp('****') - surface_bins = input('Do you want to apply an offset using surface reflection? 1/0 [1]'); - if isempty(surface_bins) - surface_bins = 1; + pressure_data = input('Do you have external pressure data? 1/0 [0]'); + if isempty(pressure_data) + pressure_data = 0; end - if surface_bins == 1 - sbins = input('Determine range of surface bins used for instrument depth correction (with aplitude plot, ie. 30:35): '); - - %% 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); - 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 - disp('****') - if strcmp(adcp.direction,'up') - [z,dpt1,offset,xnull] = adcp_surface_fit(dpt,ea,sbins,blen,tlen,lag,blnk,nbin,fbind); - elseif strcmp(adcp.direction,'dn') - z(1,:) = dpt1(1,:)+(tlen+blen+lag)*0.5+blnk; - for ii = 2:length(binmat(:,1)) - z(ii,:) = z(1,:)+(binmat(ii,:)-1.5)*blen; - end - else - error('Bin depth calculation: unknown direction!'); - end - - saveas(figure(5),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Histdiff_depth'],'fig') - saveas(figure(6),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Offset_depth'],'fig') - saveas(figure(7),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Amplitude_bins_2'],'fig') - else - dpt = sw_dpth(press,mooring.lat)'; % convert pressure to depth, press needs to have dimension (n x 1) - dpt1 = repmat(dpt,nbin,1); - z = dpt; - figure - plot(dpt) - hold on - plot(adcp.instr_depth*ones(length(dpt),1),'--r') - hold off + if pressure_data == 1 + [filename, path, ~] = uigetfile('*.cnv', 'Select external pressure file'); + press_file = fullfile(path, filename); + press_file = read_sbe(press_file); + real_depth = press_file.depSM(first:last); + dpt = sw_dpth(press,mooring.lat)'; % convert pressure to depth, press needs to have dimension (n x 1) + plot(real_depth, 'b') grid on - axis([0 length(dpt) min(min(dpt,adcp.instr_depth*ones(1,length(dpt))))-50 max(max(dpt,adcp.instr_depth*ones(1,length(dpt))))+50]) - xlabel('Ensembles') + hold on + plot(dpt, 'r') + legend('External sensor depth','ADCP sensor depth'); + xlabel('Time Index') ylabel('Depth [m]') - legend('Depth','Nominal Depth') - offset = input('Do you want to apply a manual offset? 1/0 [0]'); - if isempty(offset) - offset = 0; - elseif offset == 1 - offset = input('->Enter new offset:'); - dpt = dpt + offset; - dpt1 = dpt1 + offset; - z = z + offset; - hold on - plot(dpt,'g') - hold off - legend('Depth', 'Nominal Depth', 'Corrected depth') - elseif offset == 0 - offset = 0; + exsens_data = input('Do you want to use external sensor data? 1/0 [0]'); + if isempty(exsens_data) + exsens_data = 1; end - - binmat = repmat((1:nbin)',1,length(dpt1)); - - if strcmp(adcp.direction,'up') - z(1,:) = dpt1(1,:)-(tlen+blen+lag)*0.5-blnk; + if exsens_data + dpt = real_depth-offset_pres_sens; + dpt1 = repmat(dpt,nbin,1); + binmat = repmat((1:nbin)',1,length(dpt1)); + z(1,:) = dpt1(1,:)-fbind; for ii = 2:length(binmat(:,1)) z(ii,:) = z(1,:)-(binmat(ii,:)-1.5)*blen; end - elseif strcmp(adcp.direction,'dn') - z(1,:) = dpt1(1,:)+(tlen+blen+lag)*0.5+blnk; - for ii = 2:length(binmat(:,1)) - z(ii,:) = z(1,:)+(binmat(ii,:)-1.5)*blen; - end - else - error('Bin depth calculation: unknown direction!'); end - end -end -%% Remove bad data if ADCP is looking upward -u1=u; v1=v; w1=w; vel_err1=vel_err; ea1=ea; corr1=corr; z1=z; -if surface_bins == 1 - disp('****') - disp('Remove bad data near surface due to sidelobe'); - - if strcmp(adcp.direction,'up') - for i = 1:length(time) - sz_dpt(i) = adcp_shadowzone(dpt1(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 - 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; - corr1(sbin(i)+1:end,i) = NaN; - z1(sbin(i)+1:end,i) = NaN; + + %% If upward looking: determine range of surface bins used for instrument depth correction below! + surface_bins = 0; + if pressure_data == 0 + disp('****') + %surface_bins = input('Do you want to apply an offset using surface reflection? 1/0 [1]'); + surface_bins = 1; + if isempty(surface_bins) + surface_bins = 1; end - - if 1 - bins = nmedian(sbin)-4:nmedian(sbin)+4; - adcp_check_surface(bins,u,u1,v,v1,corr,corr1,time,bin,z,z1,sz_dpt); - % 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 - reply_sidelobe = input('Do you want to apply manual criterion? 1/0 [0]:'); - if isempty(reply_sidelobe) - reply_sidelobe = 0; + if surface_bins == 1 + sbins = input('Determine range of surface bins used for instrument depth correction (with aplitude plot, ie. 30:35): '); + + %% 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); + 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 + disp('****') + if strcmp(adcp.direction,'up') + [z,dpt1,offset,xnull] = adcp_surface_fit(dpt,ea,sbins,blen,tlen,lag,blnk,nbin,fbind); + elseif strcmp(adcp.direction,'dn') + z(1,:) = dpt1(1,:)+(tlen+blen+lag)*0.5+blnk; + for ii = 2:length(binmat(:,1)) + z(ii,:) = z(1,:)+(binmat(ii,:)-1.5)*blen; + end + else + error('Bin depth calculation: unknown direction!'); + end + + saveas(figure(5),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Histdiff_depth'],'fig') + saveas(figure(6),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Offset_depth'],'fig') + saveas(figure(7),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Amplitude_bins_2'],'fig') + else + dpt = sw_dpth(press,mooring.lat)'; % convert pressure to depth, press needs to have dimension (n x 1) + dpt1 = repmat(dpt,nbin,1); + z = dpt; + figure + plot(dpt) + hold on + plot(adcp.instr_depth*ones(length(dpt),1),'--r') + hold off + grid on + axis([0 length(dpt) min(min(dpt,adcp.instr_depth*ones(1,length(dpt))))-50 max(max(dpt,adcp.instr_depth*ones(1,length(dpt))))+50]) + xlabel('Ensembles') + ylabel('Depth [m]') + legend('Depth','Nominal Depth') + offset = input('Do you want to apply a manual offset? 1/0 [0]'); + if isempty(offset) + offset = 0; + elseif offset == 1 + offset = input('->Enter new offset:'); + dpt = dpt + offset; + dpt1 = dpt1 + offset; + z = z + offset; + hold on + plot(dpt,'g') + hold off + legend('Depth', 'Nominal Depth', 'Corrected depth') + elseif offset == 0 + offset = 0; end - if reply_sidelobe - bin_cutoff = input('->Enter new bin cutoff value:'); - u1=u; v1=v; w1=w; vel_err1=vel_err; ea1=ea; corr1=corr; z1=z; - u1(bin_cutoff+1:end,:) = NaN; - v1(bin_cutoff+1:end,:) = NaN; - w1(bin_cutoff+1:end,:) = NaN; - vel_err1(bin_cutoff+1:end,:) = NaN; - ea1(bin_cutoff+1:end,:) = NaN; - corr1(bin_cutoff+1:end,:) = NaN; - z1(bin_cutoff+1:end,:) = NaN; - adcp_check_surface(bins,u,u1,v,v1,corr,corr1,time,bin,z,z1,sz_dpt); + + binmat = repmat((1:nbin)',1,length(dpt1)); + + if strcmp(adcp.direction,'up') + z(1,:) = dpt1(1,:)-(tlen+blen+lag)*0.5-blnk; + for ii = 2:length(binmat(:,1)) + z(ii,:) = z(1,:)-(binmat(ii,:)-1.5)*blen; + end + elseif strcmp(adcp.direction,'dn') + z(1,:) = dpt1(1,:)+(tlen+blen+lag)*0.5+blnk; + for ii = 2:length(binmat(:,1)) + z(ii,:) = z(1,:)+(binmat(ii,:)-1.5)*blen; + end + else + error('Bin depth calculation: unknown direction!'); end + end end - saveas(figure(7),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Meridional_zonal_velocity'],'fig') -end -%% SAVE DATA -% More meta information -%adcp.comment=''; -adcp.config = raw.config; -if pressure_data == 0 - adcp.z_offset = offset; -end -adcp.ang = ang; -adcp.mag_dev = rot; - -% 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 = z1; -data.depth = dpt1; -data.temp = temp; -data.sspd = soundspeed; -data.lat = mooring.lat; -data.lon = mooring.lon; + %% Remove bad data if ADCP is looking upward + u1=u; v1=v; w1=w; vel_err1=vel_err; ea1=ea; corr1=corr; z1=z; + if surface_bins == 1 + disp('****') + disp('Remove bad data near surface due to sidelobe'); -disp('****') -reply_ts = input('Do you want to calculate Target Strength? (CTD profile needed) 1/0 [0]:'); -if isempty(reply_ts) - reply_ts = 0; -end -if reply_ts == 1 - % EA0: noisefloor from ADCP electronic noise; a scalar - disp('->Calculate Target Strength') - disp(['->EA0=' num2str(EA0) '(noisefloor from ADCP electronic noise) ; pH=8.1; (seawater pH)']) - file_ctd = input('->Path to CTD .nc file [''*.nc'']:'); - if exist(file_ctd, 'file') - ctd.depth = ncread(file_ctd, 'PRES'); - ctd.temp = ncread(file_ctd, 'TEMP'); - ctd.sal = ncread(file_ctd, 'PSAL'); - ctd.tempR = ones(size(z)); - ctd.salR = ones(size(z)); - % for each ADCP bin, find salinity and temperature values - for i_bin = 1:nbin - [~,ind_depth] = min(abs(z(i_bin,:)-ctd.depth(:))); - ctd.salR(i_bin,:) = ctd.sal(ind_depth); - ctd.tempR(i_bin,:) = ctd.temp(ind_depth); + if strcmp(adcp.direction,'up') + for i = 1:length(time) + sz_dpt(i) = adcp_shadowzone(dpt1(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 + 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; + corr1(sbin(i)+1:end,i) = NaN; + z1(sbin(i)+1:end,i) = NaN; + end + + if 1 + bins = nmedian(sbin)-4:nmedian(sbin)+4; + if ~any(high_att ~= 0) + adcp_check_surface(bins,u,u1,v,v1,corr,corr1,time,bin,z,z1,sz_dpt); + end + % 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 + reply_sidelobe = input('Do you want to apply manual criterion? 1/0 [0]:'); + if isempty(reply_sidelobe) + reply_sidelobe = 0; + end + if reply_sidelobe + bin_cutoff = input('->Enter new bin cutoff value:'); + u1=u; v1=v; w1=w; vel_err1=vel_err; ea1=ea; corr1=corr; z1=z; + u1(bin_cutoff+1:end,:) = NaN; + v1(bin_cutoff+1:end,:) = NaN; + w1(bin_cutoff+1:end,:) = NaN; + vel_err1(bin_cutoff+1:end,:) = NaN; + ea1(bin_cutoff+1:end,:) = NaN; + corr1(bin_cutoff+1:end,:) = NaN; + z1(bin_cutoff+1:end,:) = NaN; + adcp_check_surface(bins,u,u1,v,v1,corr,corr1,time,bin,z,z1,sz_dpt); + end + end end - % Compute TS - [sspd,~] = meshgrid(soundspeed,1:nbin); - data.ts = target_strength(data.ea,EA0,ctd.salR,ctd.tempR,sspd,temp',soundspeed',z,adcp.config.cell,adcp.config.blank,adcp.config.sysconfig.angle,freq); - figure; - colormap jet; - pcolor(data.ts.ts); - shading flat; - title('Target Strength'); colorbar; - ylabel('Bins'); - xlabel('Time index'); - saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','target_strength'],'fig') - else - disp('CTD file doesn''t exist') + saveas(figure(7),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Meridional_zonal_velocity'],'fig') + end + + %% SAVE DATA + % More meta information + %adcp.comment=''; + adcp.config = raw.config; + if pressure_data == 0 + adcp.z_offset = offset; + end + adcp.ang = ang; + adcp.mag_dev = rot; + + % 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 = z1; + data.depth = dpt1; + data.temp = temp; + data.sspd = soundspeed; + data.lat = mooring.lat; + data.lon = mooring.lon; + + disp('****') + reply_ts = input('Do you want to calculate Target Strength? (CTD profile needed) 1/0 [0]:'); + if isempty(reply_ts) reply_ts = 0; end -end + if reply_ts == 1 + % EA0: noisefloor from ADCP electronic noise; a scalar + disp('->Calculate Target Strength') + disp(['->EA0=' num2str(EA0) '(noisefloor from ADCP electronic noise) ; pH=8.1; (seawater pH)']) + file_ctd = input('->Path to CTD .nc file [''*.nc'']:'); + if exist(file_ctd, 'file') + ctd.depth = ncread(file_ctd, 'PRES'); + ctd.temp = ncread(file_ctd, 'TEMP'); + ctd.sal = ncread(file_ctd, 'PSAL'); + ctd.tempR = ones(size(z)); + ctd.salR = ones(size(z)); + % for each ADCP bin, find salinity and temperature values + for i_bin = 1:nbin + [~,ind_depth] = min(abs(z(i_bin,:)-ctd.depth(:))); + ctd.salR(i_bin,:) = ctd.sal(ind_depth); + ctd.tempR(i_bin,:) = ctd.temp(ind_depth); + end + % Compute TS + [sspd,~] = meshgrid(soundspeed,1:nbin); + data.ts = target_strength(data.ea,EA0,ctd.salR,ctd.tempR,sspd,temp',soundspeed',z,adcp.config.cell,adcp.config.blank,adcp.config.sysconfig.angle,freq); + figure; + colormap jet; + pcolor(data.ts.ts); + shading flat; + title('Target Strength'); colorbar; + ylabel('Bins'); + xlabel('Time index'); + saveas(gcf,[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','target_strength'],'fig') + else + disp('CTD file doesn''t exist') + reply_ts = 0; + end + end -%% Save raw data -disp('****') -disp('Saving raw data') -save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '.mat'],'adcp','mooring','data','-v7.3'); + %% Save raw data + disp('****') + disp('Saving raw data') + save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '.mat'],'adcp','mooring','data','z','time','-v7.3'); + +else + + load([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '_raw.mat']); + blen = raw.config.cell; + nbin = raw.config.ncells; % number of bins + bin = 1:nbin; % bin number + freq = raw.config.sysconfig.frequency; + reply_ts = 0; + clear raw + load([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '.mat'],'adcp','mooring','data','z','time'); + + +end %% Interpolate data on a regular vertical grid if strcmp(adcp.direction,'up') @@ -483,14 +518,14 @@ if strcmp(adcp.direction,'up') if reply_ts == 1 ts_interp = NaN(length(time),length(Z)); end - + 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); + u_interp(i,ind:ind+npts-1) = data.u(1:npts,i); + v_interp(i,ind:ind+npts-1) = data.v(1:npts,i); if reply_ts == 1 ts_interp(i,ind:ind+npts-1) = data.ts.ts(1:npts,i); end @@ -503,14 +538,14 @@ else if reply_ts == 1 ts_interp = NaN(length(time),length(Z)); end - + 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); + u_interp(i,ind:ind+npts-1) = data.u(1:npts,i); + v_interp(i,ind:ind+npts-1) = data.v(1:npts,i); if reply_ts == 1 ts_interp(i,ind:ind+npts-1) = data.ts.ts(1:npts,i); end @@ -520,6 +555,7 @@ end if reply_ts == 0 ts_interp = NaN; end + %% Horizontal interpolation, filtering and subsampling disp('****') horz_interp = input('Do you want to do an horizontal interpolation? 1/0 [1]:'); @@ -532,22 +568,26 @@ if horz_interp if isempty(filt) filt = 1; end - + sub = input('Do you want to subsampling data daily? 1/0 [1]:'); - if isempty(filt) - sub = 1; + if isempty(filt) + sub = 1; end if filt if sub + statut = 'tide_filtered_daily'; disp('Horizontal interpolation, Filtering & Subsampling') else + statut = 'tide_filtered'; disp('Horizontal interpolation & Filtering') end disp('Filtering using 40h Hamming filter (other disponible filter: FFT filter, tide model)') else if sub + statut = 'daily'; disp('Horizontal interpolation & Subsampling') else + statut = 'interp'; disp('Horizontal interpolation') end end @@ -558,6 +598,7 @@ end if horz_interp == 0 + statut = 'gridded'; uintfilt = u_interp'; vintfilt = v_interp'; inttim = data.time; @@ -565,14 +606,14 @@ else [uintfilt,vintfilt,tsintfilt,inttim,utid_baro,vtid_baro] = adcp_filt_sub(data,u_interp',v_interp',ts_interp',1:length(Z),reply_ts, filt, sub); end -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(1),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','data_raw_U',statut],'fig') +saveas(figure(2),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','data_raw_V',statut],'fig') %% Save interpolated data disp('****') bin_start = input('Determine first bin indice with good interpolated data: '); % bin indice where good interpolated data for the whole dataset start disp('****') -bin_end = input('Determine last bin indice with good interpolated data: '); % bin indice where good interpolated data for the whole dataset start +bin_end = input('Determine last bin indice with good interpolated data: '); % bin indice where good interpolated data for the whole dataset start data.uintfilt = uintfilt(bin_start:bin_end,:); data.vintfilt = vintfilt(bin_start:bin_end,:); if reply_ts == 1 @@ -580,7 +621,10 @@ if reply_ts == 1 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'); +save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '_' statut '.mat'],'adcp','mooring','data'); + +%% NetCDF with gridded data +extract_to_netcdf(fpath_output, mooring, data, statut, d_fillvalue) %% Figure niv_u = (-1:0.05:1); @@ -600,11 +644,8 @@ ylabel('Depth (m)'); ylim([0,round(max(max(z)))]); %change figure label in HH:MM gregtick; -if filt - title({[mooring.name, ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz (filtered from tide)']}); -else - title({[mooring.name, ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz']}); -end +title({[mooring.name, ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz (' statut ')']}, 'Interpreter', 'none'); + %v subplot(2,1,2); @@ -618,13 +659,10 @@ ylabel('Depth (m)'); ylim([0,round(max(max(z)))]); %change figure label in HH:MM gregtick; -if filt - title({[mooring.name, ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz (filtered from tide)']}); -else - title({[mooring.name, ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz']}); -end +title({[mooring.name, ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz (' statut ')']}, 'Interpreter', 'none'); + -graph_name = [fpath_output, mooring.name '_', num2str(instr), '_U_V_int_filt_sub']; +graph_name = [fpath_output, mooring.name '_', num2str(instr), '_U_V_', statut]; set(hf,'Units','Inches'); pos = get(hf,'Position'); set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); @@ -641,12 +679,9 @@ if reply_ts == 1 ylim([0,round(max(max(z,adcp.instr_depth)))]); %change figure label in HH:MM gregtick; - if filt - title({[mooring.name, ' - TARGET STRENGTH - RDI ',num2str(freq),' kHz (filtered from tide)']}); - else - title({[mooring.name, ' - TARGET STRENGTH - RDI ',num2str(freq),' kHz']}); - end - + title({[mooring.name, ' - TARGET STRENGTH - RDI ',num2str(freq),' kHz (' statut ')']}, 'Interpreter', 'none'); + + graph_name = [fpath_output, mooring.name, '_', num2str(instr), '_TS_int_filt_sub']; set(hf,'Units','Inches'); pos = get(hf,'Position'); @@ -688,41 +723,4 @@ end % 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 -disp('****') -disp('Creating .nc file') - -% Input parameters for NETCDF Global Attributes -tc_globAttFilename = fullfile('tools/input_GlobalAttrParameters.xls'); % JLL 2020/12 Il serait judicieux de remonter cette valeur en début du script template_get_adcp_data.m - -%% Prepare informations and variables required to create NETCDF file %% -[yr_start , ~, ~] = gregorian(data.inttim(1)); -[yr_end, ~, ~] = gregorian(data.inttim(length(data.inttim))); - -% Read inputs metadata required for NETCDF Global Attributes -[~,~,cell_ncAttributes] = xlsread(tc_globAttFilename); - -% Complete output path and filename -tc_ncFilenam_out = fullfile(fpath_output,['ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'_1d.nc']); - -% Assign a "4D-size" (TIME,DEPTH,LATITUDE,LONGITUDE) to the ADCP current variables : UINTFILT, VINTFILT -td_uADCP = ones(numel(data.inttim),numel(data.Z),numel(data.lat),numel(data.lon)) * d_fillvalue; -td_uADCP(:,:,1,1) = data.uintfilt'; -td_vADCP = ones(numel(data.inttim),numel(data.Z),numel(data.lat),numel(data.lon)) * d_fillvalue; -td_vADCP(:,:,1,1) = data.vintfilt'; - -% Flip for convention -data.Z = fliplr(data.Z); -td_uADCP = fliplr(td_uADCP); -td_vADCP = fliplr(td_vADCP); - -% Group general ADCP mooring informations and ADCP data to be written in NETCDF file format -struct_dataADCP = struct('mooringName', mooring.name, 'mooringLat', mooring.lat,... - 'mooringLon', mooring.lon, 'time', data.inttim, 'depth', data.Z,... - 'u', td_uADCP, 'v', td_vADCP); - -%% Write netcdf file %% -f_w_ADCP_ncOS(tc_ncFilenam_out,cell_ncAttributes,struct_dataADCP,d_fillvalue); -disp('****') \ No newline at end of file +% title({[mooring.name, ' - MERIDIONAL TIDE VELOCITY - RDI ',num2str(freq),' kHz']}); \ No newline at end of file diff --git a/tools/extract_to_netcdf.m b/tools/extract_to_netcdf.m new file mode 100644 index 0000000000000000000000000000000000000000..99f15ab0e53a9b13f34105a40ab1444f7172fa7f --- /dev/null +++ b/tools/extract_to_netcdf.m @@ -0,0 +1,40 @@ +function extract_to_netcdf(fpath_output, mooring, data, statut,d_fillvalue) + +%% Write netcdf file +disp('****') +disp('Creating .nc file') + +% Input parameters for NETCDF Global Attributes +tc_globAttFilename = fullfile('tools/input_GlobalAttrParameters.xls'); % JLL 2020/12 Il serait judicieux de remonter cette valeur en début du script template_get_adcp_data.m + +%% Prepare informations and variables required to create NETCDF file %% +[yr_start , ~, ~] = gregorian(data.inttim(1)); +[yr_end, ~, ~] = gregorian(data.inttim(length(data.inttim))); + +% Read inputs metadata required for NETCDF Global Attributes +[~,~,cell_ncAttributes] = xlsread(tc_globAttFilename); + +% Complete output path and filename +tc_ncFilenam_out = fullfile(fpath_output,['ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'_' statut '.nc']); + +% Assign a "4D-size" (TIME,DEPTH,LATITUDE,LONGITUDE) to the ADCP current variables : UINTFILT, VINTFILT +td_uADCP = ones(numel(data.inttim),numel(data.Z),numel(data.lat),numel(data.lon)) * d_fillvalue; +td_uADCP(:,:,1,1) = data.uintfilt'; +td_vADCP = ones(numel(data.inttim),numel(data.Z),numel(data.lat),numel(data.lon)) * d_fillvalue; +td_vADCP(:,:,1,1) = data.vintfilt'; + +% Flip for convention +data.Z = fliplr(data.Z); +td_uADCP = fliplr(td_uADCP); +td_vADCP = fliplr(td_vADCP); + +% Group general ADCP mooring informations and ADCP data to be written in NETCDF file format +struct_dataADCP = struct('mooringName', mooring.name, 'mooringLat', mooring.lat,... + 'mooringLon', mooring.lon, 'time', data.inttim, 'depth', data.Z,... + 'u', td_uADCP, 'v', td_vADCP); + +%% Write netcdf file %% +f_w_ADCP_ncOS(tc_ncFilenam_out,cell_ncAttributes,struct_dataADCP,d_fillvalue); +disp('****') + +end \ No newline at end of file