
load distinct_colors_20.mat

% Load observatoinal data
eir0 = xlsread('/home/nsteiger/Dropbox/data/easter_island_rain.xlsx');
eir_yr=eir0(1,:);
eir=eir0(2:end,:);

disp('Load ERA5 data')


%t2m=ncread('/d1/nsteiger/climate-data/era5/era5_t2m_swpac_daily_1979_2017.nc','t2m');
%t2m=ncread('/d1/nsteiger/climate-data/era5/era5_t2m_scpac_daily_1979_2017.nc','t2m'); % larger area
%t2m=ncread('/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc','t2m'); % larger area
%t2m=flip(permute(t2m,[2 1 3]),1);

%mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_swpac_daily_1979_2017.nc','msl');
%mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_scpac_daily_1979_2017.nc','msl'); % larger area
mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_sh_daily_1979_2017.nc','msl'); % larger area
mslp=flip(permute(mslp,[2 1 3]),1);

%pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_swpac_daily_1979_2017.nc','tp');
%pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_scpac_daily_1979_2017.nc','tp'); % larger area
pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_sh_daily_1979_2017.nc','tp'); % larger area
pcp=flip(permute(pcp,[2 1 3]),1);
pcp(pcp<0)=0; % specify 0 for missing values


%slat=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_swpac_daily_1979_2017.nc','latitude'));
%slat=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_scpac_daily_1979_2017.nc','latitude'));
slat=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc','latitude'));
slat=flip(slat);
%slon=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_swpac_daily_1979_2017.nc','longitude'));
%slon=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_scpac_daily_1979_2017.nc','longitude'));
slon=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc','longitude'));

% Where is ocean?
lat=ncread('/d1/nsteiger/climate-data/era5/landseamask_orography_era5_1deg.nc','latitude');
lat=flip(lat);
lon=ncread('/d1/nsteiger/climate-data/era5/landseamask_orography_era5_1deg.nc','longitude');
[~,islt]=intersect(lat,slat);
[~,isln]=intersect(lon,slon);

lsm=ncread('/d1/nsteiger/climate-data/era5/landseamask_orography_era5_1deg.nc','lsm');
lsm=flip(permute(lsm,[2 1]),1);
lmsk=lsm(islt,isln);
ln=find(lmsk>=0.5);oc=find(lmsk<0.5);
omsk=lmsk;omsk(ln)=0;omsk(oc)=1;

% Weight the data over the ocean by area
oids=find(omsk(:)==1);
W=cosd(repmat(slat,[1 length(slon)]));W=W(:);
wt=W(oids);

% Easter island lats/lons: -27.116667, -109.366667
eic_lat=-27.116667;
eic_lon=360-109.366667;

% Where in the model is Easter Island?
[nlm,nlat,nlon]=nearest_latlon(slat,slon,eic_lat,eic_lon);
%[nlm_h,nlat_h,nlon_h]=nearest_latlon(hlat,hlon,eic_lat,eic_lon);

% What area to analyze? 
% lat(25)=-62.24, lat(70)=-11.77, lon(185)=207, lon(252)=282 
%et=30:72;en=185:252;
%et=25:103;en=125:259; % broad tropical Pacific box
%eilat=lat;eilon=lon;

%pcpvh=reshape(pcp_hr,length(hlat)*length(hlon),size(pcp_hr,3));
pcpv=reshape(pcp,length(slat)*length(slon),size(pcp,3));
mslpv=reshape(mslp,length(slat)*length(slon),size(mslp,3));
%t2mv=reshape(t2m,length(slat)*length(slon),size(t2m,3));

ei_p0=pcpv(nlm,:);
% In the obs, 0.1 mm is the lowest measurable amount, so impose that here for the analysis
ei_p0(ei_p0<1e-4)=0; % 1e-4 m = 0.1 mm
ei_m0=mslpv(nlm,:);
%ei_t0=t2mv(nlm,:);


% Start and end dates of data
ot0=datetime(1955,1,1):datetime(2017,12,31);
ot=datenum(1955,1,1):datenum(2017,12,31);

% add NaNs for leap days since obs didn't include leap days...
leapdays=ot0(month(ot0)==2 & day(ot0)==29);
[~,ilp]=intersect(ot0,leapdays);
ida=setxor(1:length(ot0),ilp);

eirv=NaN(length(ot0),1);
eirv(ida)=eir(:);

% Remove leap days from time vector
%ot0(month(ot0)==2 & day(ot0)==29)=[];
%ot=datenum(ot0);

% ERA dates
dt=(datenum('01-Jan-1979'):datenum('31-Dec-2017'))';
dtt=(datetime('01-Jan-1979'):datetime('31-Dec-2017'))';
dsht=(datetime('01-Jan-1979'):datetime('31-Dec-1981'))';


% Use obs for the rain, but replace missing obs with ERA data
[~,idte,idto]=intersect(dtt,ot0);
% where are missing values
ei_oe=eirv(idto); % obs + era5
ei_e0=ei_p0(idte)*1e3; % scale data to mm
ei_e0=round(ei_e0,1); % obs are rounded to nearest 0.1 mm, but model data is not...
nanobs=find(isnan(ei_oe));
ei_oe(nanobs)=ei_e0(nanobs); % replace nans with scaled era

