Del pulso láser a la imagen de dosel

Cargar el CHM en el notebook

En el notebook (§0), tras el setup (que hace import tda_utils as T y define DATA), cargamos el recorte. load_chm usa rasterio y deja los nodata como NaN:

chm_full, meta = T.load_chm(DATA / "zambezi_chm_crop.tif")
print("recorte cargado:", chm_full.shape,
      "| válido:", round(float(np.isfinite(chm_full).mean()), 3),
      "| altura máx:", round(float(np.nanmax(chm_full)), 1), "m")

Salida:

recorte cargado: (1200, 1200) | válido: 0.958 | altura máx: 30.0 m

chm_full ya es una imagen (matriz de alturas); el ~4 % de celdas NaN (agua/huecos) son los píxeles blancos de la figura. Todo lo que sigue la trata como la función escalar f.

El pipeline conceptual

Significancia: ¿copa o ruido?

¿Qué es la distancia bottleneck?

Bottleneck en la práctica: el código

import gudhi as gd

def dgm(p):                                  # diagrama H0 como pares (nacimiento, muerte)
    f = np.isfinite(p.saddle)
    return np.column_stack([-p.apex[f], -p.saddle[f]])

tile_ruido = tile + rng.normal(0, 0.30, tile.shape)     # ruido del sensor, σ ≈ 0.30 m
d_B = gd.bottleneck_distance(dgm(pers), dgm(T.superlevel_persistence(tile_ruido)))

Esto es una iteración del bootstrap de \varepsilon_\text{min}: repetirlo B veces da la banda c_n (cuantil 95 %), y por estabilidad \varepsilon_\text{min} = 2\,c_n.

Bottleneck en la práctica: el resultado

Una sola perturbación mueve el diagrama d_B \approx 0.76 m. Repetido B veces, el cuantil 95 % de esas distancias es c_n \approx 0.92 m, y \varepsilon_\text{min} = 2\,c_n \approx 1.85 m.

La persistencia jerarquiza: distribución de prominencias

Sobre el CHM público de Zambezi: un pico enorme de prominencias diminutas (ruido del sensor) y una cola larga de copas reales. El umbral \varepsilon_\text{min} corta entre ambos.

Otra dimensión: H₁ = los claros del dosel

Todo lo anterior fue H₀ (copas = componentes conexas). La otra dimensión de la homología, H₁, cuenta agujeros — y en un CHM un agujero es un claro del bosque: un hueco de canopia baja encerrado por dosel alto (caída de un árbol, perturbación, entrada de luz).

Al bajar el “nivel del mar” sobre el CHM, un claro:

• nace cuando el dosel alto que lo rodea se cierra en un anillo (al nivel del paso más bajo del anillo);
• muere cuando el nivel llega al fondo del hueco y lo rellena;
• su persistencia H₁ = anillo − fondo = la profundidad del claro.

Los claros del dosel: el código

Toda la detección es una función: canopy_gaps corre H₁ y aplica el filtro Brokaw con tres umbrales ecológicos (los huecos sobre nodata/agua se descartan).

G = T.canopy_gaps(tile, ring_min=8.0, floor_max=2.0, depth_min=3.0)
print("generadores H₁:", len(G.all_depth), "| claros:", len(G.rows))

• ring_min — dosel circundante mínimo (m): un hueco en dosel alto.
• floor_max — fondo máximo (m): que el hueco sea realmente bajo.
• depth_min — profundidad mínima (m): un claro real, no ruido.

Salida:

generadores H₁: 22,503  |  claros significativos: 15

Un inventario de la porosidad del dosel

Los 15 claros sobre el CHM (izquierda, color = profundidad) y en el diagrama H₁ (derecha; gris = ruido, magenta = significativos). Es el reverso del inventario de copas: no las cimas, sino los huecos —que regulan luz, regeneración y disturbio.

Escala: ¿de qué tamaño es una “unidad”?

El TDA se corre por ventanas (tiles). El tamaño de ventana define qué estructura ves:

• Muy chica → pocas copas por tile, estimaciones inestables.
• Muy grande → mezclas rodales distintos; copas se funden en complejos amplios.
• Escala de trabajo → grande para tener suficientes copas (estimación estable), pero no tanto que mezcle rodales.

Qué obtenemos por tile

Del diagrama de persistencia de cada ventana extraemos descriptores que sí ven la estructura:

Todo sale de una función; para el tile de trabajo:

d = T.chm_descriptors(tile, eps)     # los 6 descriptores de la ventana

Salida:

n_sig=2841, prom_media=3.11, prom_p95=5.57, pers_entropy=7.88, altura_p95=24.79, cobertura=1.00

Cómo agrupamos: k-means sobre los descriptores

Regionalización: el código

