Gaussian simulations

\frametitle{Preamble}

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

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

neighU = NeighUnique_create()

ndim = 2
defineDefaultSpace(ESpaceType_RN(), ndim)

## NULL

We calculate the experimental directional variogram and fit the Model

varioparam = VarioParam_createOmniDirection(npas=40, dpas=10)
vario_raw2dir = Vario_create(varioparam, dat)
err = vario_raw2dir$compute()

fitmod = Model()
err = fitmod$fit(vario_raw2dir, 
                 types=ECov_fromKeys(c("NUGGET", "SPHERICAL", "CUBIC")))
fitmod$display()

## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 1
## Number of basic structure(s) = 1
## Number of drift function(s)  = 0
## Number of drift equation(s)  = 0
## 
## Covariance Part
## ---------------
## Spherical
## - Sill         =      1.155
## - Range        =    135.133
## Total Sill     =      1.155

## NULL

---

\frametitle{Unconditional simulation}

Use the Turning Bands method with 1 band

err = simtub(dbout=target, model=fitmod, nbtuba=1, seed=12454,
             namconv=NamingConvention("Simu1"))
p = ggDefaultGeographic()
p = p + plot(target, "Simu1")
p = p + plot.decoration(title="Simulation with 1 band")
ggPrint(p)

target$display()

## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 5
## Maximum Number of UIDs       = 5
## Total number of samples      = 11097
## Number of active samples     = 3092
## 
## Grid characteristics:
## ---------------------
## Origin :     65.000   535.000
## Mesh   :      4.938     4.963
## Number :         81       137
## 
## 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
## Column = 4 - Name = Simu1 - Locator = z1

## NULL

---

\frametitle{Unconditional simulation}

Use the Turning Bands method with 10 bands

err = simtub(dbout=target, model=fitmod, nbtuba=10, seed=12454,
             namconv=NamingConvention("Simu10"))
p = ggDefaultGeographic()
p = p + plot(target, "Simu10")
p = p + plot.decoration(title="Simulation with 10 Turning Bands")
ggPrint(p)

target$display()

## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 6
## Maximum Number of UIDs       = 6
## Total number of samples      = 11097
## Number of active samples     = 3092
## 
## Grid characteristics:
## ---------------------
## Origin :     65.000   535.000
## Mesh   :      4.938     4.963
## Number :         81       137
## 
## 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
## Column = 4 - Name = Simu1 - Locator = NA
## Column = 5 - Name = Simu10 - Locator = z1

## NULL

---

\frametitle{Unconditional simulation}

Use the Turning Bands method with 1000 bands

err = simtub(dbout=target, model=fitmod, nbtuba=1000, seed=12454,
             namconv=NamingConvention("Simu1000"))
p = ggDefaultGeographic()
p = p + plot(target,"Simu1000")
p = p + plot.decoration(title="Simulation with 1000 Turning Bands")
ggPrint(p)

target$display()

## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 7
## Maximum Number of UIDs       = 7
## Total number of samples      = 11097
## Number of active samples     = 3092
## 
## Grid characteristics:
## ---------------------
## Origin :     65.000   535.000
## Mesh   :      4.938     4.963
## Number :         81       137
## 
## 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
## Column = 4 - Name = Simu1 - Locator = NA
## Column = 5 - Name = Simu10 - Locator = NA
## Column = 6 - Name = Simu1000 - Locator = z1

## NULL

---

\frametitle{Conditional simulations}

Perform 10 conditional simulations with 1000 bands

nbsimu = 10
nbtuba = 1000
seed   = 13231

Simulations with Simple Kriging (Mean is extracted from the Data)

mean_Temperature = dbStatisticsMono(dat, c("J*temp"), 
                                    EStatOption_fromKeys(c("MEAN")))
print(mean_Temperature)

## [1] 2.81457

err = fitmod$setMeans(mean_Temperature)

err = simtub(dat, target, model=fitmod, neigh=neighU, nbsimu=nbsimu,
             nbtuba=nbtuba, seed=seed)
target$display()

## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 17
## Maximum Number of UIDs       = 17
## Total number of samples      = 11097
## Number of active samples     = 3092
## 
## Grid characteristics:
## ---------------------
## Origin :     65.000   535.000
## Mesh   :      4.938     4.963
## Number :         81       137
## 
## 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
## Column = 4 - Name = Simu1 - Locator = NA
## Column = 5 - Name = Simu10 - Locator = NA
## Column = 6 - Name = Simu1000 - Locator = NA
## Column = 7 - Name = Simu.January_temp.1 - Locator = z1
## Column = 8 - Name = Simu.January_temp.2 - Locator = z2
## Column = 9 - Name = Simu.January_temp.3 - Locator = z3
## Column = 10 - Name = Simu.January_temp.4 - Locator = z4
## Column = 11 - Name = Simu.January_temp.5 - Locator = z5
## Column = 12 - Name = Simu.January_temp.6 - Locator = z6
## Column = 13 - Name = Simu.January_temp.7 - Locator = z7
## Column = 14 - Name = Simu.January_temp.8 - Locator = z8
## Column = 15 - Name = Simu.January_temp.9 - Locator = z9
## Column = 16 - Name = Simu.January_temp.10 - Locator = z10

## NULL

---

\frametitle{Display several Simulation Outcomes}

p1 = ggDefaultGeographic()
p1 = p1 + plot(target, "Simu*temp.1")
p1 = p1 + plot(dat, flagCst=TRUE)
p2 = ggDefaultGeographic()
p2 = p2 + plot(target, "Simu*temp.2")
p2 = p2 + plot(dat, flagCst=TRUE)
p3 = ggDefaultGeographic()
p3 = p3 + plot(target, "Simu*temp.3")
p3 = p3 + plot(dat, flagCst=TRUE)
p4 = ggDefaultGeographic()
p4 = p4 + plot(target, "Simu*temp.4")
p4 = p4 + plot(dat, flagCst=TRUE)
ggarrange(p1,p2,p3,p4,nrow=2,ncol=2)

---

\frametitle{Comparing Mean of Simulations with the Kriging}

err = target$statistics(c("Simu.January_temp*"), 
                        EStatOption_fromKeys(c("MEAN")),
                        flagStoreInDb=TRUE)

err = kriging(dat, target, model=fitmod, neigh = neighU,
              namconv=NamingConvention("KS"))

p = ggplot()
p = p + plot.correlation(target, "Stats.MEAN", "KS*estim", flagDiag=TRUE, bins=100)
p = p + plot.decoration(xlab="Mean of Simulations", ylab="Simple Kriging Estimate")
ggPrint(p)

target$display()

---

\frametitle{Simulations with External Drift}

We perform the following steps:

• Load the grid containing the Elevation Map
• Prepare the Data Base: designation of External Drift
• Fit the Model of the Residuals

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

## NULL

dlfile = "https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Elevations.NF"
fileNF = "Scotland_Elevations.NF"
download.file(dlfile, fileNF)
target = DbGrid_createFromNF(fileNF)
target$setLocator("Elevation", ELoc_F())

## NULL

vario_resED = Vario_create(varioparam, dat)
model = Model()
err = model$setDriftIRF(nfex = 1)
err = vario_resED$compute(model=model)
vario_resED$display()

## 
## Variogram characteristics
## =========================
## Number of variable(s)       = 1
## Number of direction(s)      = 1
## Space dimension             = 2
## Variance-Covariance Matrix     0.363
## 
## Direction #1
## ------------
## Number of lags              = 40
## Direction coefficients      =      1.000     0.000
## Direction angles (degrees)  =      0.000     0.000
## Tolerance on direction      =     90.000 (degrees)
## Calculation lag             =     10.000
## Tolerance on distance       =     50.000 (Percent of the lag value)
## 
## For variable 1
##       Rank    Npairs  Distance     Value
##          0    19.000     3.118     0.070
##          1    89.000    10.690     0.122
##          2   168.000    20.346     0.126
##          3   244.000    30.323     0.183
##          4   316.000    40.429     0.195
##          5   385.000    50.162     0.217
##          6   399.000    60.296     0.223
##          7   463.000    70.062     0.277
##          8   450.000    79.807     0.249
##          9   473.000    90.115     0.282
##         10   549.000   100.141     0.231
##         11   484.000   109.866     0.307
##         12   491.000   119.865     0.293
##         13   458.000   130.001     0.349
##         14   474.000   139.888     0.335
##         15   472.000   149.918     0.334
##         16   390.000   159.896     0.415
##         17   410.000   169.913     0.466
##         18   385.000   179.840     0.442
##         19   425.000   189.849     0.409
##         20   354.000   200.117     0.487
##         21   341.000   209.881     0.502
##         22   325.000   219.710     0.566
##         23   296.000   229.902     0.566
##         24   277.000   239.927     0.602
##         25   236.000   250.111     0.566
##         26   201.000   259.921     0.510
##         27   177.000   270.161     0.586
##         28   166.000   279.560     0.560
##         29   177.000   289.873     0.503
##         30   148.000   299.866     0.461
##         31   116.000   310.335     0.504
##         32    95.000   319.854     0.508
##         33    92.000   329.857     0.618
##         34    64.000   339.936     0.397
##         35    57.000   350.172     0.517
##         36    61.000   359.473     0.334
##         37    51.000   369.274     0.386
##         38    33.000   380.087     0.389
##         39    31.000   389.553     0.380

