% wind_ts.m  An mfile to plot the timeseries of North Carolina wind data
%            from NDBC moored buoys, NDBC C-MAN stations, NOS water level
%            stations, and the COAMPS wind model.

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Armstrong, B.N., Warner, J.C., Voulgaris, G., List, J.H., Thieler, E.R., 
% Martini, M.M., Montgomery, E.T., 2011, Carolinas Coastal Change Processes 
% Project Data Report for Observations near Diamond Shoals, NC from January-April 
% 2009: U.S. Geological Survey Open-File Report 2011-1156.
%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                     
% Dependencies:
%       the netcdf toolbox
%       julian.m
%       timeplt.m
%       dolandscape.m
%       fixpaper.m
%       suptitle.m

clear
more off

% define data and output directories
dataDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2009DVD\supporting data';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2009DVD\supporting data';

% define data file names
dataFiles={ 'FRFwind_2009.nc','NAM_wind.nc','NDBC_wind.nc'};
nFiles=length(dataFiles);

%P-NAM-CONUS_12km-noaaport_best.ncd
% define parameters for plots
startTime=julian(2009,01,01,00);        %plot start stime
stopTime=julian(2009,05,01,00);         %plot stop time  
col=get(0,'defaultaxescolororder');     %plot color order
ylab={ {'Wind','Direction (deg)'},...   %plot ylabels
       {'Wind','Speed (m/s)'},...
       {'North Wind','Velocity (m/s)'},...
       {'East Wind','Velocity (m/s)'},...
       {'Wind','Velocity (m/s)'} };
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %past start time and stop time
ylims=[ 0 360;...   %plot ylimits
        0  25;...
      -20  20;...
      -20  20;...
      -15  15];
xt=[];
for yy=2009:2009
    for mm=1:12
        tic=julian(yy,mm,01,00); 
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; tic];   % place xticks on the first of every month
        end
    end
end

% initialize variables
jd=cell(nFiles,1);
wind_u=cell(nFiles,1);
wind_v=cell(nFiles,1);
wind_uv=cell(nFiles,1);
wind_spd=cell(nFiles,1);
wind_dir=cell(nFiles,1);
for f=1:length(dataFiles)
    station_id=dataFiles{f}(1:3);
    nc=netcdf([dataDirectory '\' dataFiles{f}]);
    jd=nc{'windtime'}(:);
    windstart=find((jd(:) >= startTime),1,'first');
    windend=find((jd(:) <= stopTime),1,'last');
    jd=jd(windstart:windend);
    wind_u=nc{'U'}(windstart:windend);
    wind_u(wind_u>=98)=NaN;
    wind_v=nc{'V'}(windstart:windend);
    wind_v(wind_v>=98)=NaN;
    close(nc);
    %FRF Winds are measured at the end of the pier at an elevation of 19.36 m
    %(NGVD) presently using an RM Young Marine anemometer. Resultant wind speeds
    %and directions are determined by vector averaging the data and are reported
    %relative to true north (Shore normal is 70 deg true). Wind directions are
    %reported as the direction the wind is blowing from, clockwise from
    %true north.
    if strcmp(station_id,'FRF')
        wind_u=-wind_u;
        wind_v=-wind_v;
    end
    %NDBC Wind direction (the direction the wind is coming from in degrees
    %clockwise from true N)
    if strcmp(station_id,'NDB')
        wind_u=-wind_u;
        wind_v=-wind_v;
        station_id='NDBC';
    end
    %calculate wind dir and spd
    wind_spd=sqrt(wind_u.^2+wind_v.^2);
    wind_dir=(atan2(wind_v,wind_u))*(180/pi());      %east component of wind, m/s
    wind_dir(wind_dir<0)=wind_dir(wind_dir<0)+360;
        
        delta_t=gmean(diff(jd(:)));     %time step of the data
        
        % nan bad data
           % for stations with hourly data, show
        % a wind stick every 3 hours 
        wind_uv=nan*ones(size(jd));
        
            interv=ceil(3/(delta_t*24));
             wind_uv=wind_u + i*wind_v;
             wind_uv(isnan(wind_uv))=0;
             
       
        % call timeplt.m to create individual plots
        % (clearing global timeplt*, the ResizeFcn,
        % and the Tag are necessary to release
        % timeplt.m's control over your figure!)
        xdata=jd(1:interv:end);
        ydata=[wind_dir(1:interv:end) wind_spd(1:interv:end) wind_v(1:interv:end) wind_u(1:interv:end) wind_uv(1:interv:end)];
        panelOrd=[1:min(size(ydata))]; 
        figure
        han=timeplt2(xdata,ydata,panelOrd);
        clear global timeplt*
        set(gcf,'ResizeFcn',[])
        set(gcf,'Tag','')
        dolandscape
        fixpaper
 
        % reformat xlims, ylims, ylabels, and line color
        for h=1:length(han)
            axes(han(h))
            set(gca,'xlim',xlims,'xtick',xt)
            set(gca,'ylim',ylims(h,:))
            ylabhan=get(gca,'Ylabel');
            set(ylabhan,'String',ylab{h})
            pos=get(ylabhan,'position');
            linehan=findobj(gca,'Type','line','color','b');
        end

        % re-do xlabels
        axes(han(1))
        oldhan=findobj(gca,'Tag','xtick labels');
        delete(oldhan);set(gca,'XTickLabel',{''})
        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)-1.5;
        text(xt(tic),xtlpos,xtl,'horizontalalignment','center',...
            'verticalalignment','top')
        end
        
        % give it a title and export to .pdf
        set(gcf,'PaperPositionMode','Auto')
        if strcmp(station_id,'NAM')
            suptitle(['Wind Calculated from the ',station_id,' Model'])
            print(fullfile(plotDirectory,['winds_',station_id]),'-dpdf')
        else
            suptitle(['Wind Observations from ',station_id,' Station']);
            print(fullfile(plotDirectory,['winds_',station_id]),'-dpdf');
        end

        % get data ready for combined plot
        % interpolate the COAMPS wind data
        % onto an hourly time base, b/c timeplt.m
        % assumes the data is the same length!
        if strcmp(station_id,'FRF')
            wind_uf=wind_u;
            wind_vf=wind_v;
            wind_spdf=wind_spd;
            wind_dirf=wind_dir;
            jdf=jd;
        end
        if strcmp(station_id,'NDBC')
            wind_un=wind_u;
            wind_vn=wind_v;
            wind_spdn=wind_spd;
            wind_dirn=wind_dir;
            jdn=jd;
        end
        
            
end
