Customize inverse distance predictions within raster classes in R

Question

I want to modify the following code so that for all values beyond the defined maximum distance I can use the average values of points which fall in the same soil order category. Any recommendations?

# load packages and data
library(gstat)
library(raster)
data(meuse)
data(meuse.grid)
################## 
meuse <- meuse[110:155,]
meuse <- na.omit(meuse)
coordinates(meuse) <- ~x+y
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
spplot(meuse.grid[4], main="soil class")
# inverse distance prediction for maximum distance of 2km.
ca.idw = idw(cadmium ~ 1, meuse, meuse.grid, omax =5, maxdist = 2000)

Answer 1

One way is to add the soil category as a (factor) predictor, and use ordinary kriging with variogram model that reaches a sill (like the Spherical model), and a small range

ca.ok = krige(cadmium ~ soil, meuse, meuse.grid, vgm(1, "Sph", 100))

which gives:

Another way would be to do separate ordinary krigings for each soil category. The difference is that in the former model, residuals are taken into account for interpolation across boundaries (the residual interpolation is global, ignores boundaries), in the latter they are not.