library(RTMB)
library(vegan)
data(mite)
data(mite.env)
data(mite.xy)
xy<-mite.xy

dat <- list(Y=mite[,31], X=mite.env$SubsDens, D=as.matrix(dist(xy)))
param <- list(mu=0, alpha=0, logPhi=0, logSigma=0, omega=numeric(nrow(mite)))

f <- function(par){
  getAll(par,dat)
  phi <- exp(logPhi)
  sigma <- exp(logSigma)
  S <- sigma^2*exp(-(D/phi)^2)
  nll <- -dmvnorm(omega,mu=0,Sigma=S,log=TRUE)
  lambda <- exp(mu+alpha*X+omega)
  nll <- nll -sum(dpois(Y, lambda, log=TRUE))
  nll  
}

obj <- MakeADFun(f, param, random = "omega", silent=TRUE)

fit <- nlminb(obj$par, obj$fn, obj$gr)

sdr <- sdreport(obj)
pl <- as.list(sdr, "Est")
plsd <- as.list(sdr, "Std")

cc <- viridis::mako(500)
lamCol<-cc[as.numeric(cut(pl$omega,breaks = 501))]
plot(xy, col=lamCol, pch=16, cex=5)