#######################################################################
#### Krige
#### Wrapper function for OK, UK and KED kriging
#### Uses geoR package
#######################################################################
### data is a dataframe with measurement coordinates in cols 1:2, data in col 3, and any covariates in the remaining columns
### predict is a dataframe with prediction location coordinates in cols 1:2, and any covariates in the remaining columns
### cov.mod is a variofit object
### trend.mod is a lm model used to provide the form for the trend
### trend.mod defaults to "cte" (OK), with options for "1st" and "2nd" for 1st and 2nd order UK
### variables in data, predict, cov.mod and trend.mod must have consistent naming
#######################################################################
### returns kriging object
### prints beta parameter estimates
#######################################################################

krige <- function(data, predict, cov.mod, trend.mod="cte"){    # default to Ordinary Kriging
  	data <- cbind(data,data[,1:2])					# rearrange data to put coordinates in covariates
		data.gd <- as.geodata(data,covar.col=4:dim(data)[2], covar.names=colnames(data[,4:dim(data)[2]]))  # put into geodata object
		predict <- cbind(predict[,1:2],data=rep(NA,dim(predict)[1]),predict)  # rearrange predict to put NA in $data and coordinates in covariates
		 # put into geodata object
		predict.gd <- as.geodata(predict, covar.col=4:dim(predict)[2], covar.names=colnames(predict[,4:dim(predict)[2]]), na.action="none")   
	if(is.character(trend.mod)) {  # cte, 1st, 2nd for default OK and UK options
		trend.d <- trend.mod
		trend.l <- trend.mod
		}
	else  {
		trend.lm <- as.formula(paste("~",as.character(formula(trend.mod))[3]))  # right side of trend formula
		trend.d <- trend.spatial(trend=trend.lm, geodata=data.gd)		# trend for locations where we have data
		trend.l <- trend.spatial(trend=trend.lm, geodata=predict.gd)	# trend for prediction locations
		}
	krige.mod <- krige.control(cov.model=cov.mod$cov.model, cov.pars=cov.mod$cov.pars, nugget=cov.mod$nugget, trend.d=trend.d, trend.l=trend.l)
	krige.out<-krige.conv(data.gd, loc=predict.gd$coords, krige=krige.mod)
	print(krige.out$beta.est)
	return(krige.out)
	}

##########################################
#### grid.poly automatically creates a regular grid matrix inside a geoR polygon
#### grid.size is the pixel or grid size or spacing
#### borders is the borders component of a geodata object
############################################################

grid.poly<-function(borders, grid.size=1){
	max.grid <- ceiling(apply(borders, MARGIN=2, FUN=max)/grid.size)*grid.size  
# find max grid cell locations for both coordinates
	min.grid <- floor(apply(borders, MARGIN=2, FUN=min)/grid.size)*grid.size
# find min grid cell locations for both coordinates
grid.out <- expand.grid(seq(min.grid[1],max.grid[1], by=grid.size),
seq(min.grid[2],max.grid[2],by=grid.size))
	colnames(grid.out) <- colnames(borders)
	return(locations.inside(grid.out, borders))     
# return clipped grid, points inside borders
	}


##########################################
#### image.poly constructs an kriged image plot for a polygon-bounded geoR data object
#### krige.out  is a kriging output object
#### locations = prediction locations (n x 2 dataframe)
#### borders is the borders component of a geodata object, not necessary
#### values contains the plotting response variable, defaulting to krige.out$predict
#### 		(can also be sqrt(krige.out$krige.var))
############################################################

image.poly<-function(krige.out, locations, borders=NA, values=krige.out$predict){
	grid.full <- expand.grid(unique(locations[,1]),unique(locations[,2])) 	# get full retangular grid
	colnames(grid.full) <- colnames(locations)			# get coordinate names		
	grid.inside <- cbind(locations,values)				
	grid.all <- merge(x=grid.full, y=grid.inside, all=T, sort=F)
	grid.all <- grid.all[order(grid.all[,2],grid.all[,1]),]
	image(krige.out, loc=grid.all[,1:2], values=grid.all[,3], col=rev(heat.colors(12)))
	if(!is.na(borders)) {polygon(borders, lwd=2)}
	}