%% PLOT, WITH POINT-WISE 90% CREDIBLE INTERVALS, ESTIMATES OF THE
% APRIL-SEPTEMBER TEMPERATURE ANOMALY FIELD, FOR FOUR YEARS. 

% file loads Temperature_Array_3pers and BARCAST_INPUT_vHNL_Jan2012

%% Set the four years, between 1400 and 2011 (inclusive)
%years=[1453, 1601, 1642, 1695];
%years=[1998, 2003, 2007, 2009];
%years=[1814, 1815, 1816, 1817];
years=[1840:30:1930];

%CHANGE NOTHING AFTER THIS. 

%% load 
load Temperature_Array_3pers %Temperature_Array_3pers

load Inputs/BARCAST_INPUT_vHNL_Jan2012 %BARCAST_INPUT
timeline=BARCAST_INPUT.Data_timeline;

%we want to extract the elements of the temperature array that correspond
%to the grid points, viz a viz. Also eliminate the first year (1399), prior
%to any observations. 

SpaceTime_3pers=Temperature_Array_3pers(BARCAST_INPUT.Inds_Central,2:1:end,:);

Locs=BARCAST_INPUT.Master_Locs(BARCAST_INPUT.Inds_Central,:);

% PLOT the field at two years, along with the uncertainty. 
Plot_Inds=years-1399;

N_Plots=length(Plot_Inds);
Full_SpaceTime_Matrix=squeeze(SpaceTime_3pers(:, Plot_Inds, 2));
Full_CI_Widths_Matrix=squeeze(SpaceTime_3pers(:, Plot_Inds, 3))-squeeze(SpaceTime_3pers(:, Plot_Inds, 1));

%figure out how to plot colours on the squares. 
Lat_pred_range=[45, 90];
Lon_pred_range=[-180, 180];

%find the locations that will be shown that will always be NaNs:
Pos_locs=[kron([Lon_pred_range(1)+2.5:5:Lon_pred_range(2)-2.5]', ones(size([Lat_pred_range(1)+2.5:5:Lat_pred_range(2)-2.5]'))) , kron(ones(size([Lon_pred_range(1)+2.5:5:Lon_pred_range(2)-2.5]')), [Lat_pred_range(1)+2.5:5:Lat_pred_range(2)-2.5]')];
Empty_locs=setdiff(Pos_locs, Locs, 'rows');

%set colormap
load HotCoolColormap_LONG
ccmap=HotCoolColormap_LONG;

%color to plot errors larger than BigErrCUT
BigErrCol=[0.85, 0.0, 0.85];
BigErrCUT=5.5;

ColLims_T=max(max(abs(SpaceTime_3pers(:, :, 2))))*[-1, 1];
ColLims_W=[min(min(squeeze(SpaceTime_3pers(:, :, 3))-squeeze(SpaceTime_3pers(:, :, 1)))),min([BigErrCUT, max(max(squeeze(SpaceTime_3pers(:, :, 3))-squeeze(SpaceTime_3pers(:, :, 1))))])];

cN=length(ccmap);
colormap(ccmap)

figure(10), clf
colormap(ccmap)

fig=figure(10);
set(fig,'DoubleBuffer','on', 'color',  [1.0 1.0 1.0], 'position', [10,10,1350,900], 'renderer', 'painters');
    
subplot(2, N_Plots, 1)
set(gca, 'fontSize', 15)
plot(-1000, -1000, 'g^',  'color', [0.85, 0.1, 0.1], 'markersize', 5.5, 'linewidth', 1.5), hold on %'markerFaceColor',[0.3, 0.8, 0.3]), hold on
plot(-1000, -1000, 'wo', 'color', [0.70, 0.70, 0.70], 'markersize', 4, 'markerFaceColor',[0.70, 0.70, 0.70]), hold on
plot(-1000, -1000, 'g+',  'color', [0.85, 0.0, 0.85], 'markersize', 5, 'linewidth', 1.5), hold on %, 'markerFaceColor',[0.85, 0.0, 0.85]), hold on
leg=legend('Proxy', 'Inst.', 'Both', 'location', [0.9 0.4725 0.0800 0.0907]);

