GEONE - DEESSE - Simulations with connectivity data

Main points addressed

• how to set hard data with connectivity data

• how to set connectivity constraints

Import what is required

import numpy as np
import matplotlib.pyplot as plt
import time
import os

# import package 'geone'
import geone as gn

# Show version of python and version of geone
import sys
print(sys.version_info)
print('geone version: ' + gn.__version__)

sys.version_info(major=3, minor=13, micro=7, releaselevel='final', serial=0)
geone version: 1.3.0

Training image (TI)

# Read file
data_dir = 'data'
filename = os.path.join(data_dir, 'ti.txt')
ti = gn.img.readImageTxt(filename)

# Values in the TI
ti.get_unique()

array([0., 1., 2.])

# Setting for categories / colors
categ_val = [0, 1, 2]
categ_col = ['lightblue', 'darkgreen', 'orange']

plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(ti, categ=True, categVal=categ_val, categCol=categ_col, title='TI')
plt.show()

Simulation grid

Define the simulation grid (number of cells in each direction, cell unit, origin).

nx, ny, nz = 100, 100, 1         # number of cells
sx, sy, sz = ti.sx, ti.sy, ti.sz # cell unit
ox, oy, oz = 0.0, 0.0, 0.0       # origin (corner of the "first" grid cell)

Hard data and connectivity data

A hard data point is defined by:

• its spatial coordinates (x, y, z) in the simulation grid,

• the value of the simulated variable, whose the name (code in the example below) identical to the name given by varname in the class geone.deesseinterface.DeesseInput (input for deesse),

• a connected component (geobody) label, whose the name (geobody_label in the example below) identical to the name given by varname in the class geone.deesseinterface.Connectivity describing the connectivity constraints.

Points with the same positive connected component label means that they should be connected in the output realizations. A negative or zero label means no connectivity constraints on that point.

npt = 5 # number of points
nv = 5   # number of variables including x, y, z coordinates
varname = ['x', 'y', 'z', 'code', 'geobody_label'] # list of variable names
v = np.array([
    [ 10.5,  40.5, 0.5,  1, 1], # x, y, z, code, geobody_label: 1st point
    [ 50.5,  35.5, 0.5,  1, 1], # 2nd point: should be connected to the 1st point
    [ 90.5,  48.5, 0.5,  1, 1], # 3rd point: should be connected to the 1st and 2nd points
    [ 75.5,  50.5, 0.5,  0, 0], # no connectivity constraint on that point
    [ 20.5,  80.5, 0.5,  2, 0], # no connectivity constraint on that point
    ]).T # variable values: (nv, npt)-array
hd = gn.img.PointSet(npt=npt, nv=nv, varname=varname, val=v)

# Get the colors for values of the variable of index 3 in the point set,
# according to the color settings used for the TI
hd_col = gn.imgplot.get_colors_from_values(hd.val[3], categ=True, categVal=categ_val, categCol=categ_col)

# Get the colors for geobody label (variable of index 4)
gb_col = gn.imgplot.get_colors_from_values(hd.val[4], categ=True, categVal=[0, 1], categCol=['gray', 'red'])

# Set an image with simulation grid geometry defined above, and no variable
im_empty = gn.img.Img(nx, ny, nz, sx, sy, sz, ox, oy, oz, nv=0)

# Plot
plt.figure(figsize=(8,5))

# Plot empty simulation grid and specify colors
gn.imgplot.drawImage2D(im_empty, categ=True, categVal=categ_val, categCol=categ_col)

# Add hard data points
plt.scatter(hd.x(), hd.y(), marker='o', s=80, color=hd_col, edgecolors=gb_col, linewidths=2)

plt.title('Hard data (and geobody label - edges)')
plt.show()

Specifiy the connectivity constraints (class geone.deesseinterface.Connectivity)

The way how to deal with the connectivity data is described by the class geone.deesseinterface.Connectivity, whose the parameters (keyword arguments) are:

• connectivityConstraintUsage: integer with the following meaning:

  • 0: no connectivity constraint

  • 1: paste connecting paths before simulation by successively binding the cells to be connected in a random order

  • 2: paste connecting paths before simulation by successively binding the cells to be connected beginning with the pair of cells with the smallest distance and then the remaining cells in increasing order according to their distance to the set of cells already connected (the distance between two cells is defined as the length (in number of cells) of the minimal path binding the two nodes in an homogeneous connected medium according to the type of connectivity connectivityType)

  • 3: check connectivity pattern during simulation (not working properly…)

• connectivityType: a character string indicating which type of connection is used:

  • connect_face: connection through the faces of the cells (“6-neighbors connection”)

  • connect_face_edge: connection through the faces or edges of the cells (“18-neighbors connection”)

  • connect_face_edge_corner: connection through the faces or edges or corner (vertices) of the cells (“26-neighbors connection”)

