% This script loads data from Emile-Geay and al 2012 NINO3.4
% reconstructions and plots them along with radiative forcing
% over the last 850 years.

% it also illustrates the incredible usefullness of a proper database
% format, which is something I'd like to see the community move towards in
% the very near future, so that it does take years of human time to
% make use of the NCDC archive. 
%  
%  Matlab code necessary to reproduce the reconstructions will soon be
%  uploaded at: http://code.google.com/p/common-climate/
%
%  Julien Emile-Geay, USC 03/12/2012

%1) Loading, manipulating and plotting reconstructions
% =====================================================
load JEG_graphics   % graphical definitions
load NINO34_ttls_3target_2012_proxy_db  % reconstructions

% organised as follows: 
    %  (1) ERSSTv3 target 
    %  (2) HadSST2i target
    %  (3) Kaplan target
    
forcings = load('millennium_trop_forcings2011.mat');  % forcing

% extract 20-year lowpass filtered dual RegEM solution with HadSST2i target
index_reg = nino34_ttls{2};
index_cps = nino34_cps{2};
% Corrrelation
R = corr(index_cps.f20,index_reg.f20);


Fv = forcings.Fv.IVI.d;  tv = forcings.Fv.IVI.t;
Fs = forcings.Fs.Steinhilber.di+1366; % add So = 1366 baseline (best match with other estimates) 
ts = forcings.Fs.Steinhilber.ti;

fig('Forcing and Recon'),clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(tv,Fv,ts,Fs); grid, hold on
ylims=[-12 0]; xlims = [1150 1995];
set(H1,'Color',skyblue,'linewidth',[1]);
set(AX(1),'YLim',ylims,'XLim',xlims,'Ytick',[-12 -10 -8 -6 -4 -2 0],'YColor',skyblue);
hYLabel1=ylabel('Volcanic forcing (W.m^{-2}) ',style_l{:}); grid ;  
hold off
%
axes(AX(2)),  yli=[1364.5 1367];
set(H2,'Color',[0 0 1],'linewidth',[2]);
set(AX(2),'YLim',yli,'XLim',xlims, 'Ytick',[1365:1:1368],'YColor',[0 0 1]);

