1. Introducción

Este notebook ilustra cómo acceder a series de tiempo de imágenes Sentinel-2 usando los paquetes rstac y gdalcubes, y cómo calcular el índice NDVI para el municipio de Tuluá, Valle del Cauca.

rstac ofrece funciones que simplifican el proceso de adquisición de datos espaciales desde un catálogo STAC (SpatioTemporal Asset Catalog). Más información aquí.

gdalcubes es un paquete de R y una librería en C++ que facilita el procesamiento de colecciones de imágenes satelitales de forma más rápida e interactiva. Más información aquí.

Usaremos el catálogo COG de Sentinel-2 disponible libremente en Amazon Web Services (AWS) y el endpoint STAC-API correspondiente en https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a.

Adicionalmente usaremos terra para leer archivos raster y tmap para la visualización.

2. Instalación de paquetes

Primero es necesario instalar los paquetes requeridos.

Usa la consola para correr el siguiente código -una línea a la vez-:

## No correr este código desde aquí - ¡Usar la consola de R!
#install.packages("rstac")
#install.packages("tidyterra")
#install.packages("gdalcubes")
#install.packages("stars")
#install.packages("mapview")
#install.packages("tmap")
#install.packages("leaflet")
#install.packages("sf")
#install.packages("geodata")

3. Configuración inicial

Limpiar la memoria y cargar las librerías:

# limpiar espacio de trabajo
rm(list = ls(all = TRUE))

# cargar librerías
library(rstac)
library(gdalcubes)
library(stars)
library(tmap)
library(tmaptools)
library(mapview)
library(leaflet)
library(sf)
library(ggplot2)
library(terra)
library(geodata)

4. Definir el área de interés (AOI)

En lugar de dibujar el polígono manualmente, obtenemos el límite oficial del municipio de Tuluá a partir de los límites administrativos de Colombia (GADM, nivel 2 - municipios), usando el paquete geodata.

# descargar límites administrativos de nivel 2 (municipios) de Colombia
col_mun <- gadm(country = "COL", level = 2, path = tempdir())
## Cached as: C:\Users\RICHIE\AppData\Local\Temp\RtmpGGcqEt/gadm/gadm41_COL_2_pk.rds
col_mun_sf <- st_as_sf(col_mun)

# filtrar el municipio de Tuluá
tulua <- col_mun_sf[col_mun_sf$NAME_2 == "Tuluá", ]

# revisar que se haya encontrado el polígono correcto
plot(st_geometry(tulua))

Guardamos el AOI como GeoJSON en el directorio de datos, para tenerlo disponible y reutilizarlo sin volver a descargar los límites administrativos:

dir.create("tulua", showWarnings = FALSE)
st_write(tulua, "tulua/map.geojson", delete_dsn = TRUE)
## Deleting source `tulua/map.geojson' using driver `GeoJSON'
## Writing layer `map' to data source `tulua/map.geojson' using driver `GeoJSON'
## Writing 1 features with 13 fields and geometry type Polygon.

En este notebook, el AOI corresponde al municipio de Tuluá, Valle del Cauca.

Luego, leer el archivo geojson guardado:

# usar la función list.files() para revisar el directorio y los nombres de archivo
tulua <- st_read("tulua/map.geojson")
## Reading layer `map' from data source `C:\Users\RICHIE\Desktop\tulua\map.geojson' using driver `GeoJSON'
## Simple feature collection with 1 feature and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.314 ymin: 3.894089 xmax: -75.81911 ymax: 4.185584
## Geodetic CRS:  WGS 84
## Reading layer `map' from data source
##   `.../tulua/map.geojson'
##   using driver `GeoJSON'
## Simple feature collection with 1 feature and n fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.28 ymin: 3.95 xmax: -76.05 ymax: 4.20
## Geodetic CRS:  WGS 84

Nota: si prefieres el límite desde OpenStreetMap en lugar de GADM, puedes usar osmdata::getbb("Tuluá, Valle del Cauca, Colombia", format_out = "sf_polygon").

5. Buscar imágenes Sentinel-2

La función stac_search() permite buscar datos Sentinel-2 dado un AOI georreferenciado en EPSG:4326 y un periodo de tiempo. Intentémoslo:

Nota sobre rendimiento: para una primera prueba usamos un rango de tiempo corto (2 meses) y menos bandas/menor resolución en las secciones siguientes. Una vez confirmes que todo el flujo funciona, puedes ampliar el rango a "2025-01-01/2025-12-31" y afinar la resolución a dx=10, dy=10.

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