• nclass and classInterval: define the number of classes, and for each class the ensemble of values that can be connected together as a (union of) interval(s) (these parameters are defined as for probability constraints, see the corresponding example)

• varname: the name of the connected component label, should be present in a conditioning data set (geobody_label in this example)

• tiAsRefFlag: boolean (used if connectivityConstraintUsage=1 or 2) indicating that the paths pasted in the simulation grid are borrowed from the TI (True) or from the image given by refConnectivityImage (False)

• refConnectivityImage and refConnectivityVarIndex : image (class geone.img.Img) and index of the variable used (integer) from which the paths pasted in the simulation grid are borrowed (if connectivityConstraintUsage=1 or 2 and tiAsRefFlag=False)

• deactivationDistance and threshold: parameters for connectivityConstraintUsage=3

# Define connectivity constraints
co = gn.deesseinterface.Connectivity(
    connectivityConstraintUsage=2,
    connectivityType='connect_face_edge_corner',
    nclass=3,
    classInterval=[[-0.5, 0.5], [0.5, 1.5], [1.5, 2.5]],
    varname='geobody_label',
    tiAsRefFlag=True
)

Fill the input structure for deesse and launch deesse

nreal = 20
deesse_input = gn.deesseinterface.DeesseInput(
    nx=nx, ny=ny, nz=nz,
    sx=sx, sy=sy, sz=sz,
    ox=ox, oy=oy, oz=oz,
    nv=1, varname='code',
    TI=ti,
    dataPointSet=hd,          # set hard data and connectivity data
    connectivity=co,          # set connectivity constraint
    outputPathIndexFlag=True, # get simulation path in output (to plot pasted paths)
    distanceType='categorical',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.25,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 19 thread(s)]
deesseRun: DeeSse run complete
deesseRun: warnings encountered (1 times in all):
#  1: WARNING 00143: conditioning variable(s) is (are) ignored from data point set (not a simulated variable!)
Elapsed time: 4.8 sec

Retrieve the results, do some statistics, and display

# Retrieve the results
sim = deesse_output['sim']
path = deesse_output['path']

# Gather the nreal realizations into one image
all_sim = gn.img.gatherImages(sim) # all_sim is one image with nreal variables

# Do statistics over all the realizations: compute the pixel-wise proportion for the given categories
all_sim_stats = gn.img.imageCategProp(all_sim, categ_val) # do statistics (pixel-wise proportion)

# Add 10 to cells coming from hard data set or from a pasted path (simulation path index is nan):
# this trick will allow to highlight these cells...
for i in range(nreal):
    sim[i].val = sim[i].val + 10 * np.isnan(path[i].val)

# Display
new_col=['blue', 'lightgreen', 'yellow'] # colors for the "new" facies: 10, 11, 12

prop_col=['blue', 'darkgreen', 'orange'] # colors for the proportion maps
cmap = [gn.customcolors.custom_cmap(['white', c]) for c in prop_col]

plt.subplots(2, 3, figsize=(17,10), sharex=True, sharey=True) # 2 x 3 sub-plots
# ... the first three realizations
for i in range(3):
    plt.subplot(2, 3, i+1) # select next sub-plot
    gn.imgplot.drawImage2D(sim[i], categ=True,
                           categVal=categ_val + [10, 11, 12], # category to be displayed (+ for concatenation)
                           categCol=categ_col + new_col,      # colors used
                           title=f'Real #{i}')
    plt.scatter(hd.x(), hd.y(), marker='o', s=80,
                color='none', edgecolors='black', linewidths=1) # add hard data points

for i in range(3):
    plt.subplot(2, 3, i+4) # select next sub-plot
    gn.imgplot.drawImage2D(all_sim_stats, iv=i, cmap=cmap[i],
                           title=f'Prop. of categ. {i} (over {nreal} real.)')
    plt.plot(hd.x(), hd.y(), '+', markersize=20, c='black') # add hard data points

plt.show()

# Do statistics on pasted paths
for i in range(nreal):
    path[i].val = 1 * np.isnan(path[i].val)

all_path = gn.img.gatherImages(path) # all_path is one image with nreal variables
all_path_stats = gn.img.imageCategProp(all_path, [1]) # do statistics (pixel-wise proportion)

# ...multiply by nreal to get the number of times a cell is in a pasted path (or coming from the hard data set)
all_path_stats.val = nreal*all_path_stats.val

# Display the number of times a cell is in a pasted path (or coming from the hard data set)
plt.figure(figsize=(8,5))
gn.imgplot.drawImage2D(all_path_stats,
                       cmap=gn.customcolors.cmap1, # color map used (this is the default one)
                       vmin=1, # ... value zero will be displayed with color 'cunder'
                               # from gn.customcolors.cmap1
                       title=f'Coming from HD or pasted paths (over {nreal} real.)')
plt.scatter(hd.x(), hd.y(), marker='o', s=80,
            color='none', edgecolors='white', linewidths=2) # add hard data points
plt.show()