%% Simple Kriging
%
% Phaedon Kyriakidis
% Department of Geography


%% Load data
load bayAreaData.mat

%% Display data
% help mappointdata
cLims = [0 15];     % color limits for display
icolCoords = [2 3]; % column #s for coordinates
bulletSize = 25;    % bullet size
icolVar = 4;        % column # for variable to display
vType = 'cont';     % continuous variable
plotPars = {'color','vbar'}; % colorscale and vertical bar
mappointdata(bayprecip,icolCoords,bulletSize,icolVar,vType,cLims,plotPars); 
hold on; 
plot(baycoast(:,1),baycoast(:,2),'k-');
hold off;
title('Rain gauge precip data');


%% Compute & display sample semivariogram
% help spassociation
% help plotspassociation
LagDistLimits = []; % empty -> automatic definition of lag distance limits
G = spassociation(bayprecip,icolCoords,icolVar,LagDistLimits,'variogram');
plotspassociation(G);

%% Define model semivariogram
% help spassociationmodel
nDim = 2; % dimensionality -> 2D
% mixed or nested model with total sill = 9.2+0.1 = 9.3 mm^2
% nugget model with partial sill 0.1
% exponential model with partial sill 9.2
% isotropic exponential model with range 1.3 degrees
ModelSpecs = {'nugget',0.1,nan;'exponential',9.2,[0 1.3 1.3]};
M = spassociationmodel(nDim,ModelSpecs);
M.gslib % model specification in format used in other programs

%% Evaluate semivariogram model
% help spassociationmodeleval
d = (0:0.01:1.5); % distance values
az = 0;           % azimuth (use any value for isotropic model)
E = spassociationmodeleval(M,'distdir',d,'variogram',az);
% distdir -> mode of computation = distances along direction
disp(E);

%% Display model fit
plotspassociation(G);
hold on
plotspassociationmodel(E)
hold off

%% Simple kriging
% help krigep2p -> p2p = point source data & point target values
% assume known constant mean = sample mean of precip data 
mData = mean(bayprecip(:,4)); % sample precip mean
sData = std(bayprecip(:,4));  % sample precip std deviation
disp(GridSpecs); % display interpolation grid specification
iMode = 1;       % mode of computation -> 1 = regular grid
iNestFilt = [ ]; % empty array = use entire variogram model
krigeType = 0;   % type of Kriging = simple kriging
precipSK = krigep2p(iMode,GridSpecs,bayprecip,icolCoords,icolVar,M.gslib,...
                    iNestFilt,krigeType,mData);
size(precipSK) % 108000 x 2 values -> kriging estimates and std deviations
% at 300*360 nodes of regular raster specified in GridSpecs  

%% Display Simple Kriging results
% You can mask out the ocean and then display the results as:
precipSK(baymask==0,:)=NaN; % Set to NaNs predictions over 0s in baymask
maprasterdata(GridSpecs,precipSK,1,vType,cLims,plotPars);
title('Simple Kriging rainfall predictions');
figure
maprasterdata(GridSpecs,precipSK,2,vType,[0 3],plotPars); % plot 2nd col
title('Simple Kriging rainfall std deviations');

%% Comparison with sample data
% We can compare the summary statistics of the SK rainfall predictions
% to those of the rain gauge data as:
mPred = nanmean(precipSK(:,1));
sPred = nanstd(precipSK(:,1));
disp('Mean of predicted vs rain gauge precip');
disp( [mPred mData] );
disp('Std deviation of predicted vs rain gauge precip');
disp( [sPred sData] );
%%
% from which we can deduce that SK yields a predicted surface with similar
% mean but considerable smaller std deviation than their sample
% counterparts computed from the rain gauge precip data.
%
% Last, we can compare the entire distribution of the SK predictions 
% versus that of the original rain gauge data as:
qqplot(bayprecip(:,4),precipSK(:,1));
xlabel('rain gauge precip');
ylabel('SK-predicted precip');
title('QQ-plot of SK-predicted versus sample precip data');
%%
% from which we can see that interpolation tends to over-estimate the 
% lower part of the sample distribution and under-estimate the upper
% part of that distribution; this is known as the smoothing effect 
% of spatial interpolation.

%% SK-derived local conditional distributions
% Recall that Kriging yields 2 output values at each prediction location:
% (i) the SK attribute prediction, and (ii) the SK prediction error 
% variance. These 2 values define, a Gaussian in this case, distribution
% of possible attribute values at that location, given the sample data
% values; this distribution is termed a *local conditional distribution*. 
% The mean of the local conditional distribution is the SK prediction, 
% and the variance of that distribution is the SK prediction error 
% variance.
% Function |postkrige| takes as input the output from the Kriging funtion
% |krigep2p| and, assuming a Gaussian local distribution (|iDistr = 0|),
% outputs simulated attribute values at each prediction location via
% Monte Carlo simulation from the corresponding local Gaussian 
% conditional distribution. This can be done with function 
% |postkrige| by setting |iMode = 4| and providing an array |simPars =
% [nSim seed]| with the # of realizations to generate and the random #
% seed:
iDistr = 0; % Guassian distribution assumed
iMode = 4;  % mode of computation -> 4 = simulation
simPars = [50 4567];
precipSKsim = postkrige(precipSK,iDistr,iMode,simPars);
disp( size(precipSKsim) ); % 108000 x 50 values
%% 
% which is an array with as many rows as prediction locations and as
% many columns as # of realizations; each column is a different realization
% which can be displayed as:
maprasterdata(GridSpecs,precipSKsim,1,'cont',[0 15],plotPars);
title('Precip realization #1');
figure
maprasterdata(GridSpecs,precipSKsim,10,'cont',[0 15],plotPars); % 10th column
title('Precip realization #10');
%%
% The mean of the simulated realizations tends towards the SK predictions
% and the variance of the realizations tends towards the SK error variance
% as the # of realizations tends to infinity...
%
% The simulated attribute values generated above at any arbitrary 
% prediction location can be used for uncertainty propagation purposes
% in single-cell (or pointwise) models with no spatial interactions. Note,
% however, that such attribute values are generated *independently*
% from one prediction location to another, hence do not account properly
% for spatial auto-correlation and lead to "speckly" realizations.

%% Expected attribute surface
% In what follows, we will consider a known “expected” precipitation 
% surface stored in GeoEAS format in array |bayprecipexpect|:
maprasterdata(GridSpecs,bayprecipexpect,1,'cont',cLims,plotPars);
title('Expected precip surface');
%%
% The values of this expected surface at the 77 grid nodes closest 
% to the 77 rain gauge locations are stored in the 5-th column of 
% array |bayprecip|:
disp(bayprecip(1:10,:));

%% Rain gauge residuals
% We can construct a new data array with residual source data as: 
bayprecipres = [bayprecip(:,1:3) bayprecip(:,4)-bayprecip(:,5)]; 
%%
% and display them using:
mappointdata(bayprecipres,icolCoords,bulletSize,4,'cont',[-6 8],plotPars);
title('Rain gauge residuals from expected precip surface');
hold on
plot(baycoast(:,1),baycoast(:,2),'k-');
hold off
%% 
% We can now compute the mean & variance of these residual data, 
% and compare them to those of the original rain gauge data:
mRes = mean(bayprecipres(:,end));
sRes = std(bayprecipres(:,end));
disp('Original data and residual mean');
disp([mData mRes]);
disp('Original data and residual std deviation');
disp([sData sRes]);
%%
% from which it is evident that the residuals have a nearly zero mean and a
% std deviation smaller than that of the original rain gauge data; this is
% due to the expected precip surface "explaining" some of the precip
% spatial variability found in the rain gauge data.

%% Semivariogram of rain gauge residuals
% We can compute and display the sample semivariogram of the rain gauge 
% precipitation residuals as: 
R = spassociation(bayprecipres,icolCoords,4,LagDistLimits,'variogram');
plotspassociation(R);
title('Sample residual semivariogram');
%%
% the resulting residual semivariogram has a smaller sill (less variance)
% and a somewhat smaller range than the semivariogram of the original
% rain gauge data.

%% Residual model semivariogram
% Assuming a residual semivariogram model given as: 
ModelSpecsRes = {'nugget',0.1,nan;'exponential',6.6,[0 0.9 0.9]};
MdlRes = spassociationmodel(nDim,ModelSpecsRes);
MdlRes.gslib % model specification in format used in other programs

%% Evaluate residual semivariogram model
% help spassociationmodeleval
Eres = spassociationmodeleval(MdlRes,'distdir',d,'variogram',az);
% distdir -> mode of computation = distances along direction
disp(Eres);

%% Display residual model fit
plotspassociation(R);
hold on
plotspassociationmodel(Eres)
hold off

%% Simple Kriging of rain gauge residuals
% We can now predict the unknown residuals at the nodes of the 
% (300 × 360) grid using: 
mRes = 0; % Known expectation of residuals
precipSKres=krigep2p(1,GridSpecs,bayprecipres,[2 3],4,MdlRes.gslib,[ ],0,mRes);

%% Display results
precipSKres(baymask==0,:)=NaN; % Set to NaNs predictions over 0s in baymask
maprasterdata(GridSpecs,precipSKres,1,'cont',[-6 8],plotPars);
title('Simple Kriging residual rainfall predictions');
figure
maprasterdata(GridSpecs,precipSKres,2,'cont',[0 3],plotPars);
title('Simple Kriging residual rainfall std deviations');
%% 
% where the SK residual std deviations have a similar spatial pattern as
% the SK std deviations based on the original data, apart from local 
% differences reflecting the smaller sill of the residual semivariogram
% model.

%% Adding the expected attribute surface to the SK residuals
% We can now add this residual surface to the expected attribute 
% surface to compute the final SK precipitation predictions as: 
precipSKexp = bayprecipexpect + precipSKres(:,1); 
maprasterdata(GridSpecs,precipSKexp,1,'cont',cLims,plotPars);
title('SK residuals + expected precip surface');
%%
% The resulting surface inherits spatial patterns from both the SK residual
% surface and the expected precip surface. There is no way one can achieve
% a similar interpolated surface only from the sparse rain gauge data; it
% is the prior expected surface that brings the interesting features into
% play in this case. In practice, that expected attribute surface would be
% constructed from a regression model linking the dependent variable 
% (here precip) to several pertinent covariates or explanatory variables.

%% Kriging with known locally varying mean (simple kriging)
% The multiple steps above (for a given residual semivariogram model) can
% be condensed as: 
precipSKexp = krigep2p(1,GridSpecs,bayprecip,[2 3],[4 5],MdlRes.gslib,[ ],0,bayprecipexpect);
%% 
maprasterdata(GridSpecs,precipSKexp,1,'cont',cLims,plotPars);
title('SK interpolated precip with spatially varying expectation');
%%
% where the *residual* semivariogram model specification |MdlRes.gslib| 
% must be passed to function |krigep2p|.

%% Local conditional distributions
% One can generate simulated attribute values at each prediction
% location via Monte Carlo simulation from the corresponding local 
% Gaussian conditional distribution. This can be done with function 
% |postkrige| by setting |iMode = 4| and providing an array |simPars =
% [nSim seed]| with the # of realizations to generate and the random #
% seed:
iMode = 4;
simPars = [50 4567];
precipSKexpSim = postkrige(precipSKexp,iDistr,iMode,simPars);
disp( size(precipSKexpSim) );
%%
maprasterdata(GridSpecs,precipSKexpSim,1,'cont',[0 15],plotPars);
title('Precip realization #1');
figure
maprasterdata(GridSpecs,precipSKexpSim,10,'cont',[0 15],plotPars);
title('Precip realization #10');
% The mean of the simulated realizations tends towards the SK predictions
% and the variance of the realizations tends towards the SK error variance,
% both stored in array |precipSKexp|, as the # of realizations tends to 
% infinity... 
%
% The simulated attribute values above are generated *independently*
% from one prediction location to another, hence do not account properly
% for spatial auto-correlation in the *residuals* and lead to "speckly" 
% attribute realizations.