ei_p=ei_oe; % use obs infilled with era5 % deseasoning makes negative rain!


%-------------------
% KMEANS CLUSTERING
%-------------------

% Load pre-computed k-means for consistency (sometimes algorithm struggles to converge)
load('kmeans_clusters_8-May-2019.mat')

% Sort by amount of precip
[cent_s,is]=sort(cent);

%rnc{1}=nrn; % no rain
rnc{1}=find(cidx==is(1)); 
rnc{2}=find(cidx==is(2)); 
rnc{3}=find(cidx==is(3)); 
rnc{4}=find(cidx==is(4)); 

pc1=NaN(size(ei_p));
pc2=NaN(size(ei_p));
pc3=NaN(size(ei_p));
pc4=NaN(size(ei_p));

pc1(rnc{1})=ei_p(rnc{1});
pc2(rnc{2})=ei_p(rnc{2});
pc3(rnc{3})=ei_p(rnc{3});
pc4(rnc{4})=ei_p(rnc{4});

% Count up the frequencies of each cluster per year
[y0,m0,~,~,~,~] = datevec(dt);
CompArray = datenum(y0,1,1);
ut = unique(CompArray);
fpc=zeros(4,length(ut));
for k = 1:length(ut)
   fpc(1,k)=sum(~isnan(pc1(CompArray==ut(k))));
   fpc(2,k)=sum(~isnan(pc2(CompArray==ut(k))));
   fpc(3,k)=sum(~isnan(pc3(CompArray==ut(k))));
   fpc(4,k)=sum(~isnan(pc4(CompArray==ut(k))));
end

% annual cluster 4 frequency
frqc4=fpc(4,:);
frqc4=frqc4(:);



% mean annual precipitation
[mp,mdt]=mean_ts(ei_p,dt,'nansum','year');
[mp_mon,~]=mean_ts(ei_p,dt,'nansum','month');

% sort into wet, dry, and intermediate years
[idx,cvs]=kmeans(mp(:),3);
[~,i]=sort(cvs);
wids=find(idx==i(end)); % wet year indices
iids=find(idx==i(2)); % intermediate year indices
dids=find(idx==i(1)); % dry year indices


% Reshape data
k=1;l=12;
mx_o=zeros(12,length(mp_mon)/12);
for n=1:length(mp_mon)/12
    mx_o(:,n)=mp_mon(k:l);
    k=k+12; l=l+12;
end

% Seasons
dseas=[1,2,10,11,12];
iseas=[3,8,9];
wseas=[4,5,6,7];

mxd=mx_o(dseas,:);
mxi=mx_o(iseas,:);
mxw=mx_o(wseas,:);

% mean seasonal data
mn_ds=mean(mxd,1);mn_ds=mn_ds(:);
mn_is=mean(mxi,1);mn_is=mn_is(:);
mn_ws=mean(mxw,1);mn_ws=mn_ws(:);

% mean cluster data
[mpc1,mdt]=mean_ts(pc1,dt,'nansum','year');
[mpc2,mdt]=mean_ts(pc2,dt,'nansum','year');
[mpc3,mdt]=mean_ts(pc3,dt,'nansum','year');
[mpc4,mdt]=mean_ts(pc4,dt,'nansum','year');
[mpc4_mon,~]=mean_ts(pc4,dt,'nansum','month');
mpc1=mpc1(:);
mpc2=mpc2(:);
mpc3=mpc3(:);
mpc4=mpc4(:);
mp=mp(:);


%-------------------------------------
% Compare against local high pressure
%-------------------------------------

% monthly mean mslp
[mslp_ds_m,mon_dt]=mean_ts(mslp,dt,'month');
%[mn_mslp,mon_dt]=mean_ts(mslp,dt,'month');

% seasonal means -----------------------------

% Reshape data
k=1;l=12;
mX_o=zeros(size(mslp_ds_m,1),size(mslp_ds_m,2),12,length(mon_dt)/12);
for n=1:length(mon_dt)/12
    mX_o(:,:,:,n)=mslp_ds_m(:,:,k:l);
    k=k+12; l=l+12;
end

mslp_dry=squeeze(mean(mX_o(:,:,dseas,:),3));
mslp_trn=squeeze(mean(mX_o(:,:,iseas,:),3));
mslp_wet=squeeze(mean(mX_o(:,:,wseas,:),3));


% Annual mean
[mn_mslp,tm_dt]=mean_ts(mslp,dt,'year');


% project pressure values onto zonal/meridional projections ----------------------

% take region of interest
%lonmin=180;lonmax=290;
lonmin=221;lonmax=288; % 221->220 deg, 288->287 deg = coastline of chile
latmin=41;latmax=76; % 41->50 degS, 76->15 degS 

mslpX=mn_mslp./100;
s1=squeeze(mean(mslpX([latmin:latmax],[lonmin:lonmax],:),1));

figure
hold on
for i=1:39
   if ismember(i,wids)
      pclr=clrs.green;
   elseif ismember(i,iids)
      pclr=clrs.grey;
   elseif ismember(i,dids)
      pclr=clrs.brown;
   end
   plot(lonmin:lonmax,s1(:,i),'-','color',pclr)
   %ylim([3.64e4 3.675e4])
