% aqdp_curr_ts.m An mfile to plot the (optionally filtered) timeseries % of current velocities from South Carolina Aquadopp data. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Use of this program is described in: % % Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., % Work, P.A., Haas, K.A., and Hanes, D.H. (2006) % South Carolina Coastal Erosion Study Data Report for Observations % October 2003 - April 2004., USGS Open-File Report 2005-1429. % % Program written in Matlab v7.1.0 SP3 % Program ran on PC with Windows XP Professional OS. % % "Although this program has been used by the USGS, no warranty, % expressed or implied, is made by the USGS or the United States % Government as to the accuracy and functioning of the program % and related program material nor shall the fact of distribution % constitute any such warranty, and no responsibility is assumed % by the USGS in connection therewith." %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % C. Sullivan 12/01/04, It is assumed that the Aquadopp data files are % stored in a directory defined by the variable % 'dataDirectory'. It assumes the user has run % spit_wad.m by G. Voulgaris in the directory to split % the .wad data file into individual burst .mat files % with pressure and velocity data. To plot low-pass % filtered current velocities set the switch 'filt' to % 'yes', otherwise (filt='no') unfiltered current % velocities are plotted. % 01/27/05, Aquadopp data has now been processed by running % READ_aquapro_beam.m which loads the raw velocites % and converts them from beam to ENU coordinates. % Spikes and bad data are nan'd by aqdp_clean.m. % Processed data is written to netcdf files by % aqdpCurr2nc.m. Timeseries plots of current velocites % are created using velocites at the bottom, 6 mbs, % 4 mbs, and the surface. % 02/24/05, Reformatting of labels and ticks, add additional % plots % 05/24/06, Update paths to files and file names. Update to % show Hs in the top panel of the plot. % 05/26/06, Add switch to show velocities in either east and % north coordinates or cross-shore and longshore % coordinates % % Dependencies: % the netcdf toolbox % rot.m % smart_interp.m % PL66TN.m % timeplt.m % dolandscape.m % suptitle.m clear more off % define data and output directories dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES'; plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\PLOTS\ADCP'; % define Aquadopp data filenames. dataFiles={'7311aqc-a.nc','7511aqc-a.nc'}; nFiles=length(dataFiles); % define switches filt='yes'; %filtering switch, 'yes' to filter data, or 'no' coord='xl'; %coordinates switch,'xl' for cross-shore, longshore; 'en' for %east,north sticks='no';%sticks switch, 'yes' to create magnitude,direction plots and %u,v plots, 'no' to create just u,v plots % set x-limits, ticks; x-limits extend 5% beyond startTime and stopTime startTime=julian(2003,10,25,00); stopTime=julian(2004,04,25,00); timediff=stopTime-startTime; xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; xt=[]; for yy=2003:2004 for mm=1:12 tic=julian(yy,mm,01,00); if tic>=xlims(1) & tic<=xlims(2) xt=[xt; tic]; end end end % set y-limits, labels ylims=[-30 30; -30 30; -30 30; -30 30; 0 2]; switch coord case 'en' yl_bott6={ {'Northward','Velocity','at the Bottom','(cm/s)'},... {'Eastward','Velocity','at the Bottom','(cm/s)'},... {'Northward','Velocity','at 6 mbs','(cm/s)'},... {'Eastward','Velocity','at 6 mbs','(cm/s)'},... {'Site 2','Significant Wave','Height (m)'} }; yl_4surf={ {'Northward','Velocity','at the 4 mbs','(cm/s)'},... {'Eastward','Velocity','at the 4 mbs','(cm/s)'},... {'Northward','Velocity','at the Surface','(cm/s)'},... {'Eastward','Velocity','at the Surface','(cm/s)'},... {'Site 2','Significant Wave','Height (m)'} }; case 'xl' yl_bott6={ {'Cross-shore','Velocity','at the Bottom','(cm/s)'},... {'Longshore','Velocity','at the Bottom','(cm/s)'},... {'Cross-shore','Velocity','at 6 mbs','(cm/s)'},... {'Longshore','Velocity','at 6 mbs','(cm/s)'},... {'Site 2','Significant Wave','Height (m)'} }; yl_4surf={ {'Cross-shore','Velocity','at the 4 mbs','(cm/s)'},... {'Longshore','Velocity','at the 4 mbs','(cm/s)'},... {'Cross-shore','Velocity','at the Surface','(cm/s)'},... {'Longshore','Velocity','at the Surface','(cm/s)'},... {'Site 2','Significant Wave','Height (m)'} }; otherwise error('Not a valid coordinates switch!') end switch sticks case 'yes' yl_sticks={ {'Velocity','at the Bottom','(cm/s)'},... {'Velocity','at 6 mbs','(cm/s)'},... {'Velocity','at 4 mbs','(cm/s)'},... {'Velocity','at the Surface','(cm/s)'},... {'Site 2','Significant Wave','Height (m)'} }; end % Interpolate all data onto the array 'jd_int' for timeplt.m. Data gaps % larger than 'maxGap' hours are not interpolated over delta_t = 900/(3600*24); %for 15 min. data jd_int = [startTime:delta_t:stopTime]'; maxGap = 1; %hr % Get the color scheme to distinguish between data sites col=get(0,'defaultaxescolororder'); % load site 2 wave data for the top panel % of the plot nc=netcdf(fullfile(dataDirectory,'7221adwp-cal.nc'), 'nowrite'); jd_waves{1} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24); waveDelta_t=gmean(diff(jd_waves{1})); hsig{1}=nc{'wh_4061'}(:); theFillValue=nc{'wh_4061'}.FillValue_(:); badInd=find(hsig{1} >= theFillValue); hsig{1}(badInd)=nan; nc=close(nc); nc=netcdf(fullfile(dataDirectory,'7421adwp-cal.nc'), 'nowrite'); jd_waves{2} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24); hsig{2}=nc{'wh_4061'}(:); theFillValue=nc{'wh_4061'}.FillValue_(:); badInd=find(hsig{2} >= theFillValue); hsig{2}(badInd)=nan; nc=close(nc); hsig_all_comb = smart_interp([jd_waves{1};jd_waves{2}], [hsig{1};hsig{2}],... jd_int, 3); switch filt case 'yes' hsig_all_comb=pl66tn(hsig_all_comb, 1, 33); end % Initialize variables jd=cell(nFiles,1); depth=cell(nFiles,1); u_bottom=cell(nFiles,1); u_6mbs=cell(nFiles,1); u_4mbs=cell(nFiles,1); u_surface=cell(nFiles,1); v_bottom=cell(nFiles,1); v_6mbs=cell(nFiles,1); v_4mbs=cell(nFiles,1); v_surface=cell(nFiles,1); varList={'u','v'}; ncVarList={'u_1205','v_1206'}; nVars=length(ncVarList); for f=1:length(dataFiles) nc=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite'); jd{f}=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24); delta_t=gmean(diff(jd{f})); depth{f}=nc{'depth'}(:); site=nc.DESCRIPT(6); for v=1:length(ncVarList) ncobj=nc{ncVarList{v}}; theFillValue=ncobj.FillValue_(:); data=ncobj(:); bads=find(data>=theFillValue); data(bads)=nan; questBins=ncobj.bins_questioned(:); %velocity in bottom bin eval([varList{v},'_bottom{f}=data(:,1);']) %velocity in surface bin where surface is upward-most good bin if isempty(questBins) eval([varList{v},'_surface{f}=data(:,end);']) else eval([varList{v},'_surface{f}=data(:,questBins(1)-1);']) %uppermost good bin at site 8 is at 5.6mbs for the second %deployment so don't show this on the plots if f==2 u_surface{f}=nan(size(data,1),1); v_surface{f}=nan(size(data,1),1); end end %velocity in bin 4 mbs [binInd,dist2bin]=value2Index(depth{f},4,0.5); eval([varList{v},'_4mbs{f}=data(:,binInd);']) %velocity in bin 6 mbs [binInd,dist2bin]=value2Index(depth{f},6,0.5); eval([varList{v},'_6mbs{f}=data(:,binInd);']) if v==length(ncVarList) % (optionally convert velocities from east, north coordinates % to longshore, cross-shore coordinates using rot.m. Assume % a shoreline orientation of 45 deg. switch coord case 'xl' disp(['Converting from north,east coords '... 'to cross,long coords']) [u_bottom{f},v_bottom{f}]=rot(u_bottom{f},v_bottom{f},-45); [u_6mbs{f},v_6mbs{f}]=rot(u_6mbs{f},v_6mbs{f},-45); [u_4mbs{f},v_4mbs{f}]=rot(u_4mbs{f},v_4mbs{f},-45); [u_surface{f},v_surface{f}]=rot(u_surface{f},v_surface{f},-45); case 'en' disp('Not converting coordinate systems') otherwise disp('Not a valid coordinate system!') end end end nc=close(nc); if mod(f,2)==0 % combine data from deployment 1 and % deployment 2 jd_all=[jd{f-1}; jd{f}]; for var=1:nVars eval([varList{var},'_bottom_all=[',... varList{var},'_bottom{f-1}; ',varList{var},'_bottom{f}];']) eval([varList{var},'_6mbs_all=[',... varList{var},'_6mbs{f-1}; ',varList{var},'_6mbs{f}];']) eval([varList{var},'_4mbs_all=[',... varList{var},'_4mbs{f-1}; ',varList{var},'_4mbs{f}];']) eval([varList{var},'_surface_all=[',... varList{var},'_surface{f-1}; ',varList{var},'_surface{f}];']) end % interpolate data onto 'jd_int' for var=1:nVars eval([varList{var},'_bottom_all_int=smart_interp(jd_all, ',... varList{var},'_bottom_all, jd_int, maxGap);']) eval([varList{var},'_6mbs_all_int=smart_interp(jd_all, ',... varList{var},'_6mbs_all, jd_int, maxGap);']) eval([varList{var},'_4mbs_all_int=smart_interp(jd_all, ',... varList{var},'_4mbs_all, jd_int, maxGap);']) eval([varList{var},'_surface_all_int=smart_interp(jd_all, ',... varList{var},'_surface_all, jd_int, maxGap);']) end % (optionally filter the velocities) switch filt case 'yes' disp('Low-pass filtering u,v current velocities') for var=1:nVars eval([varList{var},'_bottom_all_int=PL66TN(',... varList{var},'_bottom_all_int, delta_t*24, 33);']) eval([varList{var},'_6mbs_all_int=PL66TN(',... varList{var},'_6mbs_all_int, delta_t*24, 33);']) eval([varList{var},'_4mbs_all_int=PL66TN(',... varList{var},'_4mbs_all_int, delta_t*24, 33);']) eval([varList{var},'_surface_all_int=PL66TN(',... varList{var},'_surface_all_int, delta_t*24, 33);']) end otherwise disp('Not filtering current velocities') end %%%%%%%%%%%%%%%%%%%%%% STICK PLOT %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % call timeplt.m to create individual STICK plots switch sticks case 'yes' % Make vectors clear i uv_bottom_all_int = u_bottom_all_int + i.*v_bottom_all_int; uv_6mbs_all_int = u_6mbs_all_int + i.*v_6mbs_all_int; uv_4mbs_all_int = u_4mbs_all_int + i.*v_4mbs_all_int; uv_surface_all_int = u_surface_all_int + i.*v_surface_all_int; % (clearing global timeplt*, the ResizeFcn, % and the Tag are necessary to release % timeplt.m's control over your figure!) xdata=jd_int; ydata=[uv_bottom_all_int uv_6mbs_all_int,... uv_4mbs_all_int uv_surface_all_int,... hsig_all_comb]; xdata=xdata(1:8:end); %show a stick every 2 hours ydata=ydata(1:8:end,:); panelOrd=[1:min(size(ydata))]; figure han=timeplt(xdata,ydata,panelOrd,ylims); clear global timeplt* set(gcf,'ResizeFcn',[]) set(gcf,'Tag','') dolandscape fixpaper % reformat xlims,ylabels, and line color for h=1:length(han) axes(han(h)) set(gca,'xlim',xlims,'xtick',xt) ylabhan=get(gca,'Ylabel'); set(ylabhan,'String',yl_sticks{h}) pos=get(ylabhan,'position'); set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) %align ylabels nicely if h~=length(han) obj=findobj(gca,'type','line','color','b'); set(obj,'color',col(7,:)) end if h==4 xtxt=julian(2004,3,15); ytxt=15; text(xtxt,ytxt,'No Data','horizontalalignment','center') end end % re-do xlabels axes(han(1)) oldhan=findobj(gca,'Tag','xtick labels'); delete(oldhan); for tic=1:length(xt) xtg=gregorian(xt(tic)); xtd=datenum(xtg); xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label xtlpos=ylims(1,1)-2; text(xt(tic),xtlpos,xtl,'horizontalalignment','center',... 'verticalalignment','top') end % give it a title and export to .pdf switch filt case 'yes' set(gcf,'PaperPositionMode','Auto') suptitle({'Low-Passed Water Flow Observations',['Site ',site]}); print(fullfile(plotDirectory,['s',site,'_spdir_s46b_lp']),'-dpdf'); otherwise set(gcf,'PaperPositionMode','Auto') suptitle({'Water Flow Observations',['Site ',site]}); print(fullfile(plotDirectory,['s',site,'_spdir_s46b']),'-dpdf'); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%% U,V PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % call timeplt.m to create individual U, V plots at the BOTTOM % and 6 MBS xdata=jd_int; ydata=[v_bottom_all_int u_bottom_all_int,... v_6mbs_all_int u_6mbs_all_int,... hsig_all_comb]; panelOrd=[1:min(size(ydata))]; figure han=timeplt(xdata,ydata,panelOrd,ylims); clear global timeplt* set(gcf,'ResizeFcn',[]) set(gcf,'Tag','') dolandscape fixpaper % reformat xlims, ylabels, and line color for h=1:length(han) axes(han(h)) set(gca,'xlim',xlims,'xtick',xt) ylabhan=get(gca,'Ylabel'); set(ylabhan,'String',yl_bott6{h}) pos=get(ylabhan,'position'); set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) if h~=length(han) obj=findobj(gca,'type','line','color','b'); set(obj,'color',col(7,:)) end end % re-do xlabels axes(han(1)) oldhan=findobj(gca,'Tag','xtick labels'); delete(oldhan); for tic=1:length(xt) xtg=gregorian(xt(tic)); xtd=datenum(xtg); xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; xtlpos=ylims(1,1)-1.5; text(xt(tic),xtlpos,xtl,'horizontalalignment','center',... 'verticalalignment','top') end % give it a title and export to .pdf switch filt case 'yes' set(gcf,'PaperPositionMode','Auto') suptitle({'Low-Passed Water Flow Observations',['Site ',site]}); switch coord case 'en' print(fullfile(plotDirectory,['s',site,'_en_6b_lp']),'-dpdf'); case 'xl' print(fullfile(plotDirectory,['s',site,'_xl_6b_lp']),'-dpdf'); otherwise disp('Not a valid coordinate system!') end otherwise set(gcf,'PaperPositionMode','Auto') suptitle({'Water Flow Observations',['Site ',site]}); switch coord case 'en' print(fullfile(plotDirectory,['s',site,'_en_6b']),'-dpdf'); case 'xl' print(fullfile(plotDirectory,['s',site,'_xl_6b']),'-dpdf'); otherwise disp('Not a valid coordinate system!') end end % call timeplt.m to create individual U, V plots at the 4 MBS % and the SURFACE xdata=jd_int; ydata=[v_4mbs_all_int u_4mbs_all_int,... v_surface_all_int u_surface_all_int,... hsig_all_comb]; panelOrd=[1:min(size(ydata))]; figure han=timeplt(xdata,ydata,panelOrd,ylims); clear global timeplt* set(gcf,'ResizeFcn',[]) set(gcf,'Tag','') dolandscape fixpaper % reformat xlims, ylabels, and line color for h=1:length(han) axes(han(h)) set(gca,'xlim',xlims,'xtick',xt) ylabhan=get(gca,'Ylabel'); set(ylabhan,'String',yl_4surf{h}) pos=get(ylabhan,'position'); set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) if h~=length(han) obj=findobj(gca,'type','line','color','b'); set(obj,'color',col(7,:)) end if h==3 || h==4 xtxt=julian(2004,3,15); ytxt=15; text(xtxt,ytxt,'No Data','horizontalalignment','center') end end % re-do xlabels axes(han(1)) oldhan=findobj(gca,'Tag','xtick labels'); delete(oldhan); for tic=1:length(xt) xtg=gregorian(xt(tic)); xtd=datenum(xtg); xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; xtlpos=ylims(1,1)-1.5; text(xt(tic),xtlpos,xtl,'horizontalalignment','center',... 'verticalalignment','top') end % give it a title and export to .pdf switch filt case 'yes' set(gcf,'PaperPositionMode','Auto') suptitle({'Low-Passed Water Flow Observations',['Site ',site]}); switch coord case 'en' print(fullfile(plotDirectory,['s',site,'_en_s4_lp']),'-dpdf'); case 'xl' print(fullfile(plotDirectory,['s',site,'_xl_s4_lp']),'-dpdf'); otherwise disp('Not a valid coordinate system!') end otherwise set(gcf,'PaperPositionMode','Auto') suptitle({'Water Flow Observations',['Site ',site]}); switch coord case 'en' print(fullfile(plotDirectory,['s',site,'_en_s4']),'-dpdf'); case 'xl' print(fullfile(plotDirectory,['s',site,'_xl_s4']),'-dpdf'); otherwise disp('Not a valid coordinate system!') end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end end