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

Teoría de variogramas

El kriging es tan bueno como su variograma. Antes de estimar un solo punto, PyKrige necesita una función que diga cómo se pierde la similitud entre dos mediciones a medida que se alejan. Esta lección explica qué es esa función, cómo se estima con datos y cómo leer sus tres parámetros.

De la intuición a la fórmula

Tratamos la variable medida, por ejemplo la profundidad del agua en un pozo, como una función aleatoria $Z(s)$ que depende de la posición $s$. El semivariograma es la mitad de la diferencia cuadrática esperada entre dos puntos separados por un vector $h$:

$$\gamma(h) = \frac{1}{2}\,E\big[(Z(s+h) - Z(s))^2\big]$$

Para que esta cantidad dependa solo de $h$ y no de dónde estamos, se asume estacionariedad intrínseca: la media es constante y la varianza de los incrementos depende únicamente de la separación. Si además el fenómeno es isotrópico, $\gamma$ depende solo de la distancia $|h|$, no de la dirección.

El variograma experimental

Con $n$ datos hay $n(n-1)/2$ pares. Para cada par se calcula la distancia $d_{ij}$ y la semivarianza $\frac{1}{2}(z_i - z_j)^2$. Luego los pares se agrupan en intervalos de distancia (lags) y se promedia dentro de cada uno:

$$\hat{\gamma}(h_k) = \frac{1}{2N(h_k)}\sum_{(i,j)\in h_k}(z_i - z_j)^2$$

donde $N(h_k)$ es el número de pares del intervalo $k$. Ese conjunto de puntos es el variograma experimental. PyKrige lo calcula al crear el objeto y lo expone en los atributos lags y semivariance. El parámetro nlags (6 por defecto) controla cuántos intervalos se usan.

Pepita, meseta y rango

Los modelos acotados tienen tres parámetros que se leen directamente en la gráfica:

• Pepita (nugget): el valor de $\gamma$ cerca de cero. Representa error de medición y variación a escalas menores que la distancia mínima entre muestras.
• Meseta (sill): el nivel donde la curva se estabiliza. En un proceso estacionario coincide con la varianza total. La diferencia meseta − pepita es la meseta parcial (psill), la parte estructurada de la variación.
• Rango (range): la distancia a partir de la cual dos mediciones ya no están correlacionadas. Más allá del rango, conocer un vecino no ayuda a estimar.

Los modelos de PyKrige

El variograma experimental es ruidoso y no sirve para resolver el sistema de kriging; hace falta una función válida (definida positiva). PyKrige incluye estos modelos, con $p$ = psill, $r$ = rango y $n$ = pepita:

• "spherical": $\gamma(h) = p\left(\frac{3h}{2r} - \frac{h^3}{2r^3}\right) + n$ para $h \le r$, y $p + n$ después. Alcanza la meseta exactamente en $r$.
• "exponential": $\gamma(h) = p\left(1 - e^{-3h/r}\right) + n$. Se acerca a la meseta de forma asintótica; en PyKrige, $r$ es el rango práctico, donde ya se alcanzó el 95% de la meseta parcial.
• "gaussian": $\gamma(h) = p\left(1 - e^{-h^2/(4r/7)^2}\right) + n$. Arranca plano cerca del origen: describe fenómenos muy suaves.
• "linear" y "power": no tienen meseta. Sirven cuando la varianza sigue creciendo con la distancia, señal típica de una tendencia.
• "hole-effect": oscila; aparece en fenómenos con periodicidad espacial.

Ajustar y fijar parámetros

Si no le das parámetros, PyKrige ajusta el modelo al variograma experimental por mínimos cuadrados con pérdida robusta (soft_l1). También puedes fijarlos tú con un diccionario:

import numpy as np
from pykrige.ok import OrdinaryKriging

rng = np.random.default_rng(42)
x, y = rng.uniform(0, 50, (2, 80))
z = 30 + 0.2 * x + rng.normal(0, 1, 80)

auto = OrdinaryKriging(x, y, z, variogram_model="exponential", nlags=10)
print(auto.lags.round(1))
print(auto.semivariance.round(2))
print(auto.variogram_model_parameters.round(2))  # [psill, rango, pepita]

fijo = OrdinaryKriging(x, y, z, variogram_model="exponential",
                       variogram_parameters={"sill": 4.0, "range": 20, "nugget": 0.5})
print(fijo.variogram_model_parameters)  # [3.5, 20, 0.5]

Observa la última línea: con la clave "sill" PyKrige interpreta la meseta total y guarda internamente la meseta parcial (4.0 − 0.5 = 3.5). Si prefieres pasar la parcial, usa la clave "psill".

Trampas comunes

• Lags hasta la distancia máxima: PyKrige reparte los intervalos entre la distancia mínima y la máxima del conjunto. Los últimos intervalos tienen pocos pares y valores erráticos. La regla clásica es confiar solo hasta la mitad de la distancia máxima.
• Pocos pares por lag: con menos de 30 pares, un punto del variograma no es confiable. Reduce nlags si tienes pocos datos.
• Confundir rango práctico y parámetro de escala: otras bibliotecas definen el exponencial con $e^{-h/a}$. El rango de PyKrige equivale a $3a$.
• Un variograma que no se estabiliza: suele indicar tendencia. Quitarla o usar kriging universal da un variograma acotado.
• Anisotropía oculta: si el fenómeno se extiende más en una dirección, el variograma omnidireccional promedia rangos distintos. PyKrige la maneja con anisotropy_scaling y anisotropy_angle.

Cierre

Grafica siempre el variograma experimental junto con el modelo ajustado antes de interpolar. Si los primeros lags no quedan bien representados, ajusta nlags o fija los parámetros a mano: el comportamiento cerca del origen es lo que más pesa en las estimaciones del kriging.