El granizo fino pegaba contra el vidrio de la ventana mientras mi mate se enfriaba al lado del teclado. Afuera los techos de Valdivia brillaban con esa humedad espesa de siempre, y en mi pantalla una escena satelital mostraba los cerros de la cordillera de la Costa cubiertos por una especie de velo blanquecino y lechoso que nada tenía que ver con la realidad.
Llevaba un par de horas tratando de calcular un índice simple para ver el dosel cerca del río. Pero cada gráfico arrojaba números rarísimos: zonas de bosque nativo espeso marcaban como si fueran matorrales secos. Al principio creí que había importado mal las coordenadas geográficas, pero no; el problema era que estaba alimentando mis operaciones con números crudos directamente salidos del sensor.
Del número digital al tope de la atmósfera
Cuando descargas una escena de Landsat Collection 2 Level-1, las bandas vienen en archivos TIF que guardan valores digitales enteros, conocidos como DN (digital numbers). Estos números no representan luz real reflejada por las hojas ni el suelo, sino cuentas electrónicas que registró el satélite a cientos de kilómetros de altura. Es como tener los voltajes de una cámara fotográfica en vez de los colores verdaderos.
Durante un sábado lluvioso a mediados de julio me puse a buscar tutoriales y caí en un mar de frustración. Casi todo lo que encontraba en foros usaba paquetes viejos y descontinuados que ya no corren bien en versiones actuales. Yo no sé programar de forma profesional; aprendo a tropezones los fines de semana. Descubrí que la librería moderna R para esto es terra, mucho más rápida y limpia para manejar datos espaciales sin colapsar la memoria del computador.
Para convertir esos números digitales en reflectancia al tope de la atmósfera (lo que en los manuales llaman TOA, por Top of Atmosphere), necesitamos los parámetros que vienen en el archivo de texto de metadatos, ese que termina en _MTL.txt. Ahí la NASA y el USGS entregan dos valores fijos para cada banda multiespectral (que tienen una resolución espacial de 30 metros):
- Factor multiplicativo:
0.00002(REFLECTANCE_MULT_BAND) - Factor aditivo:
-0.1(REFLECTANCE_ADD_BAND)
Además, entra en juego el ángulo de elevación solar con el que se tomó la imagen. Si el sol estaba más bajo en el horizonte —típico en tomas invernales del sur— la luz recorrió un camino más largo y disperso. Para llevar las bandas a reflectancia TOA corregida por el sol, la fórmula básica que metí en mi consola fue:
library(terra)
# Cargar la banda roja (por ejemplo, banda 4 en Landsat 8)
b4 <- rast("LC08_..._B4.TIF")
# Factores de escala oficiales de metadatos
mult <- 0.00002
add <- -0.1
sun_elevation <- 34.25 # ángulo en grados según el archivo MTL
# Conversión matemática considerando el ángulo solar
sin_sun <- sin(sun_elevation * pi / 180)
toa_b4 <- ((b4 * mult) + add) / sin_sun
Ver que la imagen pasaba de miles a valores entre cero y uno fue un alivio inmenso. Pensar que un píxel con valor negativo en reflectancia no era un misterio cuántico, sino solo un ángulo solar mal interpretado en mi código, me sacó una sonrisa culpable frente a la taza vacía.
El objeto oscuro: quitando la bruma de la pantalla
Tener reflectancia al tope de la atmósfera ayuda, pero no quita el problema del aire. La luz solar viaja dos veces por la atmósfera: baja hasta la Tierra y rebota hacia el sensor. Los gases y la humedad dispersan radiación, sobre todo en las longitudes de onda cortas como el azul y el rojo. Por eso las imágenes sin corregir se ven desteñidas, como si las miraras a través de un vidrio con polvo.
Hacia fines de agosto, después de tres semanas de pruebas con scripts que daban error tras error, me crucé con la técnica de sustracción del objeto oscuro (DOS1, o Dark Object Subtraction). La lógica es hermosa en su sencillez: en cualquier escena grande existen áreas que prácticamente no reflejan luz visible. Por ejemplo, aguas muy profundas y claras, o las sombras profundas que proyectan las laderas empinadas de los cerros en un día despejado.
Si un cuerpo de agua profunda debería reflejar un valor cercano a cero pero el sensor registra 0.03, asumimos que ese 0.03 es pura interferencia atmosférica repartida de forma homogénea sobre la zona. La corrección consiste sencillamente en restarle ese valor base a toda la banda:
# Buscamos un valor mínimo representativo (evitando ceros de bordes)
# Usar un cuantil bajo suele ser más seguro que el mínimo absoluto
val_oscuro <- quantile(toa_b4[toa_b4 > 0], probs = 0.01, na.rm = TRUE)
# Sustracción del objeto oscuro
boa_b4 <- toa_b4 - val_oscuro
# Truncamos valores negativos residuales a cero
boa_b4 <- clamp(boa_b4, lower = 0, upper = 1)
El repiqueteo de la lluvia en el techo de lata mientras esperaba que la consola terminara de escribir el raster corregido me hacía dudar de si la máquina no se había pegado. Pero cuando el comando terminó y tracé el histograma de la banda roja corregida, la curva por fin se había despegado del piso y mostraba contrastes nítidos entre los pastizales ribereños y las copas tupidas de los árboles.
Una mirada honesta sobre cuándo vale la pena este trabajo
Aquí es donde las cosas toman un rumbo distinto al de la mayoría de los manuales universitarios. Después de pasarme fines de semana enteros batallando con archivos MTL y fórmulas de sustracción, empecé a notar algo curioso al comparar imágenes de distintas temporadas.
Si tu único objetivo es monitorear cambios relativos en la vegetación a lo largo del tiempo dentro de un mismo sector, aplicar correcciones atmosféricas estándar puede resultar innecesario, e incluso perjudicial. Al normalizar y forzar sustracciones lineales como DOS1, a menudo se introducen pequeñas distorsiones artificiales en píxeles mixtos o zonas de sombra parcial. Los índices de vegetación basados en proporciones de bandas (como los cocientes infrarrojo/rojo) ya cancelan de por sí una parte importante del efecto atmosférico multiplicativo.
Me di cuenta de esto comparando mis mapas de cobertura: a veces los índices calculados directamente sobre reflectancia TOA eran mucho más estables entre fechas que los que pasaban por sustracciones manuales, donde elegir mal el píxel oscuro terminaba falseando la tendencia estacional. Si quieres profundizar en cómo equilibrar datos entre temporadas, hace poco escribí sobre cómo normalizar valores de un raster en R para comparar diferentes fechas sin caer en esos sesgos de cálculo.
También me pasó tropezar con la presencia de bruma residual que parecía sombra. De hecho, en lo que aprendí al quitar nubes de imágenes Landsat en RStudio comenté cómo una mala máscara de nubes te arruina por completo la búsqueda del objeto oscuro.
Cómo dejar armado un script para no repetir el proceso
Un domingo por la tarde hace pocas semanas, cansada de volver a escribir las mismas líneas cada vez que bajaba una imagen nueva de la costa valdiviana, decidí ordenar todo en una función pequeña. La gracia de R es que una vez que desenredas la lógica, puedes encapsularla para que trabaje por ti mientras te tomas otro té:
corregir_banda <- function(ruta_tif, mult = 0.00002, add = -0.1, sun_elev) {
r <- rast(ruta_tif)
sin_s <- sin(sun_elev * pi / 180)
toa <- ((r * mult) + add) / sin_s
# Estimación del objeto oscuro en el percentil 1%
oscuro <- as.numeric(quantile(toa[toa > 0], probs = 0.01, na.rm = TRUE))
boa <- clamp(toa - oscuro, lower = 0, upper = 1)
return(boa)
}
No pretendo que esto sea un algoritmo de nivel de investigación aeroespacial ni un reemplazo para productos procesados en laboratorio. Pero para quienes aprendemos solas frente a la pantalla, entendiendo de a poco cómo interactúa la radiación con el vapor de agua y las ramas de los coigües, aplicar esta corrección básica con nuestras propias manos le quita todo el misterio a la teledetección.