close all
clear

%% Postprocessing_plotMigrationRate.m
% original by Veiko Lehsten
% modified by Deborah Zani


%% Content (see also Lehsten et al. 2019, legend of Figure 4)
% Plot spread of each species through an area of 201x201 grid cells with static climate and corridors at 200km. 
% Figure 2: the left panels display the time when each species first reached an LAI of 0.5. 
% Each species is allowed to establish freely only in the upper left corner. 
% The right panels show the distance of the grid cells with LAI>0.5 from the starting point. 
% The red line indicates the 95th percentile of the grid cells farthest away from the starting point. 
% The migration speed is calculated as slope of this line, taking only grid cells at least 5 km away 
% from the starting point into account to avoid some initial establishing effects.


%% Define variables

% LAI threshold value
threshold_value_for_plot=0.5;

% Select column with focal species (lai.out)
column=[5]; %N.B.: 4 for Abi_alb!

% Define years to plot
startyear_pre=50;
endyear=10420;
yearstep=10 ;

% Initialize
if_corner=1;
plot_lines=1;
tic


%% Read-in and order output file

% Read LPJ-GM output (lai.out)
% Define path to lai.out file
build_dir='/.../' ;
comb_tab = dlmread([build_dir  'lai.out'],'',1,0);

% Now sort into cells
clear LAI
clear LAI_year_check
cells=unique(comb_tab(:,[1 2 ]),'rows');
comb_tab(:,[3])=comb_tab(:,[3])-120;
years=unique(comb_tab(:,[3]),'rows');
LAI_sort=sortrows(comb_tab,[1,2,3]);
for ce_i=1:size(cells,1)
    for col=[4 5 6]
        try
            LAI(ce_i,col-3,:)=comb_tab(comb_tab(:,1)==cells(ce_i,1) & comb_tab(:,2)==cells(ce_i,2),col);
        catch
            'stop'
        end
    end
end


%% Calculate distance matrix to plot
wri_time=0;
exp_time=0;
down_time=0;

for iii=1:prod(size(column))
    
    % Define colorbar
    maxc=max(comb_tab(:,column(iii)))+0.2*max(comb_tab(:,column(iii)));
    minc=min(comb_tab(:,column(iii)))-0.2*min(comb_tab(:,column(iii)));
    
    % Initialize time-step for plotting
    step=0; 
    startyear=max(startyear_pre,min(comb_tab(:,3)));
    for year=startyear:yearstep:min(endyear,max(comb_tab(:,3)))
        comb_tab=sortrows(comb_tab,3);
        step=step+1;
        year
        all_years(step)=year;
        try
            thisyear=comb_tab(comb_tab(:,3)==year,:);
            before_down=toc;
            close all
            [z,R]= downscale('read_only','data_in',comb_tab,'year',year,'column',column(iii),'jemeter','plot');
            z_all(:,:,step)=z;    
            down_time=down_time+toc-before_down;
            before_exp=toc;
            caxis([minc,maxc]);
            title(num2str(year));
        catch
            'failure' 
        end
    end
end

% Find the years when LAI>0.5 
z_allflip=permute(z_all,[2 1 3]);
%for yeari=1:size(z_all,3)
%    z_allflip(:,:,yeari)=flipdim(z_all(:,:,yeari),2);
%end
thresholded=double(z_allflip>0.5);
for yeari=1:size(thresholded,3)
    thresh_years(:,:,yeari)=thresholded(:,:,yeari)*(all_years(yeari));
end
N_value=max(max(max(thresh_years)))+1000;
thresh_years(thresh_years==0)=N_value;
minyear=min(thresh_years,[],3);

% Make distance matrix
for line_i=1:size(minyear,1)
    for col_i=1:size(minyear,2)
        if if_corner
            distma(line_i,col_i)=(line_i*line_i+col_i*col_i)^(0.5);
        else
            distma(line_i,col_i)=col_i;
        end
    end
end
plo_minyear=minyear(:);
plo_minyear_dist=distma(:);
plo_minyear_dist(plo_minyear==N_value)=[];
plo_minyear(plo_minyear==N_value)=[];


%% Plot
figure
subplot (1,2,1)
minyear(minyear==N_value)=NaN;
imagesc(minyear)
colorbar
ylabel('Distance from refugia in km');
xlabel('Distance from refugia in km');
title('Year of arrival');

c=colormap(jet);
c(1,:)=[1,1,1];
colormap(c)
subplot (1,2,2)
hold on
minyear(minyear==N_value)=NaN;
plot (plo_minyear,plo_minyear_dist,'k.');

% Take out the first 5 km
plo_minyear(plo_minyear_dist<5)=[];
plo_minyear_dist(plo_minyear_dist<5)=[];
p=polyfit(plo_minyear,plo_minyear_dist,1);

try
    clear percent_max test_set_km
    pos_array=0;
    for km_i=1:max(plo_minyear_dist)
        time_km_temp=plo_minyear(plo_minyear_dist>=km_i-0.5 & plo_minyear_dist<=km_i+0.5);
        pos_array=pos_array+1;
        percent_max(pos_array)=prctile(time_km_temp,5);
        test_set_km(pos_array)=km_i;
    end
    % Remove nans
    test_set_km(isnan(percent_max))=[];
    percent_max(isnan(percent_max))=[];

    p_1=polyfit(percent_max,test_set_km,1);
    plot(percent_max,percent_max*p_1(1)+p_1(2),'r');
    title(['Migration front (speed: ' num2str(round(p_1(1)*1000)) ' m/year)']);
    set(gcf, 'Position', [100, 200, 900, 300])

    set(gcf, 'color', [1 1 1])
    xlabel('Year');
    ylabel('Distance from refugia in km');
    ylim([0 100])
    drawnow
    hold off
    %export_fig([upperfolder strrep(strrep(under_dir,'.','_'),'/','') '.png'], '-png','-r300'); % low
    %close all
catch
end


%% References
 % Lehsten, V., Mischurow, M., Lindström, E., Lehsten, D., and Lischke, H.: 
 % LPJ-GM 1.0: simulating migration efficiently in a dynamic vegetation model. 
 % Geosci. Model Dev., 12, 893–908, 
 % https://doi.org/10.5194/gmd-12-893-2019, 2019.