scallops <- read.table("scallops.txt", header=T)

hist(scallops$tcatch)

myscallops <- scallops
myscallops[,"lgcatch"] <- log(scallops$tcatch+1)
summary(myscallops$lgcatch)
hist(myscallops$lgcatch)

plot(myscallops$long, myscallops$lat)

library(akima)
int.scp <- interp.new(myscallops$long,myscallops$lat,myscallops$lgcatch)

lat.lim <- range(myscallops$lat)
lon.lim <- range(myscallops$lon)
image(int.scp, xlim=lon.lim, ylim=lat.lim, xlab="Longitude", ylab="Latitude")
contour(int.scp, add=T)
#persp(int.scp,xlim=range(myscallops$lon),ylim=range(myscallops$lat))

library(geoR)

R <- 6371
lat <- myscallops$lat*pi/180
lon <- myscallops$long*pi/180

u <- cbind(R*cos(lat)*cos(lon), R*cos(lat)*sin(lon), R*sin(lat))

#obj <- cbind(myscallops$long,myscallops$lat,myscallops$lgcatch)
obj <- cbind(u,myscallops$lgcatch)

scallops.geo <- as.geodata(obj,coords.col=1:3,data.col=4)

scallops.var <- variog(scallops.geo, estimator.type="classical")
scallops.var.robust <- variog(scallops.geo, estimator.type="modulus")

par(mfrow=c(2,1))
plot(scallops.var)
plot(scallops.var.robust)

scallops.var.fit <- variofit(scallops.var.robust,ini.cov.pars=c(1.0,2.0),cov.model="exponential", fix.nugget=FALSE,nugget=1.0)

scallops.lik.fit <- likfit(scallops.geo,ini.cov.pars=c(1.0,2.0),cov.model="exponential",trend="cte", fix.nugget=FALSE,nugget=1.0, nospatial=FALSE)

scallops.krige.bayes <- krige.bayes(scallops.geo, locations="no", borders=NULL, model=model.control(trend.d = "cte", cov.model="matern", kappa=0.5, lambda=1), prior=prior.control(beta.prior="flat",sigmasq.prior="reciprocal",tausq.rel.prior="uniform",tausq.rel.discrete=seq(from=0.0,to=1.0,by=0.01)))

out <- scallops.krige.bayes$posterior
out <- out$sample
beta.qnt <- quantile(out$beta, c(0.50,0.025,0.975))
phi.qnt <- quantile(out$phi, c(0.50,0.025,0.975))
sigmasq.qnt <- quantile(out$sigmasq, c(0.50,0.025,0.975))
tausq.rel.qnt <- quantile(out$tausq.rel, c(0.50,0.025,0.975))