Localized Uniform Conditioning¶
Overview¶
Localized Uniform conditioning (LUC) is a estimation method that models the conditional distribution of panels, which usually are blocks 2 to 4 times the size of a SMU. SMU grades are then computed through the "localization" of the panel results using kriging applied to the SMU.
Dataset¶
The 2D Walker Lake dataset will be used for demonstration purposes.
import numpy as np
import matplotlib.pyplot as plt
import geolime as geo
import plotly.io as pio
pio.renderers.default = "sphinx_gallery"
data = geo.datasets.load("walker_lake/walker_lake_sample.csv")
data.head()
| x | y | z | v | u | t | |
|---|---|---|---|---|---|---|
| 0 | 11.0 | 8.0 | 0.0 | 0.0 | NaN | 2 |
| 1 | 8.0 | 30.0 | 0.0 | 0.0 | NaN | 2 |
| 2 | 9.0 | 48.0 | 0.0 | 224.4 | NaN | 2 |
| 3 | 8.0 | 69.0 | 0.0 | 434.4 | NaN | 2 |
| 4 | 9.0 | 90.0 | 0.0 | 412.1 | NaN | 2 |
point_cloud = geo.PointCloud(name='WalkerLake', xyz=data[['x', 'y', 'z']])
point_cloud.set_property(name='V', data=data['v'])
Change of support¶
LUC is nonlinear and requires a probabilistic model for the regionalized variable of interest to be considered. It relies on Gaussian assumptions and an anamorphosis is then introduced in order to apply it through change of support. The following cells details how it is possible to compute the change of support coefficient that permits to estimate panel quantities from SMUs. We first process the non-Gaussian V variable for anamorphosis.
The V variable follows a non normal distribution
geo.histogram_plot(
[{"object": point_cloud, "property": "V"}],
width=650,
height=400
)
gauss_u = geo.gauss_score(point_cloud.data("V"))
The gaussian scores map the V values in the gaussian scpace (mean of 0 and variance of 1). This allows a unique relationship between original values and values from a normal distribution.
fig, ax = plt.subplots()
ax.scatter(gauss_u, point_cloud.data("V"), color ='b' , s = 0.2 )
ax.set_title("Gauss scores vs observed values")
plt.xlabel('Gauss scores')
plt.ylabel('Values')
plt.show()
Hermite polynomials¶
Hermite orthogonal polynomials are very conveniant tool for approximating a Gaussian random variable. In terms of computation, they are widely used in order to avoid computing a cumulative distribution function through an integral. GeoLime provides a set of API calls to expand a variable of interest (Gaussian, for instance after having transformed a non-gaussian variable with an anamorphosis through Gauss scores, as seen previously). Hermite polynomials expansion can be retrieved with the following calls.
anam = geo.Anamorphosis(
data=point_cloud.data("V"),
coefficient_number=50
)
Gaussian scores have been interpolated using Hermite polynomials so that every values in the range of the V property can be transformed into the gaussian space.
x = np.linspace(-3, 3, 100)
fig, ax = plt.subplots()
ax.plot(x, anam.theoretical_transform(x), color ='b')
ax.set_title("Hermite expansion evaluation")
plt.show()
This anamorphosis object that handles Hermite polynomials will be used later on for estimation.
Estimation¶
The data are now gaussian and a change of support can be computed, and needs 3 parameters :
- Hermite coefficients of the anamorphosis transform
- Variogram model
- A "support", which is the size of a regular block unit
Uniform conditioning requires two interpolation supports:
- a SMU support consisting in cell of the exploitation size
- a Panel support consisting in cell of the size of generally 2-4 SMU, here we take 2
grid = geo.Voxel(
name='WalkerLakeSMU',
shape=[26, 30],
axis=[[10., 0., 0.], [0., 10., 0.]],
)
panel = geo.Voxel(
name='WalkerLakePanel',
shape=[13, 15],
axis=[[20., 0., 0.], [0., 20., 0.]],
)
Here we define a punctual covariance model. The model parameters come from the Variography guide notebook.
covariance = (
7000 * geo.Nugget(dimension=2)
+ 88000 * geo.Spherical(
dimension=2,
metric=geo.Metric(angles=160,scales=[20.,15.])
)
)
neighborhood = geo.MinMaxPointsNeighborhood(
dimension=2,
min_n = 1,
max_n=25
)
uc_estimator = geo.LocalisedUniformConditioning(
covariance_model=covariance,
neighborhood_model=neighborhood,
axes = [geo.Coord.U, geo.Coord.V]
)
uc_estimator.solve(
obj=point_cloud,
obj_region=None,
obj_attribute='V',
anamorphosis=anam,
support=grid,
support_region=None,
panel_support=panel,
panel_support_region=None,
cutoff_grades=np.linspace(0, 1800, 10, dtype=np.float64),
support_attribute='test',
discr=np.array([5,5],dtype=np.int64),
panel_discr=np.array([5,5],dtype=np.int64)
)
geo.plot(grid, property="test", agg_method="mean", width=600, height=650)