for LL=1:1:N_Plots
    subplot(2, N_Plots, LL)
    set(gca, 'fontSize', 15)

    m_proj('stereographic','long',-60,'lat',90, 'radius', 45), hold on
    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on

    
    %Find the NON NaN entries for the year:
    g_inds=find(isnan(Full_SpaceTime_Matrix(:, LL))==0);    
    for kk=1:1:length(g_inds);
        indi=max(1, ceil((Full_SpaceTime_Matrix(g_inds(kk), LL)-ColLims_T(1))/diff(ColLims_T)*cN));
        lon_mat=[(+2.5+Locs(g_inds(kk),1))*ones(2,1);(-2.5+Locs(g_inds(kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Locs(g_inds(kk),2))*ones(1,1);(+2.5+Locs(g_inds(kk),2))*ones(2,1);(-2.5+Locs(g_inds(kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), ccmap(indi,:), 'edgecolor', ccmap(indi,:)), hold on
        m_plot([lon_mat, lon_mat], [lat_mat, lat_mat], 'k-', 'Color', ccmap(indi,:)), hold on
    end

    %and the NaN entries
    for kk=1:1:length(Empty_locs(:,1));
        lon_mat=[(+2.5+Empty_locs((kk),1))*ones(2,1);(-2.5+Empty_locs((kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Empty_locs((kk),2))*ones(1,1);(+2.5+Empty_locs((kk),2))*ones(2,1);(-2.5+Empty_locs((kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), [0.7, 0.7, 0.7], 'edgecolor', [0.7, 0.7, 0.7]), hold on
    end

    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on
    
    
    hh=title(['Year = ',num2str(timeline(Plot_Inds(LL)))], 'Position', [+0, +1, 15], 'fontsize', 18);
    
    if LL==1
        caxis('manual')
        caxis(ColLims_T)
        pp=colorbar('vert', 'Position', [0.93 0.61 0.0133 0.2946], 'fontsize', 14, 'XAxisLocation', 'bottom', 'YAxisLocation', 'right');
        tt=text(-1.2, -0.7,'Posterior Median', 'fontsize', 18, 'rotation', 90);
    end
    
    %for the uncertainty
    subplot(2, N_Plots, N_Plots+LL)
    set(gca, 'fontSize', 15)
    m_proj('stereographic','long',-60,'lat',90, 'radius', 45), hold on
    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on

    %Find the NON NaN entries for the year, small uncertainty:
    g_inds=find(isnan(Full_CI_Widths_Matrix(:, LL))==0 & Full_CI_Widths_Matrix(:, LL)<BigErrCUT);    
    for kk=1:1:length(g_inds);
        indi=max(1, ceil((Full_CI_Widths_Matrix(g_inds(kk), LL)-ColLims_W(1))/diff(ColLims_W)*cN));
        lon_mat=[(+2.5+Locs(g_inds(kk),1))*ones(2,1);(-2.5+Locs(g_inds(kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Locs(g_inds(kk),2))*ones(1,1);(+2.5+Locs(g_inds(kk),2))*ones(2,1);(-2.5+Locs(g_inds(kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), ccmap(indi,:), 'edgecolor', ccmap(indi,:)), hold on
    end

    %Find the NON NaN entries for the year, large uncertainty:
    g_inds=find(isnan(Full_CI_Widths_Matrix(:, LL))==0 & Full_CI_Widths_Matrix(:, LL)>=BigErrCUT);
    for kk=1:1:length(g_inds);
        indi=max(1, ceil((Full_CI_Widths_Matrix(g_inds(kk), LL)-ColLims_W(1))/diff(ColLims_W)*cN));
        lon_mat=[(+2.5+Locs(g_inds(kk),1))*ones(2,1);(-2.5+Locs(g_inds(kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Locs(g_inds(kk),2))*ones(1,1);(+2.5+Locs(g_inds(kk),2))*ones(2,1);(-2.5+Locs(g_inds(kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), [0.4, 0.7, 0.4], 'edgecolor', [0.4, 0.7, 0.4]), hold on
    end


    %and the always NaN locations
    for kk=1:1:length(Empty_locs(:,1));
        lon_mat=[(+2.5+Empty_locs((kk),1))*ones(2,1);(-2.5+Empty_locs((kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Empty_locs((kk),2))*ones(1,1);(+2.5+Empty_locs((kk),2))*ones(2,1);(-2.5+Empty_locs((kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), [0.7, 0.7, 0.7], 'edgecolor', [0.7, 0.7, 0.7]), hold on
    end
    
    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on
    
    %find where we have data 
    I_g_I=find(isnan(BARCAST_INPUT.Inst_Data(:, Plot_Inds(LL)))==0);
    I_g_P1=find(isnan(BARCAST_INPUT.Prox_Data1(:, Plot_Inds(LL)))==0);
    I_g_P2=find(isnan(BARCAST_INPUT.Prox_Data2(:, Plot_Inds(LL)))==0);
    I_g_P3=find(isnan(BARCAST_INPUT.Prox_Data3(:, Plot_Inds(LL)))==0);
    I_g_P=unique([I_g_P1;I_g_P2;I_g_P3]); 
    
     %only the Inst
     [Only_I, inds_Only_I]=setdiff(I_g_I, I_g_P);
     %only the Prox
     [Only_P, inds_Only_P]=setdiff(I_g_P, I_g_I);
     %Both
     [Both, inds_Both]=intersect(I_g_P, I_g_I);
    
     %plot
     %only the Prox
     m_plot(BARCAST_INPUT.Master_Locs(Only_P,1), BARCAST_INPUT.Master_Locs(Only_P,2), 'g^',  'color', [0.85, 0.1, 0.1], 'markersize', 5.5, 'linewidth', 1.5), hold on %'markerFaceColor',[0.3, 0.8, 0.3]), hold on
    %only the Inst
     m_plot(BARCAST_INPUT.Master_Locs(Only_I,1), BARCAST_INPUT.Master_Locs(Only_I,2), 'wo', 'color', [0.70, 0.70, 0.70], 'markersize', 4, 'markerFaceColor',[0.70, 0.70, 0.70]), hold on
    %Both
     m_plot(BARCAST_INPUT.Master_Locs(Both,1), BARCAST_INPUT.Master_Locs(Both,2), 'g+',  'color', [0.85, 0.0, 0.85], 'markersize', 5, 'linewidth', 1.5), hold on %, 'markerFaceColor',[0.85, 0.0, 0.85]), hold on
        


     if LL==1
         caxis('manual')
         caxis(ColLims_W)
%         pp=colorbar('vert', 'Position', [0.9    0.1175    0.03    0.8], 'fontsize', 14, 'XAxisLocation', 'bottom', 'YAxisLocation', 'right');
        pp=colorbar('vert', 'Position', [0.93 0.13 0.0133 0.2946], 'fontsize', 14, 'XAxisLocation', 'bottom', 'YAxisLocation', 'right');
        tt=text(-1.2, -0.85,'90% Credible Intervals', 'fontsize', 18, 'rotation', 90);
     end
    
end
