%% Compare 2014 Delft3D run with Rudy data
% AJMR, July 20, 2022

addpath(genpath('matlab/applications/delft3d_matlab'),'-begin')
    
%% Read in observation points
clear obsPtsLon obsPtsLat

obsFiles = dir("\\COASTLINES.appsci.queensu.ca\Rudy2014\Google earth objects");
obsCount = 1;

% Read in in order 1-13, L2, PE to match spreadsheets
for i = [7 13:20 8:10 12 6 4]
    [obsPtsLon(obsCount), obsPtsLat(obsCount)] = read_kml(['\\COASTLINES.appsci.queensu.ca\Rudy2014\Google earth objects\' obsFiles(i).name]);
    [obsPtsUTMx(obsCount), obsPtsUTMy(obsCount)] = wgs2utm(obsPtsLat(obsCount),obsPtsLon(obsCount),18,'N');

    obsCount = obsCount + 1;
end

%% Read in data at observation points
clear obsDat

datFiles =  dir ("D:\Alexander\Rudy2014\Point Sorted Data");

% Read in in order 1-13, L2, PE to match spreadsheets
obsCount = 1;
for i = 4:length(datFiles)
    % Read data as table
    obsDat{2, obsCount} = readtable(['\\COASTLINES.appsci.queensu.ca\Rudy2014\Point Sorted Data\' ...
        datFiles(i).name],'DataRange','X3:AQ17','VariableNamesRange',...
        'X1:AQ1', 'VariableUnitsRange', 'X2:AQ2');
    
    % Drop empty columns
    obsDat{2, obsCount} = removevars(obsDat{2, obsCount},{ ...
        'Var3', 'Var7', 'Var9', 'Var11', 'Var13', 'Var15', 'Var17', 'Var19'});
 
    % Add Date
    dateIN = readcell(['D:\Alexander\Rudy2014\Point Sorted Data\' ...
        datFiles(i).name], 'Range','C3:D3');
    try
        obsDat{1, obsCount} = datetime(dateIN{1},'InputFormat', 'MM/dd/uuuu') + ...
            days(dateIN{2});
    catch
        obsDat{1, obsCount} = datetime(dateIN{1},'InputFormat', 'MM/dd/uuuu') + ...
            hours(hour(dateIN{2})) + minutes(minute(dateIN{2}));
    end

    obsCount = obsCount + 1;

end

%% Read in Model results at observation points

%% Read in grids
Grid{1} = wlgrid('read','\\batchelor.appsci.queensu.ca\G\Fateme\Oct_7_RudyRun\G9_1.grd');
Grid{2} = wlgrid('read','\\batchelor.appsci.queensu.ca\G\Fateme\Oct_7_RudyRun\G9_2.grd');
Grid{3} = wlgrid('read','\\batchelor.appsci.queensu.ca\G\Fateme\Oct_7_RudyRun\G9_3.grd');

modelTimes = datetime(datevec(qpread(qpfopen('\\batchelor.appsci.queensu.ca\G\Fateme\Oct_7_RudyRun\trim-G9_1.dat'),...
    1,'water level','times')));


%% Match grid points
for i = 1:length(obsPtsUTMx)
    [GridMatch{1},GridDist{1}] = NearestValue([obsPtsUTMx(i), obsPtsUTMy(i)], Grid{1}.X, Grid{1}.Y);
    [GridMatch{2},GridDist{2}] = NearestValue([obsPtsUTMx(i), obsPtsUTMy(i)], Grid{2}.X, Grid{2}.Y);
    [GridMatch{3},GridDist{3}] = NearestValue([obsPtsUTMx(i), obsPtsUTMy(i)], Grid{3}.X, Grid{3}.Y);

    if sqrt(GridDist{1}(1)^2 + GridDist{1}(2)^2) < sqrt(GridDist{2}(1)^2 + GridDist{2}(2)^2) && ...
           sqrt(GridDist{1}(1)^2 + GridDist{1}(2)^2) < sqrt(GridDist{3}(1)^2 + GridDist{3}(2)^2)
        
        gridMatch(:, i) = GridMatch{1};
        gridMatchID(i)  = 1;

    elseif sqrt(GridDist{2}(1)^2 + GridDist{2}(2)^2) <= sqrt(GridDist{1}(1)^2 + GridDist{1}(2)^2) && ...
           sqrt(GridDist{2}(1)^2 + GridDist{2}(2)^2) <= sqrt(GridDist{3}(1)^2 + GridDist{3}(2)^2)
        
        gridMatch(:, i) = GridMatch{2};
        gridMatchID(i)  = 2;

    elseif sqrt(GridDist{3}(1)^2 + GridDist{3}(2)^2) <= sqrt(GridDist{1}(1)^2 + GridDist{1}(2)^2) && ...
           sqrt(GridDist{3}(1)^2 + GridDist{3}(2)^2) <= sqrt(GridDist{2}(1)^2 + GridDist{2}(2)^2)

        gridMatch(:, i) = GridMatch{3};
        gridMatchID(i)  = 3;

    end
end



%% Extract time series surface temperatures for each point
clear modelTemp

for d = 1:size(obsDat,2)
    modelTStep(d) = NearestValue(obsDat{1, d}, modelTimes);
end

for i = 1:length(gridMatch)
    modelIN = qpread(qpfopen(['\\batchelor.appsci.queensu.ca\G\Fateme\Oct_7_RudyRun\trim-G9_' ...
        num2str(gridMatchID(i)) '.dat']),1,'temperature','griddata', ...
        0,gridMatch(1, i),gridMatch(2, i),0);


    modelTemp(:,:,i) = modelIN.Val;


    disp(i)
end

%% Read in other spreadsheet
boundDat = readtable("\\COASTLINES.appsci.queensu.ca\Rudy2014\Boundary condition data.xlsx", 'Sheet', ...
    'Boundary conditions');

%% Read in WaterTrax
traxDat = readtable("\\COASTLINES.appsci.queensu.ca\Rudy2014\Copy of Alexander SheetModelling.xlsx", 'Sheet', ...
    '2014_Rudy');

%% Plotting

figure
hold on
plot(modelTimes,mean(modelTemp(:,:,13),2),'k')
plot(modelTimes,modelTemp(:,1,13),'m')
plot(modelTimes,modelTemp(:,8,13),'b')

% plot(boundDat.Date,boundDat.Temp_1)
% scatter(traxDat.Date+hours(14),traxDat.WaterTEffluent__C_,'r')
scatter(traxDat.Date(1:end-7)+hours(14),traxDat.WaterTEffluent__C_(1:end-7),'r')
for d = [1:15 17:size(obsDat,2)]
    scatter(obsDat{1, d}, obsDat{2, d}{13,"Temp"},'b')
end
legend('Model', 'Observations', 'Rudy Observations')
ylabel('Water Temperature °C')

%% Save Model Data
writematrix(squeeze(mean(modelTemp(:,:,1:13),2)),...
        'D:\Alexander\Rudy2014\modelData.xlsx')

writematrix(modelTimes,...
        'D:\Alexander\Rudy2014\modelTimes.xlsx')

%% Stats

statCount = 1;
for d = 1:length(traxDat.Date)
    traxIDXs(statCount) = NearestValue(traxDat.Date(d)+hours(14),modelTimes);
    statCount = statCount + 1;
end

% [traxR2, traxRMSE] = rsquare(traxDat.WaterTEffluent__C_, mean(modelTemp(traxIDXs,:,13),2), false)
% [traxR2, traxRMSE] = rsquare(traxDat.WaterTEffluent__C_(1:end-7), modelTemp(traxIDXs(1:end-7),8,13), false)
[traxR2, traxRMSE] = rsquare(mean(modelTemp(traxIDXs(1:end-7),:,13),2),...
    traxDat.WaterTEffluent__C_(1:end-7))

writematrix([mean(modelTemp(traxIDXs(1:end),:,13),2) traxDat.WaterTEffluent__C_(1:end)],...
    'D:\Alexander\Rudy2014\matchedWaterTrax.xlsx')


statCount = 1;
for d = [1:15 17:size(obsDat,2)]
    rudIDXs(statCount) = NearestValue(obsDat{1, d},modelTimes);
    rudDat(statCount)  = obsDat{2, d}{13,"Temp"};
    statCount = statCount + 1;
end

% [traxR2, traxRMSE] = rsquare(mean(modelTemp(traxIDXs,:,13),2), traxDat.WaterTEffluent__C_)
[rudR2, rudRMSE] = rsquare(rudDat', squeeze(modelTemp(rudIDXs,4,13)))
% [traxR2, traxRMSE] = rsquare(mean(modelTemp(traxIDXs(1:end-7),:,13),2),...
%     traxDat.WaterTEffluent__C_(1:end-7))