Pluri-Gaussian¶
In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import os
gdoc.setNoScroll()
Prepare the basic gstlearn objects¶
Initial objects are located in a specific Directory which is defined in the next Container operation. Note that his operation will modify automatically all the names of the Files retreived using Serialize / Deserialize operation (not when reading using CSV). Also note that the Container name must be ended using a "/" (as its serves as a Directory).
In [2]:
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN,ndim)
Load the data file¶
This Load is performed starting from a CSV file.
In [3]:
filename = gdoc.loadData("BRGM", "Nitrates_LANU.csv")
datCat = pd.read_csv(filename,sep=";")
datCat.head()
Out[3]:
| X | Y | LANU | HyGeo | NO3 | BG | |
|---|---|---|---|---|---|---|
| 0 | 428575.6629 | 5.541444e+06 | 2 | 1 | 0.5 | -99.0 |
| 1 | 431573.6549 | 5.542124e+06 | 2 | 1 | 150.0 | -99.0 |
| 2 | 431565.6549 | 5.542160e+06 | 2 | 1 | 150.0 | -99.0 |
| 3 | 435820.6447 | 5.543920e+06 | 3 | 1 | 78.0 | -99.0 |
| 4 | 435955.6443 | 5.543910e+06 | 3 | 1 | 64.0 | -99.0 |
Loading polygon from a file¶
The polygon is created by deserializing the Neutral Polygon File
In [4]:
filename = gdoc.loadData("BRGM", "poly_LANU.ascii")
poly = gl.Polygons.createFromNF(filename)
poly
Out[4]:
Polygons -------- Number of Polygon Sets = 1
Creation of the gstlearn data base¶
In [5]:
dat = gl.Db()
fields = ["X","Y","LANU"]
dat[fields] = datCat[fields].values
Specification of the role of each variable (named "locators" in gstlearn)¶
In [6]:
dat.setLocators(["X","Y"],gl.ELoc.X) #Coordinates
dat.setLocator("LANU",gl.ELoc.Z) #Variable of interest
dat
Out[6]:
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 = 1691 Variables --------- Column = 0 - Name = X - Locator = x1 Column = 1 - Name = Y - Locator = x2 Column = 2 - Name = LANU - Locator = z1
Creation of the output grid¶
The output grid will contain 47 x 101 nodes. It is built to cover the data file plus an extension of 10000 x 10000.
In [7]:
Result = gl.DbGrid.createCoveringDb(dat,[47,101],[],[],[50000,50000])
Add a selection (mask the cells outside the polygon)¶
The initial selection (based on the application of the Polygon to the grid data base) must be dilated in order to avoid edge effect.
In [8]:
gl.db_polygon(Result,poly)
Result.setLocator("Polygon",gl.ELoc.Z)
Result.morpho(gl.EMorpho.DILATION,0.5,1.5,option=0,verbose=False,radius=[1,1])
Result["Polygon"] = Result["Morpho.Polygon.*"]
Result.deleteColumn("Morpho.Polygon.*")
Result.setLocator("Polygon",gl.ELoc.SEL)
In [9]:
ax = Result.plot(nameRaster="Polygon",useSel=False,flagLegendRaster=False)
ax.gstpoint(dat,nameColor="LANU",size=2)
ax.polygon(poly,linewidth=1,edgecolor="r")
ax.geometry(dims=[10,10])
ax.decoration(title="Initial information")
Computation of the proportions¶
Compute global proportions (for information)¶
In [10]:
propGlob = gl.dbStatisticsFacies(dat)
ncat = len(propGlob)
for i in range(ncat):
print("Proportion of facies "+str(i+1),"=",propGlob[i])
Proportion of facies 1 = 0.2726197516262567 Proportion of facies 2 = 0.46954464813719693 Proportion of facies 3 = 0.046126552335895916 Proportion of facies 4 = 0.205204021289178 Proportion of facies 5 = 0.006505026611472502
Compute local proportions¶
The next parts will be simplified in a future dedicated API
2.2.1 Creation of the spatial regularization model for proportions
In [11]:
model = gl.Model.createFromDb(Result)
cova = gl.CovAniso.createIsotropic(model.getContext(),gl.ECov.MATERN,range=50000.,param=2.,sill=1.,)
model.addCov(cova)
In [12]:
err = gl.db_proportion_estimate(dat,Result,model)
In [13]:
dbfmt = gl.DbStringFormat()
dbfmt.setFlags(flag_stats=True)
dbfmt.setNames(["Prop.*"])
Result.display(dbfmt)
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 8 Total number of samples = 4747 Number of active samples = 752 Grid characteristics: --------------------- Origin : 306764.2835218227.694 Mesh : 4787.785 3803.924 Number : 47 101 Data Base Statistics -------------------- 4 - Name Prop.1 - Locator p1 Nb of data = 4747 Nb of active values = 752 Minimum value = 0.000 Maximum value = 0.867 Mean value = 0.218 Standard Deviation = 0.162 Variance = 0.026 5 - Name Prop.2 - Locator p2 Nb of data = 4747 Nb of active values = 752 Minimum value = 0.000 Maximum value = 0.934 Mean value = 0.415 Standard Deviation = 0.257 Variance = 0.066 6 - Name Prop.3 - Locator p3 Nb of data = 4747 Nb of active values = 752 Minimum value = 0.000 Maximum value = 0.368 Mean value = 0.042 Standard Deviation = 0.046 Variance = 0.002 7 - Name Prop.4 - Locator p4 Nb of data = 4747 Nb of active values = 752 Minimum value = 0.000 Maximum value = 0.967 Mean value = 0.161 Standard Deviation = 0.172 Variance = 0.030 8 - Name Prop.5 - Locator p5 Nb of data = 4747 Nb of active values = 752 Minimum value = 0.000 Maximum value = 0.117 Mean value = 0.005 Standard Deviation = 0.011 Variance = 0.000 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Polygon - Locator = sel Column = 3 - Name = Prop.1 - Locator = p1 Column = 4 - Name = Prop.2 - Locator = p2 Column = 5 - Name = Prop.3 - Locator = p3 Column = 6 - Name = Prop.4 - Locator = p4 Column = 7 - Name = Prop.5 - Locator = p5
Display the results¶
In [14]:
for i in range(ncat):
fig, ax = plt.subplots()
ax = Result.plot(nameRaster="Prop."+str(i+1))
ax.decoration(title="Proportion Facies #"+str(i+1))
ax.geometry(dims=[10,10], aspect=1)
ax.gstpoint(dat,nameColor="LANU",size=2,color="black")
dat.addSelectionByLimit("LANU",gl.Limits((i+1,i+1)),"SelPoint")
ax.gstpoint(dat,nameColor="LANU",size=0.8,color="red")
dat.deleteColumn("SelPoint")
ax.polygon(poly,linewidth=1,edgecolor="r")
