% PS05_TREND a m-file script to turn the images into tuplets, fit a trend
%	surface using SURFIT.M, and then produce the detrended data tuplets
%	that will ultimately be Kriged.
%
%	This script also produces a plot of the original data, the trend 
%	surface, and the detrended data. Much of this code was originally
%	in the main program PS05_MAIN.M, but moved here for modularity.
%
% Started 2000:10:06 D. Glover, WHOI
% Modif'd 2002:10:01 DMG to save trend surfaces in a .MAT file (ps04_trend.mat)
% Modif'd 2002:10:20 DMG to read in the images so this could be run separately
% Modif'd 2008:10:27 DMG changed the name from PS04 to PS05

% Start a stopwatch

t_init_trend=clock;

% Set a new colormap

load purp_map2;

purp_map2(1,:)=[1 1 1];             % Puts white at the bottom of the colorbar
cmap=colormap(purp_map2);           % so missing data comes out white (blank)

% Read in the binary satellite images

fid1=fopen('A1987130h09da_sst.bin');
fid2=fopen('A1987131h09da_sst.bin');
fid3=fopen('A1987132h09da_sst.bin');

% Assign them to X1, X2, and X3 respectively

X1=fread(fid1,[114 114],'uint8');
X2=fread(fid2,[114 114],'uint8');
X3=fread(fid3,[114 114],'uint8');

% Get some size parameters of the images

[nrows,ncols]=size(X1);

% Convert binary values into SSTs

X1sst=byte2sst(X1);
X2sst=byte2sst(X2);
X3sst=byte2sst(X3);

% Set missing data (-3) to Nan

X1sst(X1sst==-3)=nan;
X2sst(X2sst==-3)=nan;
X3sst(X3sst==-3)=nan;

% Set the values for Bermuda to missing values, but save the indices for
% later re-insertion

[iberm jberm]=find(X1sst>35);

X1sst(X1sst>35)=nan;
X2sst(X2sst>35)=nan;
X3sst(X3sst>35)=nan;

% Calculate the lat and lon of the image pixels

for ilat=1:ncols
   for ilon=1:nrows
      Alon((ilat-1)*114+ilon)=-69+(ilon-1)*(10/(nrows-1));
      Alat((ilat-1)*114+ilon)=37-(ilat-1)*(10/(ncols-1));
   end
end

% "Unroll" images into tuplets

Y1sst=X1sst(:);
Y2sst=X2sst(:);
Y3sst=X3sst(:);

Z1sst=[Alon' Alat' Y1sst];	% Note that the lon and lat had to
Z2sst=[Alon' Alat' Y2sst];	% be transposed to produce the 
Z3sst=[Alon' Alat' Y3sst];	% n x 3 data tuplet matrix

% Remove all Nan's from the tuplets (cannot trend surface fit to Nan's)
% This is a limitation of the SURFIT.M program.

Z1sste=excise(Z1sst);
Z2sste=excise(Z2sst);
Z3sste=excise(Z3sst);

% Do trend surface analysis, first fit the trend surfaces

[a1,sa1,chisqr1,covmat1]=surfit(Z1sste(:,1),Z1sste(:,2),Z1sste(:,3),1);
[a2,sa2,chisqr2,covmat2]=surfit(Z2sste(:,1),Z2sste(:,2),Z2sste(:,3),1);
[a3,sa3,chisqr3,covmat3]=surfit(Z3sste(:,1),Z3sste(:,2),Z3sste(:,3),1);

% Next calculate the trend surfaces

Z1tr=a1(1)+a1(2)*Alon+a1(3)*Alat;
Z2tr=a2(1)+a2(2)*Alon+a2(3)*Alat;
Z3tr=a3(1)+a3(2)*Alon+a3(3)*Alat;

% Now remove the trend surface from the ORIGINAL image, not the 
% "excised" version

Z1dt=[Alon' Alat' (Z1sst(:,3)-Z1tr')];
Z2dt=[Alon' Alat' (Z2sst(:,3)-Z2tr')];
Z3dt=[Alon' Alat' (Z3sst(:,3)-Z3tr')];

% Make trend surfaces into tuplets

Z1tr=[Alon' Alat' Z1tr'];
Z2tr=[Alon' Alat' Z2tr'];
Z3tr=[Alon' Alat' Z3tr'];

% Make the trend surfaces and detrended data into images for display

Z1trimg=reshape(Z1tr(:,3),nrows,ncols);
Z2trimg=reshape(Z2tr(:,3),nrows,ncols);
Z3trimg=reshape(Z3tr(:,3),nrows,ncols);

Z1dtimg=reshape(Z1dt(:,3),nrows,ncols);
Z2dtimg=reshape(Z2dt(:,3),nrows,ncols);
Z3dtimg=reshape(Z3dt(:,3),nrows,ncols);

% Display images, AS images

figure;

% NOTE: the images have to be FLIPUD'ed and transposed due to the way
% the original data were written out in binary. This is GENERALLY true
% for images that were processed via FORTRAN, although other languages
% have similar pecularities. The best thing to do is FIRST: LOOK AT 
% YOUR DATA to determine whether it is rightside up etc.

subplot(331);
imagesc([-69 -59],[27 37],flipud(X1sst'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;
title('1987 day 130','fontsize',14);ylabel('Original','fontsize',14);
subplot(332);
imagesc([-69 -59],[27 37],flipud(X2sst'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;
title('1987 day 131','fontsize',14);
subplot(333);
imagesc([-69 -59],[27 37],flipud(X3sst'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;
title('1987 day 132','fontsize',14);

subplot(334);
imagesc([-69 -59],[27 37],flipud(Z1trimg'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;
ylabel('Trend','fontsize',14);
subplot(335);
imagesc([-69 -59],[27 37],flipud(Z2trimg'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;
subplot(336);
imagesc([-69 -59],[27 37],flipud(Z3trimg'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;

subplot(337);
imagesc([-69 -59],[27 37],flipud(Z1dtimg'),[-2.5 2.5]);colorbar;axis xy;
colormap(cmap);axis square;
ylabel('Anomalies','fontsize',14);
subplot(338);
imagesc([-69 -59],[27 37],flipud(Z2dtimg'),[-2.5 2.5]);colorbar;axis xy;
colormap(cmap);axis square;
text(-80,23,'Trend Analysis of AVHRR PathFinder SST','fontsize',16);
subplot(339);
imagesc([-69 -59],[27 37],flipud(Z3dtimg'),[-2.5 2.5]);colorbar;axis xy;
text(-60,25,'PS05b 2002');
colormap(cmap);axis square;

% Save detrended images

save ps05_trend Alon Alat X1sst X2sst X3sst Z1trimg Z2trimg Z3trimg Z1dtimg Z2dtimg Z3dtimg iberm jberm

% Report trend surface analysis execution time

sprintf('%s%g\n','Trend Surface Elasped Time: ',etime(clock,t_init_trend))