end
xlabel('Longitude')
ylabel('Mean MSLP')
box on
%print(['/home/nsteiger/Dropbox/easterisle/mslp_dist_ann_bywetdryint.png'],'-dpng','-r300');

% % compute center of mass within a box off south america ----------------------
% [Xs,Ys] = meshgrid(slon,slat);
% 
% mslpX=mn_mslp./100; % convert to mb
% % just take region off South America
% mslpX([1:latmin latmax:end],:,:)=0;
% mslpX(:,[1:lonmin lonmax:end],:)=0;
% 
% cm_ann=zeros(length(tm_dt),2);
% for i=1:length(tm_dt)
%    cm_ann(i,:)=centroid(Xs(:),Ys(:),mslpX(:,:,i));
% end
% 
% mslpX=mslp_dry./100; % convert to mb
% % just take region off South America
% mslpX([1:latmin latmax:end],:,:)=0;
% mslpX(:,[1:lonmin lonmax:end],:)=0;
% 
% cm_dry=zeros(length(tm_dt),2);
% for i=1:length(tm_dt)
%    cm_dry(i,:)=centroid(Xs(:),Ys(:),mslpX(:,:,i));
% end
% 
% mslpX=mslp_wet./100; % convert to mb
% % just take region off South America
% mslpX([1:latmin latmax:end],:,:)=0;
% mslpX(:,[1:lonmin lonmax:end],:)=0;
% 
% cm_wet=zeros(length(tm_dt),2);
% for i=1:length(tm_dt)
%    cm_wet(i,:)=centroid(Xs(:),Ys(:),mslpX(:,:,i));
% end



% use edge of 1020 contour level ---------------------------------------
pcut=1020; % pressure cut off
%pcut=1018; % pressure cut off
% reset lon values for this computation
lonmin=181;lonmax=288; % 221->220 deg, 288->287 deg = coastline of chile

mslpX=mn_mslp./100; %dd convert to mb
% find values greater than threshold to make a mask of regions
mslpX(mslpX<pcut)=0;
% just take region off South America
mslpX([1:latmin latmax:end],:,:)=0;
mslpX(:,[1:lonmin lonmax:end],:)=0;

% area-weight
W=cosd(repmat(slat,[1 length(slon)]));W=W(:);
mslpXV=reshape(mslpX,length(slat)*length(slon),size(mslpX,3));
mslpXw=reshape(bsxfun(@times,mslpXV,W),length(slat),length(slon),size(mslpX,3));

% center of mass location of these 1020 pressure blobs -------------------------
[Xs,Ys] = meshgrid(slon,slat);
cm=zeros(size(mslpX,3),2);
for i=1:size(mslpX,3)
   cm(i,:)=centroid(Xs(:),Ys(:),mslpXw(:,:,i));
end

% compute distance from easter island
cmD=dist_gc(eic_lat,eic_lon,cm(:,2),cm(:,1));

disp('ann: corr(total pcp, 1020 centroid lon)')
corr(mp,cm(:,1)) % 0.4343
disp('ann: corr(cluster 4, 1020 centroid lon)')
corr(mpc4,cm(:,1)) % 0.3720
disp('ann: corr(freq cluster 4, 1020 centroid lon)')
corr(frqc4,cm(:,1)) % 0.4565
disp('ann: corr(freq cluster 4, 1020 centroid distance)')
corr(frqc4,cmD) % 0.4654


% colormap to match contour maps from NCL
hsvM=zeros(length(mdt),3);
hsvM(:,2)=0.9; % saturation
hsvM(:,3)=0.9; % value
ncolors=39;hmin=1;hmax=360;
hsvM(:,1)=(hmin:(hmax-hmin)/(ncolors-1):(hmax+1))/360; % normalized hue array
rgbM = hsv2rgb(hsvM);


figure
hold on
for i=1:length(mdt)
plot(cm(i,1),frqc4(i),'o','color',rgbM(i,:),'markerfacecolor',rgbM(i,:))
%plot(cmD(i),frqc4(i),'o','color',rgbM(i,:),'markerfacecolor',rgbM(i,:))
end
ylabel('Cluster 4 annual frequency','fontsize',16)
xlabel('Pressure centroid longitude (deg W)','fontsize',16)
%xlabel('Distance to pressure centroid (km)','fontsize',16)
% 360-(240:5:275)
xticklabels({'120','115','110','105','100','95','90','85'}) % match xticks with NCL fig
set(gca,'fontsize',14)
box on
rp=corr(frqc4,cm(:,1)); % 0.4565
%rp=corr(frqc4,cmD); % 0.465
text(0.05,0.05,['r = ' sprintf('%0.2f',rp)],'units','normalized','fontsize',14)
print(['/home/nsteiger/Dropbox/easterisle/cluster4freq_centLon_' num2str(pcut) '.png'],'-dpng','-r300');
%print(['/home/nsteiger/Dropbox/easterisle/cluster4freq_centDist_' num2str(pcut) '.png'],'-dpng','-r300');



% save as netcdf fields?
svncdf='n';
if svncdf=='y'

