%% SUB-TIDAL DO PREDICTION

%Using Victoria Clipper T-S observations

%R. Walt Deppe
%November 16, 2013

%% Import Data Sets

VC_30Jun2012 = uiimport('2012-06-30_vc_csv_export.mat');
VC_31Jul2012 = uiimport('2012-07-31_vc_csv_export.mat');
VC_31Aug2012 = uiimport('2012-08-31_vc_csv_export.mat');
VC_30Sep2012 = uiimport('2012-09-30_vc_csv_export.mat');
VC_31Oct2012 = uiimport('2012-10-31_vc_csv_export.mat');
VC_30Nov2012 = uiimport('2012-11-30_vc_csv_export.mat');
VC_31Dec2012 = uiimport('2012-12-31_vc_csv_export.mat');

%% TSNH sets

Pars_TSNH = [1 5 6 8 9 10 11 12];

VC_TSNH_raw_26 = VC_30Jun2012.TSNH;
VC_TSNH_raw_27 = VC_31Jul2012.TSNH;
VC_TSNH_raw_28 = VC_31Aug2012.TSNH;
VC_TSNH_raw_29 = VC_30Sep2012.TSNH;
VC_TSNH_raw_30 = VC_31Oct2012.TSNH;
VC_TSNH_raw_31 = VC_30Nov2012.TSNH;
VC_TSNH_raw_32 = VC_31Dec2012.TSNH;

VC_TSNH_26 = VC_TSNH_raw_26(1:size(VC_TSNH_raw_26,1),Pars_TSNH);
VC_TSNH_27 = VC_TSNH_raw_27(1:size(VC_TSNH_raw_27,1),Pars_TSNH);
VC_TSNH_28 = VC_TSNH_raw_28(1:size(VC_TSNH_raw_28,1),Pars_TSNH);
VC_TSNH_29 = VC_TSNH_raw_29(1:size(VC_TSNH_raw_29,1),Pars_TSNH);
VC_TSNH_30 = VC_TSNH_raw_30(1:size(VC_TSNH_raw_30,1),Pars_TSNH);
VC_TSNH_31 = VC_TSNH_raw_31(1:size(VC_TSNH_raw_31,1),Pars_TSNH);
VC_TSNH_32 = VC_TSNH_raw_32(1:size(VC_TSNH_raw_32,1),Pars_TSNH);

clear VC_TSNH_raw_26 VC_TSNH_raw_27 VC_TSNH_raw_28 VC_TSNH_raw_29...
   VC_TSNH_raw_30 VC_TSNH_raw_31 VC_TSNH_raw_32 VC_30Jun2012...
   VC_31Jul2012 VC_31Aug2012 VC_30Sep2012 VC_31Oct2012...
   VC_30Nov2012 VC_31Dec2012

%% combine TSNH

%Time
col = 1;
VC_TSNH_Time_sets = vertcat(...
    vertcat(VC_TSNH_26(1:size(VC_TSNH_26,1),col)),...
    vertcat(VC_TSNH_27(1:size(VC_TSNH_27,1),col)),...
    vertcat(VC_TSNH_28(1:size(VC_TSNH_28,1),col)),...
    vertcat(VC_TSNH_29(1:size(VC_TSNH_29,1),col)),...
    vertcat(VC_TSNH_30(1:size(VC_TSNH_30,1),col)),...
    vertcat(VC_TSNH_31(1:size(VC_TSNH_31,1),col)),...
    vertcat(VC_TSNH_32(1:size(VC_TSNH_32,1),col)));
%Temperature
col = 2;
VC_TSNH_Temp_sets = vertcat(...
    vertcat(VC_TSNH_26(1:size(VC_TSNH_26,1),col)),...
    vertcat(VC_TSNH_27(1:size(VC_TSNH_27,1),col)),...
    vertcat(VC_TSNH_28(1:size(VC_TSNH_28,1),col)),...
    vertcat(VC_TSNH_29(1:size(VC_TSNH_29,1),col)),...
    vertcat(VC_TSNH_30(1:size(VC_TSNH_30,1),col)),...
    vertcat(VC_TSNH_31(1:size(VC_TSNH_31,1),col)),...
    vertcat(VC_TSNH_32(1:size(VC_TSNH_32,1),col)));
