In this preamble, we load the gstlearn library, clean the workspace.

rm(list=ls())
library(gstlearn)
library(ggplot2)
library(ggpubr)

Then we download the data base dat.

dlfile = "https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Temperatures.NF"
fileNF = "Scotland_Temperatures.NF"
download.file(dlfile, fileNF)
dat = Db_createFromNF(fileNF)

Calculate the experimental variogram vario2dir (in 2 directions)

varioParamMulti = VarioParam_createMultiple(ndir=2, npas=15, dpas=15.)
vario2dir = Vario(varioParamMulti, dat)
err = vario2dir$compute()

Calculate the fitted model fitmodOK (add the Universality Condition)

fitmodOK = Model()
types = ECov_fromKeys(c("NUGGET","EXPONENTIAL","GAUSSIAN"))
err = fitmodOK$fit(vario2dir,types=types)
err = fitmodOK$addDrift(Drift1())
ggplot() + plot.varmod(vario2dir, fitmodOK)

Define the Unique Neighborhood unique.neigh:

unique.neigh = NeighUnique()

Get the extension of the Data:

dat$getExtremas()
## [[1]]
## [1]  78.2 460.7
## 
## [[2]]
## [1]  530.4 1208.9

Create the Target file grid:

grid = DbGrid_create(x0=c(65,535),dx=c(4.94, 4.96),nx=c(81,137))
dbfmt = DbStringFormat_createFromFlags(flag_resume=FALSE, flag_vars=TRUE,
                                       flag_extend=TRUE, flag_stats=FALSE,
                                       flag_array=FALSE)
grid$display(dbfmt)
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Extension
## -------------------
## Coor #1 - Min =     65.000 - Max =    460.200 - Ext = 395.2
## Coor #2 - Min =    535.000 - Max =   1209.560 - Ext = 674.56
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## NULL
err = kriging(dbin=dat, dbout=grid, model=fitmodOK, neigh=unique.neigh,
              flag_est=TRUE, flag_std=TRUE, flag_varz=FALSE,
              namconv=NamingConvention("OK"))
p = ggDefaultGeographic()
p = p + plot.grid(grid, legend.name.raster="°C", zlim=c(0,8))
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.decoration(title="Ordinary Kriging over whole Grid")
ggPrint(p)

dlfile = "https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Elevations.NF"
fileNF = "Scotland_Elevations.NF"
download.file(dlfile, fileNF)
grid = DbGrid_createFromNF(fileNF)
grid$display(dbfmt)
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Extension
## -------------------
## Coor #1 - Min =     65.000 - Max =    455.123 - Ext = 390.123
## Coor #2 - Min =    535.000 - Max =   1200.109 - Ext = 665.109
## 
## Variables
## ---------
## Column = 0 - Name = Longitude - Locator = x1
## Column = 1 - Name = Latitude - Locator = x2
## Column = 2 - Name = Elevation - Locator = f1
## Column = 3 - Name = inshore - Locator = sel
## NULL

The output grid now contains the selection inshore. Estimation is restricted to the active cells only.

err = kriging(dbin=dat, dbout=grid, model=fitmodOK, neigh=unique.neigh,
              flag_est=TRUE, flag_std=TRUE, flag_varz=FALSE,
              namconv=NamingConvention("OK"))
p = ggDefaultGeographic()
p = p + plot.grid(grid, "OK*estim", legend.name.raster="°C", zlim=c(0.,8.))
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.point(dat,name_size="January_temp",sizmax=3,color='black')
p = p + plot.decoration(title="Estimation by Ordinary Kriging")
ggPrint(p)

p = ggDefaultGeographic()
p = p + plot.grid(grid,"OK*stdev", legend.name.raster="°C", zlim=c(0.,1.))
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.point(dat,name_size="January_temp",sizmax=3,color='black')
p = p + plot.decoration(title="St. dev. by Ordinary Kriging")
ggPrint(p)

The Model fitmodOK is first duplicated into fitmodSK. Then the Universality Condition is deleted.

fitmodSK = fitmodOK$clone()
err = fitmodSK$delDrift(rank=0)
err = fitmodSK$setMean(mean=20.)

Simple Kriging is performed

err = kriging(dbin=dat, dbout=grid, model=fitmodSK, neigh=unique.neigh,
              flag_est=TRUE, flag_std=TRUE, 
              namconv=NamingConvention("SK"))
p = ggDefaultGeographic()
p = p + plot.grid(grid,"SK*estim",legend.name.raster="°C", zlim=c(0.,8.))
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.point(dat,name_size="January_temp",sizmax=3,color='black')
p = p + plot.decoration(title="Estimation by Simple Kriging")
ggPrint(p)

p = ggDefaultGeographic()
p = p + plot.grid(grid,"SK*stdev",legend.name.raster="°C", zlim=c(0.,1.))
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.point(dat,name_size="January_temp",sizmax=3,color='black')
p = p + plot.decoration(title="St. dev. by Simple Kriging")
ggPrint(p)

p = ggplot()
p = p + plot.correlation(grid,name1="OK*estim",name2="SK*estim", 
                     flagRegr=TRUE, flagSameAxes=TRUE, bins=100)
