Un raster de elevación sin proyectar entrega pendientes de 0.001 grados en un cerro que en la vida real te deja sin aliento a la primera subida. Ese es el problema más común al armar un mapa de pendientes en RStudio: los datos satelitales de elevación casi nunca vienen listos para este cálculo, y quien recién arranca en análisis espacial suele saltarse ese paso sin darse cuenta. Sigo la teledetección con R los fines de semana, sobre todo mirando el sector boscoso del Jardín Botánico de la Universidad Austral, y la pendiente es de esos cálculos que parecen simples hasta que el número que devuelve la consola no tiene ningún sentido.
Detrás de ese número absurdo hay una explicación concreta, y una vez que la entiendes el resto del proceso — elegir el paquete correcto, decidir entre grados y radianes, escoger un método de cálculo que no se vuelva ruidoso en terreno escarpado — se ordena bastante rápido.
De datos satelitales de elevación a un mapa de pendientes que sirve
El archivo de partida suele ser un Modelo digital de elevación (DEM), del tipo SRTM con resolución de 30 metros: cada cuadrito de la grilla representa esa área en el terreno real. Un raster, en el fondo, no es otra cosa que eso — una grilla de números —, así que hasta ahí no hay ningún misterio.
Ahí aparece el problema real: calcular la pendiente sin revisar qué Sistema de referencia de coordenadas trae el archivo. La mayoría de los DEM llegan en el código EPSG 4326 (WGS 84), que mide todo en grados, no en metros. Calcular una pendiente significa comparar una subida en metros contra una distancia también en metros; si mezclas grados con metros la matemática se rompe y entrega esos valores diminutos e inútiles. La regla práctica: antes de tocar terrain(), revisa en qué unidad viene tu raster. Si dice grados, hay que proyectarlo primero. Explicar a fondo qué es un sistema de referencia y por qué existen tantos da para su propio texto aparte — acá basta con la regla de dejarlo en metros antes de seguir.
raster o terra: La decisión que ahorra tiempo
Durante bastante tiempo trabajé con el paquete 'raster', que es el que aparece en casi todos los tutoriales viejos. Funciona, pero se pone lento con archivos pesados, y terra — su sucesor — resuelve lo mismo bastante más rápido. terra es apenas uno de varios paquetes espaciales que existen para R; hacer el repaso completo de esa lista es tema para otro artículo, pero para calcular pendientes, la diferencia de velocidad ya justifica el cambio.
Proyectar el raster con la función project() es como estirar la piel de una naranja sobre una mesa: pasas de una esfera medida en grados a algo plano medido en metros. En Chile lo habitual es usar el sistema UTM correspondiente a la zona. Una vez hecho ese cambio, los valores de pendiente dejan de ser 0.001 y empiezan a ser números que corresponden a lo que ves en el terreno.
Hace poco comentaba con Denise Alarcón, un contacto de un foro de R, que por qué estudiar teledetección satelital con R para ver cambios ambientales engancha justo por estos momentos: cuando dejas de ver puros errores en rojo y la forma real del terreno aparece en la pantalla. Denise domina ggplot2 pero se atora justo en el paso de convertir un raster a data frame, así que cuando comparamos apuntes terminamos completando lo que a cada una le cuesta.
Grados o radianes: La opción que cambia el mapa
Con el raster ya proyectado y sus unidades de X, Y y Z en metros, el comando terrain() calcula la pendiente. La función terrain de terra deja elegir entre "degrees" y "radians", y si eliges mal el mapa dirá que un cerro tiene 0.5 grados de inclinación cuando en realidad es un precipicio. La opción que conviene casi siempre es "degrees": un ángulo de 45 grados ya es algo que cuesta mucho subir a pie, y esa referencia es más fácil de leer que un número en radianes que hay que convertir mentalmente.
El método Horn y sus límites en terreno escarpado
Por defecto, terrain() usa el método Horn, que es el estándar en casi todos los programas de análisis espacial. Funciona bien en terreno suave, pero tiende a sobreestimar la pendiente en relieve muy escarpado, como el de buena parte del sur de Chile — el resultado se ve con textura de sal y pimienta, lleno de píxeles sueltos que no representan bien la ladera real. Ahí conviene probar una vecindad más amplia en el cálculo, que suaviza ese ruido a cambio de perder algo de detalle fino. No es que Horn esté mal: para un vistazo general basta, pero si el mapa sale demasiado picado, ese es el ajuste que hay que tocar.
Armar el mapa final: de terrain() a los colores en pantalla
El paso final es correr terrain(x, v = "slope", unit = "degrees") sobre el raster ya proyectado y ponerle una paleta de colores que tenga sentido — las zonas más empinadas en un café oscuro, casi rojizo, y las más planas en tonos claros. Ahí aparecen las sombras de las quebradas, por dónde corre el agua cuando llueve fuerte, y por qué ciertos sectores del bosque se mantienen más húmedos que otros.
Si el área que te interesa cae justo en el borde de una sola escena DEM, vas a necesitar unir varios rasters en RStudio para ver paisajes extensos, porque un solo cuadro casi nunca alcanza para cubrir toda una cuenca. Si en cambio el DEM te queda más grande de lo necesario y solo te interesa un polígono puntual — el perímetro del jardín botánico, por ejemplo —, recortarlo contra un shapefile es otro paso posible que no cubro en este texto. Y una vez calculada la pendiente, ese mismo raster sirve después para extraer el valor exacto en un punto concreto — dónde se cayó un árbol, por decir algo —, aunque ese truco de extracción también queda para otro momento.
Lo que no sirvió: Google Earth Engine y otras rutas
Antes de instalarme en R para esto, probé Google Earth Engine, pero la interfaz de JavaScript me cerró la puerta desde el principio — no venía de programación y ese entorno daba por sentado demasiado desde la primera línea. Con terra el proceso se sintió más parecido a ir armando frases: cargas el archivo, revisas el resultado, corriges lo que salió mal, repites. Cómo fui aprendiendo todo esto por mi cuenta, sin venir de programación, es una historia que da para otro texto aparte.
Mi vecino Boris Irarrázaval, que trabaja en infraestructura de servidores y no sabe nada de SIG, me preguntó una vez por qué el mapa cambiaba tanto según si elegía grados o radianes, y esa pregunta me obligó a explicarlo en palabras simples, sin código de por medio.
Nada de esto tiene que ver con calcular NDVI a partir de bandas espectrales — ese es un proceso completamente distinto, con imágenes ópticas en vez de un modelo de elevación. La vez que vi un parche rojo firme justo donde pasa el sendero de verano en esas bandas, estaba mirando vegetación, no relieve. Conseguir esas imágenes ópticas en primer lugar, decidir de dónde bajarlas y qué resolución pedir, es otro dolor de cabeza aparte del que estamos hablando acá. Y a diferencia de esas imágenes, que a veces llegan con nubes tapando medio paisaje, un DEM no tiene ese problema: el relieve no se nubla.
Antes de confiar en tu mapa de pendientes
Con todo esto ya armado, conviene revisar tres cosas antes de dar por bueno el resultado: que el raster esté proyectado en metros y no en grados, que la unidad elegida en terrain() sea la que de verdad quieres leer, y que el método de cálculo no esté generando ruido en el terreno más escarpado. Un buen chequeo final es simplemente comparar el mapa con lo que ya conoces del terreno a pie — si una ladera que sabes que es empinada aparece plana en pantalla, algo quedó mal proyectado.
La pendiente es apenas un paso dentro de un flujo más largo. Clasificar coberturas vegetales sobre este mismo terreno es una clasificación supervisada, con su propia lógica de entrenamiento, y validar esa clasificación después necesita una matriz de confusión — otro tema completo. Si el interés es ver cómo cambia el bosque de un año a otro, ahí entran las series temporales de satélite. La pendiente incluso puede servir de insumo para modelos de distribución de especies, aunque ese es un terreno que todavía no piso. Por ahora, con el mapa terminado y los colores marcando cada quebrada, ya se puede mirar el bosque del jardín botánico sabiendo por dónde corre el agua cuando llueve.