% 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:
%       julian.m

clear
more off

% define data and output directories
dataDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\files';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\PLOTS';

% define data file names
dataFiles={ 'Met_FRF_02_2010.nc','41025h2010.nc'};
nFiles=length(dataFiles);
name={'FRF waverider','NDBC 41025'};
name2={'FRF','NDBC'};
%P-NAM-CONUS_12km-noaaport_best.ncd
% define parameters for plots
startTime=julian(2010,2,1,00);
stopTime=julian(2010,2,31,00);
col=get(0,'defaultaxescolororder');     %plot color order
ylab={ {'Wind','Direction (deg)'},...   %plot ylabels
       {'Wind','Speed (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=[ -5 365;...   %plot ylimits
        0  25];
xt=[];
for yy=2010:2010
    for mm=2:2
        for dd=1:4:31
            tic=julian(yy,mm,dd,00);
            if tic>=xlims(1) & tic<=xlims(2)
                xt=[xt; tic];   % place xticks on the first of every month
            end
        end
    end
end

% initialize variables
jd=cell(nFiles,1);
wind_dir=cell(nFiles,1);
wind_speed=cell(nFiles,1);
for f=1:length(dataFiles)
    station_id=name{f};
    wind=(fullfile(dataDirectory,dataFiles{f}));
    jd=julian([ncread(wind,'Year'),ncread(wind,'Month'),...
            ncread(wind,'Day'),ncread(wind,'Hour'),ncread(wind,'Minute'),zeros(length(ncread(wind,'Minute')))]);
    windstart=find((jd(:) >= startTime),1,'first');
    windend=find((jd(:) <= stopTime),1,'last');
%     jd=jd(windstart:windend);
    wind_dir=ncread(wind,'Direction');
    wind_dir(wind_dir>=360)=NaN;
    wind_speed=ncread(wind,'Speed');
    wind_speed(wind_speed>=99)=NaN;
    
    %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.
    %NDBC Wind direction (the direction the wind is coming from in degrees
    %clockwise from true N)       
        delta_t=gmean(diff(jd(:)));     %time step of the data
        
        
        % create plot
        
        subplot(2,1,1)
        plot(jd, wind_dir,'.')
        ylabel(ylab{1});
        set(gca,'ylim',ylims(1,:))
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[])
        for tic=1:length(xt)
            xtg=gregorian(xt(tic));
            xtd=datenum(xtg);
            xtl={datestr(xtd,6)}; %2 line label
            text(xt(tic),-7,xtl,'horizontalalignment','center',...
                                 'verticalalignment','top')
        end
        title(['Wind Direction and Speed' name(f)]);

        subplot(2,1,2)
        plot(jd,wind_speed)
        ylabel(ylab{2});
        set(gca,'ylim',ylims(2,:))
        set(gca,'xlim',xlims,'xtick',xt,'xticklabel',[])
        for tic=1:length(xt)
            xtg=gregorian(xt(tic));
            xtd=datenum(xtg);
            xtl={datestr(xtd,6)}; %2 line label
            text(xt(tic),-1,xtl,'horizontalalignment','center',...
                'verticalalignment','top')
        end
        
        
        % give it a title and export to .pdf
            print(fullfile(plotDirectory,[name2{f},'_wind']),'-dpdf');            
end
