% Observations from M2 Buoy : Station ID:62091 Location:53.4800°N 05.4250°W 
% 
% More info at: http://www.marine.ie/home/publicationsdata/data/buoys/
% 
    clear ; pack;
    cd 'E:\Jon_work\PROJECTS\CURRENT PROJECTS\MICORE FP7 START June 2008\Data and data management\Irish M2 buoy'
    load Irish_M2_buoy
%   XX = char(textdata(2:64135,1));
%   time = datenum(XX(:,1:20));	
    AtmPressure= data(:,1);	WindDirection= data(:,2);	WindSpeed= data(:,3);	Gust= data(:,4);	
    DryBulb= data(:,5);	DewPoint= data(:,6);	SeaTemp= data(:,7);	WavePeriod= data(:,8);	WaveHeight= data(:,9);	
    RelativeHumidity= data(:,10);

    z1 = find(WaveHeight < 0);WaveHeight(z1) = NaN; 
    z2 = find(WavePeriod < 0);WavePeriod(z2) = NaN;
    z3 = find(WindDirection < 0);WindDirection(z3) = NaN;
    z4 = find(WindSpeed < 0);WindSpeed(z4) = NaN;
    z5 = find(Gust < 0);Gust(z5) = NaN;
    z6 = find(Gust > 100);Gust(z6) = NaN;
    z7 = find(AtmPressure < 0);AtmPressure(z7) = NaN;

    figure(1)
        subplot(5,1,1)
            plot(time,WaveHeight(1:64134),'.b');datetick('x',12);
            axis([min(time) max(time) 0 10]); xlabel('Date'); ylabel('m');title('Hs');set(gca,'fontsize',10);
        subplot(5,1,2)
            plot(time,WavePeriod(1:64134),'.k');datetick('x',12);
            axis([min(time) max(time) 0 16]); xlabel('Date'); ylabel('sec');title('Tp');set(gca,'fontsize',10);
        subplot(5,1,3)
            plot(time,WindSpeed(1:64134),'.r');hold on; %plot(time,Gust(1:64134),'.r');datetick('x',12);
            axis([min(time) max(time) 0 100]); xlabel('Date'); ylabel('m/s');title('Wind');set(gca,'fontsize',10);%legend('Mean','Gust');
        subplot(5,1,4)
            plot(time,WindDirection(1:64134),'.b');datetick('x',12);
            axis([min(time) max(time) -100 400]); xlabel('Date'); ylabel('deg');title('Wind direction');set(gca,'fontsize',10);
        subplot(5,1,5)
            plot(time,AtmPressure(1:64134),'.k');datetick('x',12);
            axis([min(time) max(time) 900 1100]); xlabel('Date'); ylabel('mb');title('Atmospheric pressure');set(gca,'fontsize',10);

    figure (2)
        subplot(3,1,1)
            plot(WindSpeed,WaveHeight,'+k');xlabel('WindSpeed (m/s)'); ylabel('WaveHeight (m)');set(gca,'fontsize',10);hold on;
            
zz1 = isnan(WindSpeed); count = 1; for index = 1:length(zz1); 
if zz1(index) ==0;WindSpeed_edit(count) = WindSpeed(index);  count = count + 1; end; 
if zz1(index) ==1;WindSpeed_edit(count) = 0;  count = count + 1; end; end;
zz2 = isnan(WaveHeight);count = 1; for index = 1:length(zz2);
if zz2(index) ==0; WaveHeight_edit(count) = WaveHeight(index); count = count + 1; end; if zz2(index) ==1; WaveHeight_edit(count) = 0; count = count + 1; end;end;
            xx = (0: 1: 50)';
            p = polyfit(WindSpeed_edit,WaveHeight_edit,2);
            fit1 = polyval(p,xx);
            plot(xx,fit1,'r-')

        subplot(3,1,2)
            plot(Gust,WaveHeight,'+k');xlabel('Gust Speed (m/s)'); ylabel('WaveHeight (m)');set(gca,'fontsize',10); hold on;
zz3 = isnan(Gust); count = 1; for index = 1:length(zz3); 
if zz3(index) ==0;Gust_edit(count) = Gust(index);  count = count + 1; end; if zz3(index) ==1;Gust_edit(count) = 0;  count = count + 1; end; end;
            xx = (0: 1: 80)';
            p = polyfit(Gust_edit,WaveHeight_edit,2);
            fit2 = polyval(p,xx);
            plot(xx,fit2,'r-')
        subplot(3,1,3)
            plot(AtmPressure,WaveHeight,'+k');xlabel('AtmPressure (mb)'); ylabel('m');title('Hs');set(gca,'fontsize',10); hold on;
zz4 = isnan(AtmPressure); count = 1; for index = 1:length(zz4); 
if zz4(index) ==0;AtmPressure_edit(count) = AtmPressure(index);  count = count + 1; end; if zz4(index) ==1;AtmPressure_edit(count) = 0;  count = count + 1; end; end;
            xx = (950: 1: 1050)';
            p = polyfit(AtmPressure_edit,WaveHeight_edit,2);
            fit3 = polyval(p,xx);
            plot(xx,fit3,'r-')
AVE=500;            
waveHT_trend = movingaverage(WaveHeight_edit,AVE);time_trend = movingaverage(time,AVE);
R1 = corrcoef(time_trend,waveHT_trend); correlation1 = R1(2)
figure(3); subplot(3,1,1); plot(time_trend,waveHT_trend,'+r');hold on;datetick('x',12);title(num2str(correlation1));
xx = (min(time_trend):1:max(time_trend));
            p = polyfit(time_trend',waveHT_trend,1);
            fit4 = polyval(p,xx);
            plot(xx,fit4,'k-')

windSP_trend = movingaverage(WindSpeed_edit,AVE);
R2 = corrcoef(time_trend,windSP_trend); correlation2 = R2(2)
figure(3); subplot(3,1,2);plot(time_trend,windSP_trend,'+r'); hold on;datetick('x',12);title(num2str(correlation2));
            p = polyfit(time_trend',windSP_trend,1);
            fit5 = polyval(p,xx);
            plot(xx,fit5,'k-')
            
AtmPressure_trend = movingaverage(AtmPressure_edit,AVE);
R3 = corrcoef(time_trend,AtmPressure_trend); correlation3 = R3(2)
figure(3); subplot(3,1,3);plot(time_trend,AtmPressure_trend,'+r'); hold on;datetick('x',12);title(num2str(correlation3));
            p = polyfit(time_trend',AtmPressure_trend,1);
            fit6 = polyval(p,xx);
            plot(xx,fit6,'k-')