#GETTING STARTED WITH geoR#

                                        #Step 1: Always set the seed.# 
set.seed(1234)

                                        #Step 2: Read a data set. We assume
                                        #that the source file (this file) is in the same directory
                                        #as the data file "scallops.txt".#
scallops <- read.table("scallops.txt", header=T)

                                        #Step 3: Draw a simple histogram#
hist(scallops$tcatch)

                                        #Step 4: Seeing the skewness, we want to use the log of the catches
                                        #Create dataframe "myscallops" with "lgcatch" as a new column
myscallops <- scallops
myscallops[,"lgcatch"] <- log(scallops$tcatch+1)
summary(myscallops$lgcatch)

                                        #Step 5: Draw another histogram of the log-catch
hist(myscallops$lgcatch)

                                        #Step 6: Plot the sites#
plot(myscallops$long, myscallops$lat)

                                        #Step 7: Exploratory tools: Image, contour and perspective plots#
library(akima)
int.scp <- interp.new(myscallops$long,myscallops$lat,myscallops$lgcatch)
#lat.lim_c(min(myscallops$lat),max(myscallops$lat))
#lon.lim_c(min(myscallops$lon),max(myscallops$lon))
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))

                                        #Step 8: Spatial analysis with geoR#
library(geoR)

                                        #Step 8a: Convert lat-lon to planar coordinates using naive scaling#
R <- 6371
lat <- myscallops$lat*pi/180
lon <- myscallops$long*pi/180

u <- cbind(R*lon, R*lat)
#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)

                                        #Step 8b: Create "geodata" object as dataframe for geoR#
scallops.geo <- as.geodata(obj,coords.col=1:2,data.col=3)

                                        #Step 8c: Create classical and robust variograms#
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)

                                        #Step 8d: Estimate variogram parameters using weighted least squares#
scallops.var.fit <- variofit(scallops.var.robust,ini.cov.pars=c(1.0,2.0),cov.model="exponential", fix.nugget=FALSE,nugget=1.0)
summary(scallops.var.fit)

                                        #Step 8e: Estimate parameters using maximum likelihood estimates#
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)
summary(scallops.lik.fit)
scallops.lik.fit$parameters.summary
scallops.lik.fit$beta.var # This produces the variance of estimated beta

                                        #Step 8f: Create locations to predict#
u.loci <- rbind(c(-73.0,39.5),c(-72.5,40.0))
u.loci <- (R*pi/180)*u.loci

                                        #Step 8g: Carry out conventional kriging#
scallops.krige.conv <- krige.conv(scallops.geo, locations=u.loci, borders=NULL, krige=krige.control(cov.pars=c(1,0.25)))

                                        #Step 8h: Carry out Bayesian kriging#
#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))