In [1]:
import matplotlib.pyplot as plt
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import numpy as np
import pandas as pd
from scipy import interpolate
gdoc.setNoScroll()
Creation of the data base for thickness¶
In [2]:
# Use either Oise_Thickness.csv, Oise_Thickness_withSides.csv or Oise_Ztop500.csv
# Data files are available here: https://soft.mines-paristech.fr/gstlearn/data-latest/Alluvial
filename = gdoc.loadData("Alluvial", "Oise_Ztop500.csv")
csv = gl.CSVformat.createStandard(gl.ECSV.FRENCH, naString="9999")
data = gl.Db.createFromCSV(filename, csv, False)
data
Out[2]:
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 1 Number of Columns = 5 Total number of samples = 2470 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = NA Column = 3 - Name = Z - Locator = NA Column = 4 - Name = Z_1 - Locator = z1
In [3]:
data.setLocator("X", gl.ELoc.X, 0)
data.setLocator("Y", gl.ELoc.X, 1)
data.setLocator("Z", gl.ELoc.L, 0)
data.setLocator("Z_1", gl.ELoc.U, 0)
sf = gl.DbStringFormat()
sf.setFlags(flag_stats=True)
thickness = data.getWithinBounds(0)
err = data.addColumns(thickness, "Elev", gl.ELoc.Z)
err = gl.DbHelper.db_duplicate(data)
data
Out[3]:
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 7 Total number of samples = 2470 Number of active samples = 2470 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Z - Locator = lower1 Column = 4 - Name = Z_1 - Locator = upper1 Column = 5 - Name = Elev - Locator = z1 Column = 6 - Name = Duplicate - Locator = sel
In [4]:
res = plt.scatter(data.getColumn("X"), data.getColumn("Y"), c=data.getColumn("Elev"))
In [5]:
res = plt.hist(thickness, bins=50)
In [6]:
# data.dumpToNF("Oise_Data_Elev500.ascii")
Creation of the grid¶
Creation of the polygon¶
In [7]:
# Data files are available here: https://soft.mines-paristech.fr/gstlearn/data-latest/Alluvial
filename = gdoc.loadData("Alluvial", "Oise_Shapefile_AlluvialPlain.csv")
poly = gl.Polygons.createFromCSV(filename, csv, False)
ax = gp.plot(poly)
In [8]:
grid = gl.DbGrid.create(
nx=[3300, 400], dx=[50.0, 50.0], x0=[630000.0, 6865000.0], angles=[40, 0]
)
In [9]:
fig, ax = gp.init()
ax.raster(grid, "x1", alpha=0.3)
ax.symbol(data)
gp.close()
In [10]:
err = gl.db_polygon(grid, poly)
grid
Out[10]:
Data Base Grid Characteristics
==============================
Data Base Summary
-----------------
File is organized as a regular grid
Space dimension = 2
Number of Columns = 4
Total number of samples = 1320000
Number of active samples = 138248
Grid characteristics:
---------------------
Origin : 630000.0006865000.000
Mesh : 50.000 50.000
Number : 3300 400
Rotation Angles = 40.000 0.000
Direct Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 -0.643
[ 1,] 0.643 0.766
Inverse Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 0.643
[ 1,] -0.643 0.766
Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = Polygon - Locator = sel
Creation of angles¶
In [11]:
# Data files are available here: https://soft.mines-paristech.fr/gstlearn/data-latest/Alluvial
filename = gdoc.loadData("Alluvial", "Oise_Shapefile_Centerline.csv")
## Loading Centerline data from file
df1 = pd.read_csv(filename, sep=";")
xc = list(df1["X"])
yc = list(df1["Y"])
x1 = xc[1:499]
y1 = yc[1:499]
# Data files are available here: https://soft.mines-paristech.fr/gstlearn/data-latest/Alluvial
filename = gdoc.loadData("Alluvial", "Oise_Shapefile_AlluvialPlain.csv")
## Loading alluvial plain contours from file & Separate into two polylines
df = pd.read_csv(filename, sep=";")
xp = list(df["X"])
yp = list(df["Y"])
x2 = xp[1 : xp.index(max(xp))]
x3 = xp[len(xp) : xp.index(max(xp)) + 1 : -1]
y2 = yp[1 : xp.index(max(xp))]
y3 = yp[len(xp) : xp.index(max(xp)) + 1 : -1]
## Adding supplementary control points at the edges
# coordinates of extremes
XA1 = 640000
YA1 = 6875000
XB1 = 740000
YB1 = 6955000
XA2 = 630000
YA2 = 6890000
XB2 = 735000
YB2 = 6980000
# Forming two supplementary vectors at the edges
n = 19
x4 = np.zeros(n)
y4 = np.zeros(n)
x5 = np.zeros(n)
y5 = np.zeros(n)
for i in range(0, len(x4)):
x4[i] = XA1 + (XB1 - XA1) / (n - 1) * i
y4[i] = YA1 + (YB1 - YA1) / (n - 1) * i
x5[i] = XA2 + (XB2 - XA2) / (n - 1) * i
y5[i] = YA2 + (YB2 - YA2) / (n - 1) * i
In [12]:
# Gradient method
# https://stackoverflow.com/questions/28269379/curve-curvature-in-numpy
dx_dt1 = np.gradient(x1)
dy_dt1 = np.gradient(y1)
velocity = np.array([[dx_dt1[i], dy_dt1[i]] for i in range(dx_dt1.size)])
ds_dt1 = np.sqrt(dx_dt1 * dx_dt1 + dy_dt1 * dy_dt1)
tangent1 = np.array([1 / ds_dt1] * 2).transpose() * velocity
dx_dt2 = np.gradient(x2)
dy_dt2 = np.gradient(y2)
velocity2 = np.array([[dx_dt2[i], dy_dt2[i]] for i in range(dx_dt2.size)])
ds_dt2 = np.sqrt(dx_dt2 * dx_dt2 + dy_dt2 * dy_dt2)
tangent2 = np.array([1 / ds_dt2] * 2).transpose() * velocity2
dx_dt3 = np.gradient(x3)
dy_dt3 = np.gradient(y3)
velocity3 = np.array([[dx_dt3[i], dy_dt3[i]] for i in range(dx_dt3.size)])
ds_dt3 = np.sqrt(dx_dt3 * dx_dt3 + dy_dt3 * dy_dt3)
tangent3 = np.array([1 / ds_dt3] * 2).transpose() * velocity3
dx_dt4 = np.gradient(x4)
dy_dt4 = np.gradient(y4)
velocity4 = np.array([[dx_dt4[i], dy_dt4[i]] for i in range(dx_dt4.size)])
ds_dt4 = np.sqrt(dx_dt4 * dx_dt4 + dy_dt4 * dy_dt4)
tangent4 = np.array([1 / ds_dt4] * 2).transpose() * velocity4
dx_dt5 = np.gradient(x5)
dy_dt5 = np.gradient(y5)
velocity5 = np.array([[dx_dt5[i], dy_dt5[i]] for i in range(dx_dt5.size)])
ds_dt5 = np.sqrt(dx_dt5 * dx_dt5 + dy_dt5 * dy_dt5)
tangent5 = np.array([1 / ds_dt5] * 2).transpose() * velocity5
tangent = np.concatenate((tangent1, tangent2, tangent3, tangent4, tangent5), axis=0)
x0 = np.concatenate((x1, x2, x3, x4, x5), axis=0)
y0 = np.concatenate((y1, y2, y3, y4, y5), axis=0)
u0 = tangent[:, 0]
v0 = tangent[:, 1]
In [13]:
# Numpy arctan2 Method
angles = np.arctan2(v0, u0) * 180 / np.pi
res = plt.hist(angles, bins=50)
In [14]:
# Interpolating the angles into angles map -> https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
xx = np.linspace(630000, 742000, 2240)
yy = np.linspace(6875000, 6982200, 2144)
xx, yy = np.meshgrid(xx, yy)
angles.shape = (np.size(angles),)
np.shape(angles)
points = np.transpose(np.vstack((x0, y0)))
angles_interp = interpolate.griddata(points, angles, (xx, yy), method="linear")
xx1 = np.concatenate(xx)
yy1 = np.concatenate(yy)
angles_interp1 = np.concatenate(angles_interp)
In [15]:
# Create a point Db point using xx, yy, angles_interp
dbv = gl.Db()
dbv.addColumns(xx1, "xx", gl.ELoc.X, 0)
dbv.addColumns(yy1, "yy", gl.ELoc.X, 1)
dbv.addColumns(angles_interp1, "angles_interp", gl.ELoc.Z, 0)
dbv.display()
grid.setLocator("Pol*")
grid.display()
Data Base Characteristics
=========================
Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension = 2
Number of Columns = 3
Total number of samples = 4802560
Variables
---------
Column = 0 - Name = xx - Locator = x1
Column = 1 - Name = yy - Locator = x2
Column = 2 - Name = angles_interp - Locator = z1
Data Base Grid Characteristics
==============================
Data Base Summary
-----------------
File is organized as a regular grid
Space dimension = 2
Number of Columns = 4
Total number of samples = 1320000
Grid characteristics:
---------------------
Origin : 630000.0006865000.000
Mesh : 50.000 50.000
Number : 3300 400
Rotation Angles = 40.000 0.000
Direct Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 -0.643
[ 1,] 0.643 0.766
Inverse Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 0.643
[ 1,] -0.643 0.766
Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = Polygon - Locator = NA
In [16]:
# Tester la valeur retournée par migrate :
grid.deleteColumn("Migr*")
gl.migrateMulti(dbv, grid, ["angles_interp"])
grid.setLocator("Migrate*", gl.ELoc.Z)
grid.display()
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 = 1320000
Grid characteristics:
---------------------
Origin : 630000.0006865000.000
Mesh : 50.000 50.000
Number : 3300 400
Rotation Angles = 40.000 0.000
Direct Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 -0.643
[ 1,] 0.643 0.766
Inverse Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 0.643
[ 1,] -0.643 0.766
Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = Polygon - Locator = NA
Column = 4 - Name = Migrate.angles_interp - Locator = z1
In [17]:
grid.useSel = False
grid.setLocator("Polygon")
uid_selection = grid.addSelectionByLimit(
"*angles_interp*", gl.Limits([-10], [10]), "angle_define"
)
In [18]:
xplt = grid.getColumnByLocator(gl.ELoc.X, 0, useSel=True)
yplt = grid.getColumnByLocator(gl.ELoc.X, 1, useSel=True)
aplt = grid.getColumnByLocator(gl.ELoc.Z, 0, useSel=True)
Creation of a selection (polygon + borders by morphological dilation)¶
In [19]:
grid.setLocator("Poly*", gl.ELoc.SEL)
In [20]:
vmin = 0.5
vmax = 1.5
nxy = grid.getNXs()
image2 = gl.BImage(nxy)
tab = grid.getColumn("Polygon", useSel=False)
image = gl.morpho_double2image(nxy, tab, vmin, vmax)
localVD = gl.VectorDouble(len(tab))
gl.morpho_dilation(0, [3, 3], image, image2)
for i in range(10):
gl.morpho_dilation(0, [1, 1], image, image2)
gl.morpho_dilation(0, [1, 1], image2, image)
# gl.morpho_dilation(0, [1,1], image, image2)
gl.morpho_image2double(image2, 0, 1.0, 0.0, localVD)
In [21]:
grid["res"] = localVD.getVector()
In [22]:
grid.setLocator("res", gl.ELoc.SEL)
grid
Out[22]:
Data Base Grid Characteristics
==============================
Data Base Summary
-----------------
File is organized as a regular grid
Space dimension = 2
Number of Columns = 7
Total number of samples = 1320000
Number of active samples = 278129
Grid characteristics:
---------------------
Origin : 630000.0006865000.000
Mesh : 50.000 50.000
Number : 3300 400
Rotation Angles = 40.000 0.000
Direct Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 -0.643
[ 1,] 0.643 0.766
Inverse Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 0.643
[ 1,] -0.643 0.766
Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = Polygon - Locator = NA
Column = 4 - Name = Migrate.angles_interp - Locator = z1
Column = 5 - Name = angle_define - Locator = NA
Column = 6 - Name = res - Locator = sel
In [23]:
res = gp.raster(grid, "Migrate.angles_interp", cmap="turbo", flagLegend=True)
Save the result¶
In [24]:
# grid.dumpToNF("Oise_GridAnglesModifFinal.ascii")
Start Kriging here once data is properly saved¶
In [25]:
# grid = gl.DbGrid.createFromNF("Oise_GridAnglesModifFinal.ascii")
grid
Out[25]:
Data Base Grid Characteristics
==============================
Data Base Summary
-----------------
File is organized as a regular grid
Space dimension = 2
Number of Columns = 7
Total number of samples = 1320000
Number of active samples = 278129
Grid characteristics:
---------------------
Origin : 630000.0006865000.000
Mesh : 50.000 50.000
Number : 3300 400
Rotation Angles = 40.000 0.000
Direct Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 -0.643
[ 1,] 0.643 0.766
Inverse Rotation Matrix
[, 0] [, 1]
[ 0,] 0.766 0.643
[ 1,] -0.643 0.766
Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = Polygon - Locator = NA
Column = 4 - Name = Migrate.angles_interp - Locator = z1
Column = 5 - Name = angle_define - Locator = NA
Column = 6 - Name = res - Locator = sel
In [26]:
# Modify the ranges to 8000, 800 if wished
ranges = [2000, 200]
# If old method is set to True, a classical SPDE kriging will be used
oldMethod = False
In [27]:
# Load data, use either with or without the sides; For the top of the alluvial plain use "Oise_Data_Elev500.ascii"
# Data files are available here: https://soft.mines-paristech.fr/gstlearn/data-latest/Alluvial
filename = gdoc.loadData("Alluvial", "Oise_Data_ThicknessSides.ascii")
data = gl.Db.createFromNF(filename)
data
Out[27]:
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 7 Total number of samples = 2334 Number of active samples = 2319 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Thickness - Locator = lower1 Column = 4 - Name = Thickness_1 - Locator = upper1 Column = 5 - Name = ThicknessSides - Locator = z1 Column = 6 - Name = Duplicate - Locator = sel
In [28]:
# Selection to calculate variogram only on non-zero values
data["selVario"] = data["Duplicate"] * (data["ThicknessSides"] > 0)
data.setLocators(["selVario"], gl.ELoc.SEL)
data
Out[28]:
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 8 Total number of samples = 2334 Number of active samples = 1010 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Thickness - Locator = lower1 Column = 4 - Name = Thickness_1 - Locator = upper1 Column = 5 - Name = ThicknessSides - Locator = z1 Column = 6 - Name = Duplicate - Locator = NA Column = 7 - Name = selVario - Locator = sel
In [29]:
myVarioParamBidir = gl.VarioParam()
mydir = gl.DirParam(40, 800.0, 0.5, 45.0, 0, 0, np.nan, np.nan, 0.0, [], [1, 1])
myVarioParamBidir.addDir(mydir)
mydir = gl.DirParam(20, 400.0, 0.5, 45.0, 0, 0, np.nan, np.nan, 0.0, [], [-1, 1])
myVarioParamBidir.addDir(mydir)
myVarioBidir = gl.Vario(myVarioParamBidir)
err = myVarioBidir.compute(data, gl.ECalcVario.VARIOGRAM)
myVarioBidir.display()
ax = gp.varmod(myVarioBidir, idir=0)
ax.decoration(title="Bi-directional Variogram for thickness(45°)")
ax = gp.varmod(myVarioBidir, idir=1)
ax.decoration(title="Bi-directional Variogram for thickness(45/135°)")
Variogram characteristics
=========================
Number of variable(s) = 1
Number of direction(s) = 2
Space dimension = 2
Variable(s) = [ThicknessSides]
Variance-Covariance Matrix 8.283
Direction #1
------------
Number of lags = 40
Direction coefficients = 1.000 1.000
Direction angles (degrees) = 45.000
Tolerance on direction = 45.000 (degrees)
Calculation lag = 800.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 1491.000 222.985 4.285
1 4135.000 817.513 6.027
2 5005.000 1597.763 5.766
3 5177.000 2403.961 7.196
4 4931.000 3202.456 7.062
5 5414.000 3996.471 6.883
6 5969.000 4804.571 6.796
7 5993.000 5607.391 6.235
8 5918.000 6392.940 6.156
9 5849.000 7198.209 7.202
10 5641.000 8004.082 6.712
11 4824.000 8800.796 7.726
12 5439.000 9591.894 6.860
13 5085.000 10419.647 6.706
14 5724.000 11192.490 6.888
15 6188.000 12005.854 6.486
16 5405.000 12774.212 6.783
17 4865.000 13605.264 6.492
18 4534.000 14393.326 6.467
19 4769.000 15203.346 6.514
20 4627.000 16001.095 5.897
21 4445.000 16792.930 6.947
22 4752.000 17616.012 7.331
23 6299.000 18422.320 7.206
24 5733.000 19188.948 7.033
25 5351.000 19996.702 6.631
26 5152.000 20791.695 7.815
27 5461.000 21593.342 7.902
28 4816.000 22408.492 7.541
29 5391.000 23190.246 9.345
30 5021.000 24009.337 7.375
31 5194.000 24797.226 7.986
32 5192.000 25607.977 7.768
33 5401.000 26409.342 8.086
34 4927.000 27190.648 7.410
35 4718.000 27990.954 7.806
36 4509.000 28793.658 7.400
37 5346.000 29648.534 6.941
38 5449.000 30397.797 7.461
39 6073.000 31199.396 7.758
Direction #2
------------
Number of lags = 20
Direction coefficients = -1.000 1.000
Direction angles (degrees) = 135.000
Tolerance on direction = 45.000 (degrees)
Calculation lag = 400.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 596.000 119.136 2.768
1 1400.000 393.701 4.532
2 1453.000 793.785 6.290
3 1224.000 1196.400 6.529
4 1032.000 1596.868 7.385
5 858.000 1983.161 7.192
6 697.000 2404.699 7.787
7 675.000 2794.542 7.751
8 419.000 3189.263 8.768
9 354.000 3597.904 11.252
10 339.000 3996.583 9.995
11 260.000 4402.766 7.896
12 232.000 4792.850 7.453
13 177.000 5178.782 7.800
14 172.000 5586.833 6.960
15 100.000 5978.383 9.960
16 89.000 6394.625 9.218
17 55.000 6800.319 5.859
18 30.000 7199.585 9.574
19 39.000 7641.196 4.828
In [30]:
model = gl.Model.createFromDb(data)
structs = [gl.ECov.NUGGET, gl.ECov.MATERN]
cons1P = gl.ConsItem.define(gl.EConsElem.PARAM, 1, type=gl.EConsType.EQUAL, value=1)
a = gl.Constraints()
a.addItem(cons1P)
err = model.fit(myVarioBidir, structs, constraints=a)
model.display()
ax = gp.varmod(myVarioBidir, model, idir=0)
ax = gp.varmod(myVarioBidir, model, idir=1)
plt.show()
# model.setMean(100)
model.display()
Model characteristics
=====================
Space dimension = 2
Number of variable(s) = 1
Number of basic structure(s) = 2
Number of drift function(s) = 0
Number of drift equation(s) = 0
Covariance Part
---------------
Nugget Effect
- Sill = 2.203
Matern (Third Parameter = 1)
- Sill = 4.715
- Ranges = 1045.034 1497.054
- Theo. Ranges = 301.675 432.162
- Angles = 45.000 0.000
- Rotation Matrix
[, 0] [, 1]
[ 0,] 0.707 -0.707
[ 1,] 0.707 0.707
Total Sill = 6.918
Known Mean(s) 0.000
Model characteristics
=====================
Space dimension = 2
Number of variable(s) = 1
Number of basic structure(s) = 2
Number of drift function(s) = 0
Number of drift equation(s) = 0
Covariance Part
---------------
Nugget Effect
- Sill = 2.203
Matern (Third Parameter = 1)
- Sill = 4.715
- Ranges = 1045.034 1497.054
- Theo. Ranges = 301.675 432.162
- Angles = 45.000 0.000
- Rotation Matrix
[, 0] [, 1]
[ 0,] 0.707 -0.707
[ 1,] 0.707 0.707
Total Sill = 6.918
Known Mean(s) 0.000
In [31]:
data.setLocators(["Duplicate"], gl.ELoc.SEL)
In [32]:
if not oldMethod:
model = gl.Model.createFromDb(data)
covMatern = gl.CovAniso.createAnisotropic(
type=gl.ECov.MATERN, ranges=ranges, sill=4.7, param=1, ctxt=model.getContext()
)
covNugget = gl.CovAniso.createIsotropic(
type=gl.ECov.NUGGET, sill=2.2, ctxt=model.getContext(), range=1
)
model.addCov(covNugget)
model.addCov(covMatern)
In [33]:
grid.setLocator("Polygon.*", gl.ELoc.SEL)
In [34]:
model.setDriftIRF(order=0)
model
Out[34]:
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 2 Number of drift function(s) = 1 Number of drift equation(s) = 1 Covariance Part --------------- Nugget Effect - Sill = 2.200 Matern (Third Parameter = 1) - Sill = 4.700 - Ranges = 2000.000 200.000 - Theo. Ranges = 577.350 57.735 Total Sill = 6.900 Drift Part ---------- Universality_Condition
In [35]:
if not oldMethod:
model.getCovAniso(1).makeAngleNoStatDb("Migrate.angles_interp", 0, grid)
In [36]:
data.setLocators(["ThicknessNosides"], gl.ELoc.Z, cleanSameLocator=True)
In [37]:
err = gl.krigingSPDE(data, grid, model)
In [38]:
# Plot kriging
fig, ax = gp.init(figsize=(20, 12))
gp.raster(grid, "KrigingSPDE.ThicknessSides.estim", cmap="turbo", flagLegend=True)
gp.close()
Save the current grid¶
The current grid can beloaded again at the start of Kriging
In [39]:
# modify the grid according to your wishes
# grid["Old_Variable"]=grid["New_variable"]
# grid.deleteColumn("TheVariableYouWant")
# grid["BedrockElevation"]=grid["PlainElevation"]-grid["AlluviumThickness"]
In [40]:
# Save the grid
# grid.dumpToNF("Oise_KrigingResults_Thickness.ascii")
Uncertainties estimation based on standard deviation¶
In [41]:
# Calculation of n simulations to estimate standard deviation (best is to change to nbsimu=100)
err = gl.simulateSPDE(data, grid, model, nbsimu=10)
In [42]:
# Calculation of Stdv by hand: gl.EStatOption.STDV in the future
grid["MeanSimu"] = grid["SimuSPDE.Thickness*"].mean(1)
grid["StdevSimu"] = np.zeros(grid["SimuSPDE.Thickness*"].shape[0])
nc = grid["SimuSPDE.Thickness*"].shape[1]
for i in range(nc):
grid["StdevSimu"] += (
(grid["SimuSPDE.Thickness*"][:, i] - grid["MeanSimu"]) ** 2
) / nc
grid["StdevSimu"] = grid["StdevSimu"] ** 0.5
In [43]:
# Plot kriging standard deviation (obtained through simulations)
fig, ax = gp.init(figsize=(20, 12))
gp.raster(grid, "StdevSimu", cmap="turbo", flagLegend=True)
gp.close()
In [ ]: