% Extract information from *_his.nc files

clear

%% Load peak discharge (180, 300) stage and depth
obs180_fnl_2010 = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012_histextent/run7009_180cms/'...
    'abq10a12bshc180cmsrst01_0000_20160104_110000_his.nc'];
obs180_fnl_hist = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqhistorical/run6009_180cms/'...
    'historicf180cmsrst03_0000_19180104_170000_his.nc'];
obs180_fnl_2010lidar = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12blhc/run8003_180cms/'...
    'abq10a12blhcf180cmsrst01_0000_20160104_010000_his.nc'];
obs180_fnl_prelidar = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqpreblhc/run9003_180cms/'...
    'abqprebhcf180cmsrst01_0000_20160104_050000_his.nc'];
obs300_fnl_2010 = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012_histextent/run7010_300cms/'...
    'abq10a12bshc300cmsrst01_0000_20160103_150000_his.nc'];
obs300_fnl_hist = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqhistorical/run6010_300cms/'...
    'historicf300cmsrst03_0000_19180104_090000_his.nc'];
obs300_fnl_2010lidar = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12blhc/run8004_300cms/'...
    'abq10a12blhcf300cmsrst01_0000_20160103_090000_his.nc'];
obs300_fnl_prelidar = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqpreblhc/run9004_300cms/'...
    'abqprebhcf300cmsrst01_0000_20160103_060000_his.nc'];

all_2010_180 = ncread(obs180_fnl_2010,'waterlevel');
all_hist_180 = ncread(obs180_fnl_hist,'waterlevel');
all_2010lidar_180 = ncread(obs180_fnl_2010lidar,'waterlevel');
all_prelidar_180 = ncread(obs180_fnl_prelidar,'waterlevel');
all_2010_300 = ncread(obs300_fnl_2010,'waterlevel');
all_hist_300 = ncread(obs300_fnl_hist,'waterlevel');
all_2010lidar_300 = ncread(obs300_fnl_2010lidar,'waterlevel');
all_prelidar_300 = ncread(obs300_fnl_prelidar,'waterlevel');

obs_2010_wse(:,1) = all_2010_180(:,end);
obs_hist_wse(:,1) = all_hist_180(:,end);
obs_2010lidar_wse(:,1) = all_2010lidar_180(:,end);
obs_prelidar_wse(:,1) = all_prelidar_180(:,end);
obs_2010_wse(:,2) = all_2010_300(:,end);
obs_hist_wse(:,2) = all_hist_300(:,end);
obs_2010lidar_wse(:,2) = all_2010lidar_300(:,end);
obs_prelidar_wse(:,2) = all_prelidar_300(:,end);

all_2010_180_d = ncread(obs180_fnl_2010,'Waterdepth');
all_hist_180_d = ncread(obs180_fnl_hist,'Waterdepth');
all_2010lidar_180_d = ncread(obs180_fnl_2010lidar,'Waterdepth');
all_prelidar_180_d = ncread(obs180_fnl_prelidar,'Waterdepth');
all_2010_300_d = ncread(obs300_fnl_2010,'Waterdepth');
all_hist_300_d = ncread(obs300_fnl_hist,'Waterdepth');
all_2010lidar_300_d = ncread(obs300_fnl_2010lidar,'Waterdepth');
all_prelidar_300_d = ncread(obs300_fnl_prelidar,'Waterdepth');

obs_2010_d(:,1) = all_2010_180_d(:,end);
obs_hist_d(:,1) = all_hist_180_d(:,end);
obs_2010lidar_d(:,1) = all_2010lidar_180_d(:,end);
obs_prelidar_d(:,1) = all_prelidar_180_d(:,end);
obs_2010_d(:,2) = all_2010_300_d(:,end);
obs_hist_d(:,2) = all_hist_300_d(:,end);
obs_2010lidar_d(:,2) = all_2010lidar_300_d(:,end);
obs_prelidar_d(:,2) = all_prelidar_300_d(:,end);

% Load local minimum elevations for contemporary and historical models
load(['/Users/Colin/Documents/MATLAB/DFlow/attenuation_chap1/data/'...
    'valley250minZcontemporary.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/attenuation_chap1/data/'...
    'rg2010valley_250m_chan_min_z.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/attenuation_chap1/data/'...
    'valley250minZhistorical.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/attenuation_chap1/data/'...
    'rg1918valley_250m_chan_min_z.mat']);

%% INPUTS for 2010 historical comparison

chan_dist_max = 18000; % km
max_Q = [180; 300; 180; 300; 180; 300]; % m^3/s
site_num = 15; % number of sites
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012_histextent/'];

% Define folder with *his.nc files
fpath = ['run7003_s180p/';'run7011_s300p/';...
    'run7005_m180p/';'run7012_m300p/';...
    'run7007_l180p/';'run7013_l300p/']; % file path
