% RTSTRUCT Recontoolbox structure utility function
%     data = rtstruct() returns an example dataset that can be used with
%     most reconstruction methods. 
%  
%     rtstruct(data) checks if the input data conforms to the structure
%     format. Throws an error if it does not.
%     
%     Structure of data: (optional fields are marked with *)
%     data.proxy                      
%         Cell array of proxy data
%     data.proxy{i}.data              
%         Matrix of proxy record. Columns equal times and rows equal proxy 
%         records. Records with identical temporal structure may be stacked.
%         Can contain NaN's.
%     data.proxy{i}.times             
%         A row vector of mid-times for each sample. Must have equal number 
%         of columns as data. Measured in years AD.
%     *data.proxy{i}.lower            
%         Lower, or older, bound in time for each proxy sample. Measured in 
%         years AD.
%     *data.proxy{i}.upper            
%         Upper, or younger, bound in time for each proxy sample. Measured in 
%         years AD.
%     *data.proxy{i}.intargetres      
%         Boolean stating if the record is in the same resolution as the 
%         target. Some methods handle such records differently.
%     *data.proxy{i}.locations        
%         Locations of proxy records. Must have equal number of rows as data.
%         First column is degrees in latitude, second in longitude.
%     data.instrumental              
%         Struct of instrumental data.
%     data.instrumental.data 
%         Instrumental values. Columns correspond to times and rows to
%         locations. Can contain NaN's.
%     data.instrumental.times        
%         Mid-times for samples in instrumental data. Measured in years AD.
%         Currently only annual instrumental data is considered.
%     *data.instrumental.locations    
%         Locations of instrumental data. Must have equal number of rows as
%         data. First column is degrees in latitude, second in longitude.
%     data.target
%         Struct of reconstruction target.
%     data.target.times
%         Mid-times for target of reconstruction. Measured in years AD.
%         Currently only annual targets are considered.


% From the paper, section 5
%The instrumental noise standard deviation for uncertainty estimation is 
%obtained by first calculating the variance of the mean estimator for each 
%year, and then calculating the square root of the average of these variances.

% Deduce lower and upper bounds or just leave a single date?
ulbnds='y';

% CREATE PROXY AND INSTRUMENTAL TIME SERIES MATCHING THIS FORMAT

latlon=xlsread('shakun_2012_minimalproxydata.xlsx',2);
plat=latlon(:,1);plon=latlon(:,2);
pdata=xlsread('shakun_2012_minimalproxydata.xlsx',1); % alternates [data dates]


% Convert all proxy data first, then take out what I want below...
dataS=struct();
j=1;k=2;
for i=1:length(latlon)
   dataS.proxy{i}.locations=[plat(i) plon(i)];
   pda=pdata(:,j)';
   pda=pda(~isnan(pda)); % Safe to do because NaNs are only on the ends...
   dataS.proxy{i}.data=pda;

   % Figure out reasonable lower and upper bounds for the proxy: assume that the
   % upper and lower bounds are just the midpoints between mean dates
   dts=1950-pdata(:,k)';
   dts=dts(~isnan(dts)); % Remove NaNs that are only on the ends
   dataS.proxy{i}.times=dts;

   % check that data are times are same length
   if length(dts)~=length(pda); error('NaN removal did not work...'); end

   % Bin high resolution data into centennial bins
   if median(abs(diff(dts)))<50
      edg0=[2000:-100:dts(end)];
      [~,dtsN,bin] = histcounts(flip(dts),flip(edg0)); % determine bins; must flip for histcount
      % Remove data on end that doesn't fit in the bins
      ridx=find(bin==0);bin(ridx)=[];fpda=flip(pda);fpda(ridx)=[];
      pdaN=accumarray(bin',fpda,[],@mean)'; % mean of each bin
      dataS.proxy{i}.data=flip(pdaN); % flip time scale back
      dts=edg0(1:end-1)-50;% mid points of edges
      dataS.proxy{i}.times=dts;    
   end

 
   if ulbnds=='y'

      % Pad time data for next step
      df=dts(1)+abs(dts(1)-dts(2));
      de=dts(end)-abs(dts(end)-dts(end-1));
      dp=[df dts de];

      lwrs=nan(1,length(dts));
      uprs=nan(1,length(dts));
      for ii=1:length(dts);
	 uprs(ii)=median(dp(ii):-1:dp(ii+1));
	 %lwrs(ii)=median(dp(ii+1):-1:dp(ii+2))+1;
	 lwrs(ii)=median(dp(ii+1):-1:dp(ii+2)); % don't worry about offsetting
      end
      
      dataS.proxy{i}.lower=lwrs;
      dataS.proxy{i}.upper=uprs;


   end

   j=j+2;k=k+2;
end


%nlns=find(plon<0); plon(nlns)=plon(nlns)+360;

% Just take proxies in the quandrant over the North Atlantic, Greenland, Arctic,
%  Europe, and Northern Africa
ti=find(-60 < plon & plon < 40 & plat > 20);
data.proxy=dataS.proxy(ti);


% %====================================================================
% % Get the instrumental time series over this region
% [X_obs,xlat,xlon,x_yrs] = load_obs_S('bearth',4,3);
% X_obs=reshape(X_obs,length(xlat),length(xlon),length(x_yrs));
% 
% lonext=[40 (360-60)];latext=[20 85];
% 
% % Find nearest lat/lon equal to or outside bounds
% ln_a=find((abs(xlon-lonext(1))-min(abs(xlon-lonext(1))))<1e-5);
% ln_b=find((abs(xlon-lonext(2))-min(abs(xlon-lonext(2))))<1e-5);
% % if two equal mins, pick one outside of bound
% ln1=min(ln_a);ln2=max(ln_b);
% 
% % Find nearest lat/lon equal to or outside bounds
% lt_a=find((abs(xlat-latext(1))-min(abs(xlat-latext(1))))<1e-5);
% lt_b=find((abs(xlat-latext(2))-min(abs(xlat-latext(2))))<1e-5);
% % if two equal mins, pick one outside of bound
% lt1=min(lt_a);lt2=max(lt_b);
% 
% % Just take over the region
% ngae=X_obs(lt1:lt2,[ln2:360 1:ln1],:);
% % area weighting
% A=cosd(repmat(xlat(lt1:lt2),[1 length([ln2:360 1:ln1])]));
% % compute index
% ngaeInd=wmean_a(ngae,A);
% 
% % estimate instrumental noise as in Hanjavari paper
% n2=reshape(ngae,length(lt1:lt2)*length([ln2:360 1:ln1]),length(x_yrs));
% nvar2=nanvar(n2); % variance of the field each year
% data.instrumental.noisestd=sqrt(mean(nvar2)); 
% 
% % Just take instrumental data up to 2000, i = 151
% data.instrumental.data=ngaeInd(1:151)';
% data.instrumental.times=x_yrs(1:151);
% %====================================================================


% Use Shakun NH recon for the instrumental scaling
snh=xlsread('shakun_nh_temp.xlsx',1); % [dates data]

data.instrumental.data=snh(:,2)';
data.instrumental.times=snh(:,1)';



save('shakun_data_04.mat','data')

























