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

xygrid <- expand.grid(x=seq(0,3, length=10), y=seq(0,10, length=11))
xy <- rbind(xy,xygrid)

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

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[!is.na(Y)], lambda[!is.na(Y)], 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[is.na(dat$Y),], col=lamCol[is.na(dat$Y)], pch=16, cex=5)