Proyecto Peridotita 4: ¿Cuándo se recuperará del todo Sierra Bermeja? Ajustando una curva matemática al NBR

Introducción

Llevamos tres entradas observando cómo se recupera Sierra Bermeja desde el espacio: primero vimos las imágenes, luego extrajimos la serie temporal completa de NBR, y después comprobamos si la lluvia acelera esa recuperación. Ahora toca una pregunta que todo el mundo hace cuando ve una gráfica que sube: ¿y si seguimos así, cuándo llegaremos arriba del todo?

Para responderla de forma rigurosa —no a ojo— vamos a ajustar una curva matemática a los datos reales de NBR, y usar esa curva para extrapolar cuántos años faltan hasta que el bosque vuelva a su nivel de antes del incendio. Esto es justo lo que hacen los científicos del clima cuando dicen “a este ritmo, el hielo ártico desaparecerá en verano hacia tal año”: no lo adivinan, ajustan un modelo matemático a los datos que tienen y proyectan la tendencia hacia el futuro.

🔬 Pregunta de investigación

A partir de la tendencia real de recuperación del NBR entre 2021 y 2026, ¿en qué año alcanzará Sierra Bermeja el 90%, el 95% y el 99% de su nivel de vegetación anterior al incendio?

Antes de nada: entender el modelo (sin fórmulas raras)

Imagina que dejas una taza de café muy caliente encima de la mesa. Al principio se enfría muy rápido; al cabo de un rato ya casi no cambia de temperatura, aunque nunca llegue a estar exactamente a la temperatura de la habitación. Esto se llama enfriamiento exponencial: la velocidad de cambio es proporcional a lo lejos que estás todavía del valor final.

Un bosque recuperándose de un incendio se comporta de forma parecida, pero al revés: en vez de perder calor, gana vida. Justo después del incendio, cuando todo está muy quemado, la recuperación es rápida (rebrota el matorral, los pinos resineros liberan semillas, aparecen las primeras plantas pioneras). Cuanto más se acerca al nivel de bosque maduro, más lenta se vuelve la recuperación (hace falta tiempo para que crezcan árboles grandes, se recupere el suelo, vuelva la fauna…). Ese patrón —rápido al principio, cada vez más lento, sin llegar nunca del todo a un límite— es exactamente lo que describe una curva exponencial de recuperación: NBR(t)=L−(L−N0)⋅e−ktNBR(t) = L – (L – N_0) \cdot e^{-kt}NBR(t)=L−(L−N0​)⋅e−kt

No hace falta memorizar la fórmula, solo entender qué significa cada letra:

Símbolo	Qué es	En nuestro caso
tt t	Tiempo transcurrido (en años) desde el punto de partida	Años desde el mínimo tras el incendio
N0N_0 N0​	El valor de NBR en el punto de partida (el peor momento)	El NBR más bajo registrado, justo después del fuego
LL L	El valor al que tiende la curva sin llegar nunca del todo (la “meta”)	El NBR del bosque sano, antes del incendio
kk k	Lo rápido que se recupera (más grande = más rápido)	Depende de la vegetación y el terreno
ee e	El número de Euler (≈2,718), la base de los logaritmos naturales	Una constante matemática, la calculadora la conoce

Requisitos previos

• El CSV de NBR generado en el Proyecto Peridotita 2, con todos los puntos de la serie (no solo los de mayo).
• Una cuenta de Google para usar Colab (ver la sección de instrucciones en el Proyecto Peridotita 3 si es la primera vez).

Paso 1: Encontrar el valor “meta” con datos reales, no supuestos

Antes de ajustar nada, necesitamos saber a qué NBR tiende la curva — el valor LL L. Lo más honesto es medirlo, no suponerlo. Y aquí hay una suerte: como la colección de imágenes empieza el 1 de septiembre de 2021 y el incendio no comenzó hasta la noche del 8 de septiembre de 2021 (dato histórico confirmado), ¡las primeras imágenes del CSV son de antes del incendio! Son nuestra línea base real.

(Pega esto en la primera celda de un cuaderno nuevo de Colab)

python

import pandas as pd
import numpy as np

# Sube tu CSV del Proyecto Peridotita 2 al panel de archivos de Colab (icono de carpeta 📁, "Subir")
df = pd.read_csv("nbr-ee-chart.csv")
df["date"] = pd.to_datetime(df["system:time_start"])
df = df.sort_values("date").reset_index(drop=True)

# Imágenes ANTES del incendio (que empezó la noche del 8 de septiembre de 2021)
prefuego = df[df["date"] <= "2021-09-08"]
print("Imágenes pre-incendio encontradas:")
print(prefuego[["date", "NBR"]])

L_real = prefuego["NBR"].mean()
print(f"\nNBR medio ANTES del incendio (nuestra meta L): {L_real:.3f}")

Para verificar este paso: deberían salir 5 imágenes entre el 3 y el 8 de septiembre de 2021, con NBR entre 0,35 y 0,40 aproximadamente — valores altos, típicos de vegetación sana.

Paso 2: Encontrar el punto de partida (el peor momento)

Ahora buscamos el valor más bajo de NBR tras el incendio — el “fondo del pozo” desde el que empieza la recuperación:

python

postfuego = df[df["date"] > "2021-09-08"].copy()
idx_min = postfuego["NBR"].idxmin()
fecha_min = postfuego.loc[idx_min, "date"]
N0_real = postfuego.loc[idx_min, "NBR"]

print(f"Mínimo NBR: {N0_real:.3f}, el {fecha_min.date()}")

Para verificar este paso: el mínimo debería caer entre finales de septiembre y noviembre de 2021, con un NBR muy bajo (por debajo de 0,10) — el momento de mayor severidad del quemado.

Paso 3: Ajustar la curva a todos los datos reales

Aquí es donde entra la magia (bien entendida) del ajuste de curvas. curve_fit, de la librería scipy, hace algo conceptualmente sencillo: prueba muchísimas combinaciones de los parámetros LL L, N0N_0 N0​ y kk k, calcula para cada combinación lo lejos que queda la curva de cada punto real (ese “lejos al cuadrado” se llama error cuadrático), y se queda con la combinación que hace ese error total lo más pequeño posible. Es el mismo principio que cuando ajustáis “a ojo” una recta a una nube de puntos en clase de Física, solo que el ordenador lo hace de forma exacta y con una curva más complicada que una recta.

python

from scipy.optimize import curve_fit

# Tiempo transcurrido en años desde el mínimo (punto de partida t=0)
recuperacion = postfuego[postfuego["date"] >= fecha_min].copy()
recuperacion["t_años"] = (recuperacion["date"] - fecha_min).dt.days / 365.25

t = recuperacion["t_años"].values
y = recuperacion["NBR"].values

# Nuestra función del modelo exponencial de recuperación
def modelo(t, L, N0, k):
    return L - (L - N0) * np.exp(-k * t)

# p0 son valores iniciales "razonables" de partida para que el ajuste no se pierda
p0 = [L_real, N0_real, 0.3]
parametros, covarianza = curve_fit(modelo, t, y, p0=p0, maxfev=10000)
L, N0, k = parametros

print(f"L (meta ajustada) = {L:.4f}")
print(f"N0 (mínimo ajustado) = {N0:.4f}")
print(f"k (velocidad de recuperación) = {k:.4f} por año")

Para verificar este paso: el valor de L que calcula el ordenador debería salir muy parecido al L_real que medisteis directamente en el Paso 1 (alrededor de 0,37). Si sale muy distinto, es señal de que algo no está bien en los datos de entrada.

💡 Esto es lo interesante: el ordenador no sabe cuál era el NBR real antes del incendio — solo ha visto los puntos de la recuperación (después del mínimo) y ha “adivinado” hacia dónde tienden, matemáticamente. Que ese valor adivinado coincida con el que medisteis de verdad en el Paso 1 es una comprobación de que el modelo tiene sentido y no es pura coincidencia.

Paso 4: ¿Qué tan bueno es el ajuste? El coeficiente R²

No toda curva que “más o menos pasa cerca” de los puntos es un buen modelo. Para medir qué tan bien explica la curva los datos reales, se usa el coeficiente de determinación, R², un número entre 0 y 1:

• R² = 1 → la curva pasa exactamente por todos los puntos (ajuste perfecto, raro con datos reales).
• R² = 0 → la curva no explica nada, es como usar simplemente la media de todos los valores.
• R² = 0,65, por ejemplo, significa que el modelo explica el 65% del comportamiento de los datos; el 35% restante se debe a ruido — en nuestro caso, sobre todo a que el NBR sube en primavera y baja en verano cada año (estrés hídrico estival), un vaivén estacional que una curva suave no puede capturar.

python

y_predicho = modelo(t, *parametros)
ss_residual = np.sum((y - y_predicho) ** 2)      # error que NO explica el modelo
ss_total = np.sum((y - np.mean(y)) ** 2)          # variación total de los datos
r2 = 1 - ss_residual / ss_total

print(f"R² = {r2:.3f}  →  el modelo explica el {r2*100:.0f}% del comportamiento observado")

Para verificar este paso: con los datos completos de Sierra Bermeja debería salir un R² alrededor de 0,6-0,7 — un ajuste razonable para datos con tanto ruido estacional, ni sospechosamente perfecto ni demasiado malo.

Paso 5: Extrapolar — ¿cuántos años faltan?

Con la curva ya ajustada, podemos calcular matemáticamente en qué momento se alcanza, por ejemplo, el 90% del camino entre el mínimo y la meta. Despejando tt t de la fórmula del modelo se llega a esta expresión (el álgebra ya está hecha, solo hay que usarla): t=−ln⁡(1−fraccioˊn)kt = \frac{-\ln(1 – \text{fracción})}{k}t=k−ln(1−fraccioˊn)​

python

from datetime import timedelta

for fraccion in [0.90, 0.95, 0.99]:
    t_necesario = -np.log(1 - fraccion) / k
    fecha_alcance = fecha_min + timedelta(days=t_necesario * 365.25)
    print(f"{int(fraccion*100)}% de recuperación: "
          f"{t_necesario:.1f} años desde el mínimo → hacia {fecha_alcance.year}")

Para verificar este paso: cuanto mayor sea el porcentaje objetivo (90% → 95% → 99%), más tarde debería salir la fecha — tiene sentido, porque acercarse del todo al valor final siempre lleva más tiempo que acercarse “casi del todo” (es la propia naturaleza de una curva que nunca toca exactamente su meta).

Resultado con los datos reales de Sierra Bermeja

Nivel de recuperación	Años desde el mínimo (sept. 2021)	Año aproximado
90%	~10 años	≈ 2032
95%	~13-14 años	≈ 2035
99% (recuperación casi completa)	~19-21 años	≈ 2042

La cifra citable del proyecto: a este ritmo, Sierra Bermeja recuperará el 90% de su vegetación anterior al incendio hacia el año 2032, pero no alcanzará un nivel prácticamente completo (99%) hasta bien entrada la década de 2040.

Paso 6: Ver la curva — porque una gráfica dice más que cien números

Con todos los cálculos ya hechos, generamos la imagen que resume todo el post: los puntos reales de NBR, la curva ajustada pasando entre ellos, y la línea de meta (LL L) a la que se acerca sin tocarla nunca.

python

import matplotlib.pyplot as plt

# Una curva suave con muchos puntos intermedios, solo para dibujar bien la línea
t_curva = np.linspace(0, 20, 400)
y_curva = modelo(t_curva, *parametros)

fig, ax = plt.subplots(figsize=(10, 6))

# Los datos reales, un punto por cada imagen de satélite
ax.scatter(t, y, s=18, alpha=0.5, color="#8a8a8a",
           label="Datos reales (cada imagen Sentinel-2)")

# La curva que ha encontrado el ajuste
ax.plot(t_curva, y_curva, color="#c0392b", linewidth=2.5,
        label="Curva ajustada (modelo exponencial)")

# La línea de meta (el valor L al que tiende sin llegar nunca del todo)
ax.axhline(L, color="#2e7d32", linestyle="--", linewidth=1.5,
           label=f"Meta L = {L:.3f} (nivel pre-incendio)")

# El punto de partida, marcado con una X
ax.scatter([0], [N0], color="#c0392b", zorder=5, s=60, marker="X",
           label=f"Mínimo post-incendio (N0={N0:.3f})")

ax.set_xlabel("Años desde el mínimo (18 sept. 2021)")
ax.set_ylabel("Índice NBR")
ax.set_title("Curva de recuperación del NBR en Sierra Bermeja (2021-2026)")
ax.legend(loc="lower right")
ax.grid(alpha=0.25)

plt.tight_layout()
plt.savefig("curva-recuperacion-nbr.png", dpi=150)
plt.show()

Para verificar este paso: la nube de puntos grises debe subir de forma irregular (por el ruido estacional) pero con tendencia clara ascendente, la línea roja debe pasar suavemente “por en medio” de esa nube sin tocar cada punto, y debe acercarse cada vez más despacio a la línea verde discontinua sin llegar nunca a cruzarla del todo — así es precisamente como se ve, visualmente, una curva exponencial de recuperación.

⚠️ Paso 7: La trampa del sobreajuste (muy importante)

Antes de cerrar el post, un aviso de los que hacen mejor científico a cualquiera: cuantos menos datos tengáis, menos os podéis fiar de un ajuste con muchos parámetros libres.

Si repetís exactamente el mismo proceso pero usando solo los 5 valores de NBR de mayo (uno por año, en vez de los más de 600 puntos completos), el ajuste “se vuelve loco”:

python

# Con solo 5 puntos (uno por año, mes de mayo)
t_mayo = np.array([0.7, 1.7, 2.7, 3.7, 4.7])  # años aprox. desde el mínimo
y_mayo = np.array([0.192, 0.171, 0.266, 0.296, 0.294])

parametros_mayo, _ = curve_fit(modelo, t_mayo, y_mayo, p0=[0.37, 0.05, 0.3], maxfev=10000)
print("Con solo 5 puntos, L sale:", round(parametros_mayo[0], 3))

Con solo 5 puntos, el valor de LL L que “inventa” el ajuste sale alrededor de 0,67 — ¡casi el doble del valor real medido (0,37)! ¿Por qué pasa esto?

Piensa en el modelo como si intentara adivinar la forma completa de una montaña a partir de solo 5 fotos, todas tomadas el mismo día del año. Con tan pocos datos, hay muchísimas curvas distintas que pasan razonablemente cerca de esos 5 puntos, y el ordenador no tiene forma de saber cuál de ellas es la “de verdad” — así que a veces elige una que se dispara hacia valores absurdos fuera del rango de lo observado. Esto se llama sobreajuste (overfitting): el modelo tiene tantos grados de libertad comparado con los datos que tiene, que puede “memorizar” esos pocos puntos sin aprender el patrón real.

La solución que usamos en el Paso 3 fue justo evitar este problema: en vez de dejar que el ajuste “inventara” el valor de LL L con pocos datos, usamos todos los puntos disponibles (más de 600, no solo 5), lo que da al modelo información suficiente para no dispararse. Esta es una de las razones más importantes por las que, en el Proyecto Peridotita 2, os pedimos sacar la serie completa de NBR y no solo las fechas de mayo.

Para verificar este paso: si repetís el ajuste con pocos puntos y os sale un valor de L mayor que el máximo NBR que jamás habéis observado en todo el estudio, es la señal inequívoca de que estáis ante un sobreajuste, no ante un resultado real.

Conclusión

Ajustar una curva matemática a datos reales no es solo “dibujar una línea que quede bonita” — es un proceso con reglas claras (mínimos cuadrados), formas de medir si el resultado es fiable (R²) y trampas conocidas en las que es fácil caer (el sobreajuste con pocos datos). Con la serie completa de NBR de Sierra Bermeja, el modelo exponencial de recuperación no solo ajusta razonablemente bien (R² ≈ 0,65), sino que “descubre” por sí solo, a partir únicamente de la tendencia de recuperación, casi el mismo valor de vegetación sana que medisteis directamente en las imágenes de antes del incendio — la mejor prueba de que el modelo tiene sentido físico, no es solo un truco matemático.

La respuesta a la pregunta con la que empezaba esta entrada es, por tanto: hacia 2032 para una recuperación notable, pero no antes de 2042 para una recuperación prácticamente completa — un plazo que encaja con lo que sabemos de los bosques mediterráneos tras un gran incendio, y un buen argumento científico para seguir vigilando Sierra Bermeja desde el espacio muchos años más.

---

Etiquetas: AIDARAC, Climate Detectives, Ciencia de Datos, curve fitting, ESA, ies monterroso, Málaga, Matemáticas, montesteam, nbr, Proyecto Peridotita, Python, Regeneración Forestal, sierra bermeja, Teledetección