%Salinity
col = 3;
VC_TSNH_Salinity_sets = vertcat(...
    vertcat(VC_TSNH_26(1:size(VC_TSNH_26,1),col)),...
    vertcat(VC_TSNH_27(1:size(VC_TSNH_27,1),col)),...
    vertcat(VC_TSNH_28(1:size(VC_TSNH_28,1),col)),...
    vertcat(VC_TSNH_29(1:size(VC_TSNH_29,1),col)),...
    vertcat(VC_TSNH_30(1:size(VC_TSNH_30,1),col)),...
    vertcat(VC_TSNH_31(1:size(VC_TSNH_31,1),col)),...
    vertcat(VC_TSNH_32(1:size(VC_TSNH_32,1),col)));
%Time_GPS
col = 4;
VC_TSNH_Time_GPS_sets = vertcat(...
    vertcat(VC_TSNH_26(1:size(VC_TSNH_26,1),col)),...
    vertcat(VC_TSNH_27(1:size(VC_TSNH_27,1),col)),...
    vertcat(VC_TSNH_28(1:size(VC_TSNH_28,1),col)),...
    vertcat(VC_TSNH_29(1:size(VC_TSNH_29,1),col)),...
    vertcat(VC_TSNH_30(1:size(VC_TSNH_30,1),col)),...
    vertcat(VC_TSNH_31(1:size(VC_TSNH_31,1),col)),...
    vertcat(VC_TSNH_32(1:size(VC_TSNH_32,1),col)));
%Lattitude
col = 5;
VC_TSNH_Latt_sets = vertcat(...
    vertcat(VC_TSNH_26(1:size(VC_TSNH_26,1),col)),...
    vertcat(VC_TSNH_27(1:size(VC_TSNH_27,1),col)),...
    vertcat(VC_TSNH_28(1:size(VC_TSNH_28,1),col)),...
    vertcat(VC_TSNH_29(1:size(VC_TSNH_29,1),col)),...
    vertcat(VC_TSNH_30(1:size(VC_TSNH_30,1),col)),...
    vertcat(VC_TSNH_31(1:size(VC_TSNH_31,1),col)),...
    vertcat(VC_TSNH_32(1:size(VC_TSNH_32,1),col)));
%Longitude
col = 6;
VC_TSNH_Long_sets = vertcat(...
    vertcat(VC_TSNH_26(1:size(VC_TSNH_26,1),col)),...
    vertcat(VC_TSNH_27(1:size(VC_TSNH_27,1),col)),...
    vertcat(VC_TSNH_28(1:size(VC_TSNH_28,1),col)),...
    vertcat(VC_TSNH_29(1:size(VC_TSNH_29,1),col)),...
    vertcat(VC_TSNH_30(1:size(VC_TSNH_30,1),col)),...
    vertcat(VC_TSNH_31(1:size(VC_TSNH_31,1),col)),...
    vertcat(VC_TSNH_32(1:size(VC_TSNH_32,1),col)));
%Speed
col = 7;
VC_TSNH_Speed_sets = vertcat(...
    vertcat(VC_TSNH_26(1:size(VC_TSNH_26,1),col)),...
    vertcat(VC_TSNH_27(1:size(VC_TSNH_27,1),col)),...
    vertcat(VC_TSNH_28(1:size(VC_TSNH_28,1),col)),...
    vertcat(VC_TSNH_29(1:size(VC_TSNH_29,1),col)),...
    vertcat(VC_TSNH_30(1:size(VC_TSNH_30,1),col)),...
    vertcat(VC_TSNH_31(1:size(VC_TSNH_31,1),col)),...
    vertcat(VC_TSNH_32(1:size(VC_TSNH_32,1),col)));