regionalize encadena todo el bloque sobre las ventanas del recorte: (1) parte el CHM en teselas y calcula sus descriptores (chm_descriptors), (2) los estandariza, (3) agrupa las ventanas con k-means en k regiones y (4) devuelve el mapa de etiquetas (una región por ventana).

R = T.regionalize(chm_full, tile=100, k=4)    # descriptores/ventana → k-means
big = T.upsample_labels(R["label_map"], R["tile"], chm_full.shape)   # etiquetas → píxel

Salida:

ventanas clasificadas: 140  |  ε_min = 1.87 m  |  regiones (k) = 4

Cada etiqueta de región es un color en el mapa: R["label_map"] guarda la región de cada ventana; R["feat"], sus descriptores.

La escala es una elección ecológica

La misma partición sobre Zambezi a 50 / 100 / 200 m. Cada color es una de las 4 regiones que agrupa k-means por sus descriptores TDA: mismo color = estructura de dosel parecida (las etiquetas son categóricas, no un orden; importa qué ventanas comparten color). Ventanas chicas resuelven copas y rodales; grandes, el paisaje.

¿Aporta la topología sobre la altura?

Agrupamos las mismas ventanas de Zambezi de dos maneras: solo con altura (p95, cobertura) y añadiendo los descriptores TDA. Las particiones difieren (ARI ≈ 0.25).

¿Qué CHM le damos al TDA?

Hasta aquí asumimos que el CHM era confiable. Pero el CHM no es un dato bruto: es una decisión de agregación — y como el TDA mide forma, esa decisión cambia la topología que medimos.

Un CHM sale de la nube de puntos aplicando un agregador por celda (máximo, media, p95…). Lo habitual es el máximo duro (f_{\max}): la altura del retorno más alto. Pero hay una pregunta metodológica que casi nadie hace:

La elección del agregador de altura no es preprocesamiento neutral: condiciona la topología que vas a medir. La pregunta no es qué tan bien estima la altura, sino qué objeto produce para el TDA.

En doseles bajos y uniformes esto puede volverse un problema real: sobre un CHM de máximo duro la topología puede resultar frágil e irreproducible. Ese es el mecanismo que motiva usar un agregador suave como LSE.

El agregador LSE: interpola entre máximo y media

Una envolvente superior suave por log-sum-exp (LSE), por celda:

f_\tau(c) = \tau \,\log\!\Big(\tfrac{1}{n_c}\textstyle\sum_{i\in c} e^{h_i/\tau}\Big)

• Familia continua entre dos extremos insatisfactorios: \tau\to 0 = máximo; \tau\to\infty = media. Siempre \bar h \le f_\tau \le \max.
• A \tau intermedio pesa los retornos altos ponderados por toda la distribución de la celda.

def chm_lse(Z, tau):                      # Z[c] = alturas de los retornos de la celda c
    m = Z.max(axis=2, keepdims=True)      # estabilización numérica
    return tau*np.log(np.mean(np.exp((Z - m)/tau), axis=2)) + m[..., 0]

El suavizado gaussiano: cómo se hace

El suavizado va dentro del estimador, no sobre un CHM ya hecho. Por celda partimos de la suma de pesos exponenciales S_\tau y el conteo n:

S_\tau(c)=\textstyle\sum_{i\in c} e^{h_i/\tau}, \qquad n(c)=\#\{\text{retornos en } c\}

Suavizamos ambos por separado con el mismo kernel gaussiano (soporte σ = 1 m) y tomamos el cociente:

f_\tau^{\text{(prod)}}(c) = \tau \log\!\Big[\frac{G_\sigma\{S_\tau\}(c)}{G_\sigma\{n\}(c)}\Big]

Es una convolución normalizada por conteo: el análogo espacial del 1/n del LSE por celda.

# por celda: S = Σ exp(h_i/τ)  y  n = nº de retornos
f_lse = T.lse_surface(S, n, tau, sigma=1.0)   # suaviza S y n (mismo kernel) → cociente → log

El suavizado en acción

A 0.5 m la densidad de retornos varía: muchas celdas tienen 0 o 1 retorno, y ahí el LSE por celda falla:

• n = 0 → hueco (indefinido) → topología espuria en los bordes del hueco.
• n = 1 → el LSE = ese único retorno = el máximo → sin robustez, justo donde más falta.

El suavizado deja que cada celda tome retornos de su vecindad (~1 m): rellena los huecos y le da suficientes retornos para que la reducción de varianza del LSE de verdad actúe.

¿Cómo se evalúa un agregador? Tres ejes, no uno

No basta con “¿es robusto?”. Un buen agregador para TDA debería cumplir tres ejes independientes:

Reproducibilidad split-half

Lo proponemos como un experimento

No es un test con receta universal: lo proponemos como un eje explícito para calibrar el agregador según la reproducibilidad de la topología que produce.

Distancia entre diagramas