% define input and output variable
inputf='era5_sh_latlon.nc';
outputf=['era5_spac_1020contours_ann.nc'];
%outputf=['era5_spac_1018contours_ann.nc'];

% grab lat/lon from model, put it into the output file
eval(['!ncks -A -v latitude,longitude ',inputf,' ',outputf]);

% Add time dimension
nccreate(outputf,'time','Dimensions',{'time',length(mdt)});
ncwrite(outputf,'time',1:length(mdt));
ncwriteatt(outputf,'time','units','Annual mean');
ncwriteatt(outputf,'time','long_name','Year');

for i=1:length(mdt)

   Xv=mslpX(:,:,i);
   cnm=['mslp1020_' num2str(1978+i)];
   %cnm=['mslp1018_' num2str(1978+i)];

   %nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)},'FillValue',1e20);
   nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
   ncwrite(outputf,cnm,permute(Xv,[2 1]));
%   ncwriteatt(outputf,cnm,'missing_value',1e20);

end

% add center of mass lat/lon values
nccreate(outputf,'cmlat','Dimensions',{'time',length(mdt)});
ncwrite(outputf,'cmlat',cm(:,2));
nccreate(outputf,'cmlon','Dimensions',{'time',length(mdt)});
ncwrite(outputf,'cmlon',cm(:,1));


end

return


% take meridional sum, pick off extreme edge
mrd_mn=squeeze(sum(mslpAnn,1));
for i=1:size(mslpAnn,3)
   mrd_znid(i)=find(mrd_mn(:,i),1);
end
% lon edge of 1020 contour
ln1020=slon(mrd_znid);

disp('ann: corr(total pcp, 1020 lon edge)')
corr(mp,ln1020) % 0.3972
disp('ann: corr(cluster 4, 1020 lon edge)')
corr(mpc4,ln1020) % 0.3511
disp('ann: corr(freq cluster 4, 1020 lon edge)')
corr(frqc4,ln1020) % 0.4347

% location of meridional max ---------------------------------

for i=1:size(mrd_mn,2)
   [mrd_mx_vl(i),mrd_mx_lon(i)]=max(mrd_mn(:,i));
end

disp('ann: corr(total pcp, meridional max)')
corr(mp,mrd_mx_lon(:)) % 0.4079
disp('ann: corr(cluster 4, meridional max)')
corr(mpc4,mrd_mx_lon(:)) % 0.3312
disp('ann: corr(freq cluster 4, meridional max)')
corr(frqc4,mrd_mx_lon(:)) % 0.3997

disp('ann: corr(total pcp, meridional max)')
corr(mp,mrd_mx_vl(:)) % -0.2617
disp('ann: corr(cluster 4, meridional max)')
corr(mpc4,mrd_mx_vl(:)) % -0.2806
disp('ann: corr(freq cluster 4, meridional max)')
corr(frqc4,mrd_mx_vl(:)) % -0.3229

% multivariate model fails: mdl_cmp_spcz=fitlm([cm(:,1) mrd_mx_vl(:)],frqc4(:))



% local pressure correlation ------------------------------------------
ei_prss=mean_ts(mslpv(nlm,:),dt,'year');

disp('ann: corr(total pcp, 1020 centroid lon)')
corr(mp,ei_prss(:)) % -0.5174
disp('ann: corr(cluster 4, 1020 centroid lon)')
corr(mpc4,ei_prss(:)) % -0.4812 
disp('ann: corr(freq cluster 4, 1020 centroid lon)')
corr(frqc4,ei_prss(:)) % -0.5180 


% lat of SPCZ and slope of SPCZ -------------------------------------------
%addpath('/home/nsteiger/MATLAB/advanced-stats/')

% spatial pcp field
[pcp_mn,mdt]=mean_ts(pcp,dt,'nansum','year');

lt1=30;i1=find((abs(slat-lt1)-min(abs(slat-lt1)))<1e-5);
lt1=-50;i2=find((abs(slat-lt1)-min(abs(slat-lt1)))<1e-5);
[spcz_mx,spcz_zm,smeta] = spcz_lon(pcp_mn,slat(i2:i1),slon,[160 240]);% lats for SPCZ

%p0=mon2ann_sum(pcp_mn(i2:i1,:,:),4,7);% try seasonally averaging first
%[spcz_wet,spcz_wzm,smeta_w] = spcz_lon(p0,slat(i2:i1),slon,[160 240]);% lats for SPCZ


disp('ann: corr(total pcp, 1020 centroid lon)')
corr(mp,spcz_zm(:)) % 0.0324
disp('ann: corr(cluster 4, 1020 centroid lon)')
corr(mpc4,spcz_zm(:)) % 0.1020 
disp('ann: corr(freq cluster 4, 1020 centroid lon)')
corr(frqc4,spcz_zm(:)) % 0.1289 


% linear regression
bd=zeros(length(mdt),2);
for i=1:length(mdt)
   mdl=fitlm(smeta.longitudes,spcz_mx(:,i));
   bd(i,:)=mdl.Coefficients.Estimate;
end

