% eda08_09
% CAC Pacific Sea Surface Temperature (SST) dataset
% read in and plot the data

% 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);

% plor data, 12 months by six years of maps per figure
for yb = [1:YBLOCKS]    
    figure(yb);
    clf;
    bw=0.9*(256-linspace(0,255,256)')/256;
    colormap([bw,bw,bw]);
    for m = [1:MONTHS]
    for y = [1:YBLOCKSIZE]
        j = m+MONTHS*(y-1)+MONTHS*YBLOCKSIZE*(yb-1);
        k = y + (m-1)*YBLOCKSIZE;
        if( j<=IMAGES)
            subplot(MONTHS, YBLOCKSIZE, k);
            hold on;
            axis equal;
            axis ij;
            set(gca,'XTick',[]); % turn off horizontal axis
            set(gca,'YTick',[]); % turn off vertical axis
            MAP=flipud(SST(:,:,j)'); % 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;
            if( y==1 )
                ylabel(sprintf('%d',themonth(j)));
            end
            if( m==12 )
                xlabel(sprintf('%d',theyear(j)));
            end
        end
    end
    end
end

