Gaussian simulations¶

In [1]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}
In [2]:
import gstlearn as gl
import gstlearn.plot as gp
import matplotlib.pyplot as plt
import numpy as np
import os
import urllib.request
In [3]:
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim);
In [4]:
url = 'https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Temperatures.NF'
temp_nf, head = urllib.request.urlretrieve(url)
dat = gl.Db.createFromNF(temp_nf)
url = 'https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Elevations.NF'
elev_nf, head = urllib.request.urlretrieve(url)
target = gl.DbGrid.createFromNF(elev_nf)
In [5]:
neighU = gl.NeighUnique.create()

We calculate the experimental directional variogram and fit the Model

In [6]:
varioparam = gl.VarioParam.createOmniDirection(npas=40, dpas=10)
vario_raw2dir = gl.Vario.create(varioparam, dat)
vario_raw2dir.compute()

fitmod = gl.Model()
err = fitmod.fit(vario_raw2dir, types=[gl.ECov.NUGGET, gl.ECov.SPHERICAL, gl.ECov.CUBIC])

Use the Turning Bands method with 1 band

In [7]:
err = gl.simtub(None, target, model=fitmod, nbtuba=1, seed=12454, namconv=gl.NamingConvention("Simu1"))

Graphic Representation

In [8]:
gp.setDefaultGeographic(dims=[8,8])
ax = target.plot("Simu1")
ax.decoration(title="Simulation with 1 Turning Band")

Use the Turning Bands method with 10 bands

In [9]:
err = gl.simtub(None, target, model=fitmod, nbtuba=10, seed=12454, namconv=gl.NamingConvention("Simu10"))

Graphic Representation

In [10]:
ax = target.plot("Simu10")
ax.decoration(title="Simulation with 10 Turning Bands")

Use the Turning Bands method with 1000 bands

In [11]:
err = gl.simtub(None, target, model=fitmod, nbtuba=1000, seed=12454, namconv=gl.NamingConvention("Simu1000"))

Graphic Representation

In [12]:
ax = target.plot("Simu1000")
ax.decoration(title="Simulation with 1000 Turning Bands")

Conditional Simulations¶

Perform 10 conditional simulations with 1000 bands

In [13]:
nbsimu = 10
nbtuba = 1000
seed   = 13231

mean_Temperature = gl.dbStatisticsMono(dat, ["J*temp"], [gl.EStatOption.MEAN])
fitmod.setMeans(mean_Temperature)

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

In [14]:
err = gl.simtub(dat, target, model=fitmod, neigh=neighU, nbsimu=nbsimu, nbtuba=nbtuba, seed=seed)

Displaying several simulation outcomes

In [15]:
fig, ax = plt.subplots(1,3, figsize=[20,8])
ax[0].gstgrid(target,"Simu*temp.1")
ax[1].gstgrid(target,"Simu*temp.2")
ax[2].gstgrid(target,"Simu*temp.3")
fig.subplots_adjust(right=0.7)
cbar_ax = fig.add_axes([0.75, 0.1, 0.02, 0.75])
im = ax[0].collections[0]
err = fig.colorbar(im, cax = cbar_ax)

Computing the Mean of the simulations and store it as a new variable

In [16]:
err = target.statistics(["Simu.January_temp*"], [gl.EStatOption.MEAN], flagStoreInDb=True)

Comparing Mean of Simulations with the Kriging

In [17]:
err = gl.kriging(dat, target, model=fitmod, neigh = neighU, namconv=gl.NamingConvention("KS"))
In [18]:
ax = gp.correlation(target, namex="Stats.MEAN", namey="KS*estim", diagLine=True, bins=100)
ax.decoration(xlabel="Mean of Simulations", ylabel="Simple Kriging Estimate")

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
In [19]:
dat = gl.Db.createFromNF(temp_nf)
dat.setLocator("Elevation", gl.ELoc.F)

target = gl.DbGrid.createFromNF(elev_nf)
target.setLocator("Elevation", gl.ELoc.F)

vario_resED = gl.Vario.create(varioparam, dat)
model = gl.Model()
model.setDriftIRF(nfex = 1)
vario_resED.compute(model=model)

model_ED = gl.Model()
err = model_ED.fit(vario_resED, types=[gl.ECov.SPHERICAL, gl.ECov.CUBIC])
model_ED.setDriftIRF(nfex=1)

Graphic representation of the experimental Variogram and fitted Model

In [20]:
ax = gp.varmod(vario_resED, model_ED)

Performing Simulations with External Drift

In [21]:
err = gl.simtub(dat, target, model=model_ED, neigh=neighU, nbsimu=nbsimu, nbtuba=nbtuba, seed=seed)

Display several simulations with Elevation as External Drift

In [22]:
fig, ax = plt.subplots(1,3,figsize=(20,8))
ax[0].gstgrid(target,"Simu*temp.1")
ax[1].gstgrid(target,"Simu*temp.2")
ax[2].gstgrid(target,"Simu*temp.3")
fig.subplots_adjust(right=0.7)
cbar_ax = fig.add_axes([0.75, 0.1, 0.02, 0.75])
im = ax[0].collections[0]
err = fig.colorbar(im, cax = cbar_ax)

Calculating the Mean of the Conditional Simulations

In [23]:
err = target.statistics(["Simu.January_temp*"], [gl.EStatOption.MEAN], flagStoreInDb=True)

Comparing Mean of Simulations with Kriging with External Drift

In [24]:
err = gl.kriging(dat, target, model=model_ED, neigh = neighU, namconv=gl.NamingConvention("KED"))
In [25]:
ax = gp.correlation(target, namex="Stats.MEAN", namey="KED*estim", diagLine=True, bins=100)
ax.decoration(xlabel="Mean of Simulations", ylabel="Kriging with External Drift")

Deriving Indicator to exceed a Cutoff per simulation outcome

In [26]:
target["Simu.January_temp*"] = target["Simu.January_temp*"] > 0

Computing the average of the Indicators which leads to the Probability to exceed the threshold

In [27]:
err = target.statistics(["Simu.January_temp*"], [gl.EStatOption.MEAN], flagStoreInDb=True,
                       namconv=gl.NamingConvention("Proba"))
In [28]:
ax = target.plot("Proba.MEAN", flagLegendRaster=True)
ax.decoration(title="Probability for positive Temperatures")