¿Cómo se mide que “dos diagramas se parecen”? Emparejando sus puntos (y con la diagonal para los que sobran):

• Bottleneck d_B = el peor emparejamiento: cuánto hay que mover el punto que más se resiste. Es el que aparece en el teorema de estabilidad de la Parte 1.
• Wasserstein-1 W_1 = la suma de todos los movimientos: sensible a muchos cambios pequeños, no solo al peor. Usamos su versión sobre los tiempos de vida (prominencias).

Bottleneck vs. Wasserstein

Mismo emparejamiento entre dos diagramas: el bottleneck solo cuenta el peor movimiento (rojo); el Wasserstein-1 suma todos (grosor ∝ longitud). Por eso W_1 siente muchos cambios pequeños que el bottleneck ignora — y por eso es el que usamos para la reproducibilidad.

Qué mide la reproducibilidad

Honestidad sobre el alcance del diagnóstico:

• No empareja los pares (nacimiento, muerte) de los diagramas. Compara la distribución 1-D de prominencias \pi_i = \text{cima}-\text{silla} con Wasserstein-1 normalizada por la mediana. Es reproducibilidad de la distribución de prominencias H₀, no de la topología completa — dos doseles con las mismas prominencias pero copas en otros sitios pasarían el test.
• Es reproducibilidad de muestreo al 50 % de densidad: un diagnóstico interno, no una réplica de un segundo vuelo (que cambiaría trayectoria, ángulo de escaneo, oclusión…).
• Se fortalece repitiendo el split 50/50 R veces —cada vez un reparto aleatorio distinto de todos los retornos (misma densidad, no menos)— y tomando la mediana de las distancias.

Un laboratorio sintético

El raster público no tiene retornos, así que no podemos calibrar τ sobre él. Construimos una nube de puntos sintética sobre un dosel cuyo número real de copas conocemos:

# 1) dosel verdadero = máximo de gaussianas → copas CONOCIDAS (ground truth)
f_true = canopy(centers, amps, widths, base)

# 2) nube de retornos por celda: tope − penetración + ruido + pico espurio ocasional
def simulate_returns(f, n=12, sigma=0.15, pen=0.8, p_out=0.03, out=5.0):
    base = f[..., None]
    return (base - rng.exponential(pen, (N, N, n))        # penetración
                 + rng.normal(0, sigma, (N, N, n))        # ruido del sensor
                 + (rng.random((N, N, n)) < p_out) * rng.uniform(1, out, (N, N, n)))  # espurios
Z = simulate_returns(f_true)

El máximo infla; el LSE recupera

fmax = chm_max(Z)            # máximo duro
flse = chm_lse(Z, tau=4.0)   # LSE, f_τ (log-sum-exp)
fig, axs = plt.subplots(1, 3, figsize=(11, 3.5))
for ax, img, ttl in [(axs[0], f_true, f"verdad — {G_true} copas"),
                     (axs[1], fmax, f"máximo duro — {n_significant(fmax)} copas"),
                     (axs[2], flse, f"LSE τ=4 — {n_significant(flse)} copas")]:
    im = ax.imshow(img, cmap="viridis", vmin=f_true.min(), vmax=f_true.max())
    ax.set_title(ttl); ax.axis("off")
fig.colorbar(im, ax=axs, shrink=0.7, label="altura (m)")
plt.show()

El dosel sintético tiene 27 copas reales. El máximo duro reporta cientos (ruido de un solo retorno espurio); el LSE las recupera. El agregador decide si el TDA ve estructura o ruido.

“Copas” = generadores H₀ con prominencia > 1.5 m (corte fijo, por ser sintético; en Zambezi el umbral es el \varepsilon_\text{min} del bootstrap).

Calibrar τ: el código

taus = np.array([0.5, 1, 2, 4, 8, 16, 32])
ncrowns, hbias = [], []
for tau in taus:
    c = chm_lse(Z, tau)
    ncrowns.append(n_significant(c))          # copas recuperadas
    hbias.append(float(np.mean(c - f_true)))  # sesgo de altura
# graficar ncrowns y hbias frente a τ, con la línea de la verdad G_true

Barremos τ y medimos, contra la verdad, el número de copas y el sesgo de altura. Al subir τ el conteo converge a la verdad (el codo → τ* ≈ 4); un τ demasiado grande sesga la altura (el LSE tiende a la media).

Calibrar τ: el resultado

El mecanismo: el código

cB = poisson_disk(24, 12, gen)                     # copas más separadas
f_deep = canopy(cB, ..., base=2.0)                 # dosel alto y variado
pA, pB = prominences(f_true), prominences(f_deep)  # sillas de cada dosel
# comparar los histogramas de pA (superficiales) y pB (profundas)

¿Cuándo importa τ? Depende de la profundidad de las sillas entre copas:

