1. Introducción

Este cuaderno muestra el procedimiento para acceder, descargar y procesar imágenes del satélite Sentinel-2 mediante el lenguaje R, tomando como área de estudio la Ciénaga Grande, ubicada en el departamento del Magdalena, Colombia.

Mediante el uso de catálogos de activos espaciotemporales (STAC) y herramientas como rstac y gdalcubes.

Las imágenes de Sentinel-2 ofrecen información multiespectral de alta resolución, con 10 bandas espectrales y resoluciones espaciales que oscilan entre 10 y 60 metros. Con estos datos se elaboraron composiciones en Color Verdadero (RGB 4-3-2) y en Falso Color Infrarrojo, permitiendo analizar y representar de mejor manera los humedales costeros, la distribución de los manglares y los cuerpos de agua presentes en la región.

2. CARGAR LAS LIBRERIAS y LIMPIAR LA MEMORIA.

rm(list = ls(all=TRUE))

library(rstac)
library(gdalcubes)
library(stars)
library(tmap)
library(tmaptools)
library(mapview)
library(leaflet)
library(sf)
library(ggplot2)

3. Definición del área de interés (AOI)

Cargo mi zona de estudio (la Ciénaga). Lo tengo en un archivo geojson que preparé antes.

union <- st_read("C:\\Users\\suare\\OneDrive\\Documentos\\GB2\\proyecto5\\cienaga.geojson")
## Reading layer `cienaga' from data source 
##   `C:\Users\suare\OneDrive\Documentos\GB2\proyecto5\cienaga.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.33083 ymin: 10.89692 xmax: -74.13611 ymax: 11.03113
## Geodetic CRS:  WGS 84

4. Búsqueda en el catálogo STAC

Nos conectamos al servidor STAC de AWS (Earth Search) para consultar la disponibilidad de imágenes Sentinel-2 L2A durante el año 2025 sobre el área de interés.

s = stac("https://earth-search.aws.element84.com/v0")

items <- s |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = c(st_bbox(union)["xmin"],
                       st_bbox(union)["ymin"],
                       st_bbox(union)["xmax"],
                       st_bbox(union)["ymax"]), 
              datetime = "2025-01-01/2025-12-31",
              limit = 60) |>
  post_request() 
items
## ###Items
## - matched feature(s): 58
## - features (58 item(s) / 0 not fetched):
##   - S2B_18PWT_20251213_0_L2A
##   - S2C_18PWT_20251208_0_L2A
##   - S2B_18PWT_20251203_0_L2A
##   - S2A_18PWT_20251130_0_L2A
##   - S2C_18PWT_20251128_0_L2A
##   - S2B_18PWT_20251123_0_L2A
##   - S2A_18PWT_20251120_0_L2A
##   - S2C_18PWT_20251118_0_L2A
##   - S2B_18PWT_20251113_0_L2A
##   - S2A_18PWT_20251110_0_L2A
##   - ... with 48 more feature(s).
## - assets: 
## AOT, B01, B02, B03, B04, B05, B06, B07, B08, B09, B11, B12, B8A, info, metadata, overview, SCL, thumbnail, visual, WVP
## - item's fields: 
## assets, bbox, collection, geometry, id, links, properties, stac_extensions, stac_version, type
range(sapply(items$features, function(x) {x$properties$datetime}))
## [1] "2025-06-01T15:30:23Z" "2025-12-13T15:30:00Z"

5. FILTRADO DE BANDAS Y NUBES.

Selecciono solo las bandas que necesito.

assets = c("B01","B02","B03","B04","B05","B06", "B07","B08","B8A","B09","B11","SCL")
#Armo la colección, pero filtro para quedarme únicamente con imágenes que tengan menos del 20% de nubes.
s2_collection = stac_image_collection(items$features, asset_names = assets, property_filter = function(x) {x[["eo:cloud_cover"]] < 20})
s2_collection
## Image collection object, referencing 12 images with 12 bands
## Images:
##                       name      left      top   bottom     right
## 1 S2C_18PWT_20251208_0_L2A -74.75914 11.75966 10.76541 -73.99249
## 2 S2B_18PWT_20251113_0_L2A -74.74231 11.75964 10.76541 -73.99249
## 3 S2C_18PWT_20251108_0_L2A -74.75529 11.75966 10.76541 -73.99249
## 4 S2C_18PWT_20250929_0_L2A -74.75612 11.75966 10.76541 -73.99249
## 5 S2A_18PWT_20250911_0_L2A -74.74295 11.75964 10.76541 -73.99249
## 6 S2C_18PWT_20250909_0_L2A -74.74587 11.75964 10.76541 -73.99249
##              datetime        srs
## 1 2025-12-08T15:30:11 EPSG:32618
## 2 2025-11-13T15:29:59 EPSG:32618
## 3 2025-11-08T15:30:14 EPSG:32618
## 4 2025-09-29T15:30:18 EPSG:32618
## 5 2025-09-11T15:30:15 EPSG:32618
## 6 2025-09-09T15:30:18 EPSG:32618
## [ omitted 6 images ] 
## 
## Bands:
##    name offset scale unit nodata image_count
## 1   B01      0     1                      12
## 2   B02      0     1                      12
## 3   B03      0     1                      12
## 4   B04      0     1                      12
## 5   B05      0     1                      12
## 6   B06      0     1                      12
## 7   B07      0     1                      12
## 8   B08      0     1                      12
## 9   B09      0     1                      12
## 10  B11      0     1                      12
## 11  B8A      0     1                      12
## 12  SCL      0     1                      12

