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
In [3]:
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim);
In [4]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Temperatures.NF
In [5]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Elevations.NF
In [6]:
fileNF = "Scotland_Temperatures.NF"
dat = gl.Db.createFromNF(fileNF)
fileNF = "Scotland_Elevations.NF"
target = gl.DbGrid.createFromNF(fileNF)
In [7]:
neighU = gl.NeighUnique.create()
We calculate the experimental directional variogram and fit the Model
In [8]:
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 [9]:
err = gl.simtub(None, target, model=fitmod, nbtuba=1, seed=12454, namconv=gl.NamingConvention("Simu1"))
Graphic Representation
In [10]:
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 [11]:
err = gl.simtub(None, target, model=fitmod, nbtuba=10, seed=12454, namconv=gl.NamingConvention("Simu10"))
Graphic Representation
In [12]:
ax = target.plot("Simu10")
ax.decoration(title="Simulation with 10 Turning Bands")
Use the Turning Bands method with 1000 bands
In [13]:
err = gl.simtub(None, target, model=fitmod, nbtuba=1000, seed=12454, namconv=gl.NamingConvention("Simu1000"))
Graphic Representation
In [14]:
ax = target.plot("Simu1000")
ax.decoration(title="Simulation with 1000 Turning Bands")
Conditional Simulations¶
Perform 10 conditional simulations with 1000 bands
In [15]:
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 [16]:
err = gl.simtub(dat, target, model=fitmod, neigh=neighU, nbsimu=nbsimu, nbtuba=nbtuba, seed=seed)
Displaying several simulation outcomes
In [17]:
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 [18]:
err = target.statistics(["Simu.January_temp*"], [gl.EStatOption.MEAN], flagStoreInDb=True)
Comparing Mean of Simulations with the Kriging
In [19]:
err = gl.kriging(dat, target, model=fitmod, neigh = neighU, namconv=gl.NamingConvention("KS"))
In [20]:
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 [21]:
fileNF = "Scotland_Temperatures.NF"
dat = gl.Db.createFromNF(fileNF)
dat.setLocator("Elevation", gl.ELoc.F)
fileNF = "Scotland_Elevations.NF"
target = gl.DbGrid.createFromNF(fileNF)
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 [22]:
ax = gp.varmod(vario_resED, model_ED)
Performing Simulations with External Drift
In [23]:
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 [24]:
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 [25]:
err = target.statistics(["Simu.January_temp*"], [gl.EStatOption.MEAN], flagStoreInDb=True)
Comparing Mean of Simulations with Kriging with External Drift
In [26]:
err = gl.kriging(dat, target, model=model_ED, neigh = neighU, namconv=gl.NamingConvention("KED"))
In [27]:
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 [28]:
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 [29]:
err = target.statistics(["Simu.January_temp*"], [gl.EStatOption.MEAN], flagStoreInDb=True,
namconv=gl.NamingConvention("Proba"))
In [30]:
ax = target.plot("Proba.MEAN", flagLegendRaster=True)
ax.decoration(title="Probability for positive Temperatures")
