% Prepare PAR data from SeaWiFS, MODIS-Aqua, and MODIS-Terra
% In response to reviews on phenology paper 2022-2023
% Goal: put all data from all three sensors into one database
% Should be size 300 x 120 x days for 1997-2022 (9248 days?)
% Will need to average MODIS data to fit coarser spatial res of SeaWiFS.

clear; format shortg; % close all hidden;
rootdir = '/Users/jsturner/satellite_data/WAP_PAR/';

%% Create date array
day_list = datenum(1997,9,6):datenum(2022,6,31);
% SeaWIFS only 1997-09-06 to 2000-02-23
% Average of SeaWiFS and Terra 2000-02-24 to 2002-07-03
% Average of SeaWiFS, Terra, and Aqua 2002-07-04 to 2010-12-10
% Average of Terra and Aqua 2010-12-11 to 2022-06-31

par_database = nan(300,120,length(day_list));

%% SeaWiFS only 1997-09-06 to 2000-02-24
s_day_list = datenum(1997,9,6):datenum(2000,02,23);
% Careful: not every day in early parts of SeaWiFS timeline
% have to pad empty days with NaN

folder_s = [rootdir 'SeaWiFS_PAR/'];

% Create filelist SeaWiFS only
dir1 = dir(folder_s);
lf=length(dir1);
filelist_s={};
for i=1:lf
    nf=dir1(i).name; tf=contains(nf,'.nc');  % is it an nc file
    if tf==1 % if this IS an nc file 
        filelist_s=[filelist_s;nf];
    end
end; clear lf nf tf

for i = 1:length(s_day_list)
    thisday = datestr(s_day_list(i),'yyyymmdd');
    % load in the file that corresponds with this date, using "contains"
    ind_thisday = find(contains(filelist_s,thisday));
    filename_s = char(filelist_s(ind_thisday));
    if isempty(filename_s)
        par_database(:,:,i) = nan;
    else
        file1 = [folder_s filename_s];
        ds = ncreadatt(file1,'/','product_name');
        datestr1 = [ds(21:24) '-' ds(25:26) '-' ds(27:28)]; % SeaWiFS (different for SeaWiFS, Aqua, and MODIS)
        par = ncread(file1,'par');
            if i==1
            lon0 = ncread(file1,'lon'); 
            lat0 = ncread(file1,'lat'); 
            lon_par_s = repmat(lon0,1,size(par,2));
            lat_par_s = repmat(lat0,1,size(par,1))';
            end
        par_database(:,:,i) = par;
    end
end
disp(' ');
disp('Done with SeaWiFS only data chunk.');

%% Average of SeaWiFS and Terra 2000-02-24 to 2002-07-03
st_day_list = datenum(2000,2,24):datenum(2002,7,3);
st_ind_start = length(s_day_list);

% Create filelist: Terra
folder_t = [rootdir 'MODIS_Terra_PAR/'];
dir2 = dir(folder_t);
lf=length(dir2);
filelist_t={};
for i=1:lf
    nf=dir2(i).name; tf=contains(nf,'.nc');  % is it an nc file
    if tf==1 % if this IS an nc file 
        filelist_t=[filelist_t;nf];
    end
end; clear lf nf tf

for i = 1:length(st_day_list)
    thisday = datestr(st_day_list(i),'yyyymmdd');
    % load in the file that corresponds with this date, using "contains"
    ind_thisday_s = find(contains(filelist_s,thisday));
    ind_thisday_t = find(contains(filelist_t,thisday));
    filename_s = char(filelist_s(ind_thisday_s));
    filename_t = char(filelist_t(ind_thisday_t));
    
    % SeaWiFS
    if isempty(filename_s)
        par_s = nan(300,120);
    else
        file1 = [folder_s filename_s];
        ds = ncreadatt(file1,'/','product_name');
        datestr1 = [ds(21:24) '-' ds(25:26) '-' ds(27:28)]; % SeaWiFS (different for SeaWiFS, Aqua, and MODIS)
        par_s = ncread(file1,'par');
    end

    % Terra
    if isempty(filename_t)
        par_coarser = nan(300,120);
    else
        file2 = [folder_t filename_t];
        ds = ncreadatt(file2,'/','product_name');
        datestr2 = [ds(13:16) '-' ds(17:18) '-' ds(19:20)]; % Terra (different for SeaWiFS, Aqua, and MODIS)
        par_t = ncread(file2,'par');
        if i ==1
            lon0 = ncread(file2,'lon'); 
            lat0 = ncread(file2,'lat'); 
            lon_par_t = repmat(lon0,1,size(par_t,2));
            lat_par_t = repmat(lat0,1,size(par_t,1))';
        end
        par_coarser = griddata(lon_par_t,lat_par_t,par_t,lon_par_s,lat_par_s);
    end
    
    data_thisday = cat(3,par_s,par_coarser);
    par_avg = nanmean(data_thisday,3);
    par_database(:,:,i + st_ind_start) = par_avg;
