%% PREDICTION METHOD (Operational Version)

%Intrusion Events & Hypoxic Intrusion Events

%R. Walt Deppe
%November 16, 2013

%% Import NOAA Tideal Elevation Data

% from 9444900 Port Townsend, WA
% [http://tidesandcurrents.noaa.gov/waterlevels.html?id=9444900]

% download data by hour in GMT in reference to MLLW
% in year-length chunks and combine
% (excel time + 695422 = matlab time)

% NOAA_Tides = uiimport('NOAA_TidalElevation.csv');
% NOAA_Time = NOAA_Tides.NOAA_Time_Matlab;
% NOAA_TE = NOAA_Tides.NOAA_TE_hr;


%% Choose time period for analysis / Set up hourly time steps

% hourly time from (05 Aug 2009) to (02 Oct 2013)
Hour = (733990:(1/24):735509)';

% Be sure that [pl66tn.m] and [dynamicDateTicks.m],
% are in the current folder

%% TEI (Tidal Elevation Index)

TEI_noaa = interp1(NOAA_Time,NOAA_TE,Hour);

% **check interpolation**
% plot(NOAA_Time,NOAA_TE,'k');
% hold on; plot(Hour,TEI_noaa,'--r');
% hold off;
% dynamicDateTicks();

%% NTI (Neap Tide Index)

% Import MoonPhase data from the Astronomical Applications Department of
% the US Naval Observatory
% (MoonTime = day from data set)
% (MoonPhase = moon phase, 1=full, 0=new)
% http://aa.usno.navy.mil/data/docs/MoonFraction.php

% Set: 1 = Peak Neap , 0 = Peak Spring
NTIraw = (-2*abs(MoonPhase-0.5))+1;
%interpolate to hourly time scale
NTI = interp1(MoonTime,NTIraw,Hour);

clear NTIraw 
% **check interpolation**
% plot(MoonTime,MoonPhase);hold on;plot(Hour,NTI,'r');hold off;

%% DII (Diurnal Inequality Index)

% alsolute value of gradient of Tidal Amplitude Index [tidal time-scale]
tgi_noaa = abs(gradient(TEI_noaa));
% match residual current time-scale [subtidal signal]
tgi66_noaa = pl66tn(tgi_noaa,1,40);
% remove parts of filtered signal biased by NaN values
tgi40_noaa = tgi66_noaa;
for i = 1:length(tgi66_noaa);
    if (i<41 || i>(length(tgi66_noaa)-40));
        tgi40_noaa(i) = NaN;
    elseif(isnan(tgi66_noaa(i)));
        tgi40_noaa((i-40):(i+40),1) = NaN(81,1);
    end
end
% normalize/reverse sign
DII = (-(tgi40_noaa-max(tgi40_noaa)))/...
            max(-(tgi40_noaa-max(tgi40_noaa)));
% 1 = Max. Diurnal Inequality ; 0 = Diurnal Equality

clear tgi tgi66 tgi40 i tgi_noaa tgi66_noaa tgi40_noaa
% **check against TEI**
% plot(Hour,TEI_noaa,'b');hold on;
% plot(Hour,DII,'r');
% hold off; dynamicDateTicks();

%% Lagged Upwelling Index (UI - 7.25 days)

% Import Raw_Upwell data from Upwelling Index from PFEL (NOAA)
% for site 48 N, 125 W
% (Raw_Ut = day from data set)
% (Raw_Upwell = upwelling index observations)
% http://www.pfeg.noaa.gov/products/PFEL/modeled/indices/upwelling/NA/data_download.html
% **Raw_Ut range currently used: [733980 735415] 26-Jul-2009 to 30-Jun-2013

% QC raw UI set
for j = 1:length(Raw_Upwell);
    if (Raw_Upwell(j,1) <= -9999);
        Raw_Upwell(j,1) = NaN;
    end
end
%due to missing values in june 2013 record
for j = 5628:5744;
    if (isnan(Raw_Upwell(j,1)));
        Raw_Upwell(j,1) = (Raw_Upwell(j-1,1)+Raw_Upwell(j+1,1))/2;
    end
end
%find start and end points (with lag) to match Hour
start = find(Raw_Ut == Hour(1)-7.25);
final = find(Raw_Ut == Hour(length(Hour))-7.25);
%find Time start and end points (withOUT lag) to match Hour
Tstart = find(Raw_Ut == Hour(1));
Tfinal = find(Raw_Ut == Hour(length(Hour)));
%select Lagged Raw Upwelling section of data corresponding to time period
Raw_Utime = Raw_Ut(Tstart:Tfinal,1);
Raw_UI = Raw_Upwell(start:final,1);
% interpolate to hourly
UI_hr = interp1(Raw_Utime,Raw_UI,Hour);

clear j start final Raw_upwell Raw_Ut
% ** check interpolation**
% plot(Raw_Utime,Raw_UI);
% hold on;
% plot(Hour,UI_hr,'r');
% hold off; grid on;
% dynamicDateTicks();

%% UPI (Upwelling Persistance Index)

% match residual current time-scale [subtidal signal]
ui66 = pl66tn(UI_hr,1,40);
% remove parts of filtered signal biased by NaN values
UI_lpf40 = ui66;
for i = 1:length(ui66);
    if (i<41 || i>(length(ui66)-40));
        UI_lpf40(i) = NaN;
    elseif(isnan(ui66(i)));
        UI_lpf40((i-40):(i+40),1) = NaN(81,1);
    end
end

clear ui66 i
% **check filter**
% plot(Hour,UI_hr);hold on;plot(Hour,UI_lpf40,'r');hold off;

% conditional integration of UI_lpf40
UI_lpf40_integral = NaN(length(UI_lpf40),1);
UI_lpf40_sign = NaN(length(UI_lpf40),1);
for i = 1:length(UI_lpf40_sign);
    if (isfinite(UI_lpf40(i)));
        if (UI_lpf40(i)>0 );
            UI_lpf40_sign(i) = 1;
        elseif (UI_lpf40(i)<0);
               UI_lpf40_sign(i) = -1;
        else 
                UI_lpf40_sign(i) = 0;
        end
    end
end
% begin integration following NaN's at beginning
Reals = find(isfinite(UI_lpf40_sign));
start = Reals(1,1);
for t=start:length(UI_lpf40_integral);
    if (t>start);
        if (isfinite(UI_lpf40_sign(t)));
            UI_lpf40_integral(t) = UI_lpf40_sign(t)+UI_lpf40_integral(t-1);
        else
            UI_lpf40_integral(t) = UI_lpf40_integral(t-1);
        end
    else
        UI_lpf40_integral(t) = UI_lpf40_sign(t);
    end
end
% detrend to create UPI
UPI = NaN(length(UI_lpf40_integral),1);
%determine 2nd order polynomial coefficients
%using time period between the First and Last
%01January in the time series (for seasonal trend consistency)
startSeason = find(Hour==datenum(2010,01,01));
endSeason = find(Hour==datenum(2013,01,01));
poly2 = polyfit(Hour(startSeason:endSeason),...
    UI_lpf40_integral(startSeason:endSeason),2);
%fit 2nd order polynomial
UItrend = polyval(poly2,Hour(start:length(UI_lpf40_integral)));
% detrend
UPI(start:length(UI_lpf40_integral)) = ...
    UI_lpf40_integral(start:length(UI_lpf40_integral))-UItrend;

% **check conditonal integration**
%  plot(Hour,UI_lpf40);hold on;
%  plot(Hour, UI_lpf40_integral,'k');
%  plot(Hour(start:length(UI_lpf40_integral)), UItrend,'c');
%  plot(Hour, UPI,'g');
%  hold off; grid on; dynamicDateTicks();
 
clear Reals start t i UI_lpf40_integral UI_lpf40_sign UItrend poly2...
        startSeason endSeason

%% FRDI (Fraser River Discharge Index)

% import Fraser Data from Hope, BC, Canada
% downloaded from the Water Survey of Canada Website
% *must download recent and older data seperately and combine
% http://www.wateroffice.ec.gc.ca/graph/graph_e.html?stn=08MF005
% http://www.wsc.ec.gc.ca/applications/H2O/graph-eng.cfm?station=08MF005&report=daily&year=2011
% rename as (Fraser_Discharge) and (Fraser_time)

% switch to UTC time
UTC_time_fraser = Fraser_time + (8/24);
% interpolate to hourly
Fdis_hr = interp1(UTC_time_fraser,Fraser_Discharge,Hour);

% **check interpolation**
% scatter(UTC_time_fraser,Fraser_Discharge,'x');
% hold on;
% plot(Hour,Fdis_hr,'r');
% hold off;

% discharge index (normalization)
dMax = max(Fdis_hr);
dMin = min(Fdis_hr);
FDI_prefilter = (Fdis_hr - dMin)/(dMax-dMin);
% match residual current time-scale [subtidal signal]
fdi66 = pl66tn(FDI_prefilter,1,40);
% remove parts of filtered signal biased by NaN values
fdi40 = fdi66;
for i = 1:length(fdi66);
    if (i<41 || i>(length(fdi66)-40));
        fdi40(i) = NaN;
    elseif(isnan(fdi66(i)));
        fdi40((i-40):(i+40),1) = NaN(81,1);
    end
end
% normalize
FRDI = (fdi40-min(fdi40))/max(fdi40-min(fdi40));

% **check filter**
% plot(Hour,FDI_prefilter);hold on;plot(Hour,FRDI,'r');hold off;

clear fdi40 fdi66 i dMax dMin UTC_time_fraser Fraser_time...
    Fraser_Discharge FDI_prefilter

%% DOAP (Dissolved Oxygen Availibility Prediction)

% use equation derived from regression
DOAP = 7.0343 - (0.0013*UPI) - (2.0535*FRDI);

%% Intrusion Event Predictions

IEconditions = find(NTI>=0.8 & DII>=0.72);

%% Hypoxic Intrusion Event Predictions

HIEconditions = find(NTI>=0.8 & DII>=0.72 & DOAP<=6.241);

%% HIFP (Hypoxic Intrusion Favorable Period)

HIFP = NaN(length(Hour),1);
% HIFP is the hour vector that is set to NAN if it is not within 8
% days of a hypoxic intrusion event prediction point
for i = 1:length(HIFP);
    if (i>192 && i<(length(HIFP)-192));
        if (any(NTI(((i-192):(i+192)),1)>=0.8 &...
                DII(((i-192):(i+192)),1)>=0.72 &...
                DOAP(((i-192):(i+192)),1)<=6.241));
            HIFP(i) = Hour(i);
        end
    elseif (i<=192);
        if (any(NTI((1:(i+192)),1)>=0.8 &...
                DII((1:(i+192)),1)>=0.72 &...
                DOAP((1:(i+192)),1)<=6.241));
            HIFP(i) = Hour(i);
        end
    else
        if (any(NTI(((i-192):length(HIFP)),1)>=0.8 &...
                DII(((i-192):length(HIFP)),1)>=0.72 &...
                DOAP(((i-192):length(HIFP)),1)<=6.241));
            HIFP(i) = Hour(i);
        end
    end
end

clear i

%% Plot Intrusion Predictions and HIFP's

IntP = NaN(size(Hour));
for i = 1:length(IntP);
    if (any(i==IEconditions));
        IntP(i) = 1;
    else
        IntP(i) = 0;
    end
end

a=subplot(2,1,1);
area(Hour,IntP,...
    'FaceColor','g','EdgeColor','none');
hold on;
area(Hour,IntP,...
    'FaceColor','none','EdgeColor','g','LineStyle','-');
plot(Hour,zeros(size(Hour)),'k');
plot(Hour,1*ones(size(Hour)),'k');
title('Intrusion Event Predictions','FontSize',22);
dynamicDateTicks();
set(gca,'YDir','normal','FontSize',18,'Layer','top');
grid on;
legend('Intrusion Event Prediction');
hold off;

HypIntP = NaN(size(Hour));
for i = 1:length(HypIntP);
    if (isfinite(HIFP(i)));
        HypIntP(i) = 10;
    else
        HypIntP(i) = 0;
    end
end

b=subplot(2,1,2);
area(Hour,HypIntP,...
    'FaceColor','c','EdgeColor','none');
hold on;
plot(Hour,DOAP,'--b');
area(Hour,HypIntP,...
    'FaceColor','none','EdgeColor','k','LineStyle','--');
plot(Hour,3*ones(size(Hour)),'k');
plot(Hour,10*ones(size(Hour)),'k');
title('Hypoxic Intrusion Event Favorable Periods','FontSize',22);
ylabel('Dissolved Oxygen (mg/L)','FontSize',20);
legend('Hypoxic Intrusion Event Favorable',...
    'Dissolved Oxygen Availibility Prediction (seasonal)');
hold off;
grid on; set(gca,'YDir','normal','FontSize',18,'Layer','top');
ylim([3 10]);

linkaxes([a b],'x');
dynamicDateTicks([a b], 'link');
xlim([Hour(1,1),Hour(length(Hour),1)]);

set(gcf, 'PaperPosition', [0.5 2 18 12]);

clear a b IntP HypIntP