6. ESTRUCTURACIÓN DEL CUBO DE DATOS

(union_box = st_bbox(union))
##      xmin      ymin      xmax      ymax 
## -74.33083  10.89692 -74.13611  11.03113

Reproyecto el polígono al Origen Nacional (EPSG:9377)

(union_b9377 <- st_bbox(st_transform(union, "EPSG:9377")))
##    xmin    ymin    xmax    ymax 
## 4854676 2762382 4875941 2777306

#Defino las reglas de mi cubo de datos.

gdalcubes_options(parallel = 8)
(v = cube_view(srs="EPSG:9377", dx=30, dy=30, dt="P1M", 
                           aggregation="median", resampling = "average",
                           extent=list(t0 = "2025-01-01", t1 = "2025-12-31",
                                       left=union_b9377["xmin"], right=union_b9377["xmax"],
                                       top=union_b9377["ymax"], bottom=union_b9377["ymin"])))
## A data cube view object
## 
## Dimensions:
##                low             high count pixel_size
## t       2025-01-01       2025-12-31    12        P1M
## y 2762374.05425008 2777314.05425008   498         30
## x  4854673.6878443  4875943.6878443   709         30
## 
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"

7. VISUALIZACIÓN: COLOR VERDADERO

Ploteo la composición RGB normal (Rojo=B04, Verde=B03, Azul=B02) para ver la Ciénaga tal cual ve el ojo humano

gdalcubes_options(threads = 2)
## Warning in gdalcubes_options(threads = 2): 'threads' option is deprecated;
## please use 'parallel' instead
raster_cube(s2_collection, v) |>
  select_bands(c("B02","B03","B04")) |>
  reduce_time(c("median(B02)", "median(B03)", "median(B04)")) |>
  plot(rgb = 3:1, zlim = c(0,2500))

Al armar la imagen en color verdadero (mezclando las bandas 4, 3 y 2), el resultado es básicamente lo que veríamos a simple vista si estuviéramos sobrevolando la zona. Observando el mapa de la Ciénaga, lo primero que noto es que los cuerpos de agua se ven súper oscuros, entre un azul muy profundo y negro, lo cual tiene sentido por la profundidad y la cantidad de sedimentos que maneja el lugar. Por otro lado, toda la zona de manglares y vegetación es muy fácil de identificar porque resalta en diferentes tonos de verde. Finalmente, esos parches más grises, blanquitos o color arena corresponden a los suelos desnudos, los playones o las áreas urbanas que están alrededor.

8. VISUALIZACIÓN: FALSO COLOR.

Ahora hago una composición de falso color usando infrarrojo (B08, B11, B04).Con esto diferencio muchísimo mejor los cuerpos de agua de los manglares y vegetación.

raster_cube(s2_collection, v) |>
  select_bands(c("B04","B11","B08")) |>
  reduce_time(c("median(B04)", "median(B11)", "median(B08)")) |>
  plot(rgb = 3:1, zlim = c(0,4500))

Cuando se cambia a la composición en falso color (usando las bandas B08, B11 y B04), la imagen dio un giro total. Como estoy usando las bandas infrarrojas, la vegetación más densa y saludable, como los bosques de manglar, ahora ‘brilla’ en tonos anaranjados. Esto pasa porque las hojas reflejan muchísimo esa luz infrarroja que el ojo huemano no ve. Lo que más me gusta de esta vista es lo fácil que es separar la tierra del agua; como el agua no refleja casi nada de infrarrojo , la ciénaga y los humedales quedan casi negros o azules oscuros, logrando un contraste súper marcado que en el color verdadero era mucho más difícil de notar.

9. LIMPIEZA FINAL Y EXPORTACIÓN

library(gdalcubes)

# 1. Definir la carpeta del proyecto y la subcarpeta de salida
tutorialDir <- path.expand("C:/Users/suare/OneDrive/Documentos/GB2/proyecto561")
outputFolder <- file.path(tutorialDir, "cienaga_10")

