đŻ Objectif gĂ©nĂ©ral¶
Comprendre comment chaque paire de points dâĂ©chantillonnage contribue Ă la nuĂ©e variographique, et observer lâeffet du nombre de points sur la qualitĂ© du variogramme obtenu.
đ§ Ce que vous allez apprendre¶
- Comment se forment les nuages de demi-variances selon la distance.
- Comment on passe dâune nuĂ©e Ă un variogramme expĂ©rimental.
- Pourquoi un petit nombre dâĂ©chantillons peut fausser lâinterprĂ©tation (variabilitĂ©, instabilitĂ©).
âïž Ce que vous pouvez faire¶
- Générer un champ spatial simulé.
- Choisir un nombre dâĂ©chantillons Ă prĂ©lever.
- Voir toutes les paires de points et leur contribution au variogramme.
- Suivre pas Ă pas la construction dâun variogramme :
- Distance entre deux points
- Valeurs mesurées
- Calcul de la demi-variance
Source
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft2, ifft2
from scipy.spatial.distance import pdist
from ipywidgets import IntText, FloatText, VBox, Button, Output, Tab
import matplotlib.gridspec as gridspec
from IPython.display import display, clear_output
class VariogramWorkshop:
def __init__(self):
# Widgets paramĂštres (recalcul)
self.n_samples_widget = IntText(value=100, description="Nb échantillons")
self.range_widget = FloatText(value=30.0, description="Portée")
self.sill_widget = FloatText(value=1.0, description="Sill")
self.nugget_widget = FloatText(value=0.1, description="Nugget")
self.bin_width_widget = FloatText(value=10.0, description="Bin (m)")
self.calc_button = Button(description="Calculer", button_style='success')
# Bouton suivant
self.next_button = Button(description="Suivant", button_style='info')
self.next_button.disabled = True
# Outputs
self.out_plot = Output()
self.out_text = Output()
# Onglets
self.tabs = Tab(children=[self.out_plot, self.out_text])
self.tabs.set_title(0, "Visualisation")
self.tabs.set_title(1, "Calculs")
# Data placeholders
self.sample_coords = None
self.sample_values = None
self.dists = None
self.gamma = None
self.n_samples = 0
self.field = None
self.pair_index = 0
# Connexions
self.calc_button.on_click(self.on_calculate)
self.next_button.on_click(self.on_next_pair)
# Affichage interface
param_box = VBox([
self.n_samples_widget, self.range_widget, self.sill_widget,
self.nugget_widget, self.bin_width_widget, self.calc_button])
full_ui = VBox([param_box, self.next_button, self.tabs])
display(full_ui)
def generate_gaussian_field_fftma(self, size, range_, sill, nugget, seed=0):
np.random.seed(seed)
extended_size = (2 * size[0], 2 * size[1])
x = np.arange(extended_size[0])
y = np.arange(extended_size[1])
X, Y = np.meshgrid(x, y, indexing='ij')
h = np.sqrt((X - size[0])**2 + (Y - size[1])**2)
cov = sill * np.exp(-h / range_)
cov[0, 0] += nugget
spectrum = np.abs(fft2(cov))**0.5
white_noise = np.random.normal(size=extended_size)
field = np.real(ifft2(fft2(white_noise) * spectrum))
field_cropped = field[:size[0], :size[1]]
std = np.std(field_cropped)
return field_cropped / std * np.sqrt(sill+nugget) if std > 0 else field_cropped
def pdist_to_indices(self, k, n):
i = n - 2 - int(np.floor(np.sqrt(-8*k + 4*n*(n-1)-7)/2.0 - 0.5))
j = k + i + 1 - n*(n-1)//2 + (n - i)*((n - i) -1)//2
return i, j
def on_calculate(self, b):
# Récupération des paramÚtres de l'interface
n_samples = self.n_samples_widget.value
range_ = self.range_widget.value
sill = self.sill_widget.value
nugget = self.nugget_widget.value
bin_width = self.bin_width_widget.value
# Taille de la grille (200x200 ici)
size = (200, 200)
# Génération d'un nouveau champ gaussien
self.field = self.generate_gaussian_field_fftma(size, range_, sill, nugget, seed=12345)
# Sélection aléatoire des échantillons à partir de la grille
sample_indices = np.random.choice(size[0]*size[1], n_samples, replace=False)
self.sample_coords = np.array(np.unravel_index(sample_indices, size)).T
self.sample_values = self.field[self.sample_coords[:, 0], self.sample_coords[:, 1]]
self.n_samples = n_samples
# Calcul des distances entre les échantillons
self.dists = pdist(self.sample_coords)
# Calcul des différences et de la demi-variance (γ)
diffs = pdist(self.sample_values[:, None])
self.gamma = 0.5 * diffs**2
# Calcul du nombre maximal de paires d'échantillons
self.max_pairs = n_samples * (n_samples - 1) // 2
self.pair_index = 0
# Déverrouiller le bouton suivant
self.next_button.disabled = False
# Mettre Ă jour les variables internes
self.range_ = range_
self.sill = sill
self.nugget = nugget
self.bin_width = bin_width
self.size = size
# Mettre Ă jour l'affichage du graphique
self.update_plot()
def update_plot(self):
pair_index = self.pair_index
n_samples = self.n_samples
i, j = self.pdist_to_indices(pair_index, n_samples)
p1 = self.sample_coords[i]
p2 = self.sample_coords[j]
z1 = self.sample_values[i]
z2 = self.sample_values[j]
dist_p = self.dists[pair_index]
gamma_p = self.gamma[pair_index]
# Binning variogramme expérimental
max_dist = 200
bins = np.arange(0, max_dist + self.bin_width, self.bin_width)
bin_centers = 0.5 * (bins[:-1] + bins[1:])
bin_indices = np.digitize(self.dists, bins) - 1
experimental = [np.mean(self.gamma[bin_indices == k]) if np.any(bin_indices == k) else np.nan for k in range(len(bin_centers))]
# ModÚle théorique
h_vals = np.linspace(0, max_dist, 200)
theoretical = self.nugget + self.sill * (1 - np.exp(-h_vals / self.range_))
# Affichage graphique
with self.out_plot:
clear_output(wait=True)
fig = plt.figure(figsize=(14, 5))
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1.2])
ax0 = fig.add_subplot(gs[0])
ax0.imshow(self.field, cmap='viridis', origin='lower')
ax0.scatter(self.sample_coords[:, 1], self.sample_coords[:, 0], c='r', s=10, label='Ăchantillons')
ax0.scatter([p1[1], p2[1]], [p1[0], p2[0]], c='yellow', s=50, edgecolors='black', label='Points choisis')
ax0.set_title("Champ Gaussien 2D")
ax0.legend(loc='upper right')
ax1 = fig.add_subplot(gs[1])
ax1.scatter(self.dists, self.gamma, alpha=0.3, label='Nuée variographique', color='cornflowerblue')
ax1.plot(h_vals, theoretical, 'r-', label='ModÚle théorique')
ax1.plot(bin_centers, experimental, 'ko-', label='Variogramme expérimental')
for b in bins:
ax1.axvline(b, color='gray', linestyle='--', linewidth=0.5, alpha=0.5)
ax1.scatter(dist_p, gamma_p, c='red', s=100, label='Couple choisi')
ax1.set_xlim(0, 200)
ax1.set_xlabel("Distance h (m)")
ax1.set_ylabel("œ(zᔹ - zⱌ)ÂČ")
ax1.set_title("Nuée variographique")
ax1.legend(loc='upper right')
plt.tight_layout()
plt.show()
# Affichage texte dans l'onglet dédié
with self.out_text:
clear_output(wait=True)
txt = (
f"Point 1 (y,x): {p1}, valeur zâ = {z1:.3f}\n"
f"Point 2 (y,x): {p2}, valeur zâ = {z2:.3f}\n"
f"Distance = ||pâ - pâ|| = {dist_p:.2f}\n"
f"Calcul demi-variance : Îł(h) = œ (zâ - zâ)ÂČ = œ * ({z1:.3f} - {z2:.3f})ÂČ = {gamma_p:.3f}"
)
print(txt)
def on_next_pair(self, b):
self.pair_index += 1
if self.pair_index >= self.max_pairs:
self.pair_index = 0
self.update_plot()
# Lancer l'interface
VariogramWorkshop()
Loading...