%Direction
col = 8;
VC_TSNH_Dir_sets = vertcat(...
    vertcat(VC_TSNH_26(1:size(VC_TSNH_26,1),col)),...
    vertcat(VC_TSNH_27(1:size(VC_TSNH_27,1),col)),...
    vertcat(VC_TSNH_28(1:size(VC_TSNH_28,1),col)),...
    vertcat(VC_TSNH_29(1:size(VC_TSNH_29,1),col)),...
    vertcat(VC_TSNH_30(1:size(VC_TSNH_30,1),col)),...
    vertcat(VC_TSNH_31(1:size(VC_TSNH_31,1),col)),...
    vertcat(VC_TSNH_32(1:size(VC_TSNH_32,1),col)));

clear col VC_TSNH_26 VC_TSNH_27 VC_TSNH_28 VC_TSNH_29...
    VC_TSNH_30 VC_TSNH_31 VC_TSNH_32

%remove blank cells
cut = [64:72];
VC_TSNH_Time_sets(cut) = [];
VC_TSNH_Temp_sets(cut) = [];
VC_TSNH_Salinity_sets(cut) = [];
VC_TSNH_Time_GPS_sets(cut) = [];
VC_TSNH_Latt_sets(cut) = [];
VC_TSNH_Long_sets(cut) = [];
VC_TSNH_Speed_sets(cut) = [];
VC_TSNH_Dir_sets(cut) = [];

clear cut

%% fix length-match problems

for i = 1:size(VC_TSNH_Time_sets,1);
    if (size(VC_TSNH_Time_sets{i,1,:},1)>...
            size(VC_TSNH_Temp_sets{i,1,:},1));
        dsize = size(VC_TSNH_Time_sets{i,1,:},1) - ...
                   size(VC_TSNH_Temp_sets{i,1,:},1);         
        Tempfix = VC_TSNH_Temp_sets{i,1,:};
        Tempfix(length(Tempfix)+1:length(Tempfix)+dsize,1) = NaN(dsize,1);
        VC_TSNH_Temp_sets{i,1,:} = Tempfix;
    end
    if (size(VC_TSNH_Time_sets{i,1,:},1)>...
            size(VC_TSNH_Salinity_sets{i,1,:},1));
        dsize = size(VC_TSNH_Time_sets{i,1,:},1) - ...
                   size(VC_TSNH_Salinity_sets{i,1,:},1);         
        Salfix = VC_TSNH_Salinity_sets{i,1,:};
        Salfix(length(Salfix)+1:length(Salfix)+dsize,1) = NaN(dsize,1);
        VC_TSNH_Salinity_sets{i,1,:} = Salfix;
    end
end

clear i Tempfix Salfix dsize

%% ready sets for interpolation, make time set

Step = 1;   %set time step in minutes

N = size(VC_TSNH_Time_sets,1);

%remove foul time value-------------
for i = 1:N;
    SetFix = VC_TSNH_Time_sets{i,1,:};
    for j = 1:length(SetFix);
        if (SetFix(j)<datenum(2012,05,31) ...
                || SetFix(j)>datenum(2013,01,02));
            SetFix(j) = NaN;
        end
    end
    VC_TSNH_Time_sets{i,1,:} = SetFix;
end
%------------------------------