% files = struct2cell(dir(fpath)); % Get files from folder
% [~,file_count] = size(files);
% fnames = files(1,3:file_count)'; % extract file names
% [fnum,~] = size(fnames); % number of files

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_2010_cell = cell(folder_num,7);
obs_wse_2010_cell = cell(folder_num,8);
% raw_cell = cell(folder_num,3);
% distance_cell = cell(folder_num,9);
% Q_cell = cell(folder_num,7);
% Ma_cell = cell(folder_num,13);
% Mo_cell = cell(folder_num,15);
% max_min_cell = cell(folder_num,37);

for fi = 1:folder_num
    
    xs_2010_cell(fi,1) = {fpath(fi,:)};
    obs_wse_2010_cell(fi,1) = {fpath(fi,:)};
%     raw_cell(fi,1) = {fpath(fi,:)};
%     distance_cell(fi,1) = {fpath(fi,:)};
%     Q_cell(fi,1) = {fpath(fi,:)};
%     Ma_cell(fi,1) = {fpath(fi,:)};
%     Mo_cell(fi,1) = {fpath(fi,:)};
%     max_min_cell(fi,1) = {fpath(fi,:)};

    if mod(fi,2) == 0
        obs_peak_wse = obs_2010_wse(:,2);
    elseif mod(fi,2) == 1
        obs_peak_wse = obs_2010_wse(:,1);
    end
    
    max_Q_fi = max_Q(fi,1);
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

    % Import channel/island areas
    load(['/Users/Colin/Documents/MATLAB/DFlow/mass_Mo_chap3/'...
        'rg_channel_island_area.mat']);

    % Import restoration site connection polygons
    rest_con = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
        'mass_Mo_chap3/abq12_restoration_feature_connections_pts']);

    rs_index = cell(1,site_num);
    % Observation ID array
    for s = 1:site_num
        rs_index(s) = {strcat('s',num2str(s))}; % create cell with site ids
    end
    obs_marker = ['l', 'r']; % array marking left or right boundaries

    %% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    
    [xs_2010_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t, ...
        max_t, max_Q_fi);
    
    %% Extract observation point WSE information
    
    [obs_wse_2010_cell(fi,2:9)] = ...
        hpeak_obs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t,...
        max_t, obs_peak_wse, rg2010valley_250m_chan_min_z);
    
    %% Extract observation point depth information

end

%% INPUTS for historical model

chan_dist_max = 18000; % km
max_Q = [180; 300; 180; 300; 180; 300]; % m^3/s
site_num = 15; % number of sites
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqhistorical/'];

% Define folder with *his.nc files
fpath = {'run6025_s180redo/';'run6026_s300redo/';...
    'run6005_abqhfm180p/';'run6006_abqhfm300p/';...
    'run6007_abqhfl180p/';'run6020_l300pfull/'}; % file path
% files = struct2cell(dir(fpath)); % Get files from folder
% [~,file_count] = size(files);
% fnames = files(1,3:file_count)'; % extract file names
% [fnum,~] = size(fnames); % number of files

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_hist_cell = cell(folder_num,7);
obs_wse_hist_cell = cell(folder_num,8);
% raw_cell = cell(folder_num,3);
% distance_cell = cell(folder_num,9);
% Q_cell = cell(folder_num,7);
% Ma_cell = cell(folder_num,13);
% Mo_cell = cell(folder_num,15);
% max_min_cell = cell(folder_num,37);

for fi = 1:folder_num
    
    xs_hist_cell(fi,1) = {fpath(fi,:)};
    obs_wse_hist_cell(fi,1) = {fpath(fi,:)};
%     raw_cell(fi,1) = {fpath(fi,:)};
%     distance_cell(fi,1) = {fpath(fi,:)};
%     Q_cell(fi,1) = {fpath(fi,:)};
%     Ma_cell(fi,1) = {fpath(fi,:)};
%     Mo_cell(fi,1) = {fpath(fi,:)};
%     max_min_cell(fi,1) = {fpath(fi,:)};
    
    if mod(fi,2) == 0
        obs_peak_wse = obs_hist_wse(:,2);
    elseif mod(fi,2) == 1
        obs_peak_wse = obs_hist_wse(:,1);
    end
    
    max_Q_fi = max_Q(fi,1);
    
    fpath_fi = strcat(fpath_root,fpath{fi,:});
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

    % Import channel/island areas
    load(['/Users/Colin/Documents/MATLAB/DFlow/mass_Mo_chap3/'...
        'rg_channel_island_area.mat']);

    % Import restoration site connection polygons
    rest_con = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
        'mass_Mo_chap3/abq12_restoration_feature_connections_pts']);

    rs_index = cell(1,site_num);
    % Observation ID array
    for s = 1:site_num
        rs_index(s) = {strcat('s',num2str(s))}; % create cell with site ids
    end
    obs_marker = ['l', 'r']; % array marking left or right boundaries

    %% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    
    [xs_hist_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t, ...
        max_t, max_Q_fi);
    
    %% Extract observation point WSE information
    
    [obs_wse_hist_cell(fi,2:9)] = ...
        hpeak_obs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t, ...
        max_t, obs_peak_wse, rg1918valley_250m_chan_min_z);
    
    %% Extract observation point depth information
    
    
    