## NULL

model_ED = Model()
err = model_ED$fit(vario_resED, 
                   types=ECov_fromKeys(c("SPHERICAL","CUBIC")))
err = model_ED$setDriftIRF(nfex=1)
model_ED$display()

## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 1
## Number of basic structure(s) = 2
## Number of drift function(s)  = 2
## Number of drift equation(s)  = 2
## 
## Covariance Part
## ---------------
## Spherical
## - Sill         =      0.178
## - Range        =     30.039
## Cubic
## - Sill         =      0.420
## - Range        =    451.351
## Total Sill     =      0.598
## 
## Drift Part
## ----------
## Universality Condition
## External Drift - Rank=0

## NULL

ggplot() + plot.varmod(vario_resED, model_ED)

---

\frametitle{Simulations with External Drift}

err = simtub(dat, target, model=model_ED, neigh=neighU, nbsimu=nbsimu,
             nbtuba=nbtuba, seed=seed)

p1 = ggDefaultGeographic() + plot(target, "Simu*temp.1")
p2 = ggDefaultGeographic() + plot(target, "Simu*temp.2")
p3 = ggDefaultGeographic() + plot(target, "Simu*temp.3")
p4 = ggDefaultGeographic() + plot(target, "Simu*temp.4")
ggarrange(p1,p2,p3,p4,ncol=2,nrow=2)

---

\frametitle{Comparing Mean of Simulations and Kriging with External Drift}

err = target$statistics(c("Simu.January_temp*"), 
                        EStatOption_fromKeys(c("MEAN")),
                        flagStoreInDb=TRUE)

err = kriging(dat, target, model=model_ED, neigh= neighU,
              namconv=NamingConvention("KED"))

p = ggplot()
p = p + plot.correlation(target, "Stats.MEAN", "KED*estim", flagDiag=TRUE,
                     bins=100)
p = p + plot.decoration(xlab="Mean of Simulations", 
                ylab="Kriging with External Drift")
ggPrint(p)

target$display()

\frametitle{Deriving probability to exceed a Cutoff}

target["Simu.January_temp*"] = target["Simu.January_temp*"] > 0

err = target$statistics(c("Simu.January_temp*"), 
                        EStatOption_fromKeys(c("MEAN")), 
                        flagStoreInDb=TRUE,
                        namconv=NamingConvention("Proba"))

p = ggDefaultGeographic()
p = p + plot(target, "Proba.MEAN", show.legend.raster=TRUE, legend.name.raster="Probability")
p = p + plot.decoration(title="Probability for positive Temperatures")
ggPrint(p)

target$display()

## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 18
## Maximum Number of UIDs       = 18
## Total number of samples      = 11097
## Number of active samples     = 3092
## 
## Grid characteristics:
## ---------------------
## Origin :     65.000   535.000
## Mesh   :      4.938     4.963
## Number :         81       137
## 
## 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
## Column = 4 - Name = Simu.January_temp.1 - Locator = NA
## Column = 5 - Name = Simu.January_temp.2 - Locator = NA
## Column = 6 - Name = Simu.January_temp.3 - Locator = NA
## Column = 7 - Name = Simu.January_temp.4 - Locator = NA
## Column = 8 - Name = Simu.January_temp.5 - Locator = NA
## Column = 9 - Name = Simu.January_temp.6 - Locator = NA
## Column = 10 - Name = Simu.January_temp.7 - Locator = NA
## Column = 11 - Name = Simu.January_temp.8 - Locator = NA
## Column = 12 - Name = Simu.January_temp.9 - Locator = NA
## Column = 13 - Name = Simu.January_temp.10 - Locator = NA
## Column = 14 - Name = Stats.MEAN - Locator = NA
## Column = 15 - Name = KED.January_temp.estim - Locator = NA
## Column = 16 - Name = KED.January_temp.stdev - Locator = NA
## Column = 17 - Name = Proba.MEAN - Locator = z1

## NULL