
tic
clear

load distinct_colors_20.mat


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');
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');
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');
pcp=flip(permute(pcp,[2 1 3]),1);


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=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'));

% 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=nearest_latlon(slat,slon,eic_lat,eic_lon);

% Grab easter island precipitation
pcpv=reshape(pcp,length(slat)*length(slon),size(pcp,3));
eip=pcpv(nlm,:);


disp('Removing the seasonal cycle...')
dt=(datenum('01-Jan-1979'):datenum('31-Dec-2017'))';
xv=reshape(mslp,length(slat)*length(slon),length(dt));
xv2=reshape(pcp,length(slat)*length(slon),length(dt));
%xv=reshape(t2m,length(slat)*length(slon),length(dt));
xv_ds=zeros(size(xv));
xv2_ds=zeros(size(xv));
parfor i=1:length(slat)*length(slon)
   xv_ds(i,:)=deseason(dt,xv(i,:));
   xv2_ds(i,:)=deseason(dt,xv2(i,:));
end

% Just take ocean values
xvo=xv_ds(oids,:);
%xvo=xv(oids,:);disp('Data is not deseasoned...')


% %----------------------------------------------------
% % Try just doing SOM analysis over large rain events
% 
% 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)); 


%------------------
%  Calculate SOMs
%------------------

% Detrend (not really much of a trend at all in the data...)
%D=detrend(xvo')';
D=xvo; disp('No detrending...')
%D=xvo(:,rnc{4}); disp('No detrending & only cluster 4 data')


toc

somRows=10;somCols=10;
disp(['SOMs = ' num2str(somRows*somCols)])
init='l'; % (l)inear or (r)andom initialization
% D[lat*lon,time]
disp('Computing SOM...')
[ptrns0,bmus,rmsErr,pfrq,sMap,sD] = som_make_njs(D,somRows,somCols,init,wt);
disp('SOM completed');toc


% Expand out the ptrns
ptrns=NaN(somRows*somCols,length(slat)*length(slon));
ptrns(:,oids)=ptrns0;

disp('Saving output...')
timestp2=strrep(datestr(now),' ','_');
save(['output_som' num2str(somRows*somCols) '_era5_mslp_' timestp2 '.mat'],'-v7.3')

toc

