Skip to article frontmatterSkip to article content

🎯 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...