Preamble

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

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

Then we download the data base dat.

fileNF = loadData("Scotland", "Scotland_Temperatures.NF")
dat = Db_createFromNF(fileNF)

We also create a Db object containing a grid covering the data points in the data base dat. To do so, we start by displaying the minimal and maximal coordinates of the points in dat using the getExtremas() method from the Db class.

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

The first (resp. second) element of the list contains the min and max coordinates in the first (resp. second) space dimension. Based on this information, we create a grid covering all the points using the DbGrid_create function. We specify the coordinates of the origin (i.e. lower left corner) of the grid (argument x0), the step sizes in each dimension (argument dx) and the number of points in each dimension (argument nx).

grid = DbGrid_create(x0=c(65,530),dx=c(4.94, 4.96),nx=c(82,138)) 

We then print a summary of the content of the grid using the display method of Db class, which we supply with a DbStringFormat object specifying that we would like information about the extent of the grid (argument flag_extend in the DbStringFormat_createFromFlags function).

dbfmt = DbStringFormat_createFromFlags(flag_extend=TRUE) 
grid$display(dbfmt)
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 3
## Total number of samples      = 11316
## 
## Grid characteristics:
## ---------------------
## Origin :     65.000   530.000
## Mesh   :      4.940     4.960
## Number :         82       138
## 
## Data Base Extension
## -------------------
## Coor #1 - Min =     65.000 - Max =    465.140 - Ext = 400.14
## Coor #2 - Min =    530.000 - Max =   1209.520 - Ext = 679.52
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## NULL

Experimental Variogram and fitted Model

We compute the experimental variogram vario2dir (in 2 directions) (cf. Variography for more details).

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

We then the fit a model fitmod

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

Simple kriging prediction

To perform simple kriging, we use the function called kriging. We specify:

Additionally, it is possible to specify whether we wish to compute, at each target point, the kriging predictor (argument flag_est, default=TRUE), the kriging standard-deviation (argument flag_std, default=TRUE) and the kriging variance (argument flag_varz, default=FALSE).

The kriging function then adds new variables to the Db entered in the dbout argument corresponding to these variables. The names of these newly created variables will start by Kriging, but this prefix can be changed using the namconv argument of the kriging function

In the next example, we perform a simple kriging prediction (with unique neighborhood) of the variable January_temp in the dat data base, on the grid defined in the grid data base. To do so, we start by selecting the variable January_temp in the dat data base (i.e. we make ensure that it is the only variable with a z locator).

dat$setLocator("January_temp",ELoc_Z())
## NULL
dat
## 
## Data Base Characteristics
## =========================
## 
## Data Base Summary
## -----------------
## File is organized as a set of isolated points
## Space dimension              = 2
## Number of Columns            = 5
## Total number of samples      = 236
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = Longitude - Locator = x1
## Column = 2 - Name = Latitude - Locator = x2
## Column = 3 - Name = Elevation - Locator = NA
## Column = 4 - Name = January_temp - Locator = z1

We then create a "neighborhood" object corresponding to the specification of a unique neighborhood: this is done using the NeighUnique function.

uniqueNeigh = NeighUnique()

We now call the kriging function to perform the kriging prediction. We use the model fitmod that we previously fitted on our data, require to compute the kriging predictor and its standard-deviation (but not its variance), and change the prefix of the newly created variables to "SK".

err = kriging(dbin=dat, dbout=grid, model=fitmod, 
              neigh=uniqueNeigh,
              flag_est=TRUE, flag_std=TRUE, flag_varz=FALSE,
              namconv=NamingConvention("SK")
              )

We see that the kriging predictor and its standard deviation have been added to the grid data base.

grid
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 5
## Total number of samples      = 11316
## 
## Grid characteristics:
## ---------------------
## Origin :     65.000   530.000
## Mesh   :      4.940     4.960
## Number :         82       138
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = SK.January_temp.estim - Locator = z1
## Column = 4 - Name = SK.January_temp.stdev - Locator = NA

Finally, we plot the kriging prediction over the grid using the plot.grid function and the data points.

p = ggDefaultGeographic()
p = p + plot.grid(grid,
                  flagLegendRaster = TRUE, palette="Spectral",legendNameRaster="°C") 
p = p + plot.point(dat,flagCst = T,pch=18,cex=1.5)
p = p + plot.decoration(title="Simple Kriging over whole Grid")
ggPrint(p)