# 2. Crear la carpeta si aún no existe
if (!dir.exists(outputFolder)) {
  dir.create(outputFolder, recursive = TRUE)
}

# 3. Configurar e integrar el cubo
s2.mask <- image_mask("SCL", values = c(3, 8, 9))
gdalcubes_options(parallel = 4, ncdf_compression_level = 0)

# 4. Generar y exportar a la carpeta (outputFolder)
raster_cube(s2_collection, v, mask = s2.mask) |> 
  select_bands(c("B04", "B08")) |> 
  reduce_time(c("median(B04)", "median(B08)")) |> 
  write_tif(outputFolder) # Se pasa la carpeta, NO una ruta final con .tif

message("¡Proceso completado exitosamente en: ", outputFolder, "!")
library(terra)
## Warning: package 'terra' was built under R version 4.5.3
## terra 1.9.34
## 
## Adjuntando el paquete: 'terra'
## The following objects are masked from 'package:gdalcubes':
## 
##     animate, crop, size
# Buscar el archivo GeoTIFF recién creado en la carpeta
f <- list.files(outputFolder, pattern = "\\.tif$", full.names = TRUE)[1]

# Cargar con terra
cienaga_s2 <- terra::rast(f)

# Extraer las medianas calculadas
# Extraer la 1a capa (Rojo / B04) y la 2a capa (NIR / B08) por su índice
red <- cienaga_s2[[1]]
nir <- cienaga_s2[[2]]

# Calcular NDVI
ndvi <- (nir - red) / (nir + red)
names(ndvi) <- "NDVI"

# Graficar para confirmar
plot(ndvi, main = "NDVI Ciénaga (Mediana 2025)", col = rev(terrain.colors(10)))

library(leaflet)
library(terra)

# 1. Reproyectar el NDVI a WGS84 (EPSG:4326) para que leaflet lo pueda mostrar
ndvi_wgs84 <- project(ndvi, "EPSG:4326")

# 2. Definir la paleta de colores para el NDVI (-1 a 1)
ndvi_palette <- colorNumeric(
  palette = c("red", "yellow", "green"), 
  domain = c(-1, 1), 
  na.color = "transparent"
)

# 3. Crear el mapa interactivo con leaflet centrado en la Ciénaga
mapa_interactivo <- leaflet() %>% 
  addProviderTiles("Esri.WorldImagery", group = "Esri Satellite") %>% 
  # Coordenadas ajustadas al centro de tu área de estudio
  setView(lng = -74.23347, lat = 10.96403, zoom = 12) %>% 
  addRasterImage(
    x = ndvi_wgs84, 
    colors = ndvi_palette, 
    opacity = 0.6, 
    group = "NDVI Layer"
  ) %>% 
  addLegend(
    pal = ndvi_palette, 
    values = c(-1, 1), 
    title = "NDVI Index",
    position = "bottomright",
    opacity = 0.8
  ) %>% 
  addLayersControl(
    baseGroups = c("Esri Satellite"),
    overlayGroups = c("NDVI Layer"),
    options = layersControlOptions(collapsed = FALSE)
  )

# 4. Desplegar el mapa
mapa_interactivo

10. ¿Qué aprendí hoy?

  1. Sobre el color verdadero (RGB 4-3-2): Me di cuenta de que al mezclar las bandas roja, verde y azul vemos el mapa tal cual como si estuviéramos ahí. Es súper útil y fácil de entender, pero también noté que a veces tiene sus límites; por ejemplo, si hay un humedal muy oscuro, se confunde fácilmente con un cuerpo de agua normal porque los tonos se ven casi iguales.

  2. Sobre el falso color: Aprendí que usar las bandas infrarrojas (como la B08 y B11) es increíble porque nos deja ver cosas que el ojo humano normalmente se perdería. Me pareció muy práctico cómo el infrarrojo hace que la vegetación sana “grite” en color naranja, y cómo nos ayuda a delimitar los bordes del agua con mucha exactitud porque absorbe la luz y queda oscura.

3.Sobre la resolución de Sentinel-2: Algo nuevo para mí fue entender que Sentinel-2 no toma toda la imagen con el mismo nivel de detalle. Maneja distintas resoluciones: tiene bandas muy nítidas a 10 metros (que son geniales para el color y el infrarrojo cercano), otras a 20 metros que sirven más para ver detalles agrícolas o de humedad, y unas menos detalladas a 60 metros que se usan más que todo para corregir cosas de la atmósfera, como las nubes.

Bibliografía.

Lizarazo, I. 2026. Accessing Sentinel-2 imagery. Disponible en: https://rpubs.com/ials2un/1436811