end           
disp('Done with Terra + SeaWiFS data chunk.');      

%% Average of SeaWiFS, Terra, and Aqua 2002-07-04 to 2010-12-10
sta_day_list = datenum(2002,7,4):datenum(2010,12,10);
sta_ind_start = length(s_day_list) + length(st_day_list);

% Create filelist: Aqua
folder_a = [rootdir 'MODIS_Aqua_PAR/'];
dir3 = dir(folder_a);
lf=length(dir3);
filelist_a={};
for i=1:lf
    nf=dir3(i).name; tf=contains(nf,'.nc');  % is it an nc file
    if tf==1 % if this IS an nc file 
        filelist_a=[filelist_a;nf];
    end
end; clear lf nf tf


for i = 1:length(sta_day_list)
    thisday = datestr(sta_day_list(i),'yyyymmdd');
    % load in the file that corresponds with this date, using "contains"
    ind_thisday_s = find(contains(filelist_s,thisday));
    ind_thisday_t = find(contains(filelist_t,thisday));
    ind_thisday_a = find(contains(filelist_a,thisday));
    filename_s = char(filelist_s(ind_thisday_s));
    filename_t = char(filelist_t(ind_thisday_t));
    filename_a = char(filelist_a(ind_thisday_a));
    
    % SeaWiFS
    if isempty(filename_s)
        par_s = nan(300,120);
    else
        file1 = [folder_s filename_s];
        ds = ncreadatt(file1,'/','product_name');
        datestr1 = [ds(21:24) '-' ds(25:26) '-' ds(27:28)]; % SeaWiFS (different for SeaWiFS, Aqua, and MODIS)
        par_s = ncread(file1,'par');
    end
    
    % Terra
    if isempty(filename_t)
        par_coarser = nan(300,120);
    else
        file2 = [folder_t filename_t];
        ds = ncreadatt(file2,'/','product_name');
        datestr2 = [ds(13:16) '-' ds(17:18) '-' ds(19:20)]; % Terra (different for SeaWiFS, Aqua, and MODIS)
        par_t = ncread(file2,'par');
        if i ==1
            lon0 = ncread(file2,'lon'); 
            lat0 = ncread(file2,'lat'); 
            lon_par_t = repmat(lon0,1,size(par_t,2));
            lat_par_t = repmat(lat0,1,size(par_t,1))';
        end
        par_coarser = griddata(lon_par_t,lat_par_t,par_t,lon_par_s,lat_par_s);
    end
    
    % Aqua
    if isempty(filename_a)
        par_coarser_a = nan(300,120);
    else
        file3 = [folder_a filename_a];
        ds = ncreadatt(file3,'/','product_name');
        datestr3 = [ds(12:15) '-' ds(16:17) '-' ds(18:19)]; % Aqua (different for SeaWiFS, Aqua, and MODIS)
        par_a = ncread(file3,'par');
        if i ==1
            lon0 = ncread(file2,'lon'); 
            lat0 = ncread(file2,'lat'); 
            lon_par_a = repmat(lon0,1,size(par_a,2));
            lat_par_a = repmat(lat0,1,size(par_a,1))';
        end
        par_coarser_a = griddata(lon_par_a,lat_par_a,par_a,lon_par_s,lat_par_s);
    end

    data_thisday = cat(3,par_s,par_coarser,par_coarser_a);
    par_avg = nanmean(data_thisday,3);
    par_database(:,:,i + sta_ind_start) = par_avg;
