%% PREDICTION METHOD (Operational Version)

%Intrusion Events & Near-hypoxic Intrusion Events

%R. Walt Deppe
%This version originally created October 31, 2013
%
%--modified on May 9, 2016 to change:
%   >"Hypoxic Intrusion Favorable Period (HIFP)"
%     to
%    "Near-hypoxic Intrusion Favorable Period (NIFP)"
%   >"Neap Tide Index (NTI)" to "Lunar Index (LI)"
%   >Window for NIFP's to 5-days on either
%    side of intrusion event (previously 8-days)

%% 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;

%% 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

%% 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

%% 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)

%% Required data inputs

display('*Be sure the following inputs are prepared and loaded:');
display('Tides: [NOAA_TE,NOAA_Time]');
display('Moon: [MoonPhase,MoonTime]');
display('Upwelling: [Raw_Upwell,Raw_Ut]');
display('River discharge: [Fraser_Discharge,Fraser_time]');
display(sprintf('*Intructions for inputs in code comments\n'));

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

% hourly time based on set of time series provided
day_start = ceil(max([Raw_Ut(1)+374 NOAA_Time(1) MoonTime(1) Fraser_time(1)]));
day_end = floor(min([Raw_Ut(end) NOAA_Time(end) MoonTime(end) Fraser_time(end)]));
Hour = (day_start:(1/24):day_end)';

display('*Upwelling data must begin at least 374 days before prediction start-time.');
display(['Prediction start-time based on inputs: ',datestr(day_start)]);
display(['Prediction end-time based on inputs: ',datestr(day_end)]);
tOK = input('Are these start and end times acceptable? (y/n) [default=y]\n','s');

if tOK=='n';
    display('New dates must be within this time period');
    sY = input('Enter new start year: ');
    sM = input('Enter new start month: ');
    sD = input('Enter new start day: ');
    eY = input('Enter new end year: ');
    eM = input('Enter new end month: ');
    eD = input('Enter new end day: ');    
    day_start = datenum(sY,sM,sD);
    day_end = datenum(eY,eM,eD);
    Hour = (day_start:(1/24):day_end)';
    display(['New prediction start-time based on inputs: ',datestr(day_start)]);
    display(['New prediction end-time based on inputs: ',datestr(day_end)]);
end

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

clear tOK eD eM eY sD sM sY 

%% 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();

%% LI (Lunar Index)

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

clear LIraw 
% **check interpolation**
% plot(MoonTime,MoonPhase);hold on;plot(Hour,LI,'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)

longHour = (Hour(1)-366:(1/24):Hour(end))';

% QC raw UI set
for j = 1:length(Raw_Upwell);
    if (Raw_Upwell(j,1) <= -9999);
        Raw_Upwell(j,1) = NaN;
    end
end
%find start and end points (with lag) to match Hour
[slH,start] = min(abs(Raw_Ut - (longHour(1)-7.25)));
[flH,final] = min(abs(Raw_Ut - (longHour(end)-7.25)));

%find Time start and end points (withOUT lag) to match Hour
[TslH,Tstart] = min(abs(Raw_Ut - longHour(1)));
[TflH,Tfinal] = min(abs(Raw_Ut - longHour(end)));

%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,longHour);

clear j start final Raw_upwell Raw_Ut slH flH TslH TflH
% ** check interpolation**
% plot(Raw_Utime,Raw_UI);
% hold on;
% plot(longHour,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(longHour,UI_hr);hold on;plot(longHour,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);
UItrend = NaN(length(UI_lpf40_integral),1);
startSeason = find(longHour==Hour(1));
endSeason = find(longHour==Hour(end));
for i = startSeason:endSeason;
UItrend(i) = nanmean(UI_lpf40_integral(find(longHour==longHour(i)-365.25):i));
% detrend by subtracting the mean of all data points in the previous year
% for each point
    UPI(i) = ...
        UI_lpf40_integral(i)-UItrend(i);
end
% **check conditonal integration**
%  plot(longHour,UI_lpf40);hold on;
%  plot(longHour, UI_lpf40_integral,'k');
%  plot(longHour, UItrend,'c');
%  plot(longHour, UPI,'g');
%  hold off; grid on; dynamicDateTicks();
 
 UPIlong = UPI;
 UPI = interp1(longHour,UPIlong,Hour);
 
clear Reals start t i UI_lpf40_integral UI_lpf40_sign UItrend poly2...
        startSeason endSeason


%% FRDI (Fraser River Discharge Index)


% 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;

% match residual current time-scale [subtidal signal]
fdi66 = pl66tn(Fdis_hr,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-691.56)/(11748-691.56);
%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.14 - (0.00127*UPI) - (2.77*FRDI);

%% Intrusion Event Predictions

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

%% Near-hypoxic Intrusion Event Predictions

NHIEconditions = find(LI>=0.8 & DII>=0.72 & DOAP<=6.2);

%% NIFP (Near-hypoxic Intrusion Favorable Period)

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

clear i

%% Plot Intrusion Predictions and NIFP's

figure;

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;

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

b=subplot(2,1,2);
area(Hour,NHypIntP,...
    'FaceColor','c','EdgeColor','none');
hold on;
plot(Hour,DOAP,'--b');
area(Hour,NHypIntP,...
    'FaceColor','none','EdgeColor','k','LineStyle','--');
plot(Hour,3*ones(size(Hour)),'k');
plot(Hour,10*ones(size(Hour)),'k');
title('Near-hypoxic Intrusion Event Favorable Periods','FontSize',22);
ylabel('Dissolved Oxygen (mg/L)','FontSize',20);
legend('Near-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 i IntP NHypIntP

