% wind_ts.m  An mfile to plot the timeseries of South 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:
%
% 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 02/22/05, all data is hourly, except the COAMPS data which is
%                       every 3 hours.  
%                     
% Dependencies:
%       the netcdf toolbox
%       julian.m
%       timeplt.m
%       dolandscape.m
%       fixpaper.m
%       suptitle.m

clear
more off

% define data and output directories
dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATAFILES';
plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\PLOTS\MET';

% define data file names
dataFiles={ '41013h1-stdMet.nc',...
            'FPSN7h1-stdMet.nc',...
            '8659897h1-stdMet.nc',...
            '8661070h1-stdMet.nc',...
            'COAMPSh3-wind.nc' };
nFiles=length(dataFiles);

% define parameters for plots
startTime=julian(2003,10,25,00);        %plot start stime
stopTime=julian(2004,04,25,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;...
      -20  20];
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];   % 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:nFiles
    
        nc=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
        jd{f}=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
        wind_u{f}=nc{'WU_422'}(:);      %east component of wind, m/s
        wind_v{f}=nc{'WV_423'}(:);      %north component of wind, m/s
        wind_spd{f}=nc{'WS_401'}(:);    %wind speed, m/s
        wind_dir{f}=nc{'WD_412'}(:);    %wind direction (FROM), degrees_T
        theFillValue=nc.VAR_FILL(:);    %netcdf fill value
        delta_t=gmean(diff(jd{f}));     %time step of the data
        inst_type=nc.INST_TYPE(:);      %station type
        inst_id=nc.instrument_id(:);    %station number
        model=nc.model_name(:);         %model name
        nc=close(nc);
        
        % isolate the subset of data between 
        % start time and stop time
        startInd=find(jd{f}(:) >= startTime,1,'first');
        stopInd=find(jd{f}(:) <= stopTime,1,'last');
        jd{f}=jd{f}(startInd:stopInd);
        wind_u{f}=wind_u{f}(startInd:stopInd);
        wind_v{f}=wind_v{f}(startInd:stopInd);
        wind_spd{f}=wind_spd{f}(startInd:stopInd);
        wind_dir{f}=wind_dir{f}(startInd:stopInd);
        
        % nan bad data
        badInd=find(wind_u{f} >= theFillValue);
        wind_u{f}(badInd)=nan;
        badInd=find(wind_v{f} >= theFillValue);
        wind_v{f}(badInd)=nan;
        badInd=find(wind_spd{f} >= theFillValue);
        wind_spd{f}(badInd)=nan;
        badInd=find(wind_dir{f} >= theFillValue);
        wind_dir{f}(badInd)=nan;
        
        % for stations with hourly data, show
        % a wind stick every 3 hours 
        wind_uv{f}=nan*ones(size(jd{f}));
        if isequal(delta_t,1/24) 
            wind_uv{f}(1:3:end)=wind_u{f}(1:3:end) + i*wind_v{f}(1:3:end);
        else
            wind_uv{f}=wind_u{f} + i*wind_v{f};
        end
               
        % 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{f};
        ydata=[wind_dir{f} wind_spd{f} wind_v{f} wind_u{f} wind_uv{f}];
        panelOrd=[1:min(size(ydata))]; 
        figure
        han=timeplt(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');
            set(get(gca,'Ylabel'),'position',[julian([2003,10,8]) pos(2) pos(3)]) %align ylabels nicely
            linehan=findobj(gca,'Type','line','color','b');
            set(linehan,'color',col(f,:));
        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
            text(xt(tic),-20,xtl,'horizontalalignment','center',...
                                 'verticalalignment','top')
        end
        
        % give it a title and export to .pdf
        set(gcf,'PaperPositionMode','Auto')
        if strcmp(model,'COAMPS')
            suptitle(['Wind Observations from the ',model,' Model'])
            print(fullfile(plotDirectory,['winds_',model]),'-dpdf')
        else
            suptitle(['Wind Observations from ',inst_type,' ',inst_id]);
            print(fullfile(plotDirectory,['winds_',inst_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(model,'COAMPS')
            wind_u{f}=smart_interp(jd{f}, wind_u{f}, jd{1}, 3);
            wind_v{f}=smart_interp(jd{f}, wind_v{f}, jd{1}, 3);
            wind_spd{f}=smart_interp(jd{f}, wind_spd{f}, jd{1}, 3);
            wind_dir{f}=smart_interp(jd{f}, wind_dir{f}, jd{1}, 3);
        end
        
        wind_u_all(:,f)=wind_u{f};
        wind_v_all(:,f)=wind_v{f};
        wind_spd_all(:,f)=wind_spd{f};
        wind_dir_all(:,f)=wind_dir{f};

end            

% now call timeplt.m to create a combined plot. this plot will not show
% sticks.
xdata=jd{1};
ydata=[wind_dir_all wind_spd_all wind_v_all wind_u_all];
panelOrd=[repmat(1,1,min(size(ydata))/4),...
          repmat(2,1,min(size(ydata))/4),...
          repmat(3,1,min(size(ydata))/4),...
          repmat(4,1,min(size(ydata))/4)];
figure
han=timeplt(xdata,ydata,panelOrd);
clear global timeplt*
set(gcf,'ResizeFcn',[])
set(gcf,'Tag','')
LEG=legend(han(4),strvcat('NDBC Moored','Buoy 41013'),...
                  strvcat('NDBC C-Man','Station FPSN7'),...
                  strvcat('NOS Water Level','Station 8659897'),...
                  strvcat('NOS Water Level','Station 8661070'),...
                  strvcat('COAMPS Wind','Model'));
set(LEG,'Location','NorthOutside')   
set(LEG,'Orientation','Horizontal')
figCol=get(gcf,'color');
set(LEG,'color',figCol)
set(LEG,'edgecolor',figCol)
dolandscape
fixpaper
pos=get(LEG,'position');
set(LEG,'position',[pos(1) pos(2)+.05 pos(3) pos(4)])

% reformat xlims, ylims, ylabels
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');
    set(get(gca,'Ylabel'),'position',[julian([2003,10,8]) pos(2) pos(3)]) %align ylabels nicely
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
    text(xt(tic),-20,xtl,'horizontalalignment','center',...
        'verticalalignment','top')
end

% give it a title and export to .pdf
set(gcf,'PaperPositionMode','Auto')
suptitle([{'Wind Observations';'All Sources'}])
print(fullfile(plotDirectory,['winds_all']),'-dpdf')