%Get info out of netcdf
file = 'thetao_surface_Omon_GFDL-ESM2M_piControl_r1i1p1_049601-050012.nc';
info = ncinfo(file);
latS=-20;latN=20;
lonW=110;lonE=300;

%%
%Extract SST, lat, lon data from netcdf
SST = ncread(file,'thetao',[1 1 1 1], [inf inf 1 inf]); %only want top level of data
LAT = ncread(file,'lat');
LON = ncread(file,'lon');
LON(LON<0) = LON(LON<0)+360;

%%
SST = squeeze(SST);

%Convert to columns
SST = reshape(SST, [size(SST,1)*size(SST,2) size(SST,3)]);

%Remove Seasonal Cycle
for i=1:12
    SSTclim(:,i) = nanmean(SST(:,i:12:end),2);
end
SST = SST-repmat(SSTclim, [1 size(SST,2)/12]);

%Keep only SST in Equatorial Pacific: 130E-280E, 15S-15N
SST_eqP = SST;
SST_eqP(LAT(:)>latN,:) = NaN; SST_eqP(LAT(:)<latS,:) = NaN;
SST_eqP(LON(:)<lonW,:) = NaN; SST_eqP(LON(:)>lonE,:) = NaN;

%remove rows with NaN
SST_nonNaN = SST_eqP(~any(isnan(SST_eqP),2),:);

%%
%calculate covariance matrix
N = size(SST_nonNaN,2);
C = SST_nonNaN*SST_nonNaN'/N;

%calculate PCs (U) and eigenvalues
[U,D] = eig(C);

%calculate time expansion coefficents
V = U'*SST_nonNaN;

%Variance explained by each EOF
trace_C=trace(C);
fraction_variance=diag(D)/trace_C;

fraction_variance(end)

%%
%Extract EOF1, and put it into same shape as LAT and LON
EOF1 = nan(size(SST,1),1);
EOF1(~any(isnan(SST_eqP),2))=U(:,end);
EOF1 = reshape(EOF1,size(LAT));

%Plot the EOF
load coast;
figure;
axesm('miller','MapLonLimit',[lonW lonE], 'MapLatLimit',[latS latN]);
pcolorm(double(LAT),double(LON),EOF1);
plotm(lat,long,'k');