% obs_ts.m An mfile to plot the (optionally filtered) % timeseries of current velocity data from South % Carolina ADV 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 03/18/05, Raw ADV data is processed and saved to netCDF files. % It is assumed that the netCDF files are stored in a % directory defined by the variable 'dataDirectory'. % To plot low-pass filtered OBS data, set the % switch 'filt' to 'yes', otherwise (filt='no') % unfiltered OBS data is plotted. % % Dependencies: % the netcdf toolbox % smart_interp.m % PL66TN.m % timeplt.m % dolandscape.m % suptitle.m clear more off % Define data and output directories dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATAFILES'; plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\PLOTS\OBS'; % Define OBS data filenames (OBS data stored w/ % ADV and PCADP data) dataFiles = {'adv7282vp-cal2.nc','adv7482vp-cal2.nc',... 'adv7261vp-cal2.nc','adv7461vp-cal2.nc',... '728pcvp-cal2.nc','748pcvp-cal2.nc',... 'adv7262vp-cal2.nc','adv7462vp-cal2.nc'}; nFiles = length(dataFiles); % Define limits and x-ticks for the plots startTime = julian(2003,10,25,00); % OBS plot start time stopTime = julian(2004,04,25,00); % OBS plot stop time timediff = stopTime-startTime; % OBS plot xlimits extend 5 percent beyond 'startTime' xlims = [(startTime-timediff*0.05) (stopTime+timediff*0.05)]; % and 'stopTime' xticks = []; for yy = 2003:2004 for mm = 1:12 tick = julian(yy,mm,01,00); if tick >= xlims(1) & tick <= xlims(2) xticks = [xticks; tick]; % fir xticks on the first of every month end end end ylims=[0 2]; % OBS plot ylimits % Define whether to filter and/or interpolate % data filtr = 'no'; intrp = 'no'; switch intrp case {'yes','y','Yes','Y'} % 'jd_int' is the time array onto which all data is % interpolated. int_delta_t = 1/24; % 1 hour time spacing obs_time_int = [startTime:int_delta_t:stopTime]'; % data gaps larger than 'maxGap' hours are not % interpolated over maxGap = 3; end % Significant wave height from site 2 is % shown in the top panel of all plots wvs1 = netcdf('7221adw-a.nc'); wvs2 = netcdf('7421adw-a.nc'); wave_time = [wvs1{'time'}(:) + wvs1{'time2'}(:)/3600/1000/24;... wvs2{'time'}(:) + wvs2{'time2'}(:)/3600/1000/24]; hsig = [wvs1{'wh_4061'}(:); wvs2{'wh_4061'}(:)]; theFillValue = wvs1{'wh_4061'}.FillValue_(:); badInd = find(hsig >= theFillValue); hsig(badInd) = nan; % Initialize variables for OBS data obs_time = cell(nFiles,1); obs1_data = cell(nFiles,1); obs2_data = cell(nFiles,1); obs_delta_t = cell(nFiles,1); obs1_serial = cell(nFiles,1); obs2_serial = cell(nFiles,1); for f = 1:length(dataFiles) % Load the OBS data if ~isempty(dataFiles{f}) obs = netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite'); obs_time{f} = obs{'time'}(:) + obs{'time2'}(:)/3600/1000/24; % 1st OBS obs1_data{f} = obs{'NEP1_56'}(:); if isempty(obs1_data{f}) obs1_data{f} = nan.*ones(size(obs_time{f})); else obs1_serial{f} = obs{'NEP1_56'}.serial(:); theFillValue = obs{'NEP1_56'}.FillValue_(:); badInd = find(obs1_data{f} >= theFillValue); obs1_data{f}(badInd) = nan; end % 2nd OBS obs2_data{f} = obs{'NEP2_56'}(:); if isempty(obs2_data{f}) obs2_data{f} = nan.*ones(size(obs_time{f})); else obs2_serial{f} = obs{'NEP2_56'}.serial(:); theFillValue = obs{'NEP2_56'}.FillValue_(:); badInd = find(obs2_data{f} >= theFillValue); obs2_data{f}(badInd) = nan; end obs_delta_t{f} = gmean(diff(obs_time{f})); %time step, julian days site=obs.DESCRIPT(6); %site # obs=close(obs); end if mod(f,2)==0 % Set the site if f==2 || f==6 site = '7A'; else site = '6'; end % Combine data from deployment 1 and % deployment 2 obs_time_all = [obs_time{f-1}; nan; obs_time{f}]; obs1_data_all = [obs1_data{f-1}; nan; obs1_data{f}]; obs2_data_all = [obs2_data{f-1}; nan; obs2_data{f}]; % Get serial numbers for OBS's if strcmp(obs1_serial{f-1},obs1_serial{f}) serial1 = obs1_serial{f}; elseif isempty(obs1_serial{f-1}) && ~isempty(obs1_serial{f}) serial1 = obs1_serial{f}; elseif ~isempty(obs1_serial{f-1}) && isempty(obs1_serial{f}) serial1 = obs1_serial{f-1}; else serial1 = 'No 1st OBS'; end if strcmp(obs2_serial{f-1},obs2_serial{f}) serial2 = obs2_serial{f}; elseif isempty(obs2_serial{f-1}) && ~isempty(obs2_serial{f}) serial2 = obs2_serial{f}; elseif ~isempty(obs2_serial{f-1}) && isempty(obs2_serial{f}) serial2 = obs2_serial{f-1}; else serial2 = 'No 2nd OBS'; end % (Optionally) Interpolate OBS data switch intrp case {'yes','y','Yes','Y'} disp('Interpolating wave data') hsig_int = smart_interp(wave_time, hsig, obs_time_int, maxGap); hsig = hsig_int; wave_time = obs_time_int; disp('Interpolating OBS data') obs1_data_all_int = smart_interp(obs_time_all, obs1_data_all, obs_time_int, maxGap); obs1_data_all = obs1_data_all_int; obs2_data_all_int = smart_interp(obs_time_all, obs2_data_all, obs_time_int, maxGap); obs2_data_all = obs2_data_all_int; obs_time_all = obs_time_int; delta_t = int_delta_t; otherwise disp('Not interpolating data') % delta_t = obs_delta_t{f}; %TODO: different delta-t's? end % % (Optionally) Filter OBS data % switch filt % case {'yes','y','Yes','Y'} % disp('Low-pass filtering OBS data'); % obs1_data_all=PL66TN(obs_data_all, delta_t*24, 33); % obs2_data_all=PL66TN(obs_data_all, delta_t*24, 33); % otherwise % disp('Not filtering data'); % end % Plot figure('pos',[100 100 1000 800]) orient landscape AX(1) = subplot(3,1,1); plot(wave_time,hsig) title('Significant Wave Height from Site 2') AX(2) = subplot(3,1,2); plot(obs_time_all, obs1_data_all) title(serial1) AX(3) = subplot(3,1,3); plot(obs_time_all, obs2_data_all) title(serial2) for a = 1:length(AX) ax = AX(a); % set ylimits and labels if a>=2 set(ax,'ylim',ylims) set(get(ax,'Ylabel'),'String','OBS voltage (v)') else set(ax,'ylim',[0 2]) set(get(ax,'Ylabel'),'String',{'Signifiant wave','height (m)'}) end % set xlimits and labels set(ax,'xlim',xlims) set(ax,'xtick',xticks) set(ax,'xticklabel',[]) if a == length(AX) for t=1:length(xticks) xtg=gregorian(xticks(t)); xtd=datenum(xtg); xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label xtlpos=-0.5; text(xticks(t),xtlpos,xtl,'horizontalalignment','center',... 'verticalalignment','top') end end set(ax,'box','on') suptitle({['Site ',site,' OBS Observations'],... [dataFiles{f-1},', ',dataFiles{f}]}) end % % combine data from all sites in an array for plotting % ind=f/2; % for var=1:nVars % eval([varList{var},'_all_comb(:,ind)=',... % varList{var},'_all_int;']) % end end end % %%%%%%%%%%%%%%%%%%%%%% PLOT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % decimate the wind % uv_wind_all_int_dec = nan*ones(size(uv_wind_all_int)); % uv_wind_all_int_dec(1:5:end) = uv_wind_all_int(1:5:end); % % call timeplt.m % xdata=jd_int; % ydata=[vstd_all_comb,... % ustd_all_comb,... % v_all_comb,... % u_all_comb,... % brange_all_comb,... % pressure_all_comb]; %uv_wind_all_int_dec]; % panelOrd=[repmat(1,1,nFiles/2),... % repmat(2,1,nFiles/2),... % repmat(3,1,nFiles/2),... % repmat(4,1,nFiles/2),... % repmat(5,1,nFiles/2),... % repmat(6,1,nFiles/2)];%repmat(6,1,1)]; % figure % han=timeplt(xdata,ydata,panelOrd,ylims); % clear global timeplt* % set(gcf,'ResizeFcn',[]) % set(gcf,'Tag','') % LEG=legend(han(5),'Site 7 ADV1','Site 7 ADV2','Site 6 ADV1','Site 6 ADV2'); % set(LEG,'Location','NorthOutside') % set(LEG,'Orientation','Horizontal') % figCol=get(gcf,'color'); % set(LEG,'color',figCol) % set(LEG,'edgecolor',figCol) % dolandscape % fixpaper % legpos=get(LEG,'position'); % pos=get(han(6),'position'); % set(LEG,'position',[pos(1) pos(2)+pos(4) pos(3) legpos(4)]) % % % reformat xlims, ylabels % for h=1:length(han) % axes(han(h)) % set(gca,'xlim',xlims,'xtick',xt) % ylabhan=get(gca,'Ylabel'); % set(ylabhan,'String',ylab{h}) % pos=get(ylabhan,'position'); % set(get(gca,'Ylabel'),'position',[julian([2003,10,6,0,0,0]) pos(2) pos(3)]) % 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-Pass Filtered Wind and Bottom Current Observations',... % 'All Sites'}]); % print(fullfile(plotDirectory,['adv_all_lp']),'-dpdf'); % otherwise % set(gcf,'PaperPositionMode','Auto') % suptitle([{'Wind and Bottom Current Observations',... % 'All Sites'}]); % print(fullfile(plotDirectory,['adv_all']),'-dpdf'); % end