From 326968488ea51ae013743adf944770cd21d8f5df Mon Sep 17 00:00:00 2001 From: rousselo <pierre.rousselot@ird.fr> Date: Sun, 28 Feb 2021 10:43:53 +0100 Subject: [PATCH] New interpolation grid if downward adcp & add specific script --- multi_adcp_mooring/merge_data.m | 182 +++++ .../merge_data_finalfile.m | 309 ++++---- .../template_combine_adcp_up_and_down.m | 101 +-- .../template_combine_adcp_up_and_up.m | 164 +++++ specific/template_get_adcp_data_10w.m | 668 +++++++++++++++++ specific/template_get_adcp_data_10w_down.m | 682 ++++++++++++++++++ specific/template_get_adcp_data_10w_olddata.m | 253 +++++++ template_get_adcp_data.m | 72 +- 8 files changed, 2215 insertions(+), 216 deletions(-) create mode 100644 multi_adcp_mooring/merge_data.m rename merge_data.m => multi_adcp_mooring/merge_data_finalfile.m (60%) create mode 100644 multi_adcp_mooring/template_combine_adcp_up_and_up.m create mode 100644 specific/template_get_adcp_data_10w.m create mode 100644 specific/template_get_adcp_data_10w_down.m create mode 100644 specific/template_get_adcp_data_10w_olddata.m diff --git a/multi_adcp_mooring/merge_data.m b/multi_adcp_mooring/merge_data.m new file mode 100644 index 0000000..06d26e9 --- /dev/null +++ b/multi_adcp_mooring/merge_data.m @@ -0,0 +1,182 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Merge subsurface ADCP mooring data % +% Autor: P. Rousselot / Date: 03/2020 % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +clear all; close all; + +%% Variables +FileToMerge1 = 'F:\Encours\PIRATA_ADCP-MOORINGS\Convert_ADCP/ADCP_23W0N_2001_2015.nc'; %.nc file to merge with +FileToMerge2 = 'F:\Encours\PIRATA_ADCP-MOORINGS\23W\0_Final_Files/ADCP_0N23W_2015_2016_1d_up_down.nc'; %.nc file to merge +% FileToMerge3 = '/media/irdcampagnes/PIRATA/PIRATA_ADCP-MOORINGS/10W/0_Final_Files/ADCP_10W0N_2004_2005_1d.nc'; %.nc file to merge +% FileToMerge1 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX/ADCP_0N130E_2013_2013_1d.nc'; %.nc file to merge with +% FileToMerge2 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX/ADCP_0N130E_2013_2016_1d.nc'; %.nc file to merge +step_subsampling = 1; % 1=daily +plot_data = 1; +mooring.name = '23W0N'; +freq = 'Variable'; +% mooring.name = '0N130E'; +% freq = '75'; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Read first .nc file +ncfile1.time = ncread(FileToMerge1,'TIME'); +ncfile1.depth = ncread(FileToMerge1,'DEPTH'); +ncfile1.u = ncread(FileToMerge1,'U')'; +ncfile1.v = ncread(FileToMerge1,'V')'; + +%% Read second .nc file +ncfile2.time = ncread(FileToMerge2,'TIME'); +ncfile2.depth = ncread(FileToMerge2,'DEPTH'); +ncfile2.u = ncread(FileToMerge2,'UCUR'); +ncfile2.v = ncread(FileToMerge2,'VCUR'); + +% %% Read third .nc file +% ncfile3.time = ncread(FileToMerge3,'TIME'); +% ncfile3.depth = ncread(FileToMerge3,'DEPTH'); +% ncfile3.u = ncread(FileToMerge3,'UCUR'); +% ncfile3.v = ncread(FileToMerge3,'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)) +depth = 0:5:1000; + +%% 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 = time + datenum(1950,1,1,0,0,0) + datenum(2020,01,01); % -0.5 is necessary to remove 12h00 introduced by datetime function starting at noon instead of midnight as performed by julian function previously used +time = time*ones(1,length(depth)); + +%% Create u/v matrix +for ii = 1:length(ncfile1.time) + u1(ii,:) = interp1(ncfile1.depth,ncfile1.u(ii,:),depth); + v1(ii,:) = interp1(ncfile1.depth,ncfile1.v(ii,:),depth); +end +for ii = 1:length(ncfile2.time) + u2(ii,:) = interp1(ncfile2.depth,ncfile2.u(ii,:),depth); + v2(ii,:) = interp1(ncfile2.depth,ncfile2.v(ii,:),depth); +end +% for ii = 1:length(ncfile3.time) +% u3(ii,:) = interp1(ncfile3.depth,ncfile3.u(ii,:),depth); +% v3(ii,:) = interp1(ncfile3.depth,ncfile3.v(ii,:),depth); +% end +u = vertcat(u1,u2); +v = vertcat(v1,v2); + +u = u'; +v = v'; + +% create a continuous series of daily data, ranging from min(d) to max(d) +final_time = ceil(min(min(time))):1:floor(max(max(time))); +ADCP_final.u = NaN(length(depth),length(final_time)); +ADCP_final.v = NaN(length(depth),length(final_time)); + +for i_time = 1:length(final_time) + for j_time = 1:length(time) + if final_time(i_time) == time(j_time,1) + ADCP_final.u(:,i_time)=u(:,j_time); + ADCP_final.v(:,i_time)=v(:,j_time); + end + end +end + +time = final_time; +time = time'*ones(1,length(depth)); +time = time'; +u = ADCP_final.u; +v = ADCP_final.v; + +%% Plot data +if plot_data + niv_u = (-1.5:0.05:1); + niv_v = (-1.8:0.05:1.5); + depth = depth.*ones(length(time),1); + depth = depth'; + %date1 = julian(ncfile1.time(1)); + %date2 = julian(ncfile2.time(end)); + date1 = datestr(ncfile1.time(1)+datenum(1950,01,01),'yyyy'); + date2 = datestr(ncfile2.time(end)+datenum(1950,01,01),'yyyy'); + + hf=figure('position', [0, 0, 1400, 1000]); + %u + subplot(2,1,1); + colormap jet + [C,h] = contourf(time,depth,u,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,max(max(depth))]); + %datetick('x','mm-yyyy'); + gregtick; + title({[mooring.name, ' - ' num2str(date1) ' to ' num2str(date2) ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz (filtered from tide)']}); + + %v + subplot(2,1,2); + [C,h] = contourf(time,depth,v,niv_v); + set(h,'LineColor','none'); + caxis(niv_v([1 end])); + h = colorbar; + ylabel(h,'V [m s^-^1]'); + set(gca,'ydir', 'reverse'); + ylabel('Depth (m)'); + ylim([0,max(max(depth))]); + gregtick; + title({[mooring.name, ' - ' num2str(date1) ' to ' num2str(date2) ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz (filtered from tide)']}); + + graph_name = ['ADCP_' mooring.name '_' num2str(date1) '_' num2str(date2) '_U_V_daily']; + set(hf,'Units','Inches'); + pos = get(hf,'Position'); + set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]) + print(hf,graph_name,'-dpdf','-r300'); +end + +%% Write netcdf file +time = time(1,:)'; +timed = gregorian(time); +timed(:,5) = 0; +timed(:,6) = 0; +time = datenum(timed)-datenum(1950,01,01); + +disp('****') +disp('Creating .nc file') +yr_start = datestr(time(1)+datenum(1950,01,01), 'yyyy'); +yr_end = datestr(time(end)+datenum(1950,01,01), 'yyyy'); + +ncid = netcdf.create(['ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'.nc'],'NC_CLOBBER'); + +%create dimension +dimidz = netcdf.defDim(ncid, 'DEPTH', size(depth,2)); +dimidt = netcdf.defDim(ncid, 'TIME', netcdf.getConstant('NC_UNLIMITED')); +%Define IDs for the dimension variables (pressure,time,latitude,...) +time_ID = netcdf.defVar(ncid,'TIME','double',dimidt); +netcdf.putAtt(ncid, time_ID, 'long_name', 'Time'); +netcdf.putAtt(ncid, time_ID, 'axis', 'T'); +netcdf.putAtt(ncid, time_ID, 'units', 'days since 1950-01-01T00:00:00Z'); +netcdf.putAtt(ncid, time_ID, 'time_origin', '1-JAN-1950:00:00:00'); +depth_ID = netcdf.defVar(ncid,'DEPTH','double',dimidz); +netcdf.putAtt(ncid, depth_ID, 'long_name', 'Depth'); +netcdf.putAtt(ncid, depth_ID, 'axis', 'Z'); +netcdf.putAtt(ncid, depth_ID, 'units', 'meters'); +netcdf.putAtt(ncid, depth_ID, 'positive', 'down'); +netcdf.putAtt(ncid, depth_ID, 'point_spacing', 'even'); +%Define the main variable () +u_ID = netcdf.defVar(ncid,'U','double',[dimidz dimidt]); +netcdf.putAtt(ncid, u_ID, 'long_name', 'Zonal velocities'); +v_ID = netcdf.defVar(ncid,'V','double',[dimidz dimidt]); +netcdf.putAtt(ncid, v_ID, 'long_name', 'Meridional velocities'); +%We are done defining the NetCdf +netcdf.endDef(ncid); +%Then store the dimension variables in +netcdf.putVar(ncid, depth_ID, depth); +netcdf.putVar(ncid, time_ID, 0, length(time), time); +%Then store my main variable +netcdf.putVar(ncid, u_ID, u); +netcdf.putVar(ncid, v_ID, v); +%We're done, close the netcdf +netcdf.close(ncid); +disp('****') \ No newline at end of file diff --git a/merge_data.m b/multi_adcp_mooring/merge_data_finalfile.m similarity index 60% rename from merge_data.m rename to multi_adcp_mooring/merge_data_finalfile.m index f4e2a4c..e1693a9 100644 --- a/merge_data.m +++ b/multi_adcp_mooring/merge_data_finalfile.m @@ -1,149 +1,162 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Merge subsurface ADCP mooring data % -% Autor: P. Rousselot / Date: 03/2020 % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -clear all; close all; - -%% Variables -FileToMerge1 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX2/ADCP_0N130E_2010_2010_1d.nc'; %.nc file to merge with -FileToMerge2 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX/ADCP_0N130E_2013_2016_1d.nc'; %.nc file to merge -FileToMerge3 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX2/ADCP_0N130E_2010_2012_1d.nc'; %.nc file to merge -% FileToMerge1 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX/ADCP_0N130E_2013_2013_1d.nc'; %.nc file to merge with -% FileToMerge2 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX/ADCP_0N130E_2013_2016_1d.nc'; %.nc file to merge -step_subsampling = 1; % 1=daily -plot_data = 1; -mooring.name = '0N130E'; -freq = '75'; -% mooring.name = '0N130E'; -% freq = '75'; -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% Read first .nc file -ncfile1.time = ncread(FileToMerge1,'time'); -ncfile1.depth = ncread(FileToMerge1,'depth'); -ncfile1.u = ncread(FileToMerge1,'u'); -ncfile1.v = ncread(FileToMerge1,'v'); - -%% Read second .nc file -ncfile2.time = ncread(FileToMerge2,'time'); -ncfile2.depth = ncread(FileToMerge2,'depth'); -ncfile2.u = ncread(FileToMerge2,'u'); -ncfile2.v = ncread(FileToMerge2,'v'); - -%% Read third .nc file -ncfile3.time = ncread(FileToMerge3,'time'); -ncfile3.depth = ncread(FileToMerge3,'depth'); -ncfile3.u = ncread(FileToMerge3,'u'); -ncfile3.v = ncread(FileToMerge3,'v'); -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% 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); - -%% Create time matrix -time = vertcat(ncfile1.time,ncfile3.time,ncfile2.time); -time = time*ones(1,length(depth)); - -%% Create u/v matrix -for ii = 1:length(ncfile1.time) - u1(ii,:) = interp1(ncfile1.depth,ncfile1.u(ii,:),depth); - v1(ii,:) = interp1(ncfile1.depth,ncfile1.v(ii,:),depth); -end -for ii = 1:length(ncfile2.time) - u2(ii,:) = interp1(ncfile2.depth,ncfile2.u(ii,:),depth); - v2(ii,:) = interp1(ncfile2.depth,ncfile2.v(ii,:),depth); -end -for ii = 1:length(ncfile3.time) - u3(ii,:) = interp1(ncfile3.depth,ncfile3.u(ii,:),depth); - v3(ii,:) = interp1(ncfile3.depth,ncfile3.v(ii,:),depth); -end -u = vertcat(u1,u3,u2); -v = vertcat(v1,v3,v2); - -%% Plot data -if plot_data - niv_u = (-1.5:0.05:0.9); - niv_v = (-1.8:0.05:1.5); - depth = depth.*ones(length(time),1); - date1 = julian(ncfile1.time(1)); - date2 = julian(ncfile2.time(end)); - - hf=figure('position', [0, 0, 1400, 1000]); - %u - subplot(2,1,1); - colormap jet - [C,h] = contourf(time,depth,u,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,max(max(depth))]); - gregtick; - title({[mooring.name, ' - ' num2str(date1(1)) ' to ' num2str(date2(1)) ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz (filtered from tide)']}); - - %v - subplot(2,1,2); - [C,h] = contourf(time,depth,v,niv_v); - set(h,'LineColor','none'); - caxis(niv_v([1 end])); - h = colorbar; - ylabel(h,'V [m s^-^1]'); - set(gca,'ydir', 'reverse'); - ylabel('Depth (m)'); - ylim([0,max(max(depth))]); - gregtick; - title({[mooring.name, ' - ' num2str(date1(1)) ' to ' num2str(date2(1)) ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz (filtered from tide)']}); - - graph_name = ['ADCP_' mooring.name '_' num2str(date1(1)) '_' num2str(date2(1)) '_U_V_daily']; - set(hf,'Units','Inches'); - pos = get(hf,'Position'); - set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]) - print(hf,graph_name,'-dpdf','-r300'); -end - -%% Write netcdf file -time = time(:,1); -timed = gregorian(time); -timed(:,5) = 0; -timed(:,6) = 0; -time = datenum(timed)-datenum(1950,01,01); - -disp('****') -disp('Creating .nc file') -yr_start = datestr(time(1)+datenum(1950,01,01), 'yyyy'); -yr_end = datestr(time(end)+datenum(1950,01,01), 'yyyy'); - -ncid = netcdf.create(['ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'.nc'],'NC_CLOBBER'); - -%create dimension -dimidz = netcdf.defDim(ncid, 'DEPTH', size(depth,2)); -dimidt = netcdf.defDim(ncid, 'TIME', netcdf.getConstant('NC_UNLIMITED')); -%Define IDs for the dimension variables (pressure,time,latitude,...) -time_ID = netcdf.defVar(ncid,'TIME','double',dimidt); -netcdf.putAtt(ncid, time_ID, 'long_name', 'Time'); -netcdf.putAtt(ncid, time_ID, 'axis', 'T'); -netcdf.putAtt(ncid, time_ID, 'units', 'days since 1950-01-01T00:00:00Z'); -netcdf.putAtt(ncid, time_ID, 'time_origin', '1-JAN-1950:00:00:00'); -depth_ID = netcdf.defVar(ncid,'DEPTH','double',dimidz); -netcdf.putAtt(ncid, depth_ID, 'long_name', 'Depth'); -netcdf.putAtt(ncid, depth_ID, 'axis', 'Z'); -netcdf.putAtt(ncid, depth_ID, 'units', 'meters'); -netcdf.putAtt(ncid, depth_ID, 'positive', 'down'); -netcdf.putAtt(ncid, depth_ID, 'point_spacing', 'even'); -%Define the main variable () -u_ID = netcdf.defVar(ncid,'U','double',[dimidz dimidt]); -netcdf.putAtt(ncid, u_ID, 'long_name', 'Zonal velocities'); -v_ID = netcdf.defVar(ncid,'V','double',[dimidz dimidt]); -netcdf.putAtt(ncid, v_ID, 'long_name', 'Meridional velocities'); -%We are done defining the NetCdf -netcdf.endDef(ncid); -%Then store the dimension variables in -netcdf.putVar(ncid, depth_ID, depth(1,:)); -netcdf.putVar(ncid, time_ID, 0, length(time), time); -%Then store my main variable -netcdf.putVar(ncid, u_ID, u'); -netcdf.putVar(ncid, v_ID, v'); -%We're done, close the netcdf -netcdf.close(ncid); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Merge subsurface ADCP mooring data % +% Autor: P. Rousselot / Date: 03/2020 % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%clear all; close all; + +%% Variables +FileToMerge1 = 'F:\Encours\PIRATA_ADCP-MOORINGS\Convert_ADCP/ADCP_23W0N_2001_2015.nc'; %.nc file to merge with +FileToMerge2 = 'F:\Encours\PIRATA_ADCP-MOORINGS\23W\0_Final_Files/ADCP_0N23W_2015_2016_1d_up_down.nc'; %.nc file to merge +%FileToMerge2 = 'F:\Encours\PIRATA_ADCP-MOORINGS\23W\v2\2008-2009\ADCP_23W0N_2008_2009_1d.nc'; +% FileToMerge3 = '/media/irdcampagnes/PIRATA/PIRATA_ADCP-MOORINGS/10W/0_Final_Files/ADCP_10W0N_2004_2005_1d.nc'; %.nc file to merge +% FileToMerge1 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX/ADCP_0N130E_2013_2013_1d.nc'; %.nc file to merge with +% FileToMerge2 = '/home/proussel/Documents/OUTILS/ADCP/ADCP_mooring_data_processing/IDMX/ADCP_0N130E_2013_2016_1d.nc'; %.nc file to merge +step_subsampling = 1; % 1=daily +plot_data = 1; +mooring.name = '23W0N'; +freq = 'Variable'; +% mooring.name = '0N130E'; +% freq = '75'; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Read first .nc file +ncfile1.time = ncread(FileToMerge1,'TIME'); +ncfile1.depth = ncread(FileToMerge1,'DEPTH'); +ncfile1.u = ncread(FileToMerge1,'U')'; +ncfile1.v = ncread(FileToMerge1,'V')'; + +%% Read second .nc file +ncfile2.time = ncread(FileToMerge2,'TIME'); +ncfile2.depth = ncread(FileToMerge2,'DEPTH'); +ncfile2.u = ncread(FileToMerge2,'UCUR'); +ncfile2.v = ncread(FileToMerge2,'VCUR'); + +% %% Read third .nc file +% ncfile3.time = ncread(FileToMerge3,'TIME'); +% ncfile3.depth = ncread(FileToMerge3,'DEPTH'); +% ncfile3.u = ncread(FileToMerge3,'UCUR'); +% ncfile3.v = ncread(FileToMerge3,'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)) +depth = 0:5:1000; + +%% 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 = time + datenum(1950,1,1,0,0,0) + datenum(2020,01,01); % -0.5 is necessary to remove 12h00 introduced by datetime function starting at noon instead of midnight as performed by julian function previously used +%time = time*ones(1,length(depth)); + +% time = time + datetime(1950,1,1,0,0,0); % -0.5 is necessary to remove 12h00 introduced by datetime function starting at noon instead of midnight as performed by julian function previously used +% time = time*ones(1,length(depth)); + +%% Create u/v matrix +for ii = 1:length(ncfile1.time) + u1(ii,:) = interp1(ncfile1.depth,ncfile1.u(ii,:),depth); + v1(ii,:) = interp1(ncfile1.depth,ncfile1.v(ii,:),depth); +end +for ii = 1:length(ncfile2.time) + u2(ii,:) = interp1(ncfile2.depth,ncfile2.u(ii,:),depth); + v2(ii,:) = interp1(ncfile2.depth,ncfile2.v(ii,:),depth); +end +% for ii = 1:length(ncfile3.time) +% u3(ii,:) = interp1(ncfile3.depth,ncfile3.u(ii,:),depth); +% v3(ii,:) = interp1(ncfile3.depth,ncfile3.v(ii,:),depth); +% end +u = vertcat(u1,u2); +v = vertcat(v1,v2); + +time = time*ones(1,length(depth)); + + +%% Plot data +if plot_data + niv_u = (-1.5:0.05:0.9); + niv_v = (-1.8:0.05:1.5); + depth = depth.*ones(length(time),1); + date1 = datestr(ncfile1.time(1)+datenum(1950,01,01),'YYYY'); + date2 = datestr(ncfile2.time(end)+datenum(1950,01,01),'YYYY'); + + hf=figure('position', [0, 0, 1400, 1000]); + %u + subplot(2,1,1); + colormap jet + [C,h] = contourf(time,depth,u,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,max(max(depth))]); + gregtick; + title({[mooring.name, ' - ' num2str(date1) ' to ' num2str(date2) ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz (filtered from tide)']}); + + %v + subplot(2,1,2); + [C,h] = contourf(time,depth,v,niv_v); + set(h,'LineColor','none'); + caxis(niv_v([1 end])); + h = colorbar; + ylabel(h,'V [m s^-^1]'); + set(gca,'ydir', 'reverse'); + ylabel('Depth (m)'); + ylim([0,max(max(depth))]); + gregtick; + title({[mooring.name, ' - ' num2str(date1) ' to ' num2str(date2) ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz (filtered from tide)']}); + + graph_name = ['ADCP_' mooring.name '_' num2str(date1) '_' num2str(date2) '_U_V_daily']; + set(hf,'Units','Inches'); + pos = get(hf,'Position'); + set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]) + print(hf,graph_name,'-dpdf','-r300'); +end + +%% Write netcdf file +time = time(:,1); +timed = gregorian(time); +timed(:,5) = 0; +timed(:,6) = 0; +time = datenum(timed)-datenum(1950,01,01); + +disp('****') +disp('Creating .nc file') +yr_start = datestr(time(1)+datenum(1950,01,01), 'yyyy'); +yr_end = datestr(time(end)+datenum(1950,01,01), 'yyyy'); + +ncid = netcdf.create(['ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'.nc'],'NC_CLOBBER'); + +%create dimension +dimidz = netcdf.defDim(ncid, 'DEPTH', size(depth,2)); +dimidt = netcdf.defDim(ncid, 'TIME', netcdf.getConstant('NC_UNLIMITED')); +%Define IDs for the dimension variables (pressure,time,latitude,...) +time_ID = netcdf.defVar(ncid,'TIME','double',dimidt); +netcdf.putAtt(ncid, time_ID, 'long_name', 'Time'); +netcdf.putAtt(ncid, time_ID, 'axis', 'T'); +netcdf.putAtt(ncid, time_ID, 'units', 'days since 1950-01-01T00:00:00Z'); +netcdf.putAtt(ncid, time_ID, 'time_origin', '1-JAN-1950:00:00:00'); +depth_ID = netcdf.defVar(ncid,'DEPTH','double',dimidz); +netcdf.putAtt(ncid, depth_ID, 'long_name', 'Depth'); +netcdf.putAtt(ncid, depth_ID, 'axis', 'Z'); +netcdf.putAtt(ncid, depth_ID, 'units', 'meters'); +netcdf.putAtt(ncid, depth_ID, 'positive', 'down'); +netcdf.putAtt(ncid, depth_ID, 'point_spacing', 'even'); +%Define the main variable () +u_ID = netcdf.defVar(ncid,'U','double',[dimidz dimidt]); +netcdf.putAtt(ncid, u_ID, 'long_name', 'Zonal velocities'); +v_ID = netcdf.defVar(ncid,'V','double',[dimidz dimidt]); +netcdf.putAtt(ncid, v_ID, 'long_name', 'Meridional velocities'); +%We are done defining the NetCdf +netcdf.endDef(ncid); +%Then store the dimension variables in +netcdf.putVar(ncid, depth_ID, depth(1,:)); +netcdf.putVar(ncid, time_ID, 0, length(time), time); +%Then store my main variable +netcdf.putVar(ncid, u_ID, u'); +netcdf.putVar(ncid, v_ID, v'); +%We're done, close the netcdf +netcdf.close(ncid); disp('****') \ No newline at end of file diff --git a/multi_adcp_mooring/template_combine_adcp_up_and_down.m b/multi_adcp_mooring/template_combine_adcp_up_and_down.m index a6e2009..01247ea 100644 --- a/multi_adcp_mooring/template_combine_adcp_up_and_down.m +++ b/multi_adcp_mooring/template_combine_adcp_up_and_down.m @@ -16,26 +16,26 @@ clear all;close all; addpath('.\moored_adcp_proc'); % Location of ADCP up and down data -fpath = '.\data_example\up_and_down\'; +fpath = 'F:\Encours\PIRATA_ADCP-MOORINGS\23W\v2\2015-2016\'; % Directory for outputs -fpath_output = '.\data_example\up_and_down\'; +fpath_output = 'F:\Encours\PIRATA_ADCP-MOORINGS\23W\v2\2015-2016\'; % Contour levels for u anv v fields niv_u = (-1.5:0.1:1.5); niv_v = (-0.5:0.1:0.5); %% combine up and down -load([fpath, 'FR22-10W0N_508_instr_01_int_filt_sub']); +load([fpath, 'up\23W0N_1140_instr_01_int_filt_sub']); adcp_up=adcp; up_u=data.uintfilt; up_v=data.vintfilt; up_z=data.Z; up_time=data.inttim; -npts_up= data.npts_up; +%npts_up= data.npts_up; freq_up=adcp_up.config.sysconfig.frequency; -load([fpath, 'FR22-10W0N_509_instr_02_int_filt_sub_down']); +load([fpath, 'dn\23W0N_2627_instr_01_int_filt_sub']); adcp_down=adcp; down_u=data.uintfilt; down_v=data.vintfilt; @@ -54,12 +54,15 @@ distance_between_instruments = 3; distance_between_up_and_down = adcp_up.config.cell/2 + adcp_up.config.blank + distance_between_instruments + adcp_down.config.blank + adcp_down.config.cell/2; %calculate real minimum down ADCP depth up_z = fliplr(up_z); -depth_down_ADCP = up_z(min(npts_up)+1) + ceil(distance_between_up_and_down); +depth_down_ADCP = up_z(end) + ceil(distance_between_up_and_down); %interval grid between up and down ADCP interval_grid = max(up_z)+adcp_up.config.cell/2-1:adcp_up.config.cell/2:min(depth_down_ADCP)-adcp_up.config.cell/2+1; %on applique l'offset sur les profondeurs de l'ADCP down offset_ADCP_down = min(down_z) - min(depth_down_ADCP); down_z = down_z - offset_ADCP_down; +% down_z = flip(down_z); +% down_u = flip(down_u); +% down_v = flip(down_v); % initialisation de la matrice combinée up + down [B,a] = size(down_u); [B1,a1] = size(up_u); @@ -79,72 +82,84 @@ for i = 1:length(down_u) end %% figure finale -z_final = [up_z interval_grid down_z]; - +Z = [up_z interval_grid down_z]; +inttim = down_time; +bin_start = 1; +bin_end = length(Z); +uintfilt = u_final; +vintfilt = v_final; + +%% Figure +niv_u = (-1:0.05:1); +niv_v = (-1:0.05:1); +close all hf=figure('position', [0, 0, 1400, 1000]); %u subplot(2,1,1); -[~,h] = contourf(down_time,z_final,u_final,niv_u); +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)'); -gregtick -title({[mooring.name, ' - ZONAL VELOCITY - ADCP UP : RDI ',num2str(freq_up),' kHz - ADCP DOWN : RDI ',num2str(freq_down),' kHz']}); +ylim([0,round(max(Z))]); +%change figure label in HH:MM +gregtick; +title({[mooring.name, ' - ZONAL VELOCITY - RDI 150 kHz and 75 kHz (filtered from tide)']}); + %v subplot(2,1,2); -[C,h] = contourf(down_time,z_final,v_final,niv_v); +[C,h] = contourf(inttim,Z(bin_start:bin_end),vintfilt(bin_start:bin_end,:),niv_v); set(h,'LineColor','none'); caxis(niv_v([1 end])); -h=colorbar; +h = colorbar; ylabel(h,'V [m s^-^1]'); set(gca,'ydir', 'reverse'); ylabel('Depth (m)'); -gregtick -title({[mooring.name, ' - MERIDIONAL VELOCITY - ADCP UP : RDI ',num2str(freq_up),' kHz - ADCP DOWN : RDI ',num2str(freq_down),' kHz']}); - -% Save combined data -data.time = down_time; -data.u_final = u_final; -data.v_final = v_final; -data.z_final = z_final; -save([fpath, mooring.name '_UP_DOWN_int_filt_sub.mat'],'adcp','mooring','data','raw'); - -graph_name = [fpath_output, mooring.name '_U_V_UP_DOWN_int_filt_sub']; +ylim([0,round(max(Z))]); +%change figure label in HH:MM +gregtick; +title({[mooring.name, ' - ZONAL VELOCITY - RDI 150 kHz and 75 kHz (filtered from tide)']}); + + +graph_name = [fpath_output, mooring.name '_up_down_U_V_int_filt_sub']; set(hf,'Units','Inches'); -pos = get(hf,'Position'); -set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]) +pos = get(hf,'Position'); +set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); print(hf,graph_name,'-dpdf','-r300'); -% rmpath -rmpath('.\moored_adcp_proc'); +%% Write netcdf file +disp('****') +disp('Creating .nc file') +[yr_start , ~, ~] = gregorian(inttim(1)); +[yr_end, ~, ~] = gregorian(inttim(length(inttim))); -%% Write netcdf file -[yr_start , ~, ~] = gregorian(down_time(1)); -[yr_end, ~, ~] = gregorian(down_time(length(down_time))); +ncid = netcdf.create([fpath_output,'ADCP_',mooring.name,'_',num2str(yr_start),'_',num2str(yr_end),'_1d.nc'],'NC_WRITE'); -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(down_time)); -dimidz = netcdf.defDim(ncid,'depth',length(z_final)); +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); +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]); +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,down_time); -netcdf.putVar(ncid,depth_ID,z_final); +netcdf.putVar(ncid,time_ID,inttim); +netcdf.putVar(ncid,depth_ID,Z); %Then store my main variable -netcdf.putVar(ncid,u_ID,u_final); -netcdf.putVar(ncid,v_ID,v_final); +netcdf.putVar(ncid,u_ID,uintfilt'); +netcdf.putVar(ncid,v_ID,vintfilt'); + %We're done, close the netcdf netcdf.close(ncid); +disp('****') +% ------------------------------------------------------------------------------------------- diff --git a/multi_adcp_mooring/template_combine_adcp_up_and_up.m b/multi_adcp_mooring/template_combine_adcp_up_and_up.m new file mode 100644 index 0000000..29990c9 --- /dev/null +++ b/multi_adcp_mooring/template_combine_adcp_up_and_up.m @@ -0,0 +1,164 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% template_combine_adcp_up_and_down.m +% ------------------------------- +% Author : Jérémie HABASQUE - IRD +% ------------------------------- +% INPUTS: +% - ADCP up data processed (output of template_get_adcp_data.m) +% - ADCP down data processed (output of template_get_adcp_data.m) +% OUTPUTS: +% - U and V fields interpolated on a regulard grid, filtered and subsampled +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +clear all;close all; + +% path +addpath('.\moored_adcp_proc'); + +% Location of ADCP up and down data +fpath = '.\data_example\up_and_down\'; + +% Directory for outputs +fpath_output = './23W/'; + +% Contour levels for u anv v fields +niv_u = (-1.5:0.1:1.5); +niv_v = (-0.5:0.1:0.5); + +%% combine up and down +load('C:\Users\proussel\Documents\OUTILS\ADCP\ADCP_mooring_data_processing\23W\2008-2009\up\23W0N_1001_instr_01_int_filt_sub'); +adcp_up=adcp; +up_u=data.uintfilt; +up_v=data.vintfilt; +up_z=data.Z; +up_time=data.inttim; +%npts_up= data.npts_up; +freq_up=adcp_up.config.sysconfig.frequency; + +load('C:\Users\proussel\Documents\OUTILS\ADCP\ADCP_mooring_data_processing\23W\2008-2009\dn\23W0N_1023_instr_02_int_filt_sub'); +adcp_down=adcp; +down_u=data.uintfilt; +down_v=data.vintfilt; +down_z=data.Z; +down_time=data.inttim; +freq_down = adcp_down.config.sysconfig.frequency; + +% distance between deepest measurement of the upward looking and +% shallowest measurement of the downward looking ADCP +% half distance of 1st bin in upward ADCP + +% blank of upward ADCP + +% distance between instruments + +% blank of downward ADCP + +% half distance of 1st bin in downward ADCP. +distance_between_instruments = 3; +distance_between_up_and_down = adcp_up.config.cell/2 + adcp_up.config.blank + distance_between_instruments + adcp_down.config.blank + adcp_down.config.cell/2; +%calculate real minimum down ADCP depth +up_z = fliplr(up_z); +depth_down_ADCP = 148 + ceil(distance_between_up_and_down); +%interval grid between up and down ADCP +interval_grid = max(up_z)+adcp_up.config.cell/2-1:adcp_up.config.cell/2:min(depth_down_ADCP)-adcp_up.config.cell/2+1; +%on applique l'offset sur les profondeurs de l'ADCP down +offset_ADCP_down = min(down_z) - min(depth_down_ADCP); +down_z = down_z - offset_ADCP_down; +down_z = flip(down_z); +down_u = flip(down_u); +down_v = flip(down_v); + +% initialisation de la matrice combinée up + down +[B,a] = size(down_u); [B1,a1] = size(up_u); +BB = B+B1+1 ; +DOWN = NaN(BB,a); +V_DOWN = NaN(BB,a); +u_final = NaN(BB,a); +v_final = NaN(BB,a); + +% on alimente la matrice avec les données up +u_final(1:B1,1:a1) = up_u(B1:-1:1,:); +v_final(1:B1,1:a1) = up_v(B1:-1:1,:); + +for i = 1:length(down_u) + u_final(B1+length(interval_grid)+1:B1+length(interval_grid)+B,i) = down_u(:,i); + v_final(B1+length(interval_grid)+1:B1+length(interval_grid)+B,i) = down_v(:,i); +end + +%% figure finale +Z = [up_z interval_grid down_z]; +inttim = down_time; +bin_start = 1; +bin_end = 52; +uintfilt = u_final; +vintfilt = v_final; + +%% Figure +niv_u = (-1:0.05:1); +niv_v = (-1:0.05:1); +close all +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,round(max(Z))]); +%change figure label in HH:MM +gregtick; +title({[mooring.name, ' - ZONAL VELOCITY - RDI 150 kHz and 75 kHz (filtered from tide)']}); + + +%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'); +caxis(niv_v([1 end])); +h = colorbar; +ylabel(h,'V [m s^-^1]'); +set(gca,'ydir', 'reverse'); +ylabel('Depth (m)'); +ylim([0,round(max(Z))]); +%change figure label in HH:MM +gregtick; +title({[mooring.name, ' - ZONAL VELOCITY - RDI 150 kHz and 75 kHz (filtered from tide)']}); + + +graph_name = [fpath_output, mooring.name '_up_down_U_V_int_filt_sub']; +set(hf,'Units','Inches'); +pos = get(hf,'Position'); +set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); +print(hf,graph_name,'-dpdf','-r300'); + +%% Write netcdf file +disp('****') +disp('Creating .nc 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(data.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,data.Z); +%Then store my main variable +netcdf.putVar(ncid,u_ID,data.uintfilt'); +netcdf.putVar(ncid,v_ID,data.vintfilt'); + +%We're done, close the netcdf +netcdf.close(ncid); +disp('****') +% ------------------------------------------------------------------------------------------- diff --git a/specific/template_get_adcp_data_10w.m b/specific/template_get_adcp_data_10w.m new file mode 100644 index 0000000..ae5109c --- /dev/null +++ b/specific/template_get_adcp_data_10w.m @@ -0,0 +1,668 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% template_get_adcp_data.m +% ------------------------------- +% Author : Pierre ROUSSELOT - IRD (pierre.rousselot@ird.fr) +% Jeremie HABASQUE - IRD (jeremie.habasque@ird.fr) +% ------------------------------- +% INPUTS: +% - binary raw file with .000 extension +% OUTPUTS: +% - U and V fields interpolated on a regulard grid, filtered and subsampled +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +close all; clear +addpath(genpath('../ADCP_mooring_data_processing')); +addpath('/home/proussel/Documents/OUTILS/TOOLS/nansuite'); % NaNSuitePath + +%% META information: +% Location rawfile +rawfile = 'FR10W000.000'; % binary file with .000 extension +fpath_output = './10W/'; % Output directory + +% Cruise/mooring info +cruise.name = 'PIRATA'; % cruise name +mooring.name = '10W0N'; % '0N10W' +mooring.lat = '00°00.000'; % latitude [°'] +mooring.lon = '-10°00.000'; % longitude [°'] +clock_drift = 0; % [seconds] + +% ADCP info +adcp.sn = 508; % ADCP serial number +adcp.type = '300 khz Quartermaster'; % Type : Quartermaster, longranger +adcp.direction = 'up'; % upward-looking 'up', downward-looking 'dn' +adcp.instr_depth = 150; % nominal instrument depth +instr = 1; % this is just for name convention and sorting of all mooring instruments + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Convert variables +latDegInd = strfind(mooring.lat,'°'); +lonDegInd = strfind(mooring.lon,'°'); +mooring.lat = str2double(mooring.lat(1:latDegInd-1))+str2double(mooring.lat(latDegInd+1:end-1))/60; +mooring.lon = str2double(mooring.lon(1:lonDegInd-1))+str2double(mooring.lon(lonDegInd+1:end-1))/60; +clock_drift = clock_drift/3600; % convert into hrs + +%% Read rawfile +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); +% load(raw_file) +% else +% fprintf('Read %s\n', rawfile); +% raw = read_os3(rawfile,'all'); +% save(raw_file,'raw','-v7.3'); +% end +load('FR22-UP-508.mat') + +%% 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 = SerEmmpersec(first:last,:)'./1000; u2(u2<-30)=NaN; +v2 = SerNmmpersec(first:last,:)'./1000; v2(v2<-30)=NaN; +vel_err = SerErmmpersec(first:last,:)'./1000; vel_err(vel_err<-30)=NaN; +pg=SerPG4(first:last,:)'; + +% 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(SerC4cnt(first:last,:)'); + +test(:,1,:) = SerEA1cnt'; test(:,2,:) = SerEA2cnt'; test(:,3,:) = SerEA3cnt'; test(:,4,:) = SerEA4cnt'; + +ea = squeeze(mean(test(:,4,first:last),2)); +%pg = squeeze(raw.pg(:,4,first:last)); +%time = raw.juliandate(first:last); +time = julian(SerYear(first:last)+2000,SerMon(first:last),SerDay(first:last),SerHour(first:last)); +time=time(1)+datenum(0,0,0,1,0,0).*[0:length(depth)-1]'; +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+(time(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+(time(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(time)); +%mag_dev = mag_dev(first:last); + +%% Correction of magnetic deviation +disp('****') +disp('Correct magnetic deviation') +% for ii = 1 : length(mag_dev) +% [u(:,ii),v(:,ii)] = uvrot(u2(:,ii), v2(:,ii), -mag_dev(ii)); +% end +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 + + + +%% 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; + 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 + +%% 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') + +%% If upward looking: determine range of surface bins used for instrument depth correction below! +disp('****') +surface_bins = input('Do you want to apply an offset using surface reflection? 1/0 [1]'); +if isempty(surface_bins) + surface_bins = 1; +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 + 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 + + 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 + +%% 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; + 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; + 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 + 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; +adcp.z_offset = offset; +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 +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'); + +%% Interpolate data on a regular vertical grid +Z = fliplr(blen/2:blen:max(z(:))+blen); +Zmax = max(Z); +u_interp = NaN(length(time),length(Z)); +v_interp = NaN(length(time),length(Z)); +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); + if reply_ts == 1 + ts_interp(i,ind:ind+npts-1) = data.ts.ts(1:npts,i); + end +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]:'); +if isempty(horz_interp) + horz_interp = 1; +end + +if horz_interp + filt = input('Do you want to filter tide? 1/0 [1]:'); + if isempty(filt) + filt = 1; + end + + sub = input('Do you want to subsampling data daily? 1/0 [1]:'); + if isempty(filt) + sub = 1; + end + if filt + if sub + disp('Horizontal interpolation, Filtering & Subsampling') + else + disp('Horizontal interpolation & Filtering') + end + disp('Filtering using 40h Hamming filter (other disponible filter: FFT filter, tide model)') + else + if sub + disp('Horizontal interpolation & Subsampling') + else + disp('Horizontal interpolation') + end + end +else + filt = 0; + sub = 0; +end + +[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); +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') + +if horz_interp == 0 + uintfilt = u_interp'; + vintfilt = v_interp'; + inttim = data.time; +end + +%% 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 +data.uintfilt = uintfilt(bin_start:bin_end,:); +data.vintfilt = vintfilt(bin_start:bin_end,:); +if reply_ts == 1 + data.tsintfilt = tsintfilt(bin_start:bin_end,:); +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'); + +%% Figure +niv_u = (-1:0.05:1); +niv_v = (-1:0.05:1); +close all +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,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 + +%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'); +caxis(niv_v([1 end])); +h = colorbar; +ylabel(h,'V [m s^-^1]'); +set(gca,'ydir', 'reverse'); +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 + +graph_name = [fpath_output, mooring.name '_', num2str(instr), '_U_V_int_filt_sub']; +set(hf,'Units','Inches'); +pos = get(hf,'Position'); +set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); +print(hf,graph_name,'-dpdf','-r300'); + +if reply_ts == 1 + hf=figure('position', [0, 0, 1400, 500]); + colormap jet + pcolor(inttim,Z(bin_start:bin_end),tsintfilt(bin_start:bin_end,:)); shading interp; + h = colorbar; + ylabel(h,'Target Strength [dB1]'); + set(gca,'ydir', 'reverse'); + ylabel('Depth (m)'); + 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 + + graph_name = [fpath_output, mooring.name, '_', num2str(instr), '_TS_int_filt_sub']; + set(hf,'Units','Inches'); + pos = get(hf,'Position'); + set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); + saveas(hf,graph_name,'jpeg'); +end + + + +% %% 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 +disp('****') +disp('Creating .nc 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(data.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]); +if reply_ts == 1 + ts_ID = netcdf.defVar(ncid,'ts','double',[dimidt dimidz]); +end +%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,data.Z); +%Then store my main variable +netcdf.putVar(ncid,u_ID,data.uintfilt'); +netcdf.putVar(ncid,v_ID,data.vintfilt'); +if reply_ts == 1 + netcdf.putVar(ncid,ts_ID,data.tsintfilt'); +end +%We're done, close the netcdf +netcdf.close(ncid); +disp('****') +% ------------------------------------------------------------------------------------------- diff --git a/specific/template_get_adcp_data_10w_down.m b/specific/template_get_adcp_data_10w_down.m new file mode 100644 index 0000000..ef18be4 --- /dev/null +++ b/specific/template_get_adcp_data_10w_down.m @@ -0,0 +1,682 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% template_get_adcp_data.m +% ------------------------------- +% Author : Pierre ROUSSELOT - IRD (pierre.rousselot@ird.fr) +% Jeremie HABASQUE - IRD (jeremie.habasque@ird.fr) +% ------------------------------- +% INPUTS: +% - binary raw file with .000 extension +% OUTPUTS: +% - U and V fields interpolated on a regulard grid, filtered and subsampled +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +close all; clear +addpath(genpath('../ADCP_mooring_data_processing')); +addpath('/home/proussel/Documents/OUTILS/TOOLS/nansuite'); % NaNSuitePath + +%% META information: +% Location rawfile +rawfile = 'FR10W000.000'; % binary file with .000 extension +fpath_output = './10W/'; % Output directory + +% Cruise/mooring info +cruise.name = 'PIRATA'; % cruise name +mooring.name = '10W0N'; % '0N10W' +mooring.lat = '00°00.000'; % latitude [°'] +mooring.lon = '-10°00.000'; % longitude [°'] +clock_drift = 0; % [seconds] + +% ADCP info +adcp.sn = 509; % ADCP serial number +adcp.type = '300 khz Quartermaster'; % Type : Quartermaster, longranger +adcp.direction = 'down'; % upward-looking 'up', downward-looking 'dn' +adcp.instr_depth = 300; % nominal instrument depth +instr = 1; % this is just for name convention and sorting of all mooring instruments + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Convert variables +latDegInd = strfind(mooring.lat,'°'); +lonDegInd = strfind(mooring.lon,'°'); +mooring.lat = str2double(mooring.lat(1:latDegInd-1))+str2double(mooring.lat(latDegInd+1:end-1))/60; +mooring.lon = str2double(mooring.lon(1:lonDegInd-1))+str2double(mooring.lon(lonDegInd+1:end-1))/60; +clock_drift = clock_drift/3600; % convert into hrs + +%% Read rawfile +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); +% load(raw_file) +% else +% fprintf('Read %s\n', rawfile); +% raw = read_os3(rawfile,'all'); +% save(raw_file,'raw','-v7.3'); +% end +load('10W0N_509_instr_02.mat') +load('FR22-DOWN-509.mat') + +%% 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): '); +first = 17; % time point of ADCP into and out of water +last = 12856; % determined by pressure time series +%% Correct velocity data with external T/S sensor + +%% Extract data +freq = raw.config.sysconfig.frequency; + +u2 = SerEmmpersec(first:last,:)'./1000; u2(u2<-30)=NaN; +v2 = SerNmmpersec(first:last,:)'./1000; v2(v2<-30)=NaN; +vel_err = SerErmmpersec(first:last,:)'./1000; vel_err(vel_err<-30)=NaN; +pg=SerPG4(first:last,:)'; + +% 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(SerC4cnt(first:last,:)'); + +test(:,1,:) = SerEA1cnt'; test(:,2,:) = SerEA2cnt'; test(:,3,:) = SerEA3cnt'; test(:,4,:) = SerEA4cnt'; + +ea = squeeze(mean(test(:,4,first:last),2)); +%pg = squeeze(raw.pg(:,4,first:last)); +% %time = raw.juliandate(first:last); +% d0=datenum(2006,06,26,01,0,0); +% d1=datenum(2008,09,28,21,0,0); +% time=d0:1/24:d1; +% [YY,MM,DD,hh,mm,ss] = datevec(time); +% time=julian(YY,MM,DD,hh,mm,ss); + +depth = AnDepthmm(first:last)./1000; +time = julian(SerYear(first:last)+2000,SerMon(first:last),SerDay(first:last),SerHour(first:last)); +time=time(1)+datenum(0,0,0,1,0,0).*[0:length(depth)-1]'; + +% 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+(time(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+(time(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(time)); +%mag_dev = mag_dev(first:last); + +%% Correction of magnetic deviation +disp('****') +disp('Correct magnetic deviation') +% for ii = 1 : length(mag_dev) +% [u(:,ii),v(:,ii)] = uvrot(u2(:,ii), v2(:,ii), -mag_dev(ii)); +% end +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 + + + +%% 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<20); % 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; +% 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 + +%% 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') + +%% If upward looking: determine range of surface bins used for instrument depth correction below! +disp('****') +surface_bins = 0;%input('Do you want to apply an offset using surface reflection? 1/0 [1]'); +if isempty(surface_bins) + surface_bins = 1; +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 +% 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 +% +% 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 + +%% 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; + 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; + 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 + saveas(figure(7),[fpath_output,mooring.name,'_',num2str(adcp.sn),'_instr_',num2str(instr),'_','Meridional_zonal_velocity'],'fig') +end + +%% SAVE DATA +z1 = data.z_bins; +dpt1 = data.depth; +z = data.z_bins; + +% More meta information +%adcp.comment=''; +adcp.config = raw.config; +%adcp.z_offset = offset; +%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 = 0;%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); + 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'); + +%% Interpolate data on a regular vertical grid +Z = (round(min(data.depth))+blen/2:blen:max(z(:))+blen); +Zmax = min(Z); +u_interp = NaN(length(time),length(Z)); +v_interp = NaN(length(time),length(Z)); +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); + if reply_ts == 1 + ts_interp(i,ind:ind+npts-1) = data.ts.ts(1:npts,i); + end +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]:'); +if isempty(horz_interp) + horz_interp = 1; +end + +if horz_interp + filt = input('Do you want to filter tide? 1/0 [1]:'); + if isempty(filt) + filt = 1; + end + + sub = input('Do you want to subsampling data daily? 1/0 [1]:'); + if isempty(filt) + sub = 1; + end + if filt + if sub + disp('Horizontal interpolation, Filtering & Subsampling') + else + disp('Horizontal interpolation & Filtering') + end + disp('Filtering using 40h Hamming filter (other disponible filter: FFT filter, tide model)') + else + if sub + disp('Horizontal interpolation & Subsampling') + else + disp('Horizontal interpolation') + end + end +else + filt = 0; + sub = 0; +end + +[uintfilt,vintfilt,tsintfilt,inttim,utid_baro,vtid_baro] = adcp_filt_sub_down(data,u_interp',v_interp',ts_interp',1:length(Z),reply_ts, filt, sub); +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') + +if horz_interp == 0 + uintfilt = u_interp'; + vintfilt = v_interp'; + inttim = data.time; +end + +%% 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 +data.uintfilt = uintfilt(bin_start:bin_end,:); +data.vintfilt = vintfilt(bin_start:bin_end,:); +if reply_ts == 1 + data.tsintfilt = tsintfilt(bin_start:bin_end,:); +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'); + +%% Figure +niv_u = (-1:0.05:1); +niv_v = (-1:0.05:1); +close all +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,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 + +%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'); +caxis(niv_v([1 end])); +h = colorbar; +ylabel(h,'V [m s^-^1]'); +set(gca,'ydir', 'reverse'); +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 + +graph_name = [fpath_output, mooring.name '_', num2str(instr), '_U_V_int_filt_sub']; +set(hf,'Units','Inches'); +pos = get(hf,'Position'); +set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); +print(hf,graph_name,'-dpdf','-r300'); + +if reply_ts == 1 + hf=figure('position', [0, 0, 1400, 500]); + colormap jet + pcolor(inttim,Z(bin_start:bin_end),tsintfilt(bin_start:bin_end,:)); shading interp; + h = colorbar; + ylabel(h,'Target Strength [dB1]'); + set(gca,'ydir', 'reverse'); + ylabel('Depth (m)'); + 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 + + graph_name = [fpath_output, mooring.name, '_', num2str(instr), '_TS_int_filt_sub']; + set(hf,'Units','Inches'); + pos = get(hf,'Position'); + set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); + saveas(hf,graph_name,'jpeg'); +end + + + +% %% 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 +disp('****') +disp('Creating .nc 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(data.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]); +if reply_ts == 1 + ts_ID = netcdf.defVar(ncid,'ts','double',[dimidt dimidz]); +end +%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,data.Z); +%Then store my main variable +netcdf.putVar(ncid,u_ID,data.uintfilt'); +netcdf.putVar(ncid,v_ID,data.vintfilt'); +if reply_ts == 1 + netcdf.putVar(ncid,ts_ID,data.tsintfilt'); +end +%We're done, close the netcdf +netcdf.close(ncid); +disp('****') +% ------------------------------------------------------------------------------------------- diff --git a/specific/template_get_adcp_data_10w_olddata.m b/specific/template_get_adcp_data_10w_olddata.m new file mode 100644 index 0000000..5c4137c --- /dev/null +++ b/specific/template_get_adcp_data_10w_olddata.m @@ -0,0 +1,253 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% template_get_adcp_data.m +% ------------------------------- +% Author : Pierre ROUSSELOT - IRD (pierre.rousselot@ird.fr) +% Jeremie HABASQUE - IRD (jeremie.habasque@ird.fr) +% ------------------------------- +% INPUTS: +% - binary raw file with .000 extension +% OUTPUTS: +% - U and V fields interpolated on a regulard grid, filtered and subsampled +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +close all; clear +addpath(genpath('../ADCP_mooring_data_processing')); +addpath('/home/proussel/Documents/OUTILS/TOOLS/nansuite'); % NaNSuitePath + +%% META information: +% Location rawfile +rawfile = 'FR10W000.000'; % binary file with .000 extension +fpath_output = './10W/'; % Output directory + +% Cruise/mooring info +cruise.name = 'PIRATA'; % cruise name +mooring.name = '10W0N'; % '0N10W' +mooring.lat = '00°00.000'; % latitude [°'] +mooring.lon = '-10°00.000'; % longitude [°'] +clock_drift = 0; % [seconds] + +% ADCP info +adcp.sn = 500; % ADCP serial number +adcp.type = '300 khz'; % Type : Quartermaster, longranger +adcp.direction = 'up'; % upward-looking 'up', downward-looking 'dn' +adcp.instr_depth = 100; % nominal instrument depth +instr = 1; % this is just for name convention and sorting of all mooring instruments + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Convert variables +latDegInd = strfind(mooring.lat,'°'); +lonDegInd = strfind(mooring.lon,'°'); +mooring.lat = str2double(mooring.lat(1:latDegInd-1))+str2double(mooring.lat(latDegInd+1:end-1))/60; +mooring.lon = str2double(mooring.lon(1:lonDegInd-1))+str2double(mooring.lon(lonDegInd+1:end-1))/60; +clock_drift = clock_drift/3600; % convert into hrs + +%% Read rawfile +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); +% load(raw_file) +% else +% fprintf('Read %s\n', rawfile); +% raw = read_os3(rawfile,'all'); +% save(raw_file,'raw','-v7.3'); +% end +load('UV_hour_10w_int10.mat') +d0=datenum(2004,2,6,0,0,0); +d1=datenum(2005,6,17,22,0,0); +adcp_2003_time=d0:1/12:d1; +adcp_2003_time = adcp_2003_time'; +u2 = uvmat(:,1:5976)/100; +v2 = uvmat(:,5977:5976*2)/100; +adcp_2003_bin_length = 4; +dpt1 = 4:adcp_2003_bin_length:104; +[YY,MM,DD,hh,mm,ss] = datevec(adcp_2003_time); +time=julian(YY,MM,DD,hh,mm,ss); + + + + +% Data structure +u_interp = u2'; +v_interp = v2'; +ts_interp = NaN; +data.time = time; +Z = dpt1; +reply_ts = 0; + +%% Horizontal interpolation, filtering and subsampling +disp('****') +horz_interp = 1; +filt = 1; +sub = 1; + + +[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); +% 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') + +if horz_interp == 0 + uintfilt = u_interp'; + vintfilt = v_interp'; + inttim = data.time; +end + +%% 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 +data.uintfilt = uintfilt(bin_start:bin_end,:); +data.vintfilt = vintfilt(bin_start:bin_end,:); +if reply_ts == 1 + data.tsintfilt = tsintfilt(bin_start:bin_end,:); +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'); + +%% Figure +niv_u = (-1:0.05:1); +niv_v = (-1:0.05:1); +close all +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,round(max(Z))]); +%change figure label in HH:MM +gregtick; +if filt + title({['10W0N - ZONAL VELOCITY - RDI 150 kHz (filtered from tide)']}); +else + title({[mooring.name, ' - ZONAL VELOCITY - RDI ',num2str(freq),' kHz']}); +end + +%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'); +caxis(niv_v([1 end])); +h = colorbar; +ylabel(h,'V [m s^-^1]'); +set(gca,'ydir', 'reverse'); +ylabel('Depth (m)'); +ylim([0,round(max(Z))]); +%change figure label in HH:MM +gregtick; +if filt + title({['10W0N - MERIDIONAL VELOCITY - RDI 150 kHz (filtered from tide)']}); +else + title({[mooring.name, ' - MERIDIONAL VELOCITY - RDI ',num2str(freq),' kHz']}); +end + +graph_name = ['10W0N_nb150khz_U_V_int_filt_sub']; +set(hf,'Units','Inches'); +pos = get(hf,'Position'); +set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); +print(hf,graph_name,'-dpdf','-r300'); + +if reply_ts == 1 + hf=figure('position', [0, 0, 1400, 500]); + colormap jet + pcolor(inttim,Z(bin_start:bin_end),tsintfilt(bin_start:bin_end,:)); shading interp; + h = colorbar; + ylabel(h,'Target Strength [dB1]'); + set(gca,'ydir', 'reverse'); + ylabel('Depth (m)'); + ylim([0,round(max(Z))]); + %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 + + graph_name = [fpath_output, mooring.name, '_nb150khz_TS_int_filt_sub']; + set(hf,'Units','Inches'); + pos = get(hf,'Position'); + set(hf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]); + saveas(hf,graph_name,'jpeg'); +end + + + +% %% 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 +disp('****') +disp('Creating .nc file') +[yr_start , ~, ~] = gregorian(inttim(1)); +[yr_end, ~, ~] = gregorian(inttim(length(inttim))); + +ncid = netcdf.create(['ADCP_10W0N_',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(data.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]); +if reply_ts == 1 + ts_ID = netcdf.defVar(ncid,'ts','double',[dimidt dimidz]); +end +%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,data.Z); +%Then store my main variable +netcdf.putVar(ncid,u_ID,data.uintfilt'); +netcdf.putVar(ncid,v_ID,data.vintfilt'); +if reply_ts == 1 + netcdf.putVar(ncid,ts_ID,data.tsintfilt'); +end +%We're done, close the netcdf +netcdf.close(ncid); +disp('****') +% ------------------------------------------------------------------------------------------- diff --git a/template_get_adcp_data.m b/template_get_adcp_data.m index 3bc09df..23ddead 100644 --- a/template_get_adcp_data.m +++ b/template_get_adcp_data.m @@ -16,21 +16,21 @@ addpath('/home/proussel/Documents/OUTILS/TOOLS/nansuite'); % NaNSuitePath %% META information: % Location rawfile -rawfile = 'IDMX2/_RDI_000.000'; % binary file with .000 extension -fpath_output = './IDMX2/'; % Output directory +rawfile = '_RDI_000.000'; % binary file with .000 extension +fpath_output = './10W/'; % Output directory % Cruise/mooring info -cruise.name = 'INDOMIX'; % cruise name -mooring.name = '0N130E'; % '0N10W' -mooring.lat = '00°04.066'; % latitude [°'] -mooring.lon = '129°12.41'; % longitude [°'] +cruise.name = 'PIRATA'; % cruise name +mooring.name = '10W0N'; % '0N10W' +mooring.lat = '00°00.000'; % latitude [°'] +mooring.lon = '-10°00.000'; % longitude [°'] clock_drift = 0; % [seconds] % ADCP info -adcp.sn = 13952; % ADCP serial number -adcp.type = '75 khz Quartermaster'; % Type : Quartermaster, longranger +adcp.sn = 509; % ADCP serial number +adcp.type = '150 khz Quartermaster'; % Type : Quartermaster, longranger adcp.direction = 'up'; % upward-looking 'up', downward-looking 'dn' -adcp.instr_depth = 632; % nominal instrument depth +adcp.instr_depth = 300; % nominal instrument depth instr = 1; % this is just for name convention and sorting of all mooring instruments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -428,23 +428,45 @@ disp('Saving raw data') save([fpath_output, mooring.name '_' num2str(adcp.sn) '_instr_' sprintf('%02d',instr) '.mat'],'adcp','mooring','data','-v7.3'); %% Interpolate data on a regular vertical grid -Z = fliplr(blen/2:blen:max(z(:))+blen); -Zmax = max(Z); -u_interp = NaN(length(time),length(Z)); -v_interp = NaN(length(time),length(Z)); -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); +if strcmp(adcp.direction,'up') + Z = fliplr(blen/2:blen:max(z(:))+blen); + Zmax = max(Z); + u_interp = NaN(length(time),length(Z)); + v_interp = NaN(length(time),length(Z)); + 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); + if reply_ts == 1 + ts_interp(i,ind:ind+npts-1) = data.ts.ts(1:npts,i); + end + end +else + Z = (round(min(data.depth))+blen/2:blen:max(z(:))+blen); + Zmax = min(Z); + u_interp = NaN(length(time),length(Z)); + v_interp = NaN(length(time),length(Z)); if reply_ts == 1 - ts_interp(i,ind:ind+npts-1) = data.ts.ts(1:npts,i); + 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); + if reply_ts == 1 + ts_interp(i,ind:ind+npts-1) = data.ts.ts(1:npts,i); + end end end -- GitLab