Spatial Data - MATLAB Recipes for Earth Sciences - page 183

S 2 GR E

1

_ –

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 fi ve 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 fi rst command computes the distance between the grid points