Example Krige interpolation

ExampleKrig.R

## Example spatial interpolation

# GIS libraries 
library(rgdal)
library(rgeos)
library(raster)

# kriging
library(gstat)

## Custom functions
vario.sum=function(vario,vario.fit.ls,max.cutoff){
  # This function extracts variogram fit information (or info that is useful)
  rslt=data.frame()
  for(i in 1:length(vario.fit.ls)){
    vario.fit=vario.fit.ls[[i]]
    fitted=variogramLine(vario.fit,maxdist=max.cutoff,dist_vector = vario$dist)        
    residual=fitted$gamma-vario$gamma      
    
    # r2 of variogram
    mss=with(vario,sum((gamma-mean(gamma))^2))
    rss=sum(residual)^2
    # mss/(mss+rss)
    
    # RMSE like lm()
    rdf=length(residual)-1
    resvar <- rss/rdf
    sigma=sqrt(resvar)
    # sigma
    
    # sigma=sqrt(mean(residual^2))
    
    # mean error, ideally 0:
    merror=mean(residual)
    # RMSE
    # rmse=sqrt(mean(residual^2))
    
   
    # partial sill to total sill
    sill.ratio=vario.fit$psill[2]/sum(vario.fit$psill)*100
    # vario.fit$psill[2] # sill
    # vario.fit$range[2] # range
    # vario.fit$psill[1] # nugget
    
    sill_nugget=as.numeric(vario.fit$psill[2])/as.numeric(vario.fit$psill[1])
    
    loglike=loglike.fun(vario$gamma,fitted$gamma)
    
    AIC=AIC.fun(loglike,p=2)
    AICc=AICc.fun(loglike,n=length(residual),p=2)
    BIC=BIC.fun(loglike,n=length(residual),p=2)
    
    tmp=data.frame(Model=as.character(vario.fit$model[2]),
                   Nugget=as.numeric(vario.fit$psill[1]),
                   Sill=as.numeric(vario.fit$psill[2]),
                   Range=as.numeric(vario.fit$range[2]),
                   sill.ratio=sill.ratio,
                   sill_nugget=sill_nugget,
                   ME=merror,
                   rmse=sigma,
                   r2=mss/(mss+rss),
                   loglike=loglike,
                   AIC=AIC,
                   AICc=AICc,
                   BIC=BIC)
    rslt=rbind(tmp,rslt)
  }
  
  if(length(vario.fit.ls)>1){
    rslt$AICcDelta=with(rslt,AICc-min(AICc))
    rslt$ModelLik=with(rslt,exp(-0.5*AICcDelta))
    rslt$AICcWt=with(rslt,ModelLik/sum(ModelLik))
    rslt=rslt[order(rslt$AICcWt),]
    rslt$cumwt=cumsum(rslt$AICcWt)
  }
  return(rslt)
  
}

krig.error.fun=function(actual,pred){
  # using kriged data and actual measured data, this function
  # determines the rmse and r2 of kriging preditions
  res=actual-pred
  mss=sum((pred-mean(pred,na.rm=T))^2,na.rm=T)
  rss=sum(res,na.rm=T)^2
  r2=mss/(mss+rss)
  rmse=sqrt(mean((res)^2,na.rm=T))
  return(data.frame(R2=r2,RMSE=rmse))
}

## I usually re-project my data into UTMs to avoid arc sec/dec deg units
## define coordinate reference systems
nad83.pro=CRS("+init=epsg:4269")
utm17=CRS("+init=epsg:26917")


# lake outline shapefile
# read data
data.folder="./GISData" # this is a folder where my data lives
## make sure the path doesn't have a "/" at the end. 
lake=readOGR(data.folder,"LakeOkeechobee_general")
# reproject the data
lake=spTransform(lake, utm17)

# sediment sample data (if a shapefile already)
sed.dat=readOGR(data.folder,"sediment_data")
sed.dat=spTransform(sed.dat, utm17)

## if not a shapefile
# sed.dat=read.csv(...)
# sed.dat=SpatialPointsDataFrame(coords=sed.dat[,c("Long","Lat")],data=sed.dat,proj4string = nad83.pro)
# sed.dat=spTransform(sed.dat, utm17)

## Make grid for kriging
## this makes a grid at the extent of the lake at 100 x 100 m resolution
grd=as.data.frame(spsample(lake,type="regular",cellsize=c(100,100)))
names(grd)       = c("UTMX", "UTMY")
coordinates(grd) = c("UTMX", "UTMY")
plot(grd)
plot(lake,add=T)
gridded(grd)=T;fullgrid(grd)=T
proj4string(grd)=utm17

## make sure you don't have double values (i.e. two values per station)
## This example first subsets the data for the 0 - 10 cm segment
## if your data isn't set up like this then diregard the subset(..)
## Ordinary kriging
surf.sed=subset(sed.dat,Depth=="0-10")
# set coordinates to match grid coordinates names
surf.sed$UTMX=coordinates(surf.sed)[,1]
surf.sed$UTMY=coordinates(surf.sed)[,2]
colnames(surf.sed@coords)=  c("UTMX", "UTMY")

## Kriging
f.1 = as.formula(TP ~ 1)
# basic model fit
v.TP=variogram(f.1 ,surf.sed ,cloud = F)
# you can adjust variogram (width[bins] and cutoff[max distance])
v.TP=variogram(f.1 ,surf.sed ,cloud = F,width=900,cutoff=10000)

# visualize variogram
plot(gamma~dist,v.TP)

## You can fit a bunch or different models
## the vgm(...) function you can specify specific sill, nugget, range, etc.

v.fit.sph = fit.variogram(v.TP,vgm("Sph"),fit.kappa=T)
v.fit.lin = fit.variogram(v.TP,vgm("Lin"),fit.kappa=T)
v.fit.exp = fit.variogram(v.TP,vgm("Exp"),fit.kappa=T)
v.fit.gau = fit.variogram(v.TP,vgm("Gau"),fit.kappa=T)
v.fit.ste = fit.variogram(v.TP,vgm("Ste"),fit.kappa=T)
v.fit.mat = fit.variogram(v.TP,vgm("Mat"),fit.kappa=T)

## adds fitted line to the plot from earlier
with(variogramLine(v.fit.sph,max(v.TP$dist)),lines(gamma~dist,col="indianred1"))
with(variogramLine(v.fit.lin,max(v.TP$dist)),lines(gamma~dist,col="blue"))
with(variogramLine(v.fit.exp,max(v.TP$dist)),lines(gamma~dist,col="green"))
with(variogramLine(v.fit.gau,max(v.TP$dist)),lines(gamma~dist,col="black"))
with(variogramLine(v.fit.ste,max(v.TP$dist)),lines(gamma~dist,col="grey"))
with(variogramLine(v.fit.mat,max(v.TP$dist)),lines(gamma~dist,col="forestgreen"))

# use the vario.sum function to evalute each model fit. 
TP.rslts=vario.sum(v.TP,
                   list(v.fit.sph,v.fit.lin,v.fit.exp,v.fit.gau,v.fit.ste,v.fit.mat),
                   max(v.TP$dist))
TP.rslts

# Krig (assume v.fit.sph is the best model)
dat.krg <- gstat::krige(f.1, surf.sed, grd, v.fit.sph)
data.r <- raster(dat.krg)
data.r.m <- mask(data.r, lake); # masks raster to lake outline
data.r.m

plot(data.r.m)

# extracts predicted values to shapefile
surf.sed$pred.TP=extract(data.r.m,surf.sed)

TP.krig.error=with(surf.sed@data,krig.error.fun(TP,pred.TP)))