Gaussian simulations
gstlearn Team
February 2023
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)## NULLWe 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## NULLUse 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## NULLUse 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## NULLUse 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## NULLPerform 10 conditional simulations with 1000 bands
nbsimu = 10
nbtuba = 1000
seed = 13231Simulations 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.81457err = 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## NULLp1 = 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)
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()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())## NULLdlfile = "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())## NULLvario_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## NULLmodel_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## NULLggplot() + plot.varmod(vario_resED, model_ED)
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)
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()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