% PS05_COKRI the m-file that runs the cokriging program COKRI.M Keep in
%   mind that PS05_TREND.M must be run at least once before this m-file
%   can be run in a standalone fashion.
%
% Started 2002:10:22 D. Glover, WHOI
% Modif'd 2004:10:18 DMG to add in the negative co-semivariogram nuggets
% Modif'd 2008:10:27 DMG changed the name from PS04 to PS05

% Start a stopwatch

t_init_cokri=clock;

% Set a colormap

load purp_map;
cmap=purp_map;
cmap(1,:)=[1 1 1];          % white for missing data
cmap(64,:)=[0.5 0.5 0.5];   % grey for land or really large values

% Load the trend surfaces and detrended data

disp('Reading in PS05_TREND.MAT');
disp(' ');
load ps05_trend;

% Unroll the detrended images

A1sst=Z1dtimg(:);
A2sst=Z2dtimg(:);
A3sst=Z3dtimg(:);

% Make the "Grand Data Matrix"

Xanom=[Alon' Alat' A1sst A2sst A3sst];

% Announce the beginning of Co-Kriging

disp('Beginning Co-Kriging');

% Declare some master cokriging variables

model=[1 1      1;          % Spherical model
       4 1.4431 1.2826];
   
% model=[4 0.6045 0.5206];    % Exponential model supports no nugget, that
                            % is to say neg. nugget => zero nugget
   
c=[0.0071 -0.0255 -0.0181;    % Spherical model
   -0.0255 0.0142 -0.0086;
   -0.0181 -0.0086 0.0370;
   0.4539 0.2743 0.1378;
   0.2743 0.2821 0.1283;
   0.1378 0.1283 0.2803];

%c=[0.4987 0.2975 0.1466;
%   0.2975 0.2895 0.1311;
%   0.1466 0.1311 0.2954];   % Exponential model

% I did the problem both ways (exponential and spherical models) and I chose
% the spherical model because the error fields were better over all. While it
% is true that the error fields were marginally better for days 131 and 132,
% the error field for day 130 was much worse for the exponential vs. the
% spherical (see ps05_cokri_exp.eps vs ps05_cokri.eps). Plus I didn't care for
% the fact that the exponential model did not support a nugget since these 
% images come from the Bermuda area where we know there is quite a bit of 
% mesoscale variability. That variability should translate into unresolved 
% "sub-grid scale" variability in the semivariograms. But the exponential model
% yielded consistently negative nuggets, an unrealistic result, so I did not 
% include a nugget model in the exponential version of PS05_COKRI.M. Other 
% interpretations are possible and your choice of theoretical semivariogram 
% model will not be held against you.

% These are some more parameters for COKRI.M

avg=[mean(Xanom(~isnan(Xanom(:,3)),3)) mean(Xanom(~isnan(Xanom(:,4)),4)) mean(Xanom(~isnan(Xanom(:,5)),5))];
block=[1 1];
nd=[1 1];
ival=0;
nk=40;	% number of nearest points to use
rad=1;	% search radius
ntok=1;

% COKRIGE the data

format compact

[x0s,s,sv,id,l]=cokri(Xanom,[Alon' Alat'],model,c,1,avg,block,nd,ival,nk,rad,ntok);

format loose

% Check for complex array in x0s, this is generally not a good sign, if 
% you give COKRI.M real numbers to krige, then your results should be real

if isreal(x0s)
   x0s_real=x0s;
else
   x0s_real=real(x0s);
   disp('There were complex numbers in x0s');
end

% Check for complex array in s

if isreal(s)
   s_real=s;
else
   s_real=real(s);
   disp('There were complex numbers in s');
end

% Convert the results back into images (i.e. add the trend back in)

X1OA=reshape(x0s_real(:,3),114,114)+Z1trimg;
X2OA=reshape(x0s_real(:,4),114,114)+Z2trimg;
X3OA=reshape(x0s_real(:,5),114,114)+Z3trimg;

% Convert the error fields into images

S1OA=reshape(s_real(:,3),114,114);
S2OA=reshape(s_real(:,4),114,114);
S3OA=reshape(s_real(:,5),114,114);

% Put Bermuda back into the plots

for k=1:length(iberm)
   X1sst(iberm(k),jberm(k))=35.2;
   X2sst(iberm(k),jberm(k))=35.2;
   X3sst(iberm(k),jberm(k))=35.2;

   X1OA(iberm(k),jberm(k))=35.2;
   X2OA(iberm(k),jberm(k))=35.2;
   X3OA(iberm(k),jberm(k))=35.2;

   S1OA(iberm(k),jberm(k))=35.2;
   S2OA(iberm(k),jberm(k))=35.2;
   S3OA(iberm(k),jberm(k))=35.2;

end%for

% Make some plots of the results (original, obj anal, error)

figure

subplot(331); % First three subplots use the original images
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(X1OA'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;
ylabel('Kriged','fontsize',14);
subplot(335);imagesc([-69 -59],[27 37],flipud(X2OA'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;
subplot(336);imagesc([-69 -59],[27 37],flipud(X3OA'),[18 26]);colorbar;axis xy;
colormap(cmap);axis square;

subplot(337);imagesc([-69 -59],[27 37],flipud(S1OA'),[-0.0075 0.42]);colorbar;axis xy;
colormap(cmap);axis square;
ylabel('Error','fontsize',14);
subplot(338);imagesc([-69 -59],[27 37],flipud(S2OA'),[-0.0075 0.42]);colorbar;axis xy;
colormap(cmap);axis square;
text(-80,23,'Objective Analysis of AVHRR PathFinder SST','fontsize',16);
subplot(339);imagesc([-69 -59],[27 37],flipud(S3OA'),[-0.0075 0.42]);colorbar;axis xy;
text(-60,25,'PS05d 2002');
colormap(cmap);axis square;

% Report cokriging execution time

sprintf('%s%g\n','Cokrig Elasped Time: ',etime(clock,t_init_cokri))