## R: Using matrix algebra to compute a Kriged value

15 October, 2010 (14:04) | R, Spatial analysis | By: Manuel Gimond

The following example was created in R and was inspired from O’Sullivan and Unwin’s Geographic Information Analysis (1st edition)

```library(maptools)
#
# Convert data to spatstat format
SP  = as(S, "SpatialPoints")# This only grabs the X,Y values
SPz = slot(S,"data") # This grabs all attribute values
P   = as(SP, "ppp")  # This creates a spatstat point pattern
marks(P) = SPz       # This adds the attribute values to P
# thus creating a marked point pattern
hist(marks(P)\$z)     # plots histogram of attribute z
```

```# Convert data to dataframe
P.frame     = as.data.frame(P)
P.frame.sub = data.frame(x=P.frame\$x,y=P.frame\$y,z=P.frame\$z)

# Get individual columns
X = P.frame\$x
Y = P.frame\$y
Z = P.frame\$z

# Now sort the vectors to match the order used in O-Sullivan & Unwinâ€™s Book
X = rev(X)
Y = rev(Y)
Z = rev(Z)

# Create distance matrix
A = dist(cbind(X,Y),method="euclidean",diag=T)

# Replace all distance values by the semivariance at the distance
P.gamma = function(X,range,sill) {P.gamma=sill*(3*X/(2*range)-0.5*
(X/range)^3);P.gamma}
A = P.gamma(A,95,100)

# Augment matrix by one row and one column
augCol = matrix(rep(1,7),nrow=7,ncol=1)
augRow = matrix(c(rep(1,7),0),1,8)
A      = cbind(as.matrix(A),augCol)
A      = rbind(as.matrix(A),augRow)

# Inverse matrix A
A.inv = solve(A)
# Create distance matrix (single column) between unknown value and all control points
# Unknown value is located at 50,50 in this example
X2 = c(50,X)
Y2 = c(50,Y)
b  = as.matrix(as.matrix(dist(cbind(X2,Y2),method="euclidean",diag=T))[,1]) # We only need dist for unk.
# Replace distance values by the semivariance
b = P.gamma(b,95,100)
# Now remove first record and augment by 1
b = as.matrix(rbind(as.matrix(b[-1,]),1))
# Finally, compute the weights
w = A.inv %*% b   # note that all but the last row sum to 1
# The unknown Z value at 50,50 can now be computed from these weights
Z.est = sum(Z * w[-8,]) # equals 10.36

# The estimated variation is computed from
est.var = sum(w*b)

# The standard error is
std.err = est.var^0.5 # equals 5.696

# To estimate a 95% confidence interval, set limits a 1.96 times std error
upper.range = 1.96 * std.err + Z.est # = 21.53
lower.range = 1.96 * std.err - Z.est # = 0.80
```