
% Define a new equal area grid using "Recursive Zonal Equal Area (EQ) Sphere Partitioning Toolbox."

% Load some data to interpolate
t2m=ncread('/d1/nsteiger/climate-data/era5/era5_1deg_hrly_2017_12.nc','t2m');
t2m=flip(permute(t2m,[2 1 3]),1);
lat=double(ncread('/d1/nsteiger/climate-data/era5/landseamask_orography_era5_1deg.nc','latitude'));
lat=flip(lat);
lon=double(ncread('/d1/nsteiger/climate-data/era5/landseamask_orography_era5_1deg.nc','longitude'));

% Give the function a 2-D sphere with the number of desired points
points_x = rad2deg(eq_point_set_polar(2,100)); % for 1 deg grid the number of points is 65,160
nlon=points_x(1,:)';
nlat=points_x(2,:)'-90;

% create grids for interpolation
%[xn,yn]=meshgrid(xlon,xlat);
[xc,yc]=meshgrid(lon,lat);
xc=xc(:);yc=yc(:);

% note masks of both data types
%nm=find(~isnan(da(:,:,end)));
%cm=find(ndm.sftlf==1);
%da(da<-99)=NaN; % put nans in fill spots but after determining mask

t2mv=t2m(:,:,1); 
t2mv=t2mv(:); 

F = scatteredInterpolant(xc,yc,t2mv,'linear','linear');
vq=F(nlat,nlon); % do the interpolation

error('THIS INTERPOLATION DOESNT SEEM TO BE WORKING AND WILL TAKE MORE TIME TO FIGURE OUT WHATS GOING ON')
figure;h=axesm('MapProjection','wagner4'); 
scatterm(nlat,nlon,10,vq)


figure;h=axesm('MapProjection','wagner4'); 
pcolorm(lat,lon,t2m(:,:,1))


%da_i=nan(length(lat),length(lon),size(da,3));
%for i=1:size(da,3)
%   % linear interpolation and extrapolation for irregular data
%   da1=da(:,:,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(lat),length(lon));
%   vd(cm)=vq; % full interpolated field
%   da_i(:,:,i)=vd;
%end