% How well does the linear model of SPCZ fit precip changes?
mdl_mp=fitlm(bd,mp);
mdl_mpc4=fitlm(bd,mpc4);
mdl_frqc4=fitlm(bd,frqc4);


% linear regression models ---------------------------------

disp('regression model: pressure + spcz')
disp('mean cluster 4 precip')
mdl_cmp_spcz=fitlm([cm(:,1) spcz_zm(:)],mpc4(:)) % p-value for spcz = 0.89
disp('frequency cluster 4 precip')
mdl_cmp_spcz=fitlm([cm(:,1) spcz_zm(:)],frqc4(:)) % p-value for spcz = 0.83

return


% pressure thresholding by season doesn't work for wet and transition season
% b/c values are often  too low to pass 1020 or 1019, while a lower threshold would
% make the mean pressure far too wide


% mslpX=mslp_dry./100; %dd convert to mb
% % find values greater than threshold to make a mask of regions
% mslpX(mslpX<pcut)=0;
% % just take region off South America
% mslpX([1:latmin latmax:end],:,:)=0;
% mslpX(:,[1:lonmin lonmax:end],:)=0;
% mslpDry=mslpX;
% 
% mslpX=mslp_trn./100; % convert to mb
% % find values greater than threshold to make a mask of regions
% mslpX(mslpX<pcut)=0;
% % just take region off South America
% mslpX([1:latmin latmax:end],:,:)=0;
% mslpX(:,[1:lonmin lonmax:end],:)=0;
% mslpTrn=mslpX;
% 
% mslpX=mslp_wet./100; % convert to mb
% % find values greater than threshold to make a mask of regions
% mslpX(mslpX<pcut)=0;
% % just take region off South America
% mslpX([1:latmin latmax:end],:,:)=0;
% mslpX(:,[1:lonmin lonmax:end],:)=0;
% mslpWet=mslpX;


figure
for i=1:39
pcolor(mslpAnn(:,:,i))
print(['/home/nsteiger/Dropbox/easterisle/mslp_ann_' num2str(pcut) '_yr' num2str(i) '.png'],'-dpng','-r300');
end

return

for i=1:39
pcolor(mslpDry(:,:,i))
print(['/home/nsteiger/Dropbox/easterisle/mslp_dry_' num2str(pcut) '_yr' num2str(i) '.png'],'-dpng','-r300');
end

for i=1:39
pcolor(mslpWet(:,:,i))
print(['/home/nsteiger/Dropbox/easterisle/mslp_wet_' num2str(pcut) '_yr' num2str(i) '.png'],'-dpng','-r300');
end

for i=1:39
pcolor(mslpTrn(:,:,i))
print(['/home/nsteiger/Dropbox/easterisle/mslp_trn_' num2str(pcut) '_yr' num2str(i) '.png'],'-dpng','-r300');
end




% find lat and lon ranges, area, and sum of pressure in the box, and also centroid lat/lon values

% zonal mean


% mean of all months, convert to mb
%mn_mslp=mean(mslp_ds_m,3)./100;



return

% find values greater than threshold to make a mask of regions
pcut=1020; % pressure cutoff
mn_mslp(mn_mslp<pcut)=NaN;
mn_mslp(mn_mslp>=pcut)=1;

% just take region off South America
mn_mslp([1:50 74:end],:)=NaN;
mn_mslp(:,[1:220 290:end],:)=NaN;

% now apply mask to mslp data
[mslp_ann,~]=mean_ts(mslp,dt,'year');
tys=size(mslp_ann,3);
mslp_sa=bsxfun(@times,mslp_ann,mn_mslp);
mslp_sa_mon=bsxfun(@times,mslp_ds_m,mn_mslp);
% take area weighted mean
mslp_mnsa=wmean(mslp_sa,slat);
mslp_mnsa_mon=wmean(mslp_sa_mon,slat);

% center of mass ----------------------
% interpolate to high resolution
% hhlat=min(slat):0.1:max(slat);
% hhlon=min(slon):0.1:max(slon);
% t=size(mslp_ann,3);
% [X,Y] = meshgrid(hhlon,hhlat);
% [X1,Y1] = meshgrid(slon,slat);
% cdata=zeros(length(hhlat),length(hhlon),t);
% for j=1:t
%    cdata(:,:,j) = interp2(X1,Y1,mslp_ann(:,:,j),X,Y);
% end

% compute center of mass
[Xs,Ys] = meshgrid(slon,slat);
cm=zeros(tys,2);
for i=1:tys
   cm(i,:)=centroid(Xs(:),Ys(:),mslp_sa(:,:,i));
end
cm_mon=zeros(length(mon_dt),2);
for i=1:length(mon_dt)
   cm_mon(i,:)=centroid(Xs(:),Ys(:),mslp_sa_mon(:,:,i));
end


disp('ann: corr(total pcp, centroid lon)')
corr(mp,cm(:,1))
disp('ann: corr(total pcp,highpres)')
corr(mp,mslp_mnsa)
disp('ann: corr(clus 4,highpres)')
corr(mpc4,mslp_mnsa)
disp('ann: corr(cluster 4, centroid lon)')
corr(mpc4,cm(:,1))

