GEONE - GEOSCLASSIC - categorical variable - Examples in 2D - non-stationary covariance model
Estimation (kriging) and simulation (Sequential Indicator Simulation, SIS)
See notebook ex_geosclassic_indicator_1d_1.ipynb for detail explanations about estimation (kriging) and simulation (Sequential Indicator Simulation, SIS) in a grid.
Non-stationary covariance model over a grid
See notebook ex_geosclassic_1d_2_non_stat_cov.ipynb for detail explanations on how to set non-stationarities in a grid.
Examples in 2D
In this notebook, examples in 2D with a non-stationary covariance model are given.
Import what is required
[1]:
import numpy as np
import matplotlib.pyplot as plt
import time
# import package 'geone'
import geone as gn
[2]:
# 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.1
Remark
The matplotlib figures can be visualized in interactive mode:
%matplotlib notebook: enable interactive mode%matplotlib inline: disable interactive mode
Category values
A list of category values (facies) must be defined. Let ncategory be the length of this list, i.e. the number of categories:
if
ncategory == 1: the unique category value given must not be equal to 0; this is used for a binary case with values (“unique category value”, 0), where 0 indicates the absence of the considered medium; conditioning data values should be “unique category value” or 0if
ncategory >= 2: this is used for a multi-category case with given values (distinct); conditioning data values should be in the list of given values
Then, set color for each category, and color maps for proportions (for further plots).
Below: select the case with ``ncategory`` greater than one or equal to one below, comment the undesired cell.
[3]:
# Case with ncategory > 1
# -----------------------
category_values = [1., 2., 3.]
ncategory = len(category_values)
# Set colors ...
categVal = category_values
categCol = ['lightblue', 'orange', 'darkgreen'] # must be of length len(categVal)
cmap_categ = [gn.customcolors.custom_cmap(['white', c]) for c in categCol]
[4]:
# # Case with ncategory = 1
# # -----------------------
# category_values = [2.] # all categories are 2. and 0.
# ncategory = len(category_values)
# # Set colors ...
# categVal = [category_values[0], 0]
# categCol = ['tab:red', 'lightblue'] # must be of length len(categVal)
# cmap_categ = [gn.customcolors.custom_cmap(['white', c]) for c in categCol]
Grid (2D)
[5]:
nx, ny = 220, 230 # number of cells
sx, sy = 1.0, 1.0 # cell unit
ox, oy = 0.0, 0.0 # origin
dimension = (nx, ny)
spacing = (sx, sy)
origin = (ox, oy)
Covariance model
In 2D, a covariance model is given by an instance of the class geone.covModel.covModel2D (or geone.covModel.covModel1D for omni-directional (isotropic) case).
Base covariance model (sationary)
The weight 'w' to every elementary contribution is set to 1.0; the method multiply_w will be used to set non-stationarities about this parameter; angle alpha is set to 0, local rotation will be set further.
[6]:
# Define the base covariance model (stationary)
cov_model = gn.covModel.CovModel2D(elem=[
('spherical', {'w':1.0, 'r':[120.0, 30.0]}), # elementary contribution
('gaussian', {'w':1.0, 'r':[120.0, 30.0]}), # elementary contribution
], alpha=0, name='model-2D example')
[7]:
# cov_model.plot_model(figsize=(15,5))
# plt.suptitle('Covariance function - base')
# plt.show()
Defining non-stationarities
[8]:
# Set an image with grid geometry defined above, and no variable
im = gn.img.Img(nx, ny, 1, sx, sy, 1., ox, oy, 0., nv=0)
# Get the x and y coordinates of the centers of grid cell (meshgrid)
xx = im.xx()[0]
yy = im.yy()[0]
# Center of the grid
x_center = 0.5*(im.xmin() + im.xmax())
y_center = 0.5*(im.ymin() + im.ymax())
# Set alpha in degrees in the grid for further estimation/simulation
# ------------------------------------------------------------------
t = 180.0/np.pi
alpha = 90.0 - np.arctan2(yy-y_center, xx-x_center)* t - 20.
# Note: `alpha` could also be a function of two parameters (x, y location)
# Define weight for gaussian model over the simulation grid
gau_w = 9. * 1. / (1. + np.exp(-(np.sqrt((xx-x_center)**2+(yy-y_center)**2)-50)/20))
# gau_w = 9. * 1. / (1. + np.exp(-(np.sqrt((xx-x_center)**2+(yy-y_center)**2)-170)/30))
sph_w = 9. - gau_w
# Set list to handle non-stationarities for further estimation/simulation in the grid
# -----------------------------------------------------------------------------------
cov_model_non_stationarity_list = [
('multiply_w', sph_w, {'elem_ind':0}), # multiply weight by `sph_w` for elem. contrib. of index 0
('multiply_w', gau_w, {'elem_ind':1}), # multiply weight by `gau_w` for elem. contrib. of index 1
]
# Note: `gau_w`, `nug_w` could also be a function of two parameters (x, y location)
[9]:
# Plot non-stationarities
# Set variable alpha, sph_w and gau_w in image im
im.append_var([alpha, sph_w, gau_w], varname=['alpha', 'sph_w', 'gau_w'])
# Plot
plt.subplots(1,3, figsize=(15, 4), sharex=True, sharey=True)
plt.subplot(1,3,1)
gn.imgplot.drawImage2D(im, iv=0, cmap='terrain', title="alpha")
len_arrow = 20.
for i in range(0, nx, 20):
for j in range(0, ny, 20):
u = xx[j, i]
v = yy[j, i]
a = -alpha[j,i]/t
plt.arrow(u, v, len_arrow*np.cos(a), len_arrow*np.sin(a), head_width=3)
plt.subplot(1,3,2)
gn.imgplot.drawImage2D(im, iv=1, cmap='terrain', title="spherical weight")
plt.subplot(1,3,3)
gn.imgplot.drawImage2D(im, iv=2, cmap='terrain', title="gaussian weight")
plt.show()
1. Example
Settings - using data (optional) and probability (constant, optional)
[10]:
if ncategory > 1:
# Case with ncategory > 1
# -----------------------
# Data
x = np.array([[ 10.42, 20.32],
[ 50.5 , 40.5],
[ 20.24, 150.7],
[200.2 , 210.1]]) # data locations (real coordinates)
v = [ 1., 2., 1., 3.] # data values
# x = None
# v = None
# Probability, proportion of each category
probability = [.1, .2, .7] # should sum to 1
# probability = None
# Type of kriging
method = 'simple_kriging'
else:
# Case with ncategory = 1
# -----------------------
# Data
x = np.array([[ 10.42, 20.32],
[ 50.5 , 40.5],
[ 20.24, 150.7],
[200.2 , 210.1]]) # data locations (real coordinates)
v = [ 0., 2., 2., 0.] # data values
# x = None
# v = None
# Probability, proportion (of non-zero category)
probability = [.7] # list of one number in the interval [0, 1]
# probability = None
# Type of kriging
method = 'simple_kriging'
Estimation of probabilities (by kriging)
[11]:
# Computational resources
nthreads = 8
t1 = time.time() # start time
geosclassic_output = gn.geosclassicinterface.estimateIndicator(
category_values, # list of categories (required)
cov_model, # covariance model(s) (required)
dimension, spacing, origin, # grid geometry (dimension is required)
x=x, v=v, # data
probability=probability, # probability
alpha=alpha, # rotation
cov_model_non_stationarity_list=cov_model_non_stationarity_list, # non-stationrities
method=method, # type of kriging
use_unique_neighborhood=False, # search neighborhood ...
searchRadius=None,
searchRadiusRelative=1.2,
nneighborMax=12,
nthreads=nthreads, # computational resources
verbose=2 # verbose mode
)
t2 = time.time() # end time
print('Elapsed time: {:.2g} sec'.format(t2-t1))
krig_img = geosclassic_output['image'] # output image
estimateIndicator: call `run_MPDSOMPGeosClassicIndicatorSim` [1 process of 8 thread(s) (OpenMP)] ...
estimateIndicator: `run_MPDSOMPGeosClassicIndicatorSim` [1 process] complete
estimateIndicator: warnings encountered (25371 times in all):
# 1: WARNING 02001: a neigbhor has been dropped (solving kriging system)
# 2: WARNING 02015: solving kriging system fails (do as if no neighbor)
Elapsed time: 9.6 sec
[12]:
# Total number of warning(s), and warning messages
geosclassic_output['nwarning'], geosclassic_output['warnings']
[12]:
(25371,
['WARNING 02001: a neigbhor has been dropped (solving kriging system)',
'WARNING 02015: solving kriging system fails (do as if no neighbor)'])
[13]:
%%script false --no-raise-error # skip this cell! (comment this line to run the cell)
# Equivalent, using other computational resources
# -----------------------------------------------
# Computational resources
nthreads = 1
t1 = time.time() # start time
geosclassic_output = gn.geosclassicinterface.estimateIndicator(
category_values, # list of categories (required)
cov_model, # covariance model(s) (required)
dimension, spacing, origin, # grid geometry (dimension is required)
x=x, v=v, # data
probability=probability, # probability
alpha=alpha, # rotation
cov_model_non_stationarity_list=cov_model_non_stationarity_list, # non-stationrities
method=method, # type of kriging
use_unique_neighborhood=False, # search neighborhood ...
searchRadius=None,
searchRadiusRelative=1.2,
nneighborMax=12,
nthreads=nthreads, # computational resources
verbose=2 # verbose mode
)
t2 = time.time() # end time
print('Elapsed time: {:.2g} sec'.format(t2-t1))
krig_img_2 = geosclassic_output['image'] # output image
print(f"Same results ? {np.allclose(krig_img.val, krig_img_2.val)}")
Simulations
[14]:
# Number of realizations
nreal = 100
# Seed
seed = 321
# Computational resources
nproc = 2
nthreads_per_proc = 4
t1 = time.time() # start time
geosclassic_output = gn.geosclassicinterface.simulateIndicator(
category_values, # list of categories (required)
cov_model, # covariance model(s) (required)
dimension, spacing, origin, # grid geometry (dimension is required)
x=x, v=v, # data
probability=probability, # probability
alpha=alpha, # rotation
cov_model_non_stationarity_list=cov_model_non_stationarity_list, # non-stationrities
method=method, # type of kriging
searchRadius=None,
searchRadiusRelative=1.2,
nneighborMax=12,
nreal=nreal, # number of realizations
seed=seed, # seed
nproc=nproc, # computational resources ...
nthreads_per_proc=nthreads_per_proc,
verbose=2 # verbose mode
)
t2 = time.time() # end time
print('Elapsed time: {:.2g} sec'.format(t2-t1))
simul_img = geosclassic_output['image'] # output image
simulateIndicator: call `run_MPDSOMPGeosClassicIndicatorSim` [2 process(es) of 4 thread(s) (OpenMP)] ...
simulateIndicator: `run_MPDSOMPGeosClassicIndicatorSim` [2 process(es)] complete
simulateIndicator: warnings encountered (324 times in all):
# 1: WARNING 02001: a neigbhor has been dropped (solving kriging system)
# 2: WARNING 02015: solving kriging system fails (do as if no neighbor)
Elapsed time: 18 sec
[15]:
# Total number of warning(s), and warning messages
geosclassic_output['nwarning'], geosclassic_output['warnings']
[15]:
(324,
[np.str_('WARNING 02001: a neigbhor has been dropped (solving kriging system)'),
np.str_('WARNING 02015: solving kriging system fails (do as if no neighbor)')])
[16]:
%%script false --no-raise-error # skip this cell! (comment this line to run the cell)
# Equivalent, using other computational resources
# -----------------------------------------------
# Computational resources
nproc = 1
nthreads_per_proc = 1
t1 = time.time() # start time
geosclassic_output = gn.geosclassicinterface.simulateIndicator(
category_values, # list of categories (required)
cov_model, # covariance model(s) (required)
dimension, spacing, origin, # grid geometry (dimension is required)
x=x, v=v, # data
probability=probability, # probability
alpha=alpha, # rotation
cov_model_non_stationarity_list=cov_model_non_stationarity_list, # non-stationrities
method=method, # type of kriging
searchRadius=None,
searchRadiusRelative=1.2,
nneighborMax=12,
nreal=nreal, # number of realizations
seed=seed, # seed
nproc=nproc, # computational resources ...
nthreads_per_proc=nthreads_per_proc,
verbose=2 # verbose mode
)
t2 = time.time() # end time
print('Elapsed time: {:.2g} sec'.format(t2-t1))
simul_img_2 = geosclassic_output['image'] # output image
print(f"Same results ? {np.allclose(simul_img.val, simul_img_2.val)}")
Plot the results
[17]:
# Compute proportion of each category (pixel-wise)
simul_img_prop = gn.img.imageCategProp(simul_img, category_values)
[18]:
# Plot
plt.subplots(1+ncategory, 2, figsize=(10, 4*(1+ncategory)), sharex=True, sharey=True)
for i in range(2):
plt.subplot(1+ncategory, 2, i+1)
gn.imgplot.drawImage2D(simul_img, iv=i, categ=True, categVal=categVal, categCol=categCol)
if x is not None:
plt.plot(*x.T, '+', c='black', markersize=5) # add conditioning point locations
plt.title(f'real #{i}')
for i in range(ncategory):
plt.subplot(1+ncategory, 2, 2*i+3)
gn.imgplot.drawImage2D(simul_img_prop, iv=i, vmin=0, vmax=1, cmap=cmap_categ[i])
if x is not None:
plt.plot(*x.T, '+', c='black', markersize=5) # add conditioning point locations
plt.title(f'Proportion of categ {categVal[i]} over {nreal} real')
plt.subplot(1+ncategory, 2, 2*i+4)
gn.imgplot.drawImage2D(krig_img, iv=i, vmin=0, vmax=1, cmap=cmap_categ[i])
if x is not None:
plt.plot(*x.T, '+', c='black', markersize=5) # add conditioning point locations
plt.title(f'Estimate probability of categ {categVal[i]} (kriging)')
plt.show()
Check results
[19]:
# Check data
# ----------
if x is not None:
# Get index of conditioning location in the grid
data_grid_index = [gn.img.pointToGridIndex(xk[0], xk[1], 0, sx, sy, 1., ox, oy, 0.) for xk in x] # (ix, iy, iz) for each data point
# Check estimation
krig_v = [krig_img.val[:, iz, iy, ix] for ix, iy, iz in data_grid_index]
if ncategory == 1:
print(f'Estimation: all data respected ? {np.all(np.asarray(krig_v).reshape(-1) == np.asarray([1 if vi == category_values[0] else 0 for vi in v]))}')
else:
print(f'Estimation: all data respected ? {np.all([np.all(krig_v[i] == np.eye(ncategory)[np.where(np.asarray(category_values) == v[i])[0][0]]) for i in range(len(x))])}')
# Check simulation
sim_v = [simul_img.val[:, iz, iy, ix] for ix, iy, iz in data_grid_index]
print(f'Simulation: all data respected ? {np.all([np.all(sim_v[i] == v[i]) for i in range(len(x))])}')
Estimation: all data respected ? True
Simulation: all data respected ? True
[20]:
# Compare probabilities
# ---------------------
if probability is not None:
print(f'Prescribed probabilities = {probability}')
print(f'Estimation: probabilities (mean over the grid) = {krig_img.val.mean(axis=(1,2,3))}')
print(f'Simulation: probabilities (mean over the grid and all real.)= {[np.mean(simul_img.val == cv) for cv in category_values]}')
Prescribed probabilities = [0.1, 0.2, 0.7]
Estimation: probabilities (mean over the grid) = [0.17822557 0.22269518 0.59907925]
Simulation: probabilities (mean over the grid and all real.)= [np.float64(0.1955203557312253), np.float64(0.24021640316205534), np.float64(0.5642632411067193)]
2. Example - using non-stationary probabilities
Setting probability (proportion) maps
[21]:
# Set an image with simulation grid geometry defined above, and no variable
im = gn.img.Img(nx, ny, 1, sx, sy, 1., ox, oy, 0., nv=0)
# Get the x and y coordinates of the centers of grid cell (meshgrid)
xx = im.xx()[0]
yy = im.yy()[0]
# Equivalent:
## xg, yg: coordinates of the centers of grid cell
#xg = ox + 0.5*sx + sx*np.arange(nx)
#yg = oy + 0.5*sy + sy*np.arange(ny)
#xx, yy = np.meshgrid(xg, yg) # create meshgrid from the center of grid cells
if ncategory > 1:
# Case with ncategory > 1
# -----------------------
# Define probability maps for each category
c = 0.9
p1 = xx - yy
p1 = c * (p1 - np.min(p1))/ (np.max(p1) - np.min(p1))
p2 = c - p1
p0 = (1. - c) * np.ones_like(p1) # 1.0 - p1 - p2 # constant map (0.1)
probability = np.array((p0, p1, p2))
else:
# Case with ncategory = 1
# -----------------------
c = 1.0
# Define probability map for non-zero category
p1 = xx - yy
p1 = c * (p1 - np.min(p1))/ (np.max(p1) - np.min(p1))
probability = p1
[22]:
# Plot
# ----
# Fill image for display
probability_img = gn.img.Img(nx, ny, 1, sx, sy, 1., ox, oy, 0., nv=ncategory, val=probability)
# Display probability maps
plt.subplots(ncategory, 1, figsize=(10,ncategory*4), sharey=True)
for i in range(ncategory):
plt.subplot(ncategory, 1, 1+i)
gn.imgplot.drawImage2D(probability_img, iv=i, title = f'Probability for categ. {categVal[i]}')
plt.show()
Settings - using data (optional)
[23]:
if ncategory > 1:
# Case with ncategory > 1
# -----------------------
# Data
x = np.array([[ 10.42, 20.32],
[ 50.5 , 40.5],
[ 20.24, 150.7],
[200.2 , 210.1]]) # data locations (real coordinates)
v = [ 1., 2., 1., 3.] # data values
# x = None
# v = None
# Probability : `probability` defined above
# Type of kriging
method = 'simple_kriging'
else:
# Case with ncategory = 1
# -----------------------
# Data
x = np.array([[ 10.42, 20.32],
[ 50.5 , 40.5],
[ 20.24, 150.7],
[200.2 , 210.1]]) # data locations (real coordinates)
v = [ 0., 2., 2., 0.] # data values
# x = None
# v = None
# Probability : `probability` defined above
# Type of kriging
method = 'simple_kriging'
Estimation of probabilities (by kriging)
[24]:
# Computational resources
nthreads = 8
t1 = time.time() # start time
geosclassic_output = gn.geosclassicinterface.estimateIndicator(
category_values, # list of categories (required)
cov_model, # covariance model(s) (required)
dimension, spacing, origin, # grid geometry (dimension is required)
x=x, v=v, # data
probability=probability, # probability
alpha=alpha, # rotation
cov_model_non_stationarity_list=cov_model_non_stationarity_list, # non-stationrities
method=method, # type of kriging
use_unique_neighborhood=False, # search neighborhood ...
searchRadius=None,
searchRadiusRelative=1.2,
nneighborMax=12,
nthreads=nthreads, # computational resources
verbose=2 # verbose mode
)
t2 = time.time() # end time
print('Elapsed time: {:.2g} sec'.format(t2-t1))
krig_img = geosclassic_output['image'] # output image
estimateIndicator: call `run_MPDSOMPGeosClassicIndicatorSim` [1 process of 8 thread(s) (OpenMP)] ...
estimateIndicator: `run_MPDSOMPGeosClassicIndicatorSim` [1 process] complete
estimateIndicator: warnings encountered (25371 times in all):
# 1: WARNING 02001: a neigbhor has been dropped (solving kriging system)
# 2: WARNING 02015: solving kriging system fails (do as if no neighbor)
Elapsed time: 11 sec
[25]:
# Total number of warning(s), and warning messages
geosclassic_output['nwarning'], geosclassic_output['warnings']
[25]:
(25371,
['WARNING 02001: a neigbhor has been dropped (solving kriging system)',
'WARNING 02015: solving kriging system fails (do as if no neighbor)'])
[26]:
%%script false --no-raise-error # skip this cell! (comment this line to run the cell)
# Equivalent, using other computational resources
# -----------------------------------------------
# Computational resources
nthreads = 1
t1 = time.time() # start time
geosclassic_output = gn.geosclassicinterface.estimateIndicator(
category_values, # list of categories (required)
cov_model, # covariance model(s) (required)
dimension, spacing, origin, # grid geometry (dimension is required)
x=x, v=v, # data
probability=probability, # probability
alpha=alpha, # rotation
cov_model_non_stationarity_list=cov_model_non_stationarity_list, # non-stationrities
method=method, # type of kriging
use_unique_neighborhood=False, # search neighborhood ...
searchRadius=None,
searchRadiusRelative=1.2,
nneighborMax=12,
nthreads=nthreads, # computational resources
verbose=2 # verbose mode
)
t2 = time.time() # end time
print('Elapsed time: {:.2g} sec'.format(t2-t1))
krig_img_2 = geosclassic_output['image'] # output image
print(f"Same results ? {np.allclose(krig_img.val, krig_img_2.val)}")
Simulations
[27]:
# Number of realizations
nreal = 100
# Seed
seed = 321
# Computational resources
nproc = 2
nthreads_per_proc = 4
t1 = time.time() # start time
geosclassic_output = gn.geosclassicinterface.simulateIndicator(
category_values, # list of categories (required)
cov_model, # covariance model(s) (required)
dimension, spacing, origin, # grid geometry (dimension is required)
x=x, v=v, # data
probability=probability, # probability
alpha=alpha, # rotation
cov_model_non_stationarity_list=cov_model_non_stationarity_list, # non-stationrities
method=method, # type of kriging
searchRadius=None,
searchRadiusRelative=1.2,
nneighborMax=12,
nreal=nreal, # number of realizations
seed=seed, # seed
nproc=nproc, # computational resources ...
nthreads_per_proc=nthreads_per_proc,
verbose=2 # verbose mode
)
t2 = time.time() # end time
print('Elapsed time: {:.2g} sec'.format(t2-t1))
simul_img = geosclassic_output['image'] # output image
simulateIndicator: call `run_MPDSOMPGeosClassicIndicatorSim` [2 process(es) of 4 thread(s) (OpenMP)] ...
simulateIndicator: `run_MPDSOMPGeosClassicIndicatorSim` [2 process(es)] complete
simulateIndicator: warnings encountered (324 times in all):
# 1: WARNING 02001: a neigbhor has been dropped (solving kriging system)
# 2: WARNING 02015: solving kriging system fails (do as if no neighbor)
Elapsed time: 19 sec
[28]:
# Total number of warning(s), and warning messages
geosclassic_output['nwarning'], geosclassic_output['warnings']
[28]:
(324,
[np.str_('WARNING 02001: a neigbhor has been dropped (solving kriging system)'),
np.str_('WARNING 02015: solving kriging system fails (do as if no neighbor)')])
[29]:
%%script false --no-raise-error # skip this cell! (comment this line to run the cell)
# Equivalent, using other computational resources
# -----------------------------------------------
# Computational resources
nproc = 1
nthreads_per_proc = 1
t1 = time.time() # start time
geosclassic_output = gn.geosclassicinterface.simulateIndicator(
category_values, # list of categories (required)
cov_model, # covariance model(s) (required)
dimension, spacing, origin, # grid geometry (dimension is required)
x=x, v=v, # data
probability=probability, # probability
alpha=alpha, # rotation
cov_model_non_stationarity_list=cov_model_non_stationarity_list, # non-stationrities
method=method, # type of kriging
searchRadius=None,
searchRadiusRelative=1.2,
nneighborMax=12,
nreal=nreal, # number of realizations
seed=seed, # seed
nproc=nproc, # computational resources ...
nthreads_per_proc=nthreads_per_proc,
verbose=2 # verbose mode
)
t2 = time.time() # end time
print('Elapsed time: {:.2g} sec'.format(t2-t1))
simul_img_2 = geosclassic_output['image'] # output image
print(f"Same results ? {np.allclose(simul_img.val, simul_img_2.val)}")
Plot the results
[30]:
# Compute proportion of each category (pixel-wise)
simul_img_prop = gn.img.imageCategProp(simul_img, category_values)
[31]:
# Plot
plt.subplots(1+ncategory, 3, figsize=(15, (1+ncategory)*4), sharex=True, sharey=True)
for i in range(3):
plt.subplot(1+ncategory, 3, i+1)
gn.imgplot.drawImage2D(simul_img, iv=i, categ=True, categVal=categVal, categCol=categCol)
if x is not None:
plt.plot(*x.T, '+', c='black', markersize=5) # add conditioning point locations
plt.title(f'real #{i}')
for i in range(ncategory):
plt.subplot(1+ncategory, 3, 3*i+4)
gn.imgplot.drawImage2D(simul_img_prop, iv=i, vmin=0, vmax=1, cmap=cmap_categ[i])
if x is not None:
plt.plot(*x.T, '+', c='black', markersize=5) # add conditioning point locations
plt.title(f'Proportion of categ {categVal[i]} over {nreal} real')
plt.subplot(1+ncategory, 3, 3*i+5)
gn.imgplot.drawImage2D(krig_img, iv=i, vmin=0, vmax=1, cmap=cmap_categ[i])
if x is not None:
plt.plot(*x.T, '+', c='black', markersize=5) # add conditioning point locations
plt.title(f'Estimate probability of categ {categVal[i]} (kriging)')
plt.subplot(1+ncategory, 3, 3*i+6)
gn.imgplot.drawImage2D(probability_img, iv=i, vmin=0, vmax=1, cmap=cmap_categ[i])
if x is not None:
plt.plot(*x.T, '+', c='black', markersize=5) # add conditioning point locations
plt.title(f'Prescribed probability for categ {categVal[i]}')
plt.show()
Check results
[32]:
# Check data
# ----------
if x is not None:
# Get index of conditioning location in the grid
data_grid_index = [gn.img.pointToGridIndex(xk[0], xk[1], 0, sx, sy, 1., ox, oy, 0.) for xk in x] # (ix, iy, iz) for each data point
# Check estimation
krig_v = [krig_img.val[:, iz, iy, ix] for ix, iy, iz in data_grid_index]
if ncategory == 1:
print(f'Estimation: all data respected ? {np.all(np.asarray(krig_v).reshape(-1) == np.asarray([1 if vi == category_values[0] else 0 for vi in v]))}')
else:
print(f'Estimation: all data respected ? {np.all([np.all(krig_v[i] == np.eye(ncategory)[np.where(np.asarray(category_values) == v[i])[0][0]]) for i in range(len(x))])}')
# Check simulation
sim_v = [simul_img.val[:, iz, iy, ix] for ix, iy, iz in data_grid_index]
print(f'Simulation: all data respected ? {np.all([np.all(sim_v[i] == v[i]) for i in range(len(x))])}')
Estimation: all data respected ? True
Simulation: all data respected ? True