clc
clear all
close all

storm_name = 'Idai';
experiment_type = 'experiment';
data_path = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' storm_name '/' experiment_type];
data_files = dir(data_path);



lat_lon_file = '/discover/nobackup/projects/smaptc/jkolassa/archive/Idai/Idai_SMAP_4D/diag/Y2019/M01/Idai_SMAP_4D.tavg1_2d_ocn_Nx.20190131_1030z.nc4';
lat_data = double(ncread(lat_lon_file,'lat'));
lon_data = double(ncread(lat_lon_file,'lon'));
frocean_data = double(ncread(lat_lon_file,'FRLAND'));

month_vec = string(['Jan';'Feb';'Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec']);

bt_data_new = [];

for f = 3:(length(data_files))

    bt_file = [data_path '/' data_files(f).name]
    bt_data = readmatrix(bt_file);
    size(bt_data)
%    if size(bt_data,2)==8 
%       bt_data(:,8) = [];
%    end

    datenum_write = datetime(bt_data(1,2),bt_data(1,3),bt_data(1,4),bt_data(1,5),0,0);

    max_traj = max(bt_data(:,1));

    for t = 1:max_traj
                
        traj_index = (bt_data(:,1)==t);
        bt_data_traj = bt_data(traj_index,:);
        bt_data_traj_new = bt_data_traj;

        tmp= (bt_data_traj(:,5)<0);
        if sum(tmp>0)
           bt_file
           t
        end

        datenum_start = datetime(bt_data_traj(1,2),bt_data_traj(1,3),bt_data_traj(1,4),bt_data_traj(1,5),0,0);

        for i = 2:size(bt_data_traj,1)
 
            datenum_new = datenum_start - (i-1)*3/24;
            datenum_str = char(datenum_new);

            bt_data_traj_new(i,2) = str2double(datenum_str(8:11));
            tmp = string(datenum_str(4:6));
            str_index = strcmp(month_vec,tmp);
            month_index = find(str_index==1);
            bt_data_traj_new(i,3) = month_index;
            bt_data_traj_new(i,4) = str2double(datenum_str(1:2));
            if length(datenum_str) > 11
               bt_data_traj_new(i,5) = str2double(datenum_str(13:14));
            end                      
            if (bt_data_traj_new(i,5)<0)
                datenum_str
            end
            bt_data_traj_new(i,8) = (i-1)*3;
        end % for i
        bt_data_new = [bt_data_new; bt_data_traj_new];
    end % for t

    if (f==3)
       start_index = 1;
    else
       start_index = stop_index + 1;
    end
    stop_index = size(bt_data_new,1);

    bt_data_tmp = bt_data_new(start_index:stop_index,:);

    ocean_mask_mat = zeros(size(frocean_data));
    ocean_mask_heat_mat = zeros(size(frocean_data));
    ocean_mask_mat_tr = zeros(size(frocean_data,1),size(frocean_data,2),max(bt_data_tmp(:,1)));
    ocean_mask_time_mat = NaN*ones(size(frocean_data));
    disp('land mask time mat 1')
    size(ocean_mask_time_mat)
    size(ocean_mask_mat)

    traj_geo = NaN*ones(max(bt_data_tmp(:,1)),2);
    traj_geo_index = 0;

    for j = 1:size(bt_data_tmp,1)

         if (j == 1)
            traj_geo_index = 1;
            traj_geo(1,:) = [bt_data_tmp(j,7) bt_data_tmp(j,6)];
         elseif (j>1) && (bt_data_tmp(j,1)>bt_data_tmp(j-1,1))
            traj_geo_index = traj_geo_index + 1;
            traj_geo(traj_geo_index,:) = [bt_data_tmp(j,7) bt_data_tmp(j,6)];
         end % if j

        if (bt_data_tmp(j,8)<= 72)
           lat_diff = sqrt((lat_data-bt_data_tmp(j,7)).^2);
           lat_index = find(lat_diff==min(lat_diff));
           lon_diff = sqrt((lon_data-bt_data_tmp(j,6)).^2);
           lon_index = find(lon_diff==min(lon_diff));

           if (frocean_data(lon_index,lat_index)<0.5)
              ocean_mask_mat(lon_index,lat_index) = 1;
              ocean_mask_mat_tr(lon_index,lat_index,bt_data_tmp(j,1)) = 1;
              ocean_mask_heat_mat(lon_index,lat_index) = ocean_mask_heat_mat(lon_index,lat_index) + 1;
              ocean_mask_time_mat(lon_index,lat_index) = ceil(bt_data_tmp(j,8)/24);
           end
        end % if bt_data_traj_new
    end % for j

    sum(ocean_mask_mat(:))

    tmp = char(datenum_start)
    tmp(12) = '_';

    save_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/Idai/land_mask/ocean_mask_2_72h_' experiment_type '_' tmp '.mat'];
    save(save_file,'ocean_mask_mat','ocean_mask_mat_tr','traj_geo','ocean_mask_heat_mat','ocean_mask_time_mat');
   
    ocean_mask_mat_new = ocean_mask_mat;
    ocean_mask_heat_mat_new = ocean_mask_heat_mat;
    ocean_mask_time_mat_new = ocean_mask_time_mat;


    for i = 1:size(ocean_mask_mat,1)
        for j = 1:size(ocean_mask_mat,2)

            if ocean_mask_mat(i,j)==1  
               n_lon_index = [i-1 i i+1];
               n_lon_index(n_lon_index<1) = [];
               n_lon_index(n_lon_index>size(ocean_mask_mat,1)) = [];

               n_lat_index = [j-1 j j+1];
               n_lat_index(n_lat_index<1) = [];
               n_lat_index(n_lat_index>size(ocean_mask_mat,2)) = [];

               for nlo = 1:length(n_lon_index)
                   for nla = 1:length(n_lat_index)
                       ocean_mask_mat_new(n_lon_index(nlo),n_lat_index(nla)) = 1;
                       ocean_mask_heat_mat_new(n_lon_index(nlo),n_lat_index(nla)) = ...
                       ocean_mask_heat_mat_new(n_lon_index(nlo),n_lat_index(nla)) + 1;
                       ocean_mask_time_mat_new(n_lon_index(nlo),n_lat_index(nla)) = ocean_mask_time_mat(i,j);
                   end
               end
             end % if
        end
     end

    ocean_mask_mat = ocean_mask_mat_new;
    ocean_mask_heat_mat = ocean_mask_heat_mat_new;
    ocean_mask_time_mat = ocean_mask_time_mat_new;

    save_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/Idai/land_mask/ocean_mask_interp_72h_' experiment_type '_' tmp '.mat'];
    save(save_file,'ocean_mask_mat','ocean_mask_heat_mat','ocean_mask_time_mat');


