%% Script to merge DUNEX-RT Archive to a "best series" dataset
% Alexander Rey, July 15, 2022

addpath(genpath('matlab/applications/delft3d_matlab'),'-begin')
addpath(genpath('D:/DUNEX_RT/Operation/Hindcast'),'-begin')
%% Get descriptive information of the hist file
runDate = datetime(2020, 07, 20);
histVars = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate,'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1);
histStations = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate,'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'water level','stations');

%% Extract data from July 20, 2020 to October 2, 2020

% Extract only the first 6 hours (36 time steps)
% Extract only stations for DUNEX (98-129)
mergedHistWL = [];
mergedHistVelX = [];
mergedHistVelY = [];
mergeHistWaveHeight = [];
mergeHistWavePer = [];
mergeHistWaveDir = [];
mergedHistTime = [];
mergeHistWindSpeed = [];
mergeHistWindDir = [];

% Loop through times 
for runDate = datetime(2020, 07, 20):hours(6):datetime(2020, 10, 02)
% for runDate = datetime(2020, 09, 16):hours(6):datetime(2020, 09, 26)
% for runDate = datetime(2020, 08, 03):hours(6):datetime(2020, 08, 06)

    % If file is missing go back and use correct time steps from later in
    % the run
    if exist(['D:\DUNEX_RT\Archive\' datestr(runDate,'yyyy-mm-dd_HH')...
        '_HistOut7.dat'],'file') == 2
        runSteps = 2:37;
        hoursBack = 0;
    elseif exist(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(6), 'yyyy-mm-dd_HH')...
        '_HistOut7.dat'],'file') == 2
        runSteps = 38:73;
        hoursBack = 6
    elseif exist(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(12),'yyyy-mm-dd_HH')...
        '_HistOut7.dat'],'file') == 2
        runSteps = 74:109;
        hoursBack = 12
    elseif exist(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(18),'yyyy-mm-dd_HH')...
        '_HistOut7.dat'],'file') == 2
        runSteps = 110:145;
        hoursBack = 18
    elseif exist(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(24),'yyyy-mm-dd_HH')...
        '_HistOut7.dat'],'file') == 2
        runSteps = 146:181;
        hoursBack = 24
    elseif exist(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(30),'yyyy-mm-dd_HH')...
        '_HistOut7.dat'],'file') == 2
        runSteps = 182:217;
        hoursBack = 30
    else
        % If missing
        mergedHistWL = [mergedHistWL nan([32,36])];
        mergedHistVelX = [mergedHistVelX nan([32,36])];
        mergedHistVelY = [mergedHistVelY nan([32,36])];
        mergeHistWaveHeight = [mergeHistWaveHeight nan([32,36])];
        mergeHistWavePer = [mergeHistWavePer nan([32,36])];
        mergeHistWaveDir = [mergeHistWaveDir nan([32,36])];       
        mergeHistWindSpeed = [mergeHistWindSpeed nan([32,36])];
        mergeHistWindDir = [mergeHistWindDir nan([32,36])];                
        mergedHistTime = [mergedHistTime datenum(runDate+minutes(10):minutes(10):runDate+hours(6))];
        continue
    end  

    % Read in model data
    modelHistWL = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(hoursBack),'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'water level','griddata',runSteps,98:129);

    modelHistVel = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(hoursBack),'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'depth averaged velocity','griddata',runSteps,98:129);

    modelHistWaveHeight = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(hoursBack),'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'significant wave height','griddata',runSteps,98:129);
    
    modelHistWavePer = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(hoursBack),'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'peak wave period','griddata',runSteps,98:129);

    modelHistWaveDir = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(hoursBack),'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'wave direction','griddata',runSteps,98:129);    

    modelHistWindSpeed = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(hoursBack),'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'wind speed','griddata',runSteps,98:129);    

    modelHistWindDir = qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(hoursBack),'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'wind direction','griddata',runSteps,98:129);    

    % Merge
    mergedHistWL = [mergedHistWL modelHistWL.Val'];
    mergedHistVelX = [mergedHistVelX modelHistVel.XComp'];
    mergedHistVelY = [mergedHistVelY modelHistVel.YComp'];
    mergeHistWaveHeight = [mergeHistWaveHeight modelHistWaveHeight.Val'];
    mergeHistWavePer = [mergeHistWavePer modelHistWavePer.Val'];
    mergeHistWaveDir = [mergeHistWaveDir modelHistWaveDir.Val'];    
    mergeHistWindSpeed = [mergeHistWindSpeed modelHistWindSpeed.Val'];
    mergeHistWindDir = [mergeHistWindDir modelHistWindDir.Val'];        
    mergedHistTime = [mergedHistTime modelHistWL.Time'];

    disp(runDate)
end

extractBed =  qpread(qpfopen(['D:\DUNEX_RT\Archive\' datestr(runDate-hours(hoursBack),'yyyy-mm-dd_HH')...
        '_HistOut7.dat']),1,'bed level','griddata',1,98:129);
extractPtX = extractBed.X(1,:);
extractPtY = extractBed.Y(1,:);
%% Plot

figure
plot(mergedHistTime, mergedHistWL(10, :))
datetick('x')
ylabel('DUNEX-RT Water Level (mNAVD88)')

figure
plot(mergedHistTime, mergeHistWaveHeight(10, :))
datetick('x')
ylabel('DUNEX-RT Significant Wave Height (m8)')


%% Export data

MergedWL_Table = array2table(mergedHistWL');
MergedWL_Table.Properties.VariableNames = histStations(98:129);
MergedWL_Table.DateTime = datetime(datevec(mergedHistTime)) ;
writetable(MergedWL_Table, 'D:\Alexander\MATLAB\DunexBounds\Dunex_RT_WL.csv')

MergedVelX_Table = array2table(mergedHistVelX');
MergedVelX_Table.Properties.VariableNames = histStations(98:129);
MergedVelX_Table.DateTime = datetime(datevec(mergedHistTime)) ;
writetable(MergedVelX_Table, 'D:\Alexander\MATLAB\DunexBounds\Dunex_RT_VelX.csv')

MergedVelY_Table = array2table(mergedHistVelY');
MergedVelY_Table.Properties.VariableNames = histStations(98:129);
MergedVelY_Table.DateTime = datetime(datevec(mergedHistTime)) ;
writetable(MergedVelY_Table, 'D:\Alexander\MATLAB\DunexBounds\Dunex_RT_VelY.csv')

MergedWaveHeight_Table = array2table(mergeHistWaveHeight');
MergedWaveHeight_Table.Properties.VariableNames = histStations(98:129);
MergedWaveHeight_Table.DateTime = datetime(datevec(mergedHistTime)) ;
writetable(MergedWaveHeight_Table, 'D:\Alexander\MATLAB\DunexBounds\Dunex_RT_SigWave.csv')

MergedWavePeriod_Table = array2table(mergeHistWavePer');
MergedWavePeriod_Table.Properties.VariableNames = histStations(98:129);
MergedWavePeriod_Table.DateTime = datetime(datevec(mergedHistTime)) ;
writetable(MergedWavePeriod_Table, 'D:\Alexander\MATLAB\DunexBounds\Dunex_RT_Period.csv')

MergedWaveDir_Table = array2table(mergeHistWaveDir');
MergedWaveDir_Table.Properties.VariableNames = histStations(98:129);
MergedWaveDir_Table.DateTime = datetime(datevec(mergedHistTime)) ;
writetable(MergedWaveDir_Table, 'D:\Alexander\MATLAB\DunexBounds\Dunex_RT_Dir.csv')

MergedWindSpeed_Table = array2table(mergeHistWindSpeed');
MergedWindSpeed_Table.Properties.VariableNames = histStations(98:129);
MergedWindSpeed_Table.DateTime = datetime(datevec(mergedHistTime)) ;
writetable(MergedWindSpeed_Table, 'D:\Alexander\MATLAB\DunexBounds\Dunex_RT_WindSpeed.csv')

MergedWindDir_Table = array2table(mergeHistWindDir');
MergedWindDir_Table.Properties.VariableNames = histStations(98:129);
MergedWindDir_Table.DateTime = datetime(datevec(mergedHistTime)) ;
writetable(MergedWindDir_Table, 'D:\Alexander\MATLAB\DunexBounds\Dunex_RT_WindDir.csv')

stationDetails = cell2table([histStations(98:129) num2cell(extractPtX') num2cell(extractPtY') num2cell(squeeze(extractBed.Val(1,1,:)))]);
stationDetails.Properties.VariableNames = {'StationName','UTMx','UTMy', 'BedLevel'};
writetable(stationDetails, 'D:\Alexander\MATLAB\DunexBounds\stationDetails.csv')