Skip to content

Instantly share code, notes, and snippets.

@SwampThingPaul
Last active September 26, 2022 16:58
Show Gist options
  • Save SwampThingPaul/504d839a28f1d76b2745397894847a2f to your computer and use it in GitHub Desktop.
Save SwampThingPaul/504d839a28f1d76b2745397894847a2f to your computer and use it in GitHub Desktop.
Example Krige interpolation
## 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)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment