#### # ICES Cooperative Research Report # Handbook of Geostatistics in R for fisheries and marine ecology # # Based on # ICES Training Course # Application of Geostatistics to analyse spatially explicit Survey data # in an Ecosystem Approach # December 2013-2014, Fontainebleau, Centre de Geosciences, Mines ParisTech # # Demo script: Ordinary kriging, cokriging, colocated cokriging and # kriging with external drift #### { library(RGeostats) liste.initiale = ls() on.exit(demo.clean(liste.initiale)) par(ask=FALSE) demo.start(FALSE) # Load packages library(mapdata) library(maps) projec.toggle(0,0) # Loading data rg.load("Demo.herring.len.scot.db.data","db.data") rg.load("Demo.herring.len.scot.poly.data","poly.data") rg.load("Demo.herring.len.scot.grid.kriging","grid.depth") # Select points inside polygon db.data <- db.polygon(db.data,poly.data) cat("Nb points: ",db.data$nech," ; outside polygon: ",db.data$nmasked,"\n") # Define the projection projec.define(projection="mean", db=db.data) # Visualizing the data set plot(db.data,title="mean length",inches=5, asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), flag.proj=FALSE) plot(grid.depth,name.contour="depth",levels=seq(100,1500,100),col=8, flag.proj=FALSE,add=T) plot(poly.data,col=4,flag.proj=FALSE,add=T) map("worldHires",add=T,fill=T,col=8) demo.pause("Display Information") ## Ordinary kriging -------------------------------------------------------- # Basic statistics print(db.data,flag.stats=TRUE,names="m.length") hist(db.extract(db.data,"m.length"),xlab="m.length",main="",nclass=10,col=8) demo.pause("Histogram") # Omni-directional variogram & Model vario.data <- vario.calc(db.data) model.vario <- model.auto(vario=vario.data, struc=c("Nugget Effect","Spherical"),wmode=2, title="Variogram Model") demo.pause("Variogram Model") # Define the neighborhood neigh.kriging <- neigh.create(ndim=2,type=0) neigh.kriging # Define grid grid.kriging <- db.locate(grid.depth,"depth",loctype=NA) # View grid plot(grid.kriging,col=1,pch="+",asp=1,title="Estimation Grid") plot(db.data,col=2,pch=19,add=T) plot(poly.data,col=4,add=T) demo.pause("Display Grid") # Perform ordinary kriging grid.kriging <- kriging(dbin=db.data,dbout=grid.kriging,model=model.vario, neigh=neigh.kriging, uc=c("1"),mean=NA,calcul="point",radix="OK") # View ordinary kriging plot(grid.kriging,name.image="z1", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Mean Length - Estimation (OK)", pos.legend=5,flag.proj=F) plot(db.data,name.prop="m.length",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Mean Length - Estimation (OK)") # Kriging variance plot(grid.kriging,name.image="z2", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Mean Length - St. Dev. (OK)", pos.legend=5,flag.proj=F) plot(db.data,name.post="m.length",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Mean Length - St. Dev. (OK)") ## Homotopic cokriging -------------------------------------------------------- # Data management db.data <- db.locate(db.data,c("depth","m.length"),loctype="z") db.data # Basic statistics print(db.data,flag.stats=TRUE,names=c("depth","m.length")) hist(db.extract(db.data,"depth"),xlab="depth",main="",nclass=30,col=8) demo.pause("Histogram of Depth") hist(db.extract(db.data,"m.length"),xlab="m.length",main="",nclass=20,col=8) demo.pause("Histogram of Mean Length") correlation(db.data,name1="depth",name2="m.length",flag.regr=TRUE) demo.pause("Corr. of Depth vs. Mean Length") # Omni-directional variogram & Model vario.data <- vario.calc(db.data) model.vario <- model.auto(vario=vario.data,wmode=2, struc=c("Nugget Effect","Spherical","Spherical"), npairpt=F,npairdw=T,inches=0.08) demo.pause("Bivariate Model") # Define grid grid.kriging <- db.locerase(grid.kriging,"z") grid.kriging # Perform ordinary cokriging grid.kriging <- kriging(dbin=db.data,dbout=grid.kriging,model=model.vario, neigh=neigh.kriging, uc=c("1"),mean=NA,calcul="point",radix="CK") # View cokriging results plot(grid.kriging,name.image="z1", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Depth - Estimation (CK)", pos.legend=5,flag.proj=F) plot(db.data,name.prop="depth",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Depth Estimation (CK)") plot(grid.kriging,name.image="z2", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Mean Length - Estimation (CK)", pos.legend=5,flag.proj=F) plot(db.data,name.prop="m.length",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Mean Length Estimation (CK)") plot(grid.kriging,name.image="z3", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Depth - St. Dev. (CK)", pos.legend=5,flag.proj=F) plot(db.data,name.post="depth",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Depth St. Dev. (CK)") plot(grid.kriging,name.image="z4", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Mean Length - St. Dev. (CK)", pos.legend=5,flag.proj=F) plot(db.data,name.post="m.length",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Mean Length St. Dev. (CK)") # Heterotopic Cokriging ------------------------------------------------------- # 45pts + 5000pts won't work because of neighborhood issues... # Better use colocated cokriging # Colocated cokriging --------------------------------------------------------- # Define grid grid.kriging <- db.locerase(grid.kriging,"z") grid.kriging # Perform colocated cokriging grid.kriging <- kriging(dbin=db.data,dbout=grid.kriging,model=model.vario, neigh=neigh.kriging, uc=c("1"),mean=NA,calcul="point",radix="CCK", rank.colcok=c(4,NA)) # View colocated cokriging results plot(grid.kriging,name.image="z1", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Depth - Estimation (CCK)", pos.legend=5,flag.proj=F) plot(db.data,name.prop="depth",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Depth Estimation (CCK)") plot(grid.kriging,name.image="z2", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Mean Length - Estimation (CCK)", pos.legend=5,flag.proj=F) plot(db.data,name.prop="m.length",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Mean Length Estimation (CCK)") plot(grid.kriging,name.image="z3", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Depth - St. Dev. (CCK)", pos.legend=5,flag.proj=F) plot(db.data,name.post="depth",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Depth St. Dev. (CCK)") plot(grid.kriging,name.image="z4", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Mean Length - St. Dev. (CCK)", pos.legend=5,flag.proj=F) plot(db.data,name.post="m.length",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Mean Length St. Dev. (CCK)") # Kriging with external drift ------------------------------------------------- # Regression db.data <- regression(db.data,names="depth",name1="m.length",flag.draw=T) # Data management db.data <- db.locerase(db.data,"z") db.data <- db.locate(db.data,"m.length","z") db.data <- db.locate(db.data,"depth","f") db.data # Omni-directional variogram & Model vario.data <- vario.calc(db.data) model.vario <- model.auto(vario=vario.data,struc=c("Nugget Effect","Spherical"), wmode=2,npairpt=F,npairdw=T,inches=0.08, title="Variogram Model") demo.pause("Variogram Model") # Define grid grid.kriging <- db.locerase(grid.kriging,"z") grid.kriging <- db.locate(grid.kriging,"depth","f") grid.kriging # View drift plot(grid.kriging,name.image="depth",title="Depth used as Drift",asp=1) plot(db.data,col=1,pch=19,add=T) plot(poly.data,col=4,add=T) demo.pause("Display Depth Drift") # Perform kriging with external drift grid.kriging <- kriging(dbin=db.data,dbout=grid.kriging,model=model.vario, neigh=neigh.kriging, uc=c("1","f1"),mean=NA,calcul="point",radix="KED") grid.kriging # View kriging with external drift plot(grid.kriging,name.image="z1", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Mean Length - Estimation (KED)", pos.legend=5,flag.proj=F) plot(db.data,name.prop="m.length",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Depth Estimation (KED)") plot(grid.kriging,name.image="z2", asp=1/cos(mean(db.extract(db=db.data,names="x2"))*pi/180), title="Mean Length - St. Dev. (CCK)", pos.legend=5,flag.proj=F) plot(db.data,name.post="m.length",col=1,pch=20,add=T,flag.proj=F) map("worldHires",add=T,fill=T,col=8) demo.pause("Mean Length St. Dev. (KED)") }