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

Proyecto final: mapa del nivel freático con kriging

En este notebook aplicas todo el curso a un caso completo. Una autoridad del agua tiene 150 pozos de monitoreo repartidos en una zona de 50 × 50 km y necesita un mapa continuo de la profundidad del nivel freático, junto con las zonas donde el agua está tan profunda que bombear deja de ser rentable. Vas a explorar los datos, detectar una tendencia, modelar el variograma de los residuos, comparar kriging ordinario, kriging universal y distancia inversa con validación cruzada, y entregar un mapa con su incertidumbre. Descarga el notebook adjunto a esta lección.

Antes de empezar

Necesitas pip install pykrige matplotlib scipy. El material se probó con PyKrige 1.7.3, NumPy 2.4, SciPy 1.17 y Matplotlib 3.10. Los datos se generan con semilla fija, así que tus números deben coincidir con los del notebook. La función generadora viene en la primera celda:

import numpy as np

def generar_pozos(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 = 25 + 0.12 * x + 0.08 * y + 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_pozos()  # x, y en km; z en metros de profundidad

Paso 1: exploración y tendencia

Empieza con un mapa de puntos coloreado por profundidad y una regresión lineal de $z$ contra las coordenadas. El agua es más profunda hacia el noreste: el plano ajustado $\hat{z} = 27.25 + 0.084x + 0.038y$ explica cerca del 31% de la varianza. Esa tendencia es la razón por la que el variograma de los datos crudos no se estabiliza bien.

Paso 2: variograma de los residuos

Calcula el variograma experimental de los residuos del plano hasta la mitad de la distancia máxima y ajusta un modelo exponencial. En los datos crudos el ajuste da un rango de unos 20 km, inflado por la tendencia; en los residuos, alrededor de 10 km con meseta parcial de 3.7.

A = np.c_[np.ones_like(x), x, y]
beta, *_ = np.linalg.lstsq(A, z, rcond=None)
residuo = z - A @ beta
print(beta.round(3))  # [27.254  0.084  0.038]

Paso 3: kriging universal con el variograma correcto

UniversalKriging estima la tendencia y el campo residual a la vez con drift_terms=["regional_linear"]. Hay un detalle clave: si no le das parámetros, PyKrige ajusta el variograma sobre los datos crudos, con tendencia incluida. Por eso conviene pasarle el variograma de los residuos que ajustaste en el paso anterior:

from pykrige.uk import UniversalKriging

vp = dict(psill=3.72, range=10.32, nugget=0.0)
uk = UniversalKriging(x, y, z, variogram_model="exponential",
                      variogram_parameters=vp, drift_terms=["regional_linear"])
malla = np.arange(0, 50.5, 1.0)
z_est, varianza = uk.execute("grid", malla, malla)

Paso 4: validación cruzada y elección del modelo

Con leave-one-out compara cinco alternativas. Los resultados que debes reproducir (RMSE en metros):

• Predecir la media: 2.22
• Distancia inversa con exponente 2: 1.61
• Kriging ordinario con variograma automático: 1.52
• Kriging universal con variograma automático (datos crudos): 1.52
• Kriging universal con variograma de los residuos: 1.43, con MSDR de 0.98

El kriging universal solo gana cuando su variograma describe los residuos. Con el variograma de los datos crudos, la tendencia se cuenta dos veces y la ventaja desaparece. El LOO del kriging universal tarda más de un minuto: es normal, porque resuelve 150 sistemas con la deriva incluida.

Paso 5: mapa, incertidumbre y zonas de decisión

Grafica tres paneles: la profundidad estimada, la desviación estándar de kriging (raíz de la varianza) con los pozos encima, y la probabilidad de que la profundidad supere un umbral $u$. Bajo el supuesto gaussiano:

$$P(Z(s) > u) = 1 - \Phi\left(\frac{u - \hat{z}(s)}{\sigma(s)}\right)$$

La desviación estándar va de unos 0.35 m junto a los pozos a casi 2 m en las esquinas sin datos. Ese contraste muestra dónde perforar el próximo pozo de monitoreo.

Trampas comunes

• Interpretar el mapa en las esquinas como si fuera tan confiable como el centro.
• Usar varianza como si fuera desviación estándar en el mapa de probabilidad.
• Ajustar el variograma con todos los datos y luego validar: reporta que el variograma no se reajustó en cada pliegue.
• Olvidar que execute("grid", ...) devuelve matrices con forma (len(y), len(x)): las filas son y.

Cierre

Cuando termines, tendrás todos los elementos para la presentación del proyecto en la siguiente lección: el modelo elegido y por qué, las métricas de validación y tres mapas que un tomador de decisiones puede leer.