items <- s |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = c(st_bbox(tulua)["xmin"],
                       st_bbox(tulua)["ymin"],
                       st_bbox(tulua)["xmax"],
                       st_bbox(tulua)["ymax"]),
              datetime = "2025-11-01/2025-12-31",
              limit = 60) |>
  post_request()

items
## ###Items
## - matched feature(s): 35
## - features (35 item(s) / 0 not fetched):
##   - S2B_18NUK_20251213_0_L2A
##   - S2C_18NUK_20251211_0_L2A
##   - S2A_18NUK_20251210_0_L2A
##   - S2A_18NVK_20251210_0_L2A
##   - S2C_18NUK_20251208_0_L2A
##   - S2C_18NVK_20251208_0_L2A
##   - S2B_18NUK_20251206_0_L2A
##   - S2B_18NUK_20251203_0_L2A
##   - S2B_18NVK_20251203_0_L2A
##   - S2C_18NUK_20251201_0_L2A
##   - ... with 25 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

Imprimir la fecha y hora de la primera y última imagen:

range(sapply(items$features, function(x) {x$properties$datetime}))
## [1] "2025-11-01T15:42:13Z" "2025-12-13T15:32:03Z"

6. Crear una colección Sentinel-2

Podemos convertir la respuesta STAC en una colección de imágenes de gdalcubes usando stac_image_collection(). Esta función espera una lista de features STAC como entrada y opcionalmente puede aplicar filtros sobre los metadatos y las bandas. A continuación creamos una colección con imágenes que tengan menos del 20% de cobertura de nubes.

# solo las bandas necesarias para RGB, falso color, NDVI y máscara de nubes
assets <- c("B02","B03","B04","B08","B11","SCL")

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 2 images with 6 bands
## Images:
##                       name      left      top   bottom     right
## 1 S2C_18NUK_20251208_0_L2A -76.34752 4.523459 3.529607 -75.81212
## 2 S2C_18NVK_20251208_0_L2A -75.90178 4.523910 3.530150 -74.91203
##              datetime        srs
## 1 2025-12-08T15:32:13 EPSG:32618
## 2 2025-12-08T15:32:10 EPSG:32618
## 
## Bands:
##   name offset scale unit nodata image_count
## 1  B02      0     1                       2
## 2  B03      0     1                       2
## 3  B04      0     1                       2
## 4  B08      0     1                       2
## 5  B11      0     1                       2
## 6  SCL      0     1                       2

7. Crear y procesar el cubo de datos

Con el objeto de colección de imágenes, ahora podemos usar gdalcubes. Definimos una vista del cubo de datos (data cube view) y algunas operaciones adicionales. En el siguiente ejemplo creamos una imagen RGB de mediana (median-composite) con resolución de 100 m a partir de un cubo de datos mensual.

(tulua_box <- st_bbox(tulua))
##       xmin       ymin       xmax       ymax 
## -76.314003   3.894089 -75.819115   4.185584
(tulua_b9377 <- st_bbox(st_transform(tulua, "EPSG:9377")))
##    xmin    ymin    xmax    ymax 
## 4632105 1988890 4687039 2021116
gdalcubes_options(parallel = 8)

(v <- cube_view(srs = "EPSG:9377", dx = 30, dy = 30, dt = "P1M",
                 aggregation = "median", resampling = "average",
                 extent = list(t0 = "2025-11-01", t1 = "2025-12-31",
                               left = tulua_b9377["xmin"],
                               right = tulua_b9377["xmax"],
                               top = tulua_b9377["ymax"],
                               bottom = tulua_b9377["ymin"])))
## A data cube view object
## 
## Dimensions:
##                low             high count pixel_size
## t       2025-11-01       2025-12-31     2        P1M
## y 1988878.00183763 2021128.00183763  1075         30
## x 4632092.19384396 4687052.19384396  1832         30
## 
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"

8. Visualizar el cubo de datos

Primero, color natural (RGB 4-3-2):

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))

Ahora, color falso (RGB 8-11-4):

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))

9. Descargar los datos

# "descargamos" los datos y los guardamos como archivos tif
tutorialDir <- path.expand("./tulua/")
outDir <- file.path(tutorialDir, "Tulua_30")

if (!file.exists(outDir)) {
  dir.create(outDir, recursive = TRUE)
  s2.mask <- image_mask("SCL", values = c(3, 8, 9))
  gdalcubes_options(parallel = 8, ncdf_compression_level = 0)
  raster_cube(s2_collection, v, mask = s2.mask) |>
    write_tif(dir = outDir)
}

Revisar los nombres de los archivos descargados:

list.files(outDir)
## [1] "cube_3600342d6c042025-11-01.tif" "cube_3600342d6c042025-12-01.tif"