end
disp('Done with Aqua + Terra + SeaWiFS data chunk.');    

%% Average of Terra and Aqua 2010-12-11 to 2022-06-31
ta_day_list = datenum(2010,12,11):datenum(2022,6,31);
ta_ind_start = length(s_day_list) + length(st_day_list) + length(sta_day_list);

for i = 1:length(ta_day_list)
    thisday = datestr(ta_day_list(i),'yyyymmdd');
    % load in the file that corresponds with this date, using "contains"
    ind_thisday_t = find(contains(filelist_t,thisday));
    ind_thisday_a = find(contains(filelist_a,thisday));
    filename_t = char(filelist_t(ind_thisday_t));
    filename_a = char(filelist_a(ind_thisday_a));
    
    % Terra
    if isempty(filename_t)
        par_coarser = nan(300,120);
    else
        file2 = [folder_t filename_t];
        ds = ncreadatt(file2,'/','product_name');
        datestr2 = [ds(13:16) '-' ds(17:18) '-' ds(19:20)]; % Terra (different for SeaWiFS, Aqua, and MODIS)
        par_t = ncread(file2,'par');
        if i ==1
            lon0 = ncread(file2,'lon'); 
            lat0 = ncread(file2,'lat'); 
            lon_par_t = repmat(lon0,1,size(par_t,2));
            lat_par_t = repmat(lat0,1,size(par_t,1))';
        end
        par_coarser = griddata(lon_par_t,lat_par_t,par_t,lon_par_s,lat_par_s);
    end
    
    % Aqua
    if isempty(filename_a)
        par_coarser_a = nan(300,120);
    else
        file3 = [folder_a filename_a];
        ds = ncreadatt(file3,'/','product_name');
        datestr3 = [ds(12:15) '-' ds(16:17) '-' ds(18:19)]; % Aqua (different for SeaWiFS, Aqua, and MODIS)
        par_a = ncread(file3,'par');
        if i ==1
            lon0 = ncread(file2,'lon'); 
            lat0 = ncread(file2,'lat'); 
            lon_par_a = repmat(lon0,1,size(par_a,2));
            lat_par_a = repmat(lat0,1,size(par_a,1))';
        end
        par_coarser_a = griddata(lon_par_a,lat_par_a,par_a,lon_par_s,lat_par_s);
    end
    
    data_thisday = cat(3,par_coarser,par_coarser_a);
    par_avg = nanmean(data_thisday,3);
    par_database(:,:,i + ta_ind_start) = par_avg;
end
disp('Done with Aqua + Terra (FINAL) data chunk.');    

%%
lon_par = lon_par_s;
lat_par = lat_par_s;

whos par_database
save('par_data_20230221.mat','par_database','day_list','lon_par','lat_par','-v7.3');

% %% Test: plot data
% mapfolder = '/Users/jsturner/satellite_data/matlab_WAP/WAP_map_SchuylerJack/';
% load([mapfolder 'bathy_and_coast.mat'],'coastlat','coastlon');
% coastcolor = [1 1 1]; 
% axislims=[-79.9792  -55.0208  -69.9792  -60.0208];
% 
% n = 6690;
% par_ex = squeeze(par_database(:,:,n));
% date_string = datestr(day_list(n),'yyyy-mm-dd');
% 
% figure 
% clf
% pcolor(lon_par_s,lat_par_s,par_ex); shading flat;
% colormap(turbo(31));
% makepretty; grid off;
% set(gca,'fontsize',12,'xtick','','ytick','','linewidth',1.5);
% caxis([0 85]);
% cw = colorbar;
% cw.LineWidth = 1;
% cw.Label.String = 'PAR (einsteins m^{-2} day^{-1})';
% cw.Location = 'southoutside';
% cw.FontSize = 16;
% title(date_string,'interpreter','none','fontsize',16);
% for i = 1:length(coastlon)
%     patch(coastlon{i},coastlat{i},coastcolor,'LineWidth',1);
% end
% 
% 


