Skip to article frontmatterSkip to article content
Source
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
import ipywidgets as widgets
from matplotlib.patches import Circle

np.random.seed(1)

# Positions fixes X des points
points_x = {'A': 30, 'B': 45, 'E': 50, 'C': 55, 'D': 70}
labels = ['A', 'B', 'E', 'C', 'D']

def exponential_cov(h, range_, sill, nugget=0):
    return nugget * (h == 0) + (sill - nugget) * np.exp(-h / range_)

def compute_simple_weights(observed_points, range_, sill):
    X = np.array(list(observed_points.values()))
    target = np.array([points_x['E'], 50]).reshape(1,2)  # Point E fixe
    
    if len(X) == 0:
        return {}
    
    d_matrix = cdist(X, X)
    d_vector = cdist(X, target).flatten()
    
    C = exponential_cov(d_matrix, range_, sill)
    c0 = exponential_cov(d_vector, range_, sill)
    
    try:
        weights = np.linalg.solve(C, c0)
    except np.linalg.LinAlgError:
        weights = np.full(len(X), np.nan)
    
    return {label: w for label, w in zip(observed_points.keys(), weights)}

def compute_ordinary_weights(observed_points, range_, sill):
    X = np.array(list(observed_points.values()))
    target = np.array([points_x['E'], 50]).reshape(1,2)  # Point E fixe
    
    if len(X) == 0:
        return {}
    
    d_matrix = cdist(X, X)
    d_vector = cdist(X, target).flatten()
    
    C = exponential_cov(d_matrix, range_, sill)
    c0 = exponential_cov(d_vector, range_, sill)
    
    n = len(X)
    C_ext = np.ones((n+1, n+1))
    C_ext[:n, :n] = C
    C_ext[-1, -1] = 0
    c0_ext = np.ones(n+1)
    c0_ext[:n] = c0
    
    try:
        sol = np.linalg.solve(C_ext, c0_ext)
        weights = sol[:-1]
    except np.linalg.LinAlgError:
        weights = np.full(n, np.nan)
    
    return {label: w for label, w in zip(observed_points.keys(), weights)}

def plot_screening_effect(yA, yB, yC, yD, show_A, show_B, show_C, show_D, range_, sill):
    observed_points = {}
    if show_A:
        observed_points['A'] = np.array([points_x['A'], yA])
    if show_B:
        observed_points['B'] = np.array([points_x['B'], yB])
    if show_C:
        observed_points['C'] = np.array([points_x['C'], yC])
    if show_D:
        observed_points['D'] = np.array([points_x['D'], yD])
    
    w_simple = compute_simple_weights(observed_points, range_, sill)
    w_ordinary = compute_ordinary_weights(observed_points, range_, sill)
    
    fig, axs = plt.subplots(1, 2, figsize=(14,6), sharey=True)
    
    def plot_weights(ax, weights, title):
        ax.set_title(title, fontsize=14)
        ax.set_xlim(15, 75)
        ax.set_ylim(30, 75)
        ax.set_aspect('equal')
        ax.grid(True)
        ax.scatter(points_x['E'], 50, c='black', marker='X', s=150)  # Point fixe E
        
        # Cercle portée autour de E
        ax.add_patch(Circle((points_x['E'], 50), radius=range_, edgecolor='gray', linestyle='--', fill=False))
        
        # Affichage des points observés
        for i, label in enumerate(['A', 'B', 'C', 'D']):
            if label in observed_points:
                w = weights.get(label, 0)
                # normalisation couleur
                norm_w = 0.5
                if weights:
                    w_vals = list(weights.values())
                    if max(w_vals) != min(w_vals):
                        norm_w = (w - min(w_vals)) / (max(w_vals) - min(w_vals))
                color = plt.cm.coolwarm(norm_w)
                ax.scatter(observed_points[label][0], observed_points[label][1], c=[color], s=300, edgecolor='k')
                # Décalage variable du texte pour éviter chevauchement
                x_text = observed_points[label][0] + 0.5
                y_text = observed_points[label][1] + 1.5 + 1.2
                ax.text(
                    x_text,
                    y_text,
                    f"{label}\n{w:.2f}",
                    fontsize=12,
                    bbox=dict(facecolor='white', alpha=0.7, edgecolor='none', pad=1)
                )
            else:
                # points cachés en gris clair petits
                x = points_x[label]
                y = {'A': yA, 'B': yB, 'C': yC, 'D': yD}[label]
                ax.scatter(x, y, c='lightgray', s=60, alpha=0.3)
                ax.text(x+0.3, y+0.3, label, fontsize=12, color='gray')
    
    plot_weights(axs[0], w_simple, "Krigeage simple")
    plot_weights(axs[1], w_ordinary, "Krigeage ordinaire")
    
    plt.tight_layout()
    plt.show()


# Sliders
yA_slider = widgets.FloatSlider(value=50, min=35, max=70, step=1, description='A')
yB_slider = widgets.FloatSlider(value=50, min=35, max=70, step=1, description='B')
yC_slider = widgets.FloatSlider(value=50, min=35, max=70, step=1, description='C')
yD_slider = widgets.FloatSlider(value=50, min=35, max=70, step=1, description='D')

show_A = widgets.Checkbox(value=True, description='Afficher A')
show_B = widgets.Checkbox(value=True, description='Afficher B')
show_C = widgets.Checkbox(value=True, description='Afficher C')
show_D = widgets.Checkbox(value=True, description='Afficher D')

range_slider = widgets.FloatSlider(value=15, min=5, max=40, step=1, description='Portée')
sill_slider = widgets.FloatSlider(value=1.0, min=0.5, max=2.0, step=0.1, description='Sill')

widgets.interact(
    plot_screening_effect,
    yA=yA_slider,
    yB=yB_slider,
    yC=yC_slider,
    yD=yD_slider,
    show_A=show_A,
    show_B=show_B,
    show_C=show_C,
    show_D=show_D,
    range_=range_slider,
    sill=sill_slider,
)
Loading...