#this code will guide you through median polish kriging,
#detrending using regression-like methods followed by ordinary kriging,
#detecting anisotropy

#as usual, load geoR and akima
library(geoR)
library(akima)

#one of the data sets used as illustration is the Wheat data
wheat<-read.table("http://www.public.iastate.edu/~pcaragea/S40608/Compute/wheat.txt",header=TRUE)
attach(wheat)


#image plot
#you can control the smoothness of the interpolation by 
#changing xo and yo. Here I wanted to use the original data grid
#if data is not on a regular lattice, this will not work as nice.
wheat.raw.i=interp(wheat[,1],wheat[,2],grain,xo=unique(u),yo=unique(v))
image(wheat.raw.i,col=gray(seq(0,1,length=20)))
title(main="Raw grain data")

wheatg<-as.geodata(wheat[,c(1,2,3)])#first two columns are u and v, the third is the yield measured by "grain"

plot(wheatg)#reality check, plot the data and get some summaries
points(wheatg,cex.min=0.3,cex.max=2,pch=20)
title("Wheat Data")
summary(wheatg)

#empirical robust variogram
wheat.variog1<-variog(wheatg,max.dist=10,method="robust")
plot(wheat.variog1,pch=20)
title("Omnidirectional variogram for Original Data",cex=1.5)

#directional variograms, default directions
wheat.variog2<-variog4(wheatg,max.dist=10,method="robust")
plot(wheat.variog2,lty=7,lwd=3)
title("Directional Variograms for Original Data",cex=1.5)

#here is how you would compute directional variograms, one direction at a time
par(mfrow=c(2,2))#plots a 2 by 2 martix of graphs
wheat.variogNS<-variog(wheatg,max.dist=15,method="robust",dir=0,unit.angle="deg")
plot(wheat.variogNS,pch=20)
title("N-S")
wheat.variogEW<-variog(wheatg,max.dist=15,method="robust",dir=90,unit.angle="deg")
plot(wheat.variogEW,pch=20)
title("E-W")
wheat.variogSWNE<-variog(wheatg,max.dist=15,method="robust",dir=45,unit.angle="deg")
plot(wheat.variogSWNE,pch=20)
title("SW-NE")
wheat.variogSENW<-variog(wheatg,max.dist=15,method="robust",dir=135,unit.angle="deg")
plot(wheat.variogSENW,pch=20)
title("SE-NW")

par(mfrow=c(1,1))#returns to usual plotting, one at a time

#would like to remove the obvious trend in the data, 
#but would not like to try to explain it. Median polish.

#To do this in R, I need to first arrange the data in matrix format. To repeat this for another data set, you will need to replace "wheatg" by the new data name
wheat.mat<-tapply(wheatg$data,list(factor(wheatg$coords[,1]),factor(wheatg$coords[,2])),function(x)x)

#apply median polish. The function is medpolish.
#you need to specify the data set in matrix format (by locations)
#if you do not remove missing values you will run into problems
wheat.med<-medpolish(wheat.mat,na.rm=T)
#this function returns several quantities: the overall effect ($overall), thee row effect ($row), the column effect ($col), the residuals ($residuals). We will need the "trend" part later on. One way to get to that is to extract the residuals from the original data, by typing:
wheat.trend.m<-wheatg$data-wheat.med$resid

#rearange the trend into a vector and plot
wheat.trend<-na.omit(as.vector(wheat.trend.m))
#to plot, interpolate first
wheat.trend.interp=interp(wheat[,1],wheat[,2],wheat.trend,xo=unique(u),yo=unique(v))
#now a plot of the trend surface, finally!!
image(wheat.trend.interp,col=gray(seq(0,1,length=20)))
title(main="Median polish trend")


#arranging the residuals back into geostat data format
#first need to put them in a vector
wheat.res<-na.omit(as.vector(wheat.med$res))
wheat.res.g<-as.geodata(cbind(wheatg$coords,wheat.res))
plot(wheat.res)
#or interpolate for an image plot
wheat.res.interp=interp(wheat[,1],wheat[,2],wheat.res,xo=unique(u),yo=unique(v))
#now plot, finally!!
image(wheat.res.interp,col=gray(seq(0,1,length=20)))
title(main="Median polish residuals")