TIME_sets = cell(N,1);
Raw_TimeSets = VC_TSNH_Time_sets;
Raw_TempSets = VC_TSNH_Temp_sets;
Raw_DirSets = VC_TSNH_Dir_sets;
Raw_LattSets = VC_TSNH_Latt_sets;
Raw_LongSets = VC_TSNH_Long_sets;
Raw_SpeedSets = VC_TSNH_Speed_sets;
Raw_SalinitySets = VC_TSNH_Salinity_sets;
for i = 1:N;
    %create time set
    Set = [VC_TSNH_Time_sets{i,:,:}];
    TIME_sets{i,:,:} = (floor(min(Set)):((1/24)/(60/Step)):ceil(max(Set)))';
    %remove time-repeats in data
    RawSetTime = [Raw_TimeSets{i,:,:}];
    RawSetTemp = [Raw_TempSets{i,:,:}];
    RawSetDir = [Raw_DirSets{i,:,:}];
    RawSetLatt = [Raw_LattSets{i,:,:}];
    RawSetLong = [Raw_LongSets{i,:,:}];
    RawSetSpeed = [Raw_SpeedSets{i,:,:}];
    RawSetSalinity = [Raw_SalinitySets{i,:,:}];
    for j = 2:length(Set);
        if (Set(j,1)==Set(j-1,1));
            RawSetTime(j,1)=NaN;
        end
    end
    cut = find(isnan(RawSetTime));
    RawSetTime(cut) = [];
    RawSetTemp(cut) = [];
    RawSetDir(cut) = [];
    RawSetLatt(cut) = [];
    RawSetLong(cut) = [];
    RawSetSpeed(cut) = [];
    RawSetSalinity(cut) = [];
    Raw_TimeSets{i,:,:} = RawSetTime;
    Raw_TempSets{i,:,:} = RawSetTemp;
    Raw_DirSets{i,:,:} = RawSetDir;
    Raw_LattSets{i,:,:} = RawSetLatt;
    Raw_LongSets{i,:,:} = RawSetLong;
    Raw_SpeedSets{i,:,:} = RawSetSpeed;
    Raw_SalinitySets{i,:,:} = RawSetSalinity;
end

clear N i j Set cut RawSetTime RawSetTemp RawSetDir RawSetLatt...
    RawSetLong RawSetSpeed RawSetSalinity VC_TSNH_Time_sets VC_TSNH_Temp_sets...
    VC_TSNH_Dir_sets VC_TSNH_Latt_sets VC_TSNH_Long_sets VC_TSNH_Speed_sets...
    VC_TSNH_Salinity_sets VC_TSNH_Time_GPS_sets SetFix

%% interpolate to time step specified

N = size(TIME_sets,1);
%Salinity
Salinity_sets = cell(N,1);
for i = 1:N;
    Set = [Raw_SalinitySets{i,:,:}];
    RawTimeSet = [Raw_TimeSets{i,:,:}];
    timeSet = [TIME_sets{i,:,:}];
    Salinity_sets{i,:,:} = interp1(RawTimeSet,Set,timeSet);
end
%Temp
Temp_sets = cell(N,1);
for i = 1:N;
    Set = [Raw_TempSets{i,:,:}];
    RawTimeSet = [Raw_TimeSets{i,:,:}];
    timeSet = [TIME_sets{i,:,:}];
    Temp_sets{i,:,:} = interp1(RawTimeSet,Set,timeSet);
end
%Dir
Dir_sets = cell(N,1);
for i = 1:N;
    Set = [Raw_DirSets{i,:,:}];
    RawTimeSet = [Raw_TimeSets{i,:,:}];
    timeSet = [TIME_sets{i,:,:}];
    Dir_sets{i,:,:} = interp1(RawTimeSet,Set,timeSet);
end
%Latt
Latt_sets = cell(N,1);
for i = 1:N;
    Set = [Raw_LattSets{i,:,:}];
    RawTimeSet = [Raw_TimeSets{i,:,:}];
    timeSet = [TIME_sets{i,:,:}];
    Latt_sets{i,:,:} = interp1(RawTimeSet,Set,timeSet);
end
%Long
Long_sets = cell(N,1);
for i = 1:N;
    Set = [Raw_LongSets{i,:,:}];
    RawTimeSet = [Raw_TimeSets{i,:,:}];
    timeSet = [TIME_sets{i,:,:}];
    Long_sets{i,:,:} = interp1(RawTimeSet,Set,timeSet);
end
%Speed
Speed_sets = cell(N,1);
for i = 1:N;
    Set = [Raw_SpeedSets{i,:,:}];
    RawTimeSet = [Raw_TimeSets{i,:,:}];
    timeSet = [TIME_sets{i,:,:}];
    Speed_sets{i,:,:} = interp1(RawTimeSet,Set,timeSet);
end

clear Set RawTimeSet timeSet i N Raw_TimeSets Raw_TempSets...
    Raw_SalinitySets Raw_DirSets Raw_LattSets Raw_LongSets...
    Raw_SpeedSets
    