%  Font styles
hYLabel2=ylabel('Solar forcing (W.m^{-2})',style_l{:}); grid ; hold off
hTitle=title('Volcanic and Solar Forcing over the Tropics',style_t{:});
set(gca,style_a{:});
% Axes
set(AX(2),'Box', 'off', 'TickDir','out','TickLength',[.02 .02],'XMinorTick','on','YMinorTick','on', 'YGrid','on')
set(AX(2),'XColor' , [.3 .3 .3], 'LineWidth', 1);
set(AX(1),'Box', 'off', 'TickDir','out','TickLength',[.02 .02],'XMinorTick','on','YMinorTick','on', 'YGrid','off')
set(AX(1),'XColor' , [.3 .3 .3], 'LineWidth', 1);
set(AX(1),'Xticklabel',[]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NINO3 (uistack did not work)
subplot(2,1,2)
clear AX H1 H2
H2 = plot(t,index_reg.tot,'color',dkgr); % plot raw RegEM solution
set(H2,'color',dkgr,'linewidth',[1]);

set(gca,'xlim',xlims,'ylim',[-3, 3],'YAxisLocation','left','YColor',[.3 .3 .3]);
set(gca,'Box', 'off', 'TickDir','out','TickLength',[.02 .02],'XMinorTick','on','YMinorTick','on', 'YGrid','on')
set(gca,'XColor' , [.3 .3 .3], 'LineWidth', 1);
hold on
% plot 20-year lowpass filtered RegEM TTLS solution 
h1=plot(gca,t,index_cps.f20,'g-','linewidth',[2]); xlim(xlims);
h2=plot(gca,t,index_reg.f20,'r-','linewidth',[2]); xlim(xlims);
% plot regem-based confidence intervals
hs=area_fill(t',index_reg.f20p',index_reg.f20m',red,red,0,0.3); hold off;
lab{1}=['NINO3.4 CPS (20-year lowpass)'];
lab{2}=['NINO3.4 RegEM TTLS (raw)'  ];	
lab{3}=['NINO3.4 RegEM TTLS (20-year lowpass)'];
lab{4}=['95% C.I.'] ;
hold off
hl=legend([h1 H2 h2 hs],lab{:},'location','northwest'); legend boxoff  % may be H2 instead
set(hl,'FontName','Times','FontSize',9), hYLabel2=ylabel('NINO3.4(K)',style_t{:});
hTitle=title('NINO3.4 reconstructions',style_t{:});
set(gca,style_a{:}); hold off
text(1830,-2.2,['R^2 =',sprintf('%2.2g',R^2)]);


% print out figure as high-resolution postscript (vector) format 
%   Matlab's Painters renderer is the best for that, but unfortunately does
%   not respect transparency, so figures must be adjusted manually in
%   Illustrator (sample attached). 
% For best results, the transparency setting should be set to 30% in
% Illustrator (or some other postscrip-editing software of your choice)
hepta_figprint('nino34_recon_and_forcings')


% ===============================================================
% 2) The amazings things one can do with a proper database syntax
% ===============================================================

% load database
load NINO34_proxy_database_annual_2012ds 
% S  = raw proxy structure
% Sd = proxy structure with corals detrended by a 150-yr spline 

% extract relevant vectors
tmin = [S.Year_i]; tmax = [S.Year_f];
dt =   [S.Resolution];
lat =  [S.Lat]; 

% select a subset of the proxy database (this may be changed at will)
select = find(tmin <= 1850 & tmax >= 1975 & dt <= 5 & abs(lat <= 35)); 

% final structure
Sf = S(select); nr = length(Sf);
sitesf = {Sf.Site}';
type = {Sf.Type}';
for r = 1:nr
   if strcmpi(type(r),'\delta^{18}O') | strcmpi(type(r),'km^3');
      plab{r} = ['$', char(type(r)), '$'];
   else
      plab{r} = char(type(r));
   end
end

% put all records in one matrix 
tm = [1995:-1:1001];
nr = len(Sf); nm = length(tm); % number of years
proxy = nan(nm,nr);

for k=1:nr
	t=Sf(k).Chronology;
	d=Sf(k).Data;	
   ti = intersect(t,tm);
	proxy(ismember(tm,ti),k)=d(ismember(t,ti));
	clear d t 	
end
v = [1:nr]; ind = find(tm >= 1700);
fig('Proxy matrix'),clf
pcolor(tm(ind), v, standardize(proxy(ind,:))')
set(gca,'Ytick',v,'YtickLabel',plab,'FontSize',13); plotTickLatex2D('xlabeldy',0.00);
hTitle=title('Proxy Matrix View');
hYLabel=ylabel([]); hXLabel=xlabel('Time');
FontSize = 16;
set(gca,'FontName','Times','FontSize', round(FontSize*0.71));
set([hTitle, hXLabel, hYLabel], 'FontName' , 'Palatino');
set([hXLabel, hYLabel],'FontSize', round(FontSize*0.86));
set(hTitle, 'FontSize', FontSize, 'FontWeight' , 'bold');

% plot site name
lim = get(gca,'Xlim'); offset=0.02
text(repmat(lim(2)+offset*(lim(2)-lim(1)),length(v),1),v, ...
   sitesf,'VerticalAlignment','middle',...
   'HorizontalAlignment','left','rotation',45);
colorbar2('horiz')

hepta_figprint('proxy_matrix');

%==========================================================
%	Plot spatio-temporal locations 
% (requires the m_map package : http://www.eos.ubc.ca/~rich/map.html)
%==========================================================
Obj   = {Sf.Object};
class =  unique(Obj);
nclass = numel(class);
plon = [Sf.Lon];
plat = [Sf.Lat];
pcode = zeros(nr,1);
% Define proxy symbols for each type
icon{1,1}=ornj;         icon{1,2}= 'o'; icon{1,3} = 'Coral'; icon{1,4} = bck;
i_coral = strmatch('CORAL',Obj); pcode(i_coral)   = 1;
icon{2,1}= skyblue;      icon{2,2}= '*'; icon{2,3} = 'Ice core'; icon{2,4} = icon{2,1};
i_cryo = strmatch('ICE',Obj); pcode(i_cryo)   = 2;
icon{3,1}= Lawn_green; icon{3,2}= '^'; icon{3,3} = 'Tree'; icon{3,4} = bck;
i_cryo = strmatch('TREE',Obj); pcode(i_cryo)   = 3;
icon{4,1}= maroon;       icon{4,2}= 's'; icon{4,3} = 'Sediment'; icon{4,4} = bck;
i_sed  = strmatch('SEDI',Obj); pcode(i_sed)   = 4;
icon{5,1}= hotpink;      icon{5,2}= 'd'; icon{5,3} = 'Speleothem'; icon{5,4} = icon{5,1};
i_stal  = strmatch('SPELEO',Obj); pcode(i_stal)   = 5;

% Temporal pattern
proxycolor = zeros(nclass,3);
for m = 1:nclass
   nproxy(:,m) = sum(~isnan(proxy(:,pcode == m)),2);
   proxycolor(m,:) = icon{m,1};
end

%  FIGURE
fig('Database in color'),clf
subplot(3,1,1:2) 
% prepare canvas
m_proj('Robinson','clong',170); 
m_coast('patch',ligr);
m_grid('box','off','xtick',6,'ytick',9,'xlabeldir','middle', 'fontsize',5,'fontname','Times');
for n = 1:nr
	hl(n)=m_line(plon(n),plat(n),'marker',icon{pcode(n),2},'color',icon{pcode(n),4},'MarkerFaceColor',icon{pcode(n),1},'linewidth',[1],'MarkerSize',[8],'linestyle','none');
end
text(-2,1.85,['a) Proxy Network for NINO3.4 reconstruction (',num2str(nr) , ' records)'],style{:});
hn = m_line([190 190 240 240 190], [-5 5 5 -5 -5],'color',dkgr)
%
subplot(3,1,3)
hb = bar(tm,sq(nproxy),'stacked')
axis([1000 2000 0 ceil(nr/10)*10])
fancyplot_deco('b) Proxy availability','Time','# proxies');
colormap(proxycolor)
[LEGH,OBJH,OUTH,OUTM]=legend(hb,icon{:,3},'West'); % 4 is for Lower right-hand corner
set(LEGH,'Fontsize',[8],'FontName','Palatino','Position',[0.1574    0.1859    0.1674    0.1268]); 
legend boxoff

hepta_figprint('proxy_database');


% Plot individual timeseries

for n=1:6:nr
	fig(['ANNUAL records from TP database ' num2str(n:n+5)]), clf
	orient landscape
	k = n;
   while k <= nr && k < n+6
		%S(k).Site	
		subplot(3,2,k-n+1)
		t= Sf(k).Chronology;
		X= Sf(k).Data;
		plot(t,X,'r-'), grid; axis tight;
		tmin = max(1000,Sf(k).Year_i);
		xlim([tmin Sf(k).Year_f]);
		ylabel(char(Sf(k).Type),'FontName' ,'Times', 'FontSize', [11])
		xlabel('Time',style_l{:})
		% title
		rec=[num2str(k) ')'];
		obj=char(Sf(k).Object)
		type=char(Sf(k).Type);
		lon_str=[num2str(round(Sf(k).Lon)),'{}^{\circ}E'];
		lat=Sf(k).Lat;                                               
		if (lat>=0) 
			lat_str=[num2str(round(lat)),'{}^{\circ}N']
		else
			lat_str=[num2str(-round(lat)),'{}^{\circ}S']
		end
		str=[rec,': ', obj,' ',type,' in ', deblank(Sf(k).Site),' (',lon_str,',',lat_str,')'];
		title(str,'FontName','Times','FontSize',[14])
      k = k + 1;
   end
%orient(gcf,'landscape','PaperPositionMode','auto');   
filen=['proxyfigs_', num2str(n),'to',num2str(n+5)]; 
hepta_figprint(filen)
end

% This concludes our introduction to the power of a carefully designed proxy
% database. The possibilities are limitless, but they all derive from the 
% basic operations and metadata illustrated above. 