end

%% INPUTS for 2010 lidar model

chan_dist_max = 25000; % km
max_Q = [180; 300; 180; 300; 180; 300]; % m^3/s
site_num = 15; % number of sites
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12blhc/'];

% Define folder with *his.nc files
fpath = {'run8012_s180pfull/';'run8013_s300pfull/';...
    'run8014_m180pfull/';'run8015_m300pfull/';...
    'run8016_l180pfull/';'run8017_l300pfull/'}; % file path
% files = struct2cell(dir(fpath)); % Get files from folder
% [~,file_count] = size(files);
% fnames = files(1,3:file_count)'; % extract file names
% [fnum,~] = size(fnames); % number of files

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_2010lidar_cell = cell(folder_num,7);
obs_wse_2010lidar_cell = cell(folder_num,8);
% raw_cell = cell(folder_num,3);
% distance_cell = cell(folder_num,9);
% Q_cell = cell(folder_num,7);
% Ma_cell = cell(folder_num,13);
% Mo_cell = cell(folder_num,15);
% max_min_cell = cell(folder_num,37);

for fi = 1:folder_num
    
    xs_2010lidar_cell(fi,1) = {fpath(fi,:)};
    obs_wse_2010lidar_cell(fi,1) = {fpath(fi,:)};
%     raw_cell(fi,1) = {fpath(fi,:)};
%     distance_cell(fi,1) = {fpath(fi,:)};
%     Q_cell(fi,1) = {fpath(fi,:)};
%     Ma_cell(fi,1) = {fpath(fi,:)};
%     Mo_cell(fi,1) = {fpath(fi,:)};
%     max_min_cell(fi,1) = {fpath(fi,:)};
    
    if mod(fi,2) == 0
        obs_peak_wse = obs_2010lidar_wse(:,2);
    elseif mod(fi,2) == 1
        obs_peak_wse = obs_2010lidar_wse(:,1);
    end
    
    max_Q_fi = max_Q(fi,1);
    
    fpath_fi = strcat(fpath_root,fpath{fi,:});
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

    %% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    
    [xs_2010lidar_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t, ...
        max_t, max_Q_fi);
    
    %% Extract observation point WSE information
    
    [obs_wse_2010lidar_cell(fi,2:9)] = ...
        hpeak_obs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t, ...
        max_t, obs_peak_wse, rg2010valley_250m_chan_min_z);
    
    %% Extract observation point depth information
    
    
    
end

%% INPUTS for pre-restoration lidar model

chan_dist_max = 25000; % km
max_Q = [180; 300; 180; 300; 180; 300]; % m^3/s
site_num = 15; % number of sites
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqpreblhc/'];

% Define folder with *his.nc files
fpath = ['run9012_s180pfull/';'run9013_s300pfull/';...
    'run9014_m180pfull/';'run9015_m300pfull/';...
    'run9016_l180pfull/';'run9017_l300pfull/']; % file path
% files = struct2cell(dir(fpath)); % Get files from folder
% [~,file_count] = size(files);
% fnames = files(1,3:file_count)'; % extract file names
% [fnum,~] = size(fnames); % number of files

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_prelidar_cell = cell(folder_num,7);
obs_wse_prelidar_cell = cell(folder_num,8);
% raw_cell = cell(folder_num,3);
% distance_cell = cell(folder_num,9);
% Q_cell = cell(folder_num,7);
% Ma_cell = cell(folder_num,13);
% Mo_cell = cell(folder_num,15);
% max_min_cell = cell(folder_num,37);

for fi = 1:folder_num
    
    xs_prelidar_cell(fi,1) = {fpath(fi,:)};
    obs_wse_prelidar_cell(fi,1) = {fpath(fi,:)};
%     raw_cell(fi,1) = {fpath(fi,:)};
%     distance_cell(fi,1) = {fpath(fi,:)};
%     Q_cell(fi,1) = {fpath(fi,:)};
%     Ma_cell(fi,1) = {fpath(fi,:)};
%     Mo_cell(fi,1) = {fpath(fi,:)};
%     max_min_cell(fi,1) = {fpath(fi,:)};
    
    if mod(fi,2) == 0
        obs_peak_wse = obs_prelidar_wse(:,2);
    elseif mod(fi,2) == 1
        obs_peak_wse = obs_prelidar_wse(:,1);
    end
    
    max_Q_fi = max_Q(fi,1);
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

    %% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    
    [xs_prelidar_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t, ...
        max_t, max_Q_fi);
    
    %% Extract observation point WSE information
    
    [obs_wse_prelidar_cell(fi,2:9)] = ...
        hpeak_obs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t, ...
        max_t, obs_peak_wse, rg2010valley_250m_chan_min_z);
    
    %% Extract observation point depth information
    
    
    
end