%% TSNH combine full series
TIME_VC_TSNH = cell2mat(TIME_sets);

Salinity_VC_TSNH = cell2mat(Salinity_sets);
 for i=1:length(Salinity_VC_TSNH);
     if (Salinity_VC_TSNH(i)<=0 || Salinity_VC_TSNH(i)>40);
         Salinity_VC_TSNH(i)=NaN;
     end
 end
Temp_VC_TSNH = cell2mat(Temp_sets);
 for i=1:length(Temp_VC_TSNH);
     if (Temp_VC_TSNH(i)<=0 || Temp_VC_TSNH(i)>30);
         Temp_VC_TSNH(i)=NaN;
     end
 end
Dir_VC_TSNH = cell2mat(Dir_sets);
Latt_VC_TSNH = cell2mat(Latt_sets);
Long_VC_TSNH = cell2mat(Long_sets);
Speed_VC_TSNH = cell2mat(Speed_sets);

clear i TIME_sets Salinity_sets Temp_sets Dir_sets Latt_sets...
        Long_sets Speed_sets

%% GPS Range Limit

GPS_ADM = find(Long_VC_TSNH<=-122.6 & Long_VC_TSNH>=-122.8 ...
            & Latt_VC_TSNH>=48.1 & Latt_VC_TSNH<=48.2);

%** check location **
% scatter(Long_VC_TSNH,Latt_VC_TSNH);hold on;
% scatter(Long_VC_TSNH(GPS_ADM),Latt_VC_TSNH(GPS_ADM),'x');
% scatter(-122.6862,48.1529,'or','filled');
% hold off;
% grid on;

%% Limit Data to Admiralty

VCsal_Adm = NaN(length(Salinity_VC_TSNH),1);
VCtemp_Adm = NaN(length(Temp_VC_TSNH),1);
VClatt_Adm = NaN(length(Latt_VC_TSNH),1);
VClong_Adm = NaN(length(Long_VC_TSNH),1);
T_Adm = NaN(length(TIME_VC_TSNH),1);
for i = GPS_ADM;
    VCsal_Adm(i) = Salinity_VC_TSNH(i);
    VCtemp_Adm(i) = Temp_VC_TSNH(i);
    VClatt_Adm(i) = Latt_VC_TSNH(i);
    VClong_Adm(i) = Long_VC_TSNH(i);
    T_Adm(i) = TIME_VC_TSNH(i);
end

%** check location **
% scatter(VClong_Adm,VClatt_Adm,'x');hold on;
% scatter(-122.6862,48.1529,'or','filled');
% hold off;
% grid on;

clear i


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

% hourly time from (16 May 2009) to (16 Jan 2013)
Hour = (datenum(2012,05,16):(1/24):datenum(2013,01,16))';

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

%% Interpolate VC Salinity data at Admiralty to hourly time scale

VCsal_raw = VCsal_Adm;
VCtime_raw = T_Adm;
cut = find(isnan(T_Adm));
VCtime_raw(cut) = [];
VCsal_raw(cut) = [];

%set up distinct value vector to eliminate NaN periods through rounding
Sal_VC_Finite = NaN(length(Salinity_VC_TSNH),1);
count = 0;
for i = 1:length(Sal_VC_Finite);
    if (isfinite(Salinity_VC_TSNH(i)));
        count = count+1;
        Sal_VC_Finite(i) = count;
    end
end

VCsal_raw_Nan = Sal_VC_Finite;
VCtime_raw_Nan = TIME_VC_TSNH;
cutNan = find(isnan(Sal_VC_Finite));
VCtime_raw_Nan(cutNan) = [];
VCsal_raw_Nan(cutNan) = [];

VCsal_hr = interp1(VCtime_raw,VCsal_raw,Hour);
VCsal_hr_NAN = interp1(VCtime_raw_Nan,VCsal_raw_Nan,Hour);
 for i = 1:length(VCsal_hr);
     if (VCsal_hr_NAN(i)~=round(VCsal_hr_NAN(i)));
         VCsal_hr(i) = NaN;
     end
 end