end % for f


datenum_vec = datenum(bt_data_new(:,2), bt_data_new(:,3),bt_data_new(:,4),bt_data_new(:,5),0,0);
bt_data_new = [bt_data_new datenum_vec];

%bt_data_interp = [];
%
%for i = 1:(size(bt_data_new,1)-1)
%
%    if (bt_data_new(i,1) ~= bt_data_new(i+1,1))
%       continue
%    end
%
%    datenum_start = datetime(bt_data_new(i,2),bt_data_new(i,3),bt_data_new(i,4),bt_data_new(i,5),0,0);
%    datenum_end = datetime(bt_data_new(i+1,2),bt_data_new(i+1,3),bt_data_new(i+1,4),bt_data_new(i+1,5),0,0);
%    datenum_new = datenum_start:(-0.5/24):datenum_end;    
%
%    date_vec_orig = bt_data_new(i:(i+1),8);
%    date_vec_new = bt_data_new(i,8):(-0.5/24):(bt_data_new(i+1,8)+(0.5/24));
%
%    lon_interp = interp1(date_vec_orig,bt_data_new(i:(i+1),6),date_vec_new,'linear');
%    lat_interp = interp1(date_vec_orig,bt_data_new(i:(i+1),7),date_vec_new,'linear');
%
%    bt_interp = repmat(bt_data_new(i,:),(length(datenum_new)-1),1);
%
%    datenum_str = char(datenum_new);
%
%    for j = 1:(length(datenum_new)-1)
%        bt_interp(j,2) = str2double(datenum_str(j,8:11));
%        bt_interp(j,4) = str2double(datenum_str(j,1:2));
%        bt_interp(j,5) = str2double(datenum_str(j,13:14));
%        tmp = string(datenum_str(j,4:6));
%        str_index = strcmp(month_vec,tmp);
%        month_index = find(str_index==1);
%        bt_interp(j,3) = month_index;
%    end
%
%    bt_interp(:,6) = lon_interp;
%    bt_interp(:,7) = lat_interp;
%    bt_interp(:,8) = date_vec_new;
%
%    bt_data_interp = [bt_data_interp; bt_interp];
%end
%
%lat_lon_file = '/discover/nobackup/projects/smaptc/jkolassa/archive/Idai/Idai_SMAP_4D/diag/Y2019/M01/Idai_SMAP_4D.tavg1_2d_ocn_Nx.20190131_1030z.nc4';
%lat_data = double(ncread(lat_lon_file,'lat'));
%lon_data = double(ncread(lat_lon_file,'lon'));
%frocean_data = double(ncread(lat_lon_file,'FRLAND'));
%
%unique_datenum = unique(bt_data_interp(:,8));
%
%for i = 1:length(unique_datenum)
%
%    date_index = (bt_data_interp(:,8) == unique_datenum(i));
%    bt_data_datenum = bt_data_interp(date_index,:);
%
%    land_mask_mat = zeros(size(frocean_data));
% 
%    for j = 1:size(bt_data_datenum)
%
%        lat_diff = sqrt((lat_data-bt_data_datenum(j,7)).^2);
%        lat_index = find(lat_diff==min(lat_diff));
%        lon_diff = sqrt((lon_data-bt_data_datenum(j,6)).^2);
%        lon_index = find(lon_diff==min(lon_diff));
%
%        if (frocean_data(lon_index,lat_index)>=0.5)
%           land_mask_mat(lon_index,lat_index) = 1;
%        end
%    end
%
%    sum(land_mask_mat(:))
% 
%    tmp = char(datestr(unique_datenum(i)))
%    tmp(12) = '_'; 
%
%    save_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/Idai/land_mask/land_mask_' tmp '.mat'];
%    save(save_file,'land_mask_mat');
%end
