library(gstat)
library(sp)
library(raster)
data(meuse)
data(meuse.grid)
coordinates(meuse) <- ~x + y
coordinates(meuse.grid) <- ~x + y

v <- fit.variogram(variogram(log(elev) ~ 1, meuse), vgm(.1, "Sph", 1000, .6))
elev <- krige(elev ~ 1, meuse, meuse.grid, v, nmax = 30)
  gridded(elev) <- TRUE
  names(elev@data) <- c("elev","var")
elev2 <- elev
  elev2@data$elev <- elev2@data$elev + runif(nrow(elev), -3, 3)

plot(stack(raster(elev,1),raster(elev2,1)))

rmse <- function(x, y) sqrt(mean((x - y)^2, na.rm=TRUE))

rmse(elev$elev, elev2$elev  ) # global RMSE

d <- 200
nb <- spdep::dnearneigh(sp::coordinates(elev), 0, d)
  xy.rmse <- rep(NA,nrow(elev))
    for (i in 1:length(nb)) {
      x.var <- elev@data[nb[[i]], ][1][, 1]
      y.var <- elev2@data[nb[[i]], ][1][, 1]
      xy.rmse[i] <- rmse(x.var,y.var)
    }

elev$rmse <- xy.rmse
plot(raster(elev,3))