##########################
## plot.gd
###########################
## plot.gd is a function to plot a geodata object
## default (col=T) is to plot in reverse heat colors(red is high)
## otherwise plots in grey scale
## size of circles proportional to Z values
##########################################

plot.gd <- function(data.gd, col=T){
	if(col){
		points(data.gd, col=rev(heat.colors(12)), xlab=NA, ylab=NA)
		}
	if(!col){
		points(data.gd, col="grey", xlab=NA, ylab=NA)
		}
	}


###########################################
## variog.508 is a function to compute an empirical variogram
## data.gd is a geodata object
## max.dist defaults to half of max pairwise distance
## number of bins (nbins) defaults to 20
## Computes classical variogram by default (type = "class")
## set type = "mod" for Cressie's robust empirical variogram
## pairs.min = 10 by default
##########################################

variog.508 <- function(data.gd, nbins=20, max.dist=NA, type="class", pairs.min = 10, plots=F){
	library(geoR)
	if(is.na(max.dist)){
		max.dist <- 0.5*max(dist(data.gd$coords)) # max.dist = half maximum distance
		}
	if(type=="class"){
		type = "classical"
		}
	if(type=="mod"){
		type = "modulus"
		}
	data.vg <- variog(data.gd, max.dist=max.dist, uvec = nbins, pairs.min=pairs.min, option="bin", estimator.type=type, messages=F)

	if(plots){
		plot.vg(data.vg, err=T)
		}
	return(data.vg)
	}


##########################
## plot.vg
###########################
## plot.vg is a function to plot the empirical variogram
## data.vg is a calculated empirical variogram
## err indicated whether to plot error bars (95% CI), default = F
## pairs indicates whether to plot number of pairs in each bin, default = T,
##     but pairs only plotted if err = T
##########################################

plot.vg <- function(data.vg, err=F, pairs=T) {
	library(geoR)
	plot(data.vg,pch=16, col="darkblue",cex=1.2, xlab="lag distance (h)",ylab="gamma")

	if(err){
		library(Hmisc)
		yup <- data.vg$v + 1.96*sqrt(2)*data.vg$v/sqrt(data.vg$n)
		ydn <- data.vg$v - 1.96*sqrt(2)*data.vg$v/sqrt(data.vg$n)
		errbar(data.vg$u, data.vg$v,yup,ydn, pch=16, col="darkblue", cex=1.2, errbar.col="darkblue", add=T)
		if(pairs){
			text(data.vg$u, data.vg$v,labels=data.vg$n,pos=1,offset=.5, cex=0.8, col="darkred", font=2)
			}
		}
	}


###########################################
## rand.vg
###########################################
## rand.vg is a function to plot the empirical variogram with non-free randomization envelope
## Test for significance of autocorrelation
## data.gd is the geodata object used to create the variogram
## data.vg is the computed empirical variogram
## nsim = number of simulations (default is 99)
##########################################

rand.vg <- function(data.gd, data.vg, nsim=99){
	data.env<-variog.mc.env(data.gd, obj.variog=data.vg, nsim=99, messages=F)
	plot(data.vg,pch=16, col="darkblue",cex=1.2, envelope=data.env)
	}


##########################
## Parametric Bootstrap
##########################
bootstrap.vg <- function(data.gd, data.vg, data.fit, nsim=19){

	theta <- data.frame(matrix(NA, nrow=nsim, ncol=3) )   
	colnames(theta) <- c("range", "nugget/sill", "sill")

	i <- 1
	j <- 1
	while(i < (nsim+1))  {   # run nsim simulations
		j <- j + 1
		simdat <- grf(grid = data.gd$coords, cov.model=data.fit$cov.model, cov.pars=data.fit$cov.pars, nugget=data.fit$nugget,messages=F, RF=F)
		simdat.lm <- lm(simdat$data ~ simdat$coords[,1]*simdat$coords[,2])
		if(summary(simdat.lm)$r.squared > 0.1) {next}

		simdat.vg <- variog(simdat,uvec=data.vg$uvec,messages=F)
		simdat.fit <- variofit(simdat.vg, ini.cov.pars=data.fit$cov.pars, nugget=data.fit$nugget, cov.model=data.fit$cov.model, weights="cressie", messages=F)
		theta[i,] <- round(c(simdat.fit$cov.pars[2],simdat.fit$nugget/simdat.fit$cov.pars[1],simdat.fit$cov.pars[1]),3)
		i <- i+1
		}
	print(paste("Fraction of simulations with 1st order stationarity: ",round(i/j,2),sep=""))
	return(rbind(apply(theta, MARGIN=2, FUN=min), apply(theta, MARGIN=2, FUN=max)))
	}


###########################################
## check.vg
###########################################
## provides information to evaluate empirical variogram assumptions 
## Input geodata object
## output points plot
## output qq-normal plot
## output directional variogram to half of maximum distance
## output omnidirectional variogram with non-free randomization envelope (nsim = 99)
## output adjusted r-squared value from linear trend model fit
##      (parametric assumptions may not be met, so r-squared only approximate)
##########################################

check.vg <- function(data.gd) {
   library(geoR)
   mar.default<-par("mar")  # save defaults
   par(mfrow=c(2,2), mar=c(4,4,1,0))	# sets up to plot in 2 by 2 panels with no margins  

   plot.gd(data.gd)   
   qqnorm(data.gd$data)
   data.vg <- variog.508(data.gd, pairs.min=1, plots=F)
   plot(variog4(data.gd, max.dist=data.vg$max.dist, messages=F), lwd=3, legend=F)
   legend(x="bottomright", legend=c("0","45","90","135"),lty = 1:4, col=1:4, lwd=3)
   rand.vg(data.gd,data.vg,nsim-99)
   par(mfrow=c(1,1), mar=mar.default)	# returns plot in 1 panel and regular margins
   data.lm<-lm(data.gd$data ~ data.gd$coords[,1]*data.gd$coords[,2])
   cat("linear trend surface approximate r-squared = ",round(summary(data.lm)$adj.r.squared,2),"\n",sep="")
   }