load SeaWiFS_chl_NewEngland

% Animation of data fields: Noted the missing values which would cause trouble for EOF analysis.
% We need to exclude them.
figure
for t=1:size(chl,1)
  pcolor(lon,lat,squeeze(chl(t,:,:))), shading interp
  plotnewengland
  caxis([0 5])
  colorbar
  title(date_str{t})
  pause(0.2)
end

% look at all cloudy pixels
test_clouds = sum(chl);
figure
pcolor(lon,lat,squeeze(test_clouds)), shading interp
plotnewengland
title('show cloudy pixels')

% EOF analysis:
% 1) transform 3D chl array into 2D data matrix and remove NaNs
% 1.1) determine land and cloud points (NaNs)
bad_ind = find(isnan(test_clouds));
good_ind = find(~isnan(test_clouds));
% 1.2) create 2D data matrix without NaNs
for t=1:size(chl,1)
  x = chl(t,:);
  x(bad_ind) = [];
  chl2D(t,:) = x;
end
% 1.3) dimensions of our data matrix are:
[M N] = size(chl2D);

% 2) calculate anomalies
chl2D_mn = sum(chl2D)/M; % means
chl2D_an = chl2D - ones(M,1)*chl2D_mn; % anomalies

% 3) calculate eigenvalues and eigenvectors of cov-mat
% 3.1) calculate covavriance matrix
C = chl2D_an'*chl2D_an/(M-1);
% 3.2) calculate SVD (and take the time)
tic
[U,S,V] = svd(C,0);
toc
disp('(for SVD only)')

% 4) look at the first 4 EOFs, variance explained and scores
% 4.1) turn back into 2D matrices and plot EOFs
aux = NaN*ones(size(chl,2)*size(chl,3),1); % initialize
for i=1:4
  aux(good_ind) = V(:,i);
  EOF{i} = reshape(aux,size(chl,2),size(chl,3));
end
% plot EOFs
figure
for i=1:4
  subplot(2,2,i)
  pcolor(lon,lat,EOF{i}), shading interp
  plotnewengland, colorbar
  title(['EOF ' num2str(i)])
end
% 4.2) look at variance explained
var_explained = diag(S)/trace(S);
figure
plot(var_explained(1:10)*100,'o-')
ylabel('variance explained [%]')
xlabel('EOF #')
% 4.3) look at EOF amplitude through time (PC scores)
for i=1:4
  score{i} = chl2D_an*V(:,i);
end
figure
plot(score{1},'-')
hold on
plot(score{2},'r-')
legend('1','2','Location','SouthWest')
xlabel(['month since ' date_str{1}])
ylabel('Principal Component Score')