• Dosel bajo y uniforme → sillas superficiales, cerca del piso de ruido → un mal agregador las cruza → topología frágil → hay que calibrar τ.
• Dosel alto y variado → sillas profundas → ningún ruido razonable las borra → τ importa poco.

El mecanismo: el resultado

Del inventario topológico al carbono

De altura a carbono

Para Zambezi usamos la alometría calibrada en campo de Stringer et al. 2015 —la misma del análisis de carbono del proyecto—: una relación stand-level directa de altura a carbono, sin densidad de tallos.

C_\text{AGB} = 59.6\,H^{0.562}\ \ [\text{Mg C/ha}] \;\xrightarrow{\times\,44/12}\; \text{CO}_2

• Ajustada a los estratos de campo de Stringer (p. ej. 2.5 m → 99 Mg C/ha; 22.5 m → 341 Mg C/ha).
• Alternativa por árbol (Komiyama 2005: H\to D\to AGB) cuando sí se conoce la densidad de tallos del sitio.

Cada sitio necesita su propia calibración alométrica; el pipeline TDA es idéntico. Corre en vivo en el notebook.

La cifra de carbono, con honestidad

El agregador también condiciona el carbono

Como el máximo duro siempre es \ge el log-sum-exp (\text{máx} \ge \text{LSE}), el máximo tiende a sobreestimar la altura y, con ella, el carbono estimado tile a tile.

• El sesgo es sistemático (siempre en el mismo sentido), así que no se promedia a cero al agregar sobre muchas hectáreas.
• Su magnitud depende de \tau y del sitio; cuantificarla sobre un CHM real es un ejercicio abierto.

Lo que corremos en vivo

Un solo notebook autocontenido (notebook/tda_forestal.qmd; en Google Colab, un clic):

1. Descargar el CHM de Zambezi; elegir una ventana y mirarlo como imagen.
2. Superniveles con gudhi → diagrama H₀ y distribución de prominencias.
3. Estimar \varepsilon_\text{min} por bootstrap → copas significativas.
4. H₁ → los claros del dosel (agujeros) con su profundidad.
5. Descriptores por ventana → regiones, sensibilidad a la escala y clásico vs. TDA.
6. Altura → carbono con la cadena alométrica, e IC topológico por bootstrap.
7. Laboratorio de τ (nube sintética): la evidencia del Bloque 4, con ground truth.

Trasladar el método a otros bosques

• La metodología es idéntica; solo cambian los datos.
• En doseles bajos y uniformes, las sillas de copa son superficiales → la topología es más sensible a cómo se construye el CHM (de ahí el LSE calibrado).
• En doseles altos y variados (como Zambezi), la altura ya lleva mucha señal → el marco sigue funcionando, con menos dependencia del agregador.

Los límites del método

Para no sobre-vender —qué no hace esto:

• Con solo un raster no puedes auditar el agregador original: la calibración de τ necesita los retornos.
• \varepsilon_\text{min} depende del modelo de ruido vertical que asumas (\sigma).
• H₀ detecta prominencias, no “árboles individuales” garantizados: una copa ancha o dos fundidas pueden contar distinto.
• H₁ detecta huecos cerrados; los claros abiertos al borde del tile requieren cuidado.
• La regionalización depende del tamaño de ventana y de las variables incluidas.

Para practicar

El eslabón que Zambezi no permite (raster sin retornos) sí se puede hacer donde hay nube de puntos:

1 · En el laboratorio sintético (público, reproducible). Construye el LSE (f_\tau) desde los retornos y corre H₀ + \varepsilon_\text{min} + H₁ sobre él; compáralo con el máximo duro.

flse = chm_lse(Z, tau=4.0)     # LSE, f_τ, desde los retornos
fmax = chm_max(Z)              # máximo duro
# corre significancia / claros sobre cada uno y compara:
# ¿cómo cambian nº de copas, ε_min y los claros?

2 · Reto (tu propia nube). Si tienes un LAS, reconstruye el CHM con chm_lse y corre el pipeline completo, comparándolo con el máximo duro sobre tu sitio.

Recursos

• Notebook (todo en uno): notebook/tda_forestal.qmd — de la descarga al carbono, más el laboratorio sintético de τ. Corre en Google Colab con un clic.
• Funciones: scripts/tda_utils.py (persistencia, umbral, descriptores, regiones, claros).
• Datos: recorte público de Zambezi en data/zambezi/ (NASA CMS #1357).

Lecturas mínimas

• Edelsbrunner & Harer (2010), Computational Topology: An Introduction.
• Cohen-Steiner, Edelsbrunner & Harer (2007), Stability of persistence diagrams.
• Fasy et al. (2014), Confidence sets for persistence diagrams, Ann. Statist.
• GUDHI (2023), manual de referencia — gudhi.inria.fr.