disp('mon: corr(total pcp, centroid lon)')
corr(mp_mon(:),cm_mon(:,1))
disp('mon: corr(total pcp,highpres)')
corr(mp_mon(:),mslp_mnsa_mon)
disp('mon: corr(cluster 4 pcp, centroid lon)')
corr(pc4(:),cm_mon(:,1))

% use this example
% mdl=fitlm(X4(rsp,xm{i}),y4(rsp));
% rsqrd(i,ii,k)=mdl.Rsquared.Adjusted;
% beta_bs(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.Estimate;
% p_bs(i,[1 (xm{i}+1)],ii,k)=mdl.Coefficients.pValue;
% aic_bs(i,ii,k)=mdl.ModelCriterion.AICc;   
% bic_bs(i,ii,k)=mdl.ModelCriterion.BIC;   

% linear regression models ---------------------------------
disp('mean pressure reg model')
mdl_mp=fitlm(mslp_mnsa,mp(:));% mean pressure model
%r2_mp=mdl_mp.Rsquared.Adjusted;

disp('centroid lon reg model')
mdl_clon=fitlm(cm(:,1),mp(:))
%r2_clon=mdl_mp.Rsquared.Adjusted;

disp('mean pressure + centroid lon reg model')
mdl_mpclon=fitlm([cm(:,1) mslp_mnsa],mp(:))
%r2_mpclon=mdl_mp.Rsquared.Adjusted;



figure
plot(mp(:),cm(:,1),'o','color',clrs.purple)
xlabel('Annual precipitation (mm)','fontsize',16)
ylabel('Pressure centroid longitude (deg E)','fontsize',16)
set(gca,'fontsize',14)
rp=corr(cm(:,1),mp(:));
disp(['total pcp vs centroid lon, corr =' num2str(rp) ])
text(0.05,0.05,['r = ' sprintf('%0.2f',rp)],'units','normalized','fontsize',14)
print(['/home/nsteiger/Dropbox/easterisle/totalPcp_centLon_' num2str(pcut) '.png'],'-dpng','-r300');


figure
plot(mp_mon(:),cm_mon(:,1),'.','color',clrs.purple)
xlabel('Monthly precipitation (mm)','fontsize',16)
ylabel('Pressure centroid longitude (deg E)','fontsize',16)
set(gca,'fontsize',14)
rp=corr(cm_mon(:,1),mp_mon(:));
disp(['mon: total pcp vs centroid lon, corr =' num2str(rp) ])
text(0.05,0.05,['r = ' sprintf('%0.2f',rp)],'units','normalized','fontsize',14)
print(['/home/nsteiger/Dropbox/easterisle/totalPcp_centLon_' num2str(pcut) '_mon.png'],'-dpng','-r300');

% figure
% plot(mpc4(:),cm(:,1),'o')
% xlabel('cluster4 precip','fontsize',14)
% ylabel('centroid lon','fontsize',14)
% title('cluster4 pcp vs centroid lon')
% rp=corr(cm(:,1),mpc4(:));
% disp(['cluster 4 vs centroid lon, corr =' num2str(rp) ])
% text(0.05,0.05,['r = ' sprintf('%0.2f',rp)],'units','normalized','fontsize',14)
% print(['/home/nsteiger/Dropbox/easterisle/clus4Pcp_centLon_' num2str(pcut) '.png'],'-dpng','-r300');
% 
% figure
% plot(mp(:,1),mslp_mnsa(:),'o')
% xlabel('total precip','fontsize',14)
% ylabel('mean pressure (mb)','fontsize',14)
% title('total pcp vs mean pressure')
% rp=corr(mslp_mnsa(:),mp(:));
% text(0.05,0.05,['r = ' sprintf('%0.2f',rp)],'units','normalized','fontsize',14)
% print(['/home/nsteiger/Dropbox/easterisle/totalPcp_meanPress_' num2str(pcut) '.png'],'-dpng','-r300');


return

%-------------------------------
% Compare against ITCZ and SPCZ
%-------------------------------
addpath('/home/nsteiger/MATLAB/advanced-stats/')

[pcp_mn,mdt]=mean_ts(pcp,dt,'nansum','month');

lt1=30;i1=find((abs(slat-lt1)-min(abs(slat-lt1)))<1e-5);
lt1=-50;i2=find((abs(slat-lt1)-min(abs(slat-lt1)))<1e-5);
p0=mon2ann_sum(pcp_mn(i2:i1,:,:),4,7);% try seasonally averaging first
[spcz_wet,spcz_wzm,smeta_w] = spcz_lon(p0,slat(i2:i1),slon,[160 240]);% lats for SPCZ
p0=mon2ann_sum(pcp_mn(i2:i1,:,:),10,2);% try seasonally averaging first
[spcz_dry,spcz_dzm,smeta_d] = spcz_lon(p0,slat(i2:i1),slon,[160 240]);% lats for SPCZ

p0=pcp_mn(i2:i1,:,:); % MONTHLY
[spcz_mon,spcz_monzm,smeta_d] = spcz_lon(p0,slat(i2:i1),slon,[160 240]);% lats for SPCZ

[itcz_mx,~] = itcz_lon(pcp_mn(i2:i1,:,:),slat(i2:i1),slon,[130 170]);% lats for SPCZ

mp_w=mon2ann_sum(mp_mon,4,7);
mp_d=mon2ann_sum(mp_mon,10,2);

%corr(spcz_wet,repmat(mp_w,[1 81])');


% linear regression
bd=zeros(length(mp_d),2);
for i=1:length(mp_d)
   mdl=fitlm(smeta_d.longitudes,spcz_dry(:,i));
   bd(i,:)=mdl.Coefficients.Estimate;
end

% How well does the linear SPCZ fit precip changes?
mdl_d=fitlm(bd,mp_d);

bw=zeros(length(mp_w),2);
for i=1:length(mp_w)
   mdl=fitlm(smeta_w.longitudes,spcz_wet(:,i));
   bw(i,:)=mdl.Coefficients.Estimate;
end

% How well does the linear SPCZ fit precip changes?
mdl_w=fitlm(bw,mp_w);


%p = polyfit(smeta_d.longitudes,spcz_dry(:,1),1)
%f = polyval(p,smeta_d.longitudes);
%x=smeta_d.longitudes;y=spcz_dry(:,1);figure;plot(x,y,'o',x,f,'-')

% 
% % Reshape data
% k=1;l=12;
% mx_o=zeros(12,length(mp_mon)/12);
% sp=zeros(12,length(mp_mon)/12);
% it=zeros(12,length(mp_mon)/12);
% for n=1:length(mp_mon)/12
%     mx_o(:,n)=mp_mon(k:l);
%     sp(:,n)=spcz_monzm(k:l);
%     it(:,n)=itcz_mx(k:l);
%     k=k+12; l=l+12;
% end
% 
% % Seasons
% dseas=[1,2,10,11,12];
% iseas=[3,8,9];
% wseas=[4,5,6,7];
% 
% mxd=mx_o(dseas,:);
% mxi=mx_o(iseas,:);
% mxw=mx_o(wseas,:);
% 
% spd=sp(dseas,:);
% spi=sp(iseas,:);
% spw=sp(wseas,:);
% 
% itd=it(dseas,:);
% iti=it(iseas,:);
% itw=it(wseas,:);
% 
% figure
% hold on
% plot(spd(:),mxd(:),'.')
% plot(spi(:),mxi(:),'.')
% plot(spw(:),mxw(:),'.')
% title('SPCZ')
% legend('d','t','w')
% 
% figure
% hold on
% plot(itd(:),mxd(:),'.')
% plot(iti(:),mxi(:),'.')
% plot(itw(:),mxw(:),'.')
% title('ITCZ')
% legend('d','t','w')
% 

ltt1=150;ltt2=210;
p0=pcp_mn(i2:i1,:,:); % MONTHLY
[spcz_mon,spcz_monzm,smeta_d] = spcz_lon(p0,slat(i2:i1),slon,[ltt1 ltt2]);% lats for SPCZ


% Reshape data
k=1;l=12;
sp=zeros(12,length(mp_mon)/12);
for n=1:length(mp_mon)/12
    sp(:,n)=spcz_monzm(k:l);
    k=k+12; l=l+12;
end

spd=sp(dseas,:);
spi=sp(iseas,:);
spw=sp(wseas,:);

disp('corr(spcz dry,pcp dry)')
corr(spd(:),mxd(:))
disp('corr(spcz trans,pcp trans)')
corr(spi(:),mxi(:))
disp('corr(spcz wet,pcp wet)')
corr(spw(:),mxw(:))



% figure
% hold on
% plot(spd(:),mxd(:),'.')
% plot(spi(:),mxi(:),'.')
% plot(spw(:),mxw(:),'.')
% box on
% title(['SPCZ [' num2str(ltt1) ',' num2str(ltt2) ']' ])
% legend('dry','trans','wet','location','northwest')
% 



% % march through different lat bands
% 
% for i=150:10:240
% 
% ltt1=i-30;ltt2=i;
% p0=pcp_mn(i2:i1,:,:); % MONTHLY
% [spcz_mon,spcz_monzm,smeta_d] = spcz_lon(p0,slat(i2:i1),slon,[ltt1 ltt2]);% lats for SPCZ
% 
% % Reshape data
% k=1;l=12;
% mx_o=zeros(12,length(mp_mon)/12);
% sp=zeros(12,length(mp_mon)/12);
% for n=1:length(mp_mon)/12
%     mx_o(:,n)=mp_mon(k:l);
%     sp(:,n)=spcz_monzm(k:l);
%     k=k+12; l=l+12;
% end
% 
% % Seasons
% dseas=[1,2,10,11,12];
% iseas=[3,8,9];
% wseas=[4,5,6,7];
% 
% mxd=mx_o(dseas,:);
% mxi=mx_o(iseas,:);
% mxw=mx_o(wseas,:);
% 
% spd=sp(dseas,:);
% spi=sp(iseas,:);
% spw=sp(wseas,:);
% 
% 
% figure
% hold on
% plot(spd(:),mxd(:),'.')
% plot(spi(:),mxi(:),'.')
% plot(spw(:),mxw(:),'.')
% title(['SPCZ [' num2str(ltt1) ',' num2str(ltt2) ']'])
% legend('d','t','w')
% 
% end



% %[mpcp,mdt]=mean_ts(pcp,dt,'nansum','year');
% [mpcp,mdt]=mean_ts(pcp,dt,'nansum','month');
% lt1=30;i1=find((abs(slat-lt1)-min(abs(slat-lt1)))<1e-5);
% lt1=-30;i2=find((abs(slat-lt1)-min(abs(slat-lt1)))<1e-5);
% [spcz_mx,~] = spcz_lon(mpcp(i2:i1,:,:),slat(i2:i1),slon,[160 240]);% lats for SPCZ
% %[itcz_mx,~] = itcz_lon(mpcp(i2:i1,:,:),slat(i2:i1),slon,[130 170]);% lats for SPCZ
% disp('corr(itcz, cluster 4)')
% disp('1-12')
% corr(mon2ann(zscore(itcz_mx(:)),1,12),mon2ann(zscore(mp_mon(:)),1,12))
% disp('10-2')
% corr(mon2ann(zscore(itcz_mx(:)),10,2),mon2ann(zscore(mp_mon(:)),10,2)) %dry season
% disp('4-7')
% corr(mon2ann(zscore(itcz_mx(:)),4,7),mon2ann(zscore(mp_mon(:)),4,7)) %wet season




%---------------------
% Compare against MJO
%---------------------
mjo0dt=(datetime('01-Jan-1979'):datetime('28-Aug-2018'))';
mjo=load('/home/nsteiger/Dropbox/data/OMI_MJO_index_NOAA.txt');
[~,~,midt]=intersect(dtt,mjo0dt); % use just up to 2017

mjoi=mjo(midt,7);
[mjoi_mon,~]=mean_ts(mjoi,midt,'nanmean','month');
%figure;scatter(mon2ann_sum(mp_mon(:),10,2),mon2ann(mjoi_mon(:),10,2))
disp('corr(MJO mon,precip mon)')
corr(mon2ann_sum(mp_mon(:),10,2),mon2ann(mjoi_mon(:),10,2))



%----------------------
% Compare against ENSO
%----------------------

[X_obs,xlat,xlon,time0] = load_obs_Smon('bearth');
X_obs=reshape(X_obs,length(xlat),length(xlon),length(time0)*12);
[t_nino0]=nino(X_obs,xlat,xlon,3,'n');
[~,ier,ibe] = intersect(1979:2017,time0);
tninoann=mon2ann(t_nino0,1,12);
disp('corr(nino34, cluster 4)')
corr(tninoann(ibe),mpc4(ier))




%---------------------
% COMPARE AGAINST SAM
%---------------------


% Load SH MSLP
%mslp_sh=ncread('/d1/nsteiger/climate-data/era5/era5_msl_sh_daily_1979_2017.nc','msl'); % larger area
%mslp_sh=flip(permute(mslp_sh,[2 1 3]),1);

%shlat=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc','latitude'));
%shlat=flip(shlat);
%shlon=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc','longitude'));

% extract mslp at 45 and 60 degrees south
lt1=-65;idl_65=find((abs(slat-lt1)-min(abs(slat-lt1)))<1e-5);
lt1=-40;idl_40=find((abs(slat-lt1)-min(abs(slat-lt1)))<1e-5);

% potentially extract lon region of interest
ln1=130;iln1=find((abs(slon-ln1)-min(abs(slon-ln1)))<1e-5);
ln2=170;iln2=find((abs(slon-ln2)-min(abs(slon-ln2)))<1e-5);

% monthly mean mslp
[mslp_ds_m,mon_dt]=mean_ts(reshape(mslp,length(slat),length(slon),length(dt)),dt,'month');

% Full SAM region
zm_mslp_40=mean(squeeze(mslp_ds_m(idl_40,:,:)),1);
zm_mslp_65=mean(squeeze(mslp_ds_m(idl_65,:,:)),1);

% Partial sam region
%zm_mslp_40=mean(squeeze(mslp_ds_m(idl_40,iln1:iln2,:)),1);
%zm_mslp_65=mean(squeeze(mslp_ds_m(idl_65,iln1:iln2,:)),1);

sami=zscore(zm_mslp_40)-zscore(zm_mslp_65);

disp('corr(sam, cluster 4)')
disp('1-12')
corr(mon2ann(zscore(sami(:)),1,12),mon2ann(zscore(mp_mon(:)),1,12)) 
disp('10-2')
corr(mon2ann(zscore(sami(:)),10,2),mon2ann(zscore(mp_mon(:)),10,2))
disp('4-7')
corr(mon2ann(zscore(sami(:)),4,7),mon2ann(zscore(mp_mon(:)),4,7))

sami_ann=mon2ann(sami(:),1,12);

figure
plot(mp(:),sami_ann(:),'o')
xlabel('total precip','fontsize',14)
ylabel('SAM index','fontsize',14)
title('total pcp vs SAM index')
rp=corr(sami_ann(:),mp(:));
text(0.05,0.05,['r = ' sprintf('%0.2f',rp)],'units','normalized','fontsize',14)
print(['/home/nsteiger/Dropbox/easterisle/totalPcp_SAMindx.png'],'-dpng','-r300');



