% eda08_11
% CAC Pacific Sea Surface Temperature (SST) dataset
% plot EOF's and timeseries of corresponding amplitude timeseries

% load the data
D=load('cac_sst.txt');

% set up sizes of arrays
[I, J] = size(D);
NBLOCK=85;
NLON = NBLOCK-1;
NLAT = J-1;
IMAGES = floor(I/NBLOCK);
Dt = 1/12;

% arrays for SST, date, lats, lons
SST = zeros(NLON,NLAT,IMAGES);
thedate = zeros(IMAGES,1); % in month.year format
theyear = zeros(IMAGES,1); % year extracted from thedate
themonth = zeros(IMAGES,1); % month extracted from thedate
lats=zeros(NLAT,1);
lons=zeros(NLON,1);

% cut up data into SST arrays
lats = D(1,2:J)';
lons = D(2:NBLOCK,1);
for j=[1:IMAGES]
    k1 = 2+(j-1)*NBLOCK;
    k2 = k1+NLON-1;
    thedate(j)=D(k1-1,1);
    themonth(j)=floor(thedate(j));
    theyear(j)=floor(10000*(thedate(j)-themonth(j))+0.1);
    SST(:,:,j) = D(k1:k2,2:J);
end

% some definitions related to plotting
MONTHS=12;
YEARS=ceil(IMAGES/12);
YBLOCKSIZE = 6;
YBLOCKS = ceil(YEARS/YBLOCKSIZE);

% fold out into array
SSTV = zeros( IMAGES, NLAT*NLON );
for j=[1:IMAGES]
    for p = [1:NLAT]
    for q = [1:NLON]
        k = p + (q-1)*NLAT;
        SSTV(j,k) = SST( q, p, j );
    end
    end
end

% singular value decomposition
[U,S,V] = svd(SSTV,0);
s=diag(S);
Ns=length(s);
C=U*S;
F=V';

% fold EOF's back into maps
Nf = 12;
SSTF = zeros(NLON,NLAT,Nf);
for j=[1:Nf]
    for p = [1:NLAT]
    for q = [1:NLON]
        k = p + (q-1)*NLAT;
        SSTF( q, p, j ) = F(j,k);
    end
    end
end

% plot first Nf EOF's
figure(YBLOCKS+3);
clf;
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
for f = [1:Nf]
    subplot(3, 4, f);
    axis equal;
    hold on;
    axis ij;
    set(gca,'XTick',[]); % turn off horizontal axis
    set(gca,'YTick',[]); % turn off vertical axis
    MAP=flipud(SSTF(:,:,f)'); % reorient so it plots as map
    range=max(max(MAP))-min(min(MAP));
    if( range==0 )
        range=1;
    end
    imagesc( (MAP-min(min(MAP)))/range );
    axis tight;
    title(sprintf('EOF %d',f));
end

% plot corresponding amplitude timeseries
figure(YBLOCKS+4);
clf;
for f = [1:Nf]
    subplot(3, 4, f);
    plot( theyear(1)+[1:IMAGES]'/12, C(:,f), 'k-' );
    hold on;
    axis( [theyear(1), theyear(IMAGES)+1, min(C(:,1)), max(C(:,1)) ] );
    title(sprintf('C(%d,t)',f));
end