10. Calcular un índice de vegetación (NDVI) y visualizarlo

Primero, leer un archivo GTiff (por ejemplo, la composición de diciembre):

# buscar automáticamente el archivo correspondiente a diciembre de 2025
archivos <- list.files(outDir, pattern = "2025-12-01\\.tif$", full.names = TRUE)
f <- archivos[1]
tulua_s2 <- terra::rast(f)

Luego, calcular el índice NDVI:

# revisar los nombres de las bandas disponibles en el archivo
names(tulua_s2)
## [1] "B02" "B03" "B04" "B08" "B11" "SCL"
# 1. Definir las bandas por nombre (más robusto que por índice)
# si los nombres exactos "B04"/"B08" no existen, se buscan por patrón
banda_red <- grep("B04", names(tulua_s2), value = TRUE)[1]
banda_nir <- grep("B08", names(tulua_s2), value = TRUE)[1]

red <- tulua_s2[[banda_red]]
nir <- tulua_s2[[banda_nir]]

# 2. Calcular el NDVI
ndvi <- (nir - red) / (nir + red)

ndvi
## class       : SpatRaster
## size        : 1075, 1832, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4632092, 4687052, 1988878, 2021128  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s)   : memory
## varname     : cube_3600342d6c042025-12-01
## name        :       B08
## min value   : -0.997622
## max value   :  0.998863

Visualización con tmap:

tm_shape(ndvi) +
  tm_raster(
    col_palette = "RdYlGn",
    style = "cont",
    title = "Valor NDVI"
  ) +
  tm_layout(main.title = "NDVI Sentinel-2 - Tuluá, dic. 2025")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
## `tm_scale_continuous()`.
## [v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
## map variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
## [tm_raster()] Argument `col_palette` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [`tm_scale_continous()`] Variable(s) "col" contains positive and negative
## values, so midpoint is set to 0. Set midpoint = NA to show the full range of
## visual values.

Visualización interactiva con leaflet:

ndvi_palette <- colorNumeric(
  palette = c("red", "yellow", "green"),
  domain = c(-1, 1),
  na.color = "transparent"
)

map <- leaflet() %>%
  addProviderTiles("Esri.WorldImagery", group = "Esri Satellite") %>%
  setView(lng = -76.20, lat = 4.09, zoom = 12) %>%
  addRasterImage(
    x = ndvi,
    colors = ndvi_palette,
    opacity = 0.6,
    group = "NDVI Layer"
  ) %>%
  addLegend(
    pal = ndvi_palette,
    values = c(-1, 1),
    title = "Índice NDVI",
    position = "bottomright",
    opacity = 0.8
  ) %>%
  addLayersControl(
    baseGroups = c("Esri Satellite"),
    overlayGroups = c("NDVI Layer"),
    options = layersControlOptions(collapsed = FALSE)
  )

map

11. Referencias

Appel, M., & Pebesma, E. (2019). On-demand processing of data cubes from satellite image collections with the gdalcubes library. Data, 4(3), 92. https://doi.org/10.3390/data4030092

Cheng, J., Schloerke, B., Karambelkar, B., Xie, Y., & Aden-Buie, G. (2025). leaflet: Create interactive web maps with the JavaScript ‘Leaflet’ library (R package). https://rstudio.github.io/leaflet/

Element84. (n.d.). Earth Search STAC API [Data catalog]. Amazon Web Services. https://earth-search.aws.element84.com/v0

Global Administrative Areas (GADM). (2024). GADM database of global administrative areas (Version 4.1) [Data set]. https://gadm.org

Hijmans, R. J. (2025a). terra: Spatial data analysis (R package version 1.8-76). https://rspatial.org/

Hijmans, R. J. (2025b). geodata: Access geographic data (R package version 0.6-8). https://github.com/rspatial/geodata

Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. The R Journal, 10(1), 439–446. https://doi.org/10.32614/RJ-2018-009

Rouse, J. W., Haas, R. H., Schell, J. A., & Deering, D. W. (1974). Monitoring vegetation systems in the Great Plains with ERTS. NASA Special Publication, 351, 309–317.

Simoes, R., Souza, F., Zaglia, M., Queiroz, G. R. de, Santos, R. D. C. dos, & Ferreira, K. R. (2021). Rstac: An R package to access Spatiotemporal Asset Catalog satellite imagery. In 2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS (pp. 7674–7677). https://doi.org/10.1109/IGARSS47720.2021.9553518

Tennekes, M. (2018). tmap: Thematic maps in R. Journal of Statistical Software, 84(6), 1–39. https://doi.org/10.18637/jss.v084.i06