Geoestadística con PyKrige: Variogramas e Interpolación Espacial

Creación de variogramas con PyKrige

En este notebook construyes un variograma de principio a fin: generas un campo espacial con parámetros conocidos, calculas el variograma experimental, ajustas modelos y comparas lo que obtienes con PyKrige contra lo que obtienes a mano. Como conoces la verdad (meseta parcial 4, rango 15 km, pepita 0.25), puedes juzgar cada ajuste. Descarga el notebook adjunto a esta lección y ejecútalo celda por celda.

Antes de empezar

Instala las dependencias con pip install pykrige matplotlib. El material se probó con PyKrige 1.7.3, NumPy 2.4, SciPy 1.17 y Matplotlib 3.10. No se descarga nada: los datos se simulan con una semilla fija, así que tus números deben coincidir con los del texto.

Paso 1: simular un campo con estructura conocida

Generamos 150 puntos en un cuadrado de 50 × 50 km. El campo tiene covarianza exponencial y se le suma ruido de medición, que es lo que el variograma verá como pepita:

import numpy as np

def generar_campo(n=150, semilla=42):
    rng = np.random.default_rng(semilla)
    x, y = rng.uniform(0, 50, (2, n))
    d = np.hypot(x[:, None] - x, y[:, None] - y)
    L = np.linalg.cholesky(4.0 * np.exp(-3 * d / 15) + 1e-8 * np.eye(n))
    z = 10 + L @ rng.standard_normal(n) + rng.normal(0, 0.5, n)
    return x.round(3), y.round(3), z.round(2)

x, y, z = generar_campo()
print(z.mean().round(2), z.var().round(2))  # 10.38 3.97

La varianza muestral (3.97) es cercana a la meseta total teórica, 4 + 0.25. Esa coincidencia es la que esperas en un proceso estacionario.

Paso 2: el variograma experimental a mano

Calcular el variograma tú mismo deja claro qué hace la biblioteca y te permite elegir el corte de distancia. Usamos pdist para obtener todas las distancias y semivarianzas por par, y promediamos por intervalos hasta la mitad de la distancia máxima:

from scipy.spatial.distance import pdist

def variograma_experimental(x, y, z, corte, n_lags=10):
    d = pdist(np.c_[x, y])
    g = 0.5 * pdist(z[:, None], "sqeuclidean")
    bordes = np.linspace(0, corte, n_lags + 1)
    k = np.digitize(d, bordes) - 1
    h = np.array([d[k == i].mean() for i in range(n_lags)])
    gam = np.array([g[k == i].mean() for i in range(n_lags)])
    pares = np.array([(k == i).sum() for i in range(n_lags)])
    return h, gam, pares

d_max = pdist(np.c_[x, y]).max()
h, gam, pares = variograma_experimental(x, y, z, d_max / 2)
print(gam.round(2), pares.min())

La semivarianza sube de 1.66 en el primer intervalo a cerca de 4 hacia los 15 km y luego se estabiliza. El intervalo con menos pares tiene 166, muy por encima del mínimo recomendado de 30.

Paso 3: ajustar un modelo y pasarlo a PyKrige

Ajustamos el modelo exponencial de PyKrige con curve_fit y le entregamos los parámetros con la clave "psill":

from scipy.optimize import curve_fit
from pykrige.ok import OrdinaryKriging

def exponencial(h, psill, rango, pepita):
    return psill * (1 - np.exp(-3 * h / rango)) + pepita

p, _ = curve_fit(exponencial, h, gam, p0=[gam.max(), d_max / 4, gam.min()],
                 bounds=([0, 0.1, 0], [10 * gam.max(), d_max / 2, gam.max()]))
print(p.round(2))  # [ 4.02 15.67  0.26]
ok = OrdinaryKriging(x, y, z, variogram_model="exponential",
                     variogram_parameters=dict(psill=p[0], range=p[1], nugget=p[2]))

El ajuste recupera la estructura real casi sin error: meseta parcial 4.02, rango 15.67 km y pepita 0.26.

Paso 4: el ajuste automático y el papel de nlags

Ahora deja que PyKrige ajuste solo, con distintos nlags y modelos. Los parámetros se leen en variogram_model_parameters y los puntos experimentales en lags y semivariance:

for nlags in (6, 10, 20):
    for modelo in ("exponential", "spherical", "gaussian"):
        ok = OrdinaryKriging(x, y, z, variogram_model=modelo, nlags=nlags)
        print(nlags, modelo, ok.variogram_model_parameters.round(2))

Con nlags=10 el exponencial da [4.13, 14.18, 0.0]: buena meseta y buen rango, pero pierde la pepita. Con nlags=6, el esférico asigna 1.84 a la pepita y solo 2.13 a la meseta parcial. La razón es la de la lección de teoría: PyKrige reparte los intervalos hasta la distancia máxima (66 km aquí), así que con pocos lags el primer punto ya promedia pares de 0 a 11 km y la forma cerca del origen se pierde.

Paso 5: graficar experimental contra modelo

El notebook superpone en una sola figura los puntos experimentales (con tamaño proporcional al número de pares), el modelo ajustado a mano y el modelo automático. También muestra enable_plotting=True, la gráfica rápida que trae PyKrige, y get_variogram_points(), que devuelve el modelo evaluado en los lags.

Paso 6: anisotropía

Por último se simula un campo con rango 25 km en dirección este-oeste y 10 km en norte-sur. Calculando el variograma solo con pares en cada dirección, verás dos curvas con rangos distintos. En PyKrige esto se corrige con anisotropy_angle, que rota el sistema de coordenadas (en grados, en sentido antihorario), y anisotropy_scaling, que estira el eje y del sistema rotado. Aquí basta con un ángulo de 0 y una escala de 2.5, el cociente entre el rango mayor y el menor: al estirar el eje norte-sur, las dos direcciones quedan con el mismo rango.

Trampas comunes

• Pasar "sill" creyendo que es la meseta parcial: PyKrige le resta la pepita.
• Confiar en los últimos lags del ajuste automático: tienen pocos pares y suelen bajar o dispararse.
• Usar el modelo gaussiano con pepita cero: el sistema de kriging queda mal condicionado y las estimaciones se vuelven inestables.
• Elegir el modelo por cómo se ve la gráfica sin validar. En la sección de validación verás cómo decidir con números.

Cierre

Repite el notebook con semilla=7: con los mismos parámetros verdaderos, el ajuste automático se va al rango máximo permitido. Ver cuánto cambia el variograma de una muestra a otra es la mejor vacuna contra la confianza excesiva en un solo ajuste.