Page 193 - MATLAB Recipes for Earth Sciences

188 7 Spatial Data

S 2 _ 1 – GR E

For our calculations with MATLAB we need the matrix of coeffi cients de-
rived from the distance matrix D and a variogram model. D was calculated
in the variography section above and we use the exponential variogram
model with nugget, sill and range from the previous section:
G_mod = (nugget + sill*(1 - exp(-3*D/range))).*(D>0);

Then we get the number of observations and add a column and row vector of
all ones to the G_mod matrix and a zero at the lower left corner:

n = length(x);
G_mod(:,n+1) = 1;
G_mod(n+1,:) = 1;
G_mod(n+1,n+1) = 0;
Now the G_mod matrix has to be inverted:

G_inv = inv(G_mod);
A grid with the locations of the unknown values is needed. Here we use a
grid cell size of five within a quadratic area ranging from 0 to 200 in x and y

direction, respectively. The coordinates are created in matrix form by:
R = 0:5:200;
[Xg1,Xg2] = meshgrid(R,R);
and converted to vectors by:

Xg = reshape(Xg1,[],1);
Yg = reshape(Xg2,[],1);
Then we allocate memory for the kriging estimates Zg and the kriging vari-
ance s2_k by:
Zg = Xg*NaN;
s2_k = Xg*NaN;
Now we are kriging the unknown at each grid point:
for k = 1:length(Xg)
DOR = ((x - Xg(k)).^2+(y - Yg(k)).^2).^0.5;
G_R = (nugget + sill*(1 - exp(-3*DOR/range))).*(DOR>0);
G_R(n+1) = 1;
E = G_inv*G_R;
Zg(k) = sum(E(1:n,1).*z);
s2_k(k) = sum(E(1:n,1).*G_R(1:n,1))+E(n+1,1);
end
Here, the first command computes the distance between the grid points