% 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

disp('Loading proxy data...')

efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';
load([pth,efl],'proxy','xmeta')

clat=xmeta{1}.lat;
clon=xmeta{1}.lon;

% Find proxies over southern South America
p_lat=zeros(length(proxy),1);
p_lon=zeros(length(proxy),1);
for i=1:length(proxy)
   p_lat(i)=proxy{i}.location(1);
   p_lon(i)=proxy{i}.location(2);
end
% Find which proxies are over south america land areas
%pi=find(p_lat<0 & p_lat>-57 & p_lon>270 & p_lon<320);
% Find which proxies are in the north-central Chile region
pi=find(p_lat<-26 & p_lat>-45 & p_lon>270 & p_lon<293);


% Which proxies to use?
%prxytyp='p'; % PDSI-like trees
prxytyp='a'; % all proxies

j=1;
data=struct();
for i=1:length(pi)
   if strcmp(prxytyp,'p')
      if strcmp(proxy{pi(i)}.ptype,'Tree Rings')
	 if strcmp(proxy{pi(i)}.psm_mthd,'p')
	 %if strcmp(proxy{pi(i)}.psm_mthd,'t')

	    % Proxy locations
	    data.proxy{j}.locations=[proxy{pi(i)}.location(1) proxy{pi(i)}.location(2)];
            % Proxy data, make sure the sign of the proxy data lines up with positive correlation with PDSI
	    %data.proxy{j}.data=(proxy{pi(i)}.data)*sign(proxy{pi(i)}.pcr);
            data.proxy{j}.data=proxy{pi(i)}.data; % don't manually adjust sign
	    % Proxy years
	    data.proxy{j}.times=proxy{pi(i)}.age_yng';

	    j=j+1;
	 end
      end
   elseif strcmp(prxytyp,'a')

      % Proxy locations
      data.proxy{j}.locations=[proxy{pi(i)}.location(1) proxy{pi(i)}.location(2)];
      % Proxy data
      data.proxy{j}.data=proxy{pi(i)}.data;
      % Proxy years
      data.proxy{j}.times=proxy{pi(i)}.age_yng';
      j=j+1;
   end
end



%====================================================================
% Get the instrumental time series over this region

disp('Interpolating obs to CESM grid')

% LOAD OBSERVATIONAL DATA
mon_avg_o=12;mon_avg_f=2;
obs_dtp1='vanderS_pdsi';
%mon_avg_o=4;mon_avg_f=3;
%obs_dtp1='bearth'; % particular data type specified by input
[O_t,xlat,xlon,x_yrs] = load_obs_S(obs_dtp1,mon_avg_o,mon_avg_f);
O_t=reshape(O_t,length(xlat),length(xlon),length(x_yrs));



% land mask
msk=load('chile_mask.mat');
lndmsk=msk.chl;

% create grids for interpolation
[xn,yn]=meshgrid(xlon,xlat);
[xc,yc]=meshgrid(clon,clat);
% note masks of both data types
nm=find(~isnan(O_t(:,:,end)));
cm=find(~isnan(lndmsk));

o_i=nan(length(clat),length(clon),length(x_yrs));
for i=1:length(x_yrs)
   % linear interpolation and extrapolation for irregular data
   da1=O_t(:,:,i);
   F = scatteredInterpolant(xn(nm),yn(nm),da1(nm),'linear','linear');
   vq=F(xc(cm),yc(cm)); % do the interpolation
   % expand out the interpolated data points onto the climate model mask
   vd=nan(length(clat),length(clon));
   vd(cm)=vq; % full interpolated field
   o_i(:,:,i)=vd;
end


% Latitude weighting
A=cosd(repmat(clat,[1 length(clon)]));
pdsi_vds_chl=wmean_a(o_i,A); 

% estimate instrumental noise as in Hanjavari paper
oiv=reshape(o_i,length(clat)*length(clon),length(x_yrs));
nvar2=nanvar(oiv(:,3:end)); % variance of the field each year
data.instrumental.noisestd=sqrt(mean(nvar2)); 

%instrumental data without first 3 years for PDSI spin-up
data.instrumental.data=pdsi_vds_chl(3:end)';
data.instrumental.times=x_yrs(3:end);
%====================================================================


disp('Saving data...')
%save('phyda_data_chile.mat','data')
%save('phyda_data_chile_nosignadj.mat','data')
save('phyda_data_chile_allproxies.mat','data')

