p = p + plot.decoration(title="Estimation Simple vs. Ordinary", 
                    xlab="Ordinary Kriging", ylab="Simple Kriging")
ggPrint(p)

p = ggplot()
p = p + plot.correlation(grid,name1 = "OK*stdev",name2="SK*stdev", 
                     flagRegr=TRUE, flagSameAxes=TRUE, bins=100)
p = p + plot.decoration(title="St. dev. Simple vs. Ordinary", 
                    xlab="Ordinary Kriging", ylab="Simple Kriging")
ggPrint(p)

err = xvalid(db=dat, model=fitmodOK, neigh=unique.neigh, 
             flag_xvalid_est=1, flag_xvalid_std=1,  
             namconv=NamingConvention_create("Xvalid", flag_locator = FALSE))
p = ggplot()
p = p + plot.hist(dat,name="*esterr*",bins=30,fill="blue")
p = p + plot.decoration(xlab="Estimation Errors", title="Cross-Validation")
ggPrint(p)

p = ggplot()
p = p + plot.hist(dat,name="*stderr*",bins=30,fill="blue")
p = p + plot.decoration(xlab="Standardized Errors", title="Cross-Validation")
ggPrint(p)

mean(dat$getColumn("*esterr*"),na.rm=TRUE)
## [1] -0.004167872
mean(dat$getColumn("*esterr*")^2,na.rm=TRUE)
## [1] 0.2393726
mean(dat$getColumn("*stderr*")^2,na.rm=TRUE)
## [1] 0.9117794
p = ggDefaultGeographic()
p = p + plot.grid(grid,"inshore")
p = p + plot.point(dat,name_size="*esterr",sizmax=3)
p = p + plot.decoration(title="Cross-Validation scores")
ggPrint(p)

p = ggDefaultGeographic()
p = p + plot.grid(grid,"inshore")
p = p + plot.point(dat,name_size="*esterr",sizmax=3,flagAbsSize=TRUE)
p = p + plot.decoration(title="Cross-Validation scores (abs. value)")
ggPrint(p)

We design a small Moving Neighborhood small.neigh with only 1 sample per neighborhood.

small.neigh = NeighMoving_create(nmini=1, nmaxi=1, radius=1000000)

We perform Ordinary Kriging

err = kriging(dbin=dat, dbout=grid, model=fitmodOK, neigh=small.neigh,
              flag_est=TRUE, flag_std=TRUE, 
              namconv=NamingConvention("Small"))
p = ggDefaultGeographic()
p = p + plot.grid(grid,"Small*estim",zlim=c(0.,8.), legend.name.raster="°C")
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.point(dat,name_size="January_temp")
p = p + plot.decoration(title="Estimation by Ordinary Kriging (Small Neigh.)")
ggPrint(p)

Building a reasonable Moving Neighborhood, although with a limited extension (radius)

moving.neigh = NeighMoving_create(nmini=1, nmaxi=10, radius=20)

Running the Ordinary Kriging

err = kriging(dat,grid,fitmodOK,moving.neigh,
              flag_est=TRUE, flag_std=TRUE, 
              namconv=NamingConvention("Reduced"))
p = ggDefaultGeographic()
p = p + plot.grid(grid,"Reduced*estim", zlim=c(0.,8.), legend.name.raster="°C")
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.point(dat,name_size="January_temp")
p = p + plot.decoration(title="Estimation by Ordinary Kriging (Reduced Moving Neigh.)")
ggPrint(p)

Lots of target sites are not estimated as no sample is found within the neighborhood.

Building a reasonable Moving Neighborhood correctly tuned: 10 samples (maximum) selected in a radius of 150 around the target site.

moving.neigh = NeighMoving_create(nmini=1, nmaxi=10, radius=150)

Running the Ordinary Kriging

err = kriging(dat,grid,fitmodOK,moving.neigh,
              flag_est=TRUE, flag_std=TRUE, 
              namconv=NamingConvention("Moving"))
p = ggDefaultGeographic()
p = p + plot.grid(grid,"Moving*estim",zlim=c(0.,8.), legend.name.raster="°C")
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.point(dat,name_size="January_temp")
p = p + plot.decoration(title="Estimation by Ordinary Kriging (Moving Neigh.)")
ggPrint(p)

p = ggDefaultGeographic()
p = p + plot.grid(grid,"Moving*stdev", zlim=c(0.,1.), legend.name.raster="°C")
## Warning in geom_tile(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.point(dat,name_size="January_temp")
p = p + plot.decoration(title="St. dev. by Ordinary Kriging (Moving Neigh.)")
ggPrint(p)

p = ggplot()
p = p + plot.correlation(grid,name1 = "OK*estim",name2="Moving*estim", 
                         bins=100, flagDiag=TRUE)
p = p + plot.decoration(title="Ordinary Kriging Estimation", 
                    xlab="Unique", ylab="Moving")
ggPrint(p)

p = ggplot()
p = p + plot.correlation(grid,name1 = "OK*stdev",name2="Moving*stdev", 
                         bins=100, flagDiag=TRUE)
p = p + plot.decoration(title="Ordinary Kriging St. Dev.", 
                    xlab="Unique", ylab="Moving")
ggPrint(p)