clc
clear all
close all

storm_name = 'Idai';
save_path = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/' storm_name '/'];
control_dir = ['/discover/nobackup/projects/smaptc/jkolassa/archive/' storm_name '/'];
control_name = 'LADAS_control_4D_new';
bt_time_str = '72h';   % string specifying how far back in time backtrajectories are computed that are used to create land_mask

lat_range = [-30 -10];
lon_range = [24 52];

start_year = 2019;
start_month = 3;
start_day = 8;

end_year = 2019;
end_month = 3;
end_day = 15;

leapYears = 1980:4:2022;
monthLength_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
monthLength_ly = [31 29 31 30 31 30 31 31 30 31 30 31];
month_vec = ['Jan';'Feb';'Mar';'Apr';'May';,'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec'];

% define latitude and longitude index for domain

coord_file = [control_dir control_name '/diag/Y' num2str(start_year,'%04i') '/M' num2str(start_month,'%02i') ...
             '/' control_name '.tavg1_2d_lnd_Nx.' num2str(start_year,'%04i') num2str(start_month,'%02i') ...
              num2str(start_day,'%02i') '_2130z.nc4'];

lat_vec = double(ncread(coord_file,'lat'));
lon_vec = double(ncread(coord_file,'lon'));

lat_index = (lat_vec<lat_range(1)) + (lat_vec>lat_range(2));
lat_index = lat_index == 0;
lon_index = (lon_vec<lon_range(1)) + (lon_vec>lon_range(2));
lon_index = lon_index == 0;

% read track data

track_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/track_files/' storm_name '_track.ascii'];
track_data = dlmread(track_file);
track_datenum = datenum(start_year, track_data(:,1), track_data(:,2), track_data(:,3), 0, 0);

% define color maps

cmap_mat_LM = cmocean('ice',12);
cmap_mat_LM = flipud(cmap_mat_LM);
caxis_min_LM = 0;
caxis_max_LM = 1;

cmap_mat_LMH = cmocean('thermal',12);
cmap_mat_LMH = flipud(cmap_mat_LMH);
caxis_min_LMH = 0;
caxis_max_LMH = 10;

% define plot domain
[lon_mat,lat_mat] = meshgrid(lon_vec,lat_vec);
load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);
latlim = lat_range;
lonlim = lon_range;

% loop over years
for y = start_year:end_year

    disp('read data')

    y_str = num2str(y,'%04i');

    % define month length based on leap year status
    if sum(y==leapYears)>0
        monthLength = monthLength_ly;
    else
        monthLength = monthLength_nly;
    end

    % loop over months
    for m = start_month:end_month

        m_str = num2str(m,'%02i');

        % define start and end day for current month
        day_first = 1;
        day_last = monthLength(m);

        if (m == start_month)
           day_first = start_day;
        end

        if (m == end_month)
           day_last = end_day;
        end

        % loop over days
        for d = day_first:day_last

            d_str = num2str(d,'%02i');
            datenum_vec = datenum(y,m,d);

            land_mask_file_control = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' ...
                                     storm_name '/land_mask/land_mask_2_' bt_time_str '_control_' d_str '-' month_vec(m,:) '-' y_str '_' ...
                                     '.mat'];

            land_mask_file_control_interp = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' ...
                                     storm_name '/land_mask/land_mask_interp_' bt_time_str '_control_' d_str '-' month_vec(m,:) '-' y_str '_' ...
                                     '.mat'];

            land_mask_file_exp = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' ...
                                 storm_name '/land_mask/land_mask_2_' bt_time_str '_experiment_' d_str '-' month_vec(m,:) '-' y_str '_' ...
                                 '.mat'];

            land_mask_file_exp_interp = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' ...
                                 storm_name '/land_mask/land_mask_interp_' bt_time_str '_experiment_' d_str '-' month_vec(m,:) '-' y_str '_' ...
                                 '.mat'];

            land_mask_exp = load(land_mask_file_exp);
            land_mask_control = load(land_mask_file_control);
            land_mask_exp_interp = load(land_mask_file_exp_interp);
            land_mask_control_interp = load(land_mask_file_control_interp);

            datenum_dist = sqrt((track_datenum-datenum_vec).^2);
            date_index = find(datenum_dist==min(datenum_dist));

            land_mask_exp.land_mask_mat(land_mask_exp.land_mask_mat==0) = NaN;
            land_mask_control.land_mask_mat(land_mask_control.land_mask_mat==0) = NaN;
            land_mask_exp.land_mask_heat_mat(land_mask_exp.land_mask_heat_mat==0) = NaN;
            land_mask_control.land_mask_heat_mat(land_mask_control.land_mask_heat_mat==0) = NaN;

            max(land_mask_exp.land_mask_heat_mat(:))
            max(land_mask_control.land_mask_heat_mat(:))
            figure
            ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
            framem
            gridm on;
            setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
            'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
            'MLabelParallel','south','FontSize',12)
            geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
            pcolorm(lat_vec,lon_vec,land_mask_exp.land_mask_mat')
            caxis([caxis_min_LM caxis_max_LM]);
            set(gcf,'Colormap',cmap_mat_LM);
            plotm(coastlat, coastlon,'k')
            hold on
            plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-m','LineWidth',2)
            hold off
            set(gca,'XLim',lonlim);
            set(gca,'YLim',latlim);
            tightmap
            title(['land mask EXP ' datestr(datenum_vec,'yyyy-mm-dd')])
            h = colorbar;
            set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
            set(gcf,'paperunits','centimeters')
            set(gcf,'paperposition',[1 1 17 15])
            saveas(gcf,[save_path storm_name '_land_mask_EXP_' datestr(datenum_vec,'yyyy-mm-dd')],'png')

            figure
            ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
            framem
            gridm on;
            setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
            'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
            'MLabelParallel','south','FontSize',12)
            geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
            pcolorm(lat_vec,lon_vec,land_mask_control.land_mask_mat')
            caxis([caxis_min_LM caxis_max_LM]);
            set(gcf,'Colormap',cmap_mat_LM);
            plotm(coastlat, coastlon,'k')
            hold on
            plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-m','LineWidth',2)
            hold off
            set(gca,'XLim',lonlim);
            set(gca,'YLim',latlim);
            tightmap
            title(['land mask CNTRL ' datestr(datenum_vec,'yyyy-mm-dd')])
            h = colorbar;
            set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
            set(gcf,'paperunits','centimeters')
            set(gcf,'paperposition',[1 1 17 15])
            saveas(gcf,[save_path storm_name '_land_mask_CNTRL_' datestr(datenum_vec,'yyyy-mm-dd')],'png')

            figure
            ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
            framem
            gridm on;
            setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
            'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
            'MLabelParallel','south','FontSize',12)
            geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
            pcolorm(lat_vec,lon_vec,land_mask_exp.land_mask_heat_mat')
            caxis([caxis_min_LMH caxis_max_LMH]);
            set(gcf,'Colormap',cmap_mat_LMH);
            plotm(coastlat, coastlon,'k')
            hold on
            plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-m','LineWidth',2)
            hold off
            set(gca,'XLim',lonlim);
            set(gca,'YLim',latlim);
            tightmap
            title(['land mask heat map EXP ' datestr(datenum_vec,'yyyy-mm-dd')])
            h = colorbar;
            set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
            set(gcf,'paperunits','centimeters')
            set(gcf,'paperposition',[1 1 17 15])
            saveas(gcf,[save_path storm_name '_land_mask_heat_EXP_' datestr(datenum_vec,'yyyy-mm-dd')],'png')

            figure
            ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
            framem
            gridm on;
            setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
            'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
            'MLabelParallel','south','FontSize',12)
            geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
            pcolorm(lat_vec,lon_vec,land_mask_control.land_mask_heat_mat')
            caxis([caxis_min_LMH caxis_max_LMH]);
            set(gcf,'Colormap',cmap_mat_LMH);
            plotm(coastlat, coastlon,'k')
            hold on
            plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-m','LineWidth',2)
            hold off
            set(gca,'XLim',lonlim);
            set(gca,'YLim',latlim);     
            tightmap
            title(['land mask heat map CNTRL ' datestr(datenum_vec,'yyyy-mm-dd')])
            h = colorbar;
            set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
            set(gcf,'paperunits','centimeters')
            set(gcf,'paperposition',[1 1 17 15])
            saveas(gcf,[save_path storm_name '_land_mask_heat_CNTRL_' datestr(datenum_vec,'yyyy-mm-dd')],'png')

            land_mask_exp_interp.land_mask_heat_mat(land_mask_exp_interp.land_mask_heat_mat==0) = NaN;

            figure
            ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
            framem
            gridm on;
            setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
            'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
            'MLabelParallel','south','FontSize',12)
            geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
            pcolorm(lat_vec,lon_vec,land_mask_exp_interp.land_mask_heat_mat')
            caxis([caxis_min_LMH caxis_max_LMH]);
            set(gcf,'Colormap',cmap_mat_LMH);
            plotm(coastlat, coastlon,'k')
            hold on
            plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-m','LineWidth',2)
            hold off
            set(gca,'XLim',lonlim);
            set(gca,'YLim',latlim);
            tightmap
            title(['land mask heat map EXP ' datestr(datenum_vec,'yyyy-mm-dd')])
            h = colorbar;
            set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
            set(gcf,'paperunits','centimeters')
            set(gcf,'paperposition',[1 1 17 15])
            saveas(gcf,[save_path storm_name '_land_mask_heat_interp_EXP_' datestr(datenum_vec,'yyyy-mm-dd')],'png')

            land_mask_control_interp.land_mask_heat_mat(land_mask_control_interp.land_mask_heat_mat==0) = NaN;

            figure
            ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
            framem
            gridm on;
            setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
            'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
            'MLabelParallel','south','FontSize',12)
            geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
            pcolorm(lat_vec,lon_vec,land_mask_control_interp.land_mask_heat_mat')
            caxis([caxis_min_LMH caxis_max_LMH]);
            set(gcf,'Colormap',cmap_mat_LMH);
            plotm(coastlat, coastlon,'k')
            hold on
            plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-m','LineWidth',2)
            hold off
            set(gca,'XLim',lonlim);
            set(gca,'YLim',latlim);
            tightmap
            title(['land mask heat map CNTRL ' datestr(datenum_vec,'yyyy-mm-dd')])
            h = colorbar;
            set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
            set(gcf,'paperunits','centimeters')
            set(gcf,'paperposition',[1 1 17 15])
            saveas(gcf,[save_path storm_name '_land_mask_heat_interp_CNTRL_' datestr(datenum_vec,'yyyy-mm-dd')],'png')


        end % d
     end % m
 end % y