%***check interp
% plot(TIME_VC_TSNH,VCsal_Adm,'r');hold on;
% scatter(Hour,VCsal_hr,'xk');
% plot(Hour,VCsal_hr);hold off;grid on;dynamicDateTicks();
% xlim([TIME_VC_TSNH(1) TIME_VC_TSNH(length(TIME_VC_TSNH))]);

clear cut VCsal_raw VCtime_raw cutNan VCtime_raw_Nan VCsal_raw_Nan...
        VCsal_hr_NAN i count

%% Interpolate VC Temp data at Admiralty to hourly time scale

VCtemp_raw = VCtemp_Adm;
VCtime_raw = T_Adm;
cut = find(isnan(T_Adm));
VCtime_raw(cut) = [];
VCtemp_raw(cut) = [];

Temp_VC_Finite = NaN(length(Temp_VC_TSNH),1);
count = 0;
for i = 1:length(Temp_VC_Finite);
    if (isfinite(Temp_VC_TSNH(i)));
        count = count+1;
        Temp_VC_Finite(i) = count;
    end
end

VCtemp_raw_Nan = Temp_VC_Finite;
VCtime_raw_Nan = TIME_VC_TSNH;
cutNan = find(isnan(Temp_VC_Finite));
VCtime_raw_Nan(cutNan) = [];
VCtemp_raw_Nan(cutNan) = [];

VCtemp_hr = interp1(VCtime_raw,VCtemp_raw,Hour);
VCtemp_hr_NAN = interp1(VCtime_raw_Nan,VCtemp_raw_Nan,Hour);
 for i = 1:length(VCtemp_hr);
     if (VCtemp_hr_NAN(i)~=round(VCtemp_hr_NAN(i)));
         VCtemp_hr(i) = NaN;
     end
 end

%***check interp
% plot(TIME_VC_TSNH,VCtemp_Adm,'r');hold on;
% scatter(Hour,VCtemp_hr,'xk');
% plot(Hour,VCtemp_hr);hold off;grid on;dynamicDateTicks();
% xlim([TIME_VC_TSNH(1) TIME_VC_TSNH(length(TIME_VC_TSNH))]);

clear cut VCtemp_raw VCtime_raw cutNan VCtime_raw_Nan VCtemp_raw_Nan...
        VCtemp_hr_NAN i count
    
%% 40-hour LowPassFilter & Seasonal Filter
% match residual current time-scale [subtidal signal]

VCsal_lpf40 = pl66tn(VCsal_hr,1,40);
% **check filter**
% plot(TIME_VC_TSNH,VCsal_Adm,'g');hold on;
% scatter(Hour,VCsal_hr,'xk');
% plot(Hour,VCsal_lpf40,'b');
% hold off;grid on;dynamicDateTicks();
% xlim([TIME_VC_TSNH(1) TIME_VC_TSNH(length(TIME_VC_TSNH))]);

VCtemp_lpf40 = pl66tn(VCtemp_hr,1,40);
% **check filter**
% plot(TIME_VC_TSNH,VCtemp_Adm,'g');hold on;
% scatter(Hour,VCtemp_hr,'xk');
% plot(Hour,VCtemp_lpf40,'b');
% hold off;grid on;dynamicDateTicks();
% xlim([TIME_VC_TSNH(1) TIME_VC_TSNH(length(TIME_VC_TSNH))]);

%% ST_DOP (Sub-Tidal Dissolved Oxygen Prediction)

% use equation derived from regression
ST_DOP = 32.6253 - (0.6094*VCsal_lpf40) - (0.8084*VCtemp_lpf40);

%% ST_DOP plot

plot(Hour,ST_DOP,'--b');
title('Empirical Sub-Tidal Dissolved Oxygen Prediction form Victoria Clipper Data'...
    ,'FontSize',22);
ylabel('Dissolved Oxygen (mg/L)','FontSize',22);
set(gca,'YDir','normal','FontSize',18);
legend('Sub-tidal DO Prediction');
grid on;
ylim([4 9]);
dynamicDateTicks();

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

clear ticks
