#load packages
require(spatial)
require(MASS)
library(gstat)
require(rgdal)
require(sp)
#grid
xy <- expand.grid(1:100, 1:100)
names(xy) <- c("x","y")

# create spatially autocorrelated data
g <- gstat(formula=z~1, locations=~x+y, dummy=T, beta=1,
           model=vgm(psill=0.1, model="Exp", range=10), nmax=10)

# simulate
yy <- predict(g, newdata=xy, nsim=1)
coordinates(yy)=~x+y
spplot(yy, col.regions = rainbow(7)) #plot

#derive autocorrelation structure
a<-data.frame(yy)
colnames(a)<-c("x","y","z")
a <- surf.ls(2, a) #fits a trend surface by least-squares
cc<-correlogram(a, 100, xlim=c(0,100), ylab="Autocorr",
                xlab="Distance") #plot correlogram
cc #details

library(sp)
library(spdep)
data(meuse)
coordinates(meuse) <- ~x+y

k.nn <- spdep::knn2nb(spdep::knearneigh(coordinates(meuse)))
  bw <- max(unlist(spdep::nbdists(k.nn, coordinates(meuse))))
( AutoCov <- spdep::autocov_dist(meuse@data[,"copper"], xy = coordinates(meuse),
                                 nbs = bw, style = "W", type = "inverse") )

nbd <- nbdists(k.nn, coordinates(meuse))
gl <- lapply(nbd, function(x) 1/x)
lw <- nb2listw(k.nn, glist=gl)
lags <- lag(lw, meuse$copper)
all.equal(lags, AutoCov)