% Use of this program is described in:
%
% Armstrong, Brandy N., Warner, John C., List, Jeff H., Martini, 
% Marinna A., Montgomery, Ellyn T., Voulgaris, George 
% 2014, Coastal Change Processes Project data report for oceanographic 
% observations near Fire Island, January to April 
% 2012 : U.S. Geological Survey Open-File Report 2014-xxxx, available 
% online at http://pubs.usgs.gov/of/2014/xxxx/.
%
% Program written in Matlab v8.0.0.755
% Program ran on PC with Windows 7 64bit

%
% "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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% B. Armstrong 6/20/14, Raw met station data is processed and saved to netCDF files.
%                       It is assumed that the netCDF files are stored in the current  
%                       directory. NOAA tide station data is acquired using
%                       THREDDS.
%                       
% Dependencies:  
%       replace.m
%       tight_subplot.m
%       gregorian.m
%       julian.m
%       draw.m

%NOAA tide stations information
%http://dods.ndbc.noaa.gov:8080/thredds/dodsC/data/stdmet/idid#/idid#hyyyy.nc
%http://dods.ndbc.noaa.gov:8080/thredds/dodsC/data/stdmet/

plotDirectory = 'I:\FireIsland\FI_12_DataDVD\PLOTS';

%% NDBC station ids
stat_id={'44065','44025'};

for n=1:length(stat_id)
    %get NDBC data
    eval(['nc=''http://dods.ndbc.noaa.gov:8080/thredds/dodsC/data/stdmet/',stat_id{n},'/',stat_id{n},'h2012.nc'';']);
    eval(['b',stat_id{n},'.hsig=ncread(nc,''wave_height'');']);
    eval(['b',stat_id{n},'.windspd=ncread(nc,''wind_spd'');']);
    eval(['b',stat_id{n},'.winddir=ncread(nc,''wind_dir'');']);
    eval(['b',stat_id{n},'.gust=ncread(nc,''gust'');']);
    eval(['b',stat_id{n},'.airp=ncread(nc,''air_pressure'');']);
    eval(['b',stat_id{n},'.airt=((double(ncread(nc,''air_temperature''))));']);
%     eval(['b',stat_id{n},'.airt=((double(ncread(nc,''air_temperature'')))*1.8)+32;']);
    eval(['b',stat_id{n},'.lon=ncread(nc,''longitude'');']);
    eval(['b',stat_id{n},'.lat=ncread(nc,''latitude'');']);
    eval(['b',stat_id{n},'.ndbctime=ncread(nc,''time'');']);
    eval(['b',stat_id{n},'.ndbct=datenum(gregorian(double(b',stat_id{n},'.ndbctime)/86400+julian(1970,1,1,0)));']);
end
    clf
%%
startTime=datenum(gregorian(julian(2012,1,27,00)));                   
stopTime=datenum(gregorian(julian(2012,05,19,00)));         
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %past start time and stop time
xt=[];
for yy=2012:2012
    for mm=1:12
        tic=datenum(gregorian(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
ylab={ {'Air Temperature','(\circC)'}...   %plot ylabels
    {'Barometric Pressure','(hPa)'},...
    {'Wind Speed','(m/s)'},...
    {'Wind Gust','(m/s)'},...
    {'Wind Direction','(degrees)'}};
%%
met_st={'9191met'};
met9191='9191met-prfix1.cdf';

buoyvariable={'airt' 'airp' 'windspd' 'gust' 'winddir'};
metvariable={'T_21' 'BP_915' 'WS_401' 'WG_402' 'WD_410'};
titles={'Air Temperature (deg C)' 'Barometric Pressure (hPa)' 'Wind Speed (m/s)' 'Wind Gust (m/s)' 'Wind Direction degrees'};
ylims=[ -10 40;970 1070;0 30;0 30;0 360;];
 clf
set(gcf,'position',[100 25 850 1100])
h = tight_subplot(5,1,[.04 .3],[.1 .1],[.1 .1]) ;

for p=1:length(metvariable)
    h(p)=subplot(length(metvariable),1,p);hold on;
    var2plot=squeeze(ncread(met9191,metvariable{p}));
    cols={'r','b','k','m'};
    for n=1:length(stat_id) %plot ndbc stations
        if strmatch(buoyvariable{p},'winddir')
        eval(['plot(h(p),b',stat_id{n},'.ndbct,squeeze(b',stat_id{n},'.',buoyvariable{p},'),''',cols{n},'.'');']);
        else
        eval(['plot(h(p),b',stat_id{n},'.ndbct,squeeze(b',stat_id{n},'.',buoyvariable{p},'),''',cols{n},''');']);
        end
        hold on;
    end
    if strmatch(metvariable{p},'WS_401')
        var2plot(var2plot>20)=nan;
    elseif  strmatch(metvariable{p},'WG_402')
        var2plot(var2plot>30)=nan;
    end
    if strmatch(metvariable{p},'WD_410')
        plot(datenum(gregorian(double(ncread(met9191,'time')) + (double(ncread(met9191,'time2'))/3600/1000/24))),var2plot-13.5,'g.');
        legend([stat_id,met_st],'Location','NorthEast');
    else
        plot(datenum(gregorian(double(ncread(met9191,'time')) + (double(ncread(met9191,'time2'))/3600/1000/24))),var2plot,'g');
        legend([stat_id,met_st],'Location','NorthEast','Orientation','Horizontal');
    end
    title(titles{p});datetick('x',20);ylim([ylims(p,:)]);
    xlim([datenum(2012,01,01) datenum(2012,06,01)]);
end
for hind=1:length(h)
        set(h(hind),'xlim',xlims,'xtick',xt)
        set(h(hind),'xticklabel',{[datestr(xt,'mmm dd yyyy')]})
        ylabhan=get(h(hind),'Ylabel');
        set(ylabhan,'String',ylab{hind})
        pos=get(ylabhan,'position');
        
end
set(gcf,'PaperPositionMode','Auto','PaperOrientation','portrait')
title(h(1),{['Fire island 2012: Meteorological Buoy and NDBC Stations Comparison']})
print(fullfile(plotDirectory,['FI12_met_sites']),'-dpdf');          