%% Plot global map and latitudinal cross-section of growing season averaged LRU
%
%  v02 adds the LRU histogram of Whelan et al. (2018) to the latitudinal
%  plot; in addition we mask C4 pixels in the map and omit them from the
%  calculation of the longitudinal means
%
%  last edit: 21.12.2022, Georg


%% clean up worspace and start timer
clear
close all
clc
tic


%% constants and settings
fs = 24; % font size


%% load colormap
disp('Loading data ...')
load batlow.mat % colormap


%% load coastlines for global map
load coastlines % coastlines


%% load pre-calculated vegetation period averaged LRU data
load global_vegperiod_lru_v02.mat


%% load pft map from Megan model
data = ncread('mksrf_landuse_rc2000_c110913.nc', 'PCT_PFT'); % percentage plant functional type
% PFTs are as in Guenther et al. (2012): O for bare ground, 1-15 plant PFTs
% and #16 is zero everywhere
% warm C4 grass is pft # 14
c4fraccov = flipud(data(:,:,15)'); % C4 fractional cover
clear data


%% load Whelan et al. (2018) histogram data
load whelan_lru_bin.mat
edges = (0:0.5:10)'; 
idx = discretize(edges_lru(1:end-1), edges);
bins_w = grpstats(bins_lru, idx, 'sum');
idx = unique(idx); 
edges_w = edges(idx); 
clear idx edges edges_lr bins_lru


%% plot global map and latitudinal cross-section 
clc; disp('Plotting data ...')
figure('Units','Normalized', 'Position', [0 0 1 1]);
intm = gp_lru;
idx = lm .* bm == 0;
sp(1) = axes('Position', [0.075 0.15 0.7 0.8]); 
    intm(idx) = NaN; % filter for land-water and bare-ground mask
    ih = imagesc(intm);
    c4fraccov(idx) = NaN;
    set(ih, 'AlphaData', 1 - c4fraccov .* 0.01)
    cb = colorbar; cb.Title.String = 'LRU (-)'; cb.Title.HorizontalAlignment = 'left';
    hold on
    plot((coastlon + 180) .* 20, (-1 .* coastlat + 90) .* 20, '-k')
    grid on; box on;
    xlabel('Longitude (°)'); ylabel('Latitude (°)')
    axis equal
    set(gca, 'Xlim', [0 7200], 'Xtick', 7200/360*0:7200/360*30:7200/360*360, ...
        'Xticklabels', -180:30:180, ...
        'Ylim', [0 3600], 'Ytick', 3600/180*0:3600/180*30:3600/180*180, ...
        'Yticklabels', 90:-30:-90, 'Fontsize', fs)
    caxis([0.6 2.2])
sp(2) = axes('Position', [0.835 0.27 0.15 0.56]);
    idx = c4fraccov > 50;
    intm(idx) = NaN; % filter for C4 fractional cover
    lru_avg = mean(intm, 2, 'omitnan'); 
    lru_std = std(intm,0 , 2, 'omitnan'); 
    shadedErrorBar(lat, lru_avg, lru_std, 'k');
    hold on
    base = 90;
    for i = 1:numel(bins_w)-1
        patch([base -bins_w(i)*2+base -bins_w(i)*2+base base], ...
            [edges_w(i) edges_w(i) edges_w(i+1) edges_w(i+1)], ...
            [0.5 0.5 0.5])
    end % for i
    plot([-90 90], [1.68 1.68], '--k', 'LineWidth', 2)
    ylabel('LRU (-)');
    grid on; box on;
    set(gca, 'Xlim', [-90 90], 'Xtick', -90:30:90, 'XtickLabels', 90:-30:-90, ...
        'Ylim', [0 4], 'Ytick', 0:1:4, 'view', [90 90], 'Fontsize', fs)
colormap(batlow)
pos1 = plotboxpos(sp(1));
pos2 = plotboxpos(sp(2));
set(sp(2), 'Position', [pos2(1), pos1(2), pos2(3), pos1(4)])

% 95 % percentile range and median
prctile(intm,[2.5 50 97.5], 'all')


%% done
toc