summary(wheat.res)#can't see any trend
wheat.finalvar<-variog(wheat.res.g,max.dist=25,estimator.type="modulus")
plot(wheat.finalvar,xlab="Distance",ylab="Gamma",pch=20)
title("MP Residuals Omnidirectional Variogram",cex=1.5)

wheat.variog.d<-variog4(wheat.res.g,max.dist=10,estimator.type="modulus")
plot(wheat.variog.d,lty=7,lwd=2)
title(main="MP Residuals Directional Variograms")
#looks like a pure nugget variogram now

#if residuals were spatially distributed, you could
#krig residuals and add back to trend surface to get final map


#removing trend using rergression
reg1<-lm(grain~u+v)#regression of grain on both u and v coordinates
summary(reg1)#looks like v should not be in the model

reg2<-lm(grain~u)
summary(reg2)

#extract trend and plot
grain.fit=reg2$fit
grain.fit.interp=interp(wheat[,1],wheat[,2],grain.fit,xo=unique(u),yo=unique(v))
#now a plot of the trend surface, finally!!
image(grain.fit.interp,col=gray(seq(0,1,length=20)))
title(main="Linear trend")

grain.res<-reg2$res#extract residuals and plot
grain.res.interp=interp(wheat[,1],wheat[,2],grain.res,xo=unique(u),yo=unique(v))
#now a plot of the trend surface, finally!!
image(grain.res.interp,col=gray(seq(0,1,length=20)))
title(main="Residuals from Linear trend")

#make it geodata
wheatres.reg.g<-as.geodata(cbind(wheat[,1],wheat[,2],grain.res))
plot(wheatres.reg.g)

#look at directional variograms
grain.res.reg.variog2<-variog4(wheatres.reg.g,max.dist=25,estimator.type="modulus")
plot(grain.res.reg.variog2,lty=7,lwd=2)
title(main="Linear trend Residuals Directional Variograms")

#omnidirectional variogram
wheatres.reg.finalvar <- variog(wheatres.reg.g, max.dist=25, method="robust")
plot(wheatres.reg.finalvar,pch=20)
title(main="Linear trend residuals Omnidirectional Variogram")
#fit variogram to residuals
vg.res=variofit(wheatres.reg.finalvar,ini.cov.pars=c(1,0),cov.model="sph",weights="npairs")
#krig residuals and add back to trend surface to get final map


#universal kriging: model trend surface and spatial dependence simultaneously
#assume you wanted to krig on a fine grid
#define the grid in the usual way
xg=seq(min(wheatg$coords[,1]),max(wheatg$coords[,1]),length=40)
yg=seq(min(wheatg$coords[,2]),max(wheatg$coords[,2]),length=40)
new.locs<-expand.grid(xg,yg)
#it's easiest to use ksline function here
#when you specify "kt" under m0 you are asking for Universal kriging. If you don't specify anything for m0, it will do ordinary kriging
#for "kt" you need to specify a trend. Possible values are 1 for a linear trend on the coordinates, or 2 (quadratic trend)
#you could also specify "kte" if you wanted to fit a model with an external covariate.
#in this case, you would need to also specify ktedata (a vector with the values of the covariate at the original locations) and ktelocations (a vector with the values of the external trend (covariates) at the prediction locations)
wheat.UK=ksline(wheatg,locations=new.locs,cov.pars=vg.res$cov.pars,nugget=vg.res$nugget,cov.model="sph",m0="kt",trend=1,ktedata=NULL,ktelocations=NULL)

image(wheat.UK,xlab="X",ylab="Y",col=gray(seq(0,1,l=40)))
image(wheat.UK,val=sqrt(wheat.UK$krige.var),xlab="X",ylab="Y",col=gray(seq(0,1,l=40)))
#note here that the surface of standard errors is basically flat with several peaks on the edges. I guess this has something to do with the rather smooth trend surface and the pure nugget form of the variogram for the residuals. This is more of an exception than the rule.