Acceso a series de tiempo de imágenes Sentinel-2
Estudio de caso: Puerto Leguízamo, Putumayo
Victor Zarate
2026-07-24
1 Introducción
Este cuaderno explica cómo acceder a series de tiempo de imágenes Sentinel-2 usando los paquetes rstac y gdalcubes, y calcular el índice NDVI para dos fechas distintas sobre un área de interés (AOI) propia. Se replica la metodología del profesor Ivan Lizarazo, aplicada aquí al municipio de Puerto Leguízamo, Putumayo, en la Amazonía colombiana.
2 Configuración del entorno
limpia cualquier objeto que quede en la memoria de R de sesiones anteriores, y carga todas las librerías necesarias: rstac/gdalcubes para buscar y procesar las imágenes satelitales, sf para manejar el polígono del AOI, terra para leer los GeoTIFF descargados, y tmap/leaflet para la visualización.
3 Área de interés (AOI)
leemos el polígono que dibujamos en geojson.io sobre Puerto Leguízamo y confirma sus límites geográficos (bounding box), que se usarán en el siguiente paso para acotar la búsqueda de imágenes.
## Reading layer `map' from data source
## `C:\Users\ASUS\OneDrive\Documents\unal\GB2\proyecto 4\Datos\map.geojson'
## using driver `GeoJSON'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.77973 ymin: -0.200972 xmax: -74.77194 ymax: -0.192455
## Geodetic CRS: WGS 84## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.77973 ymin: -0.200972 xmax: -74.77194 ymax: -0.192455
## Geodetic CRS: WGS 84
## geometry
## 1 POLYGON ((-74.77973 -0.1950...## xmin ymin xmax ymax
## -74.779726 -0.200972 -74.771941 -0.192455Resultado: 1 feature, tipo POLYGON, bbox aproximado xmin: -74.78, ymin: -0.20, xmax: -74.77, ymax: -0.19.
4 Búsqueda de imágenes Sentinel-2 vía STAC
se conecta al catálogo STAC de Sentinel-2 alojado en AWS y busca todas las escenas que intersectan el bbox del AOI dentro del rango de fechas indicado. stac_search() no descarga imágenes todavía — solo consulta el catálogo de metadatos.
s <- stac("https://earth-search.aws.element84.com/v0")
items <- s |>
stac_search(
collections = "sentinel-s2-l2a-cogs",
bbox = c(
st_bbox(aoi)["xmin"],
st_bbox(aoi)["ymin"],
st_bbox(aoi)["xmax"],
st_bbox(aoi)["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_18MWE_20251213_0_L2A
## - S2A_18MWE_20251210_0_L2A
## - S2C_18MWE_20251208_0_L2A
## - S2B_18MWE_20251203_0_L2A
## - S2A_18MWE_20251130_0_L2A
## - S2C_18MWE_20251128_0_L2A
## - S2B_18MWE_20251123_0_L2A
## - S2A_18MWE_20251120_1_L2A
## - S2C_18MWE_20251118_0_L2A
## - S2B_18MWE_20251113_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, typese revisan las fechas de la primera y última escena encontrada, para confirmar el rango temporal real cubierto por los resultados.
## [1] "2025-06-03T15:33:28Z" "2025-12-13T15:33:09Z"Resultado obtenido: 58 escenas encontradas, entre el 2025-06-03 y el 2025-12-13.
5 Creación de la colección de imágenes (filtro de nubosidad)
Qué hace: convierte la lista de escenas STAC en una colección de imágenes de gdalcubes, seleccionando las bandas espectrales de interés (assets) y filtrando solo las escenas con menos de 20% de nubosidad (eo:cloud_cover < 20), usando la propiedad de metadatos de cada escena.
assets <- c("B01","B02","B03","B04","B05","B06",
"B07","B08","B8A","B09","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 10 images with 12 bands
## Images:
## name left top bottom right
## 1 S2B_18MWE_20251213_0_L2A -75.00017 -9e-06 -0.99248 -74.41824
## 2 S2C_18MWE_20251208_0_L2A -75.00017 -9e-06 -0.99248 -74.42956
## 3 S2B_18MWE_20251203_0_L2A -75.00017 -9e-06 -0.99248 -74.41123
## 4 S2A_18MWE_20250921_1_L2A -75.00017 -9e-06 -0.99248 -74.41986
## 5 S2B_18MWE_20250914_0_L2A -75.00017 -9e-06 -0.99248 -74.41285
## 6 S2C_18MWE_20250830_0_L2A -75.00017 -9e-06 -0.99248 -74.42525
## datetime srs
## 1 2025-12-13T15:33:09 EPSG:32718
## 2 2025-12-08T15:33:19 EPSG:32718
## 3 2025-12-03T15:33:07 EPSG:32718
## 4 2025-09-21T15:33:27 EPSG:32718
## 5 2025-09-14T15:33:10 EPSG:32718
## 6 2025-08-30T15:33:30 EPSG:32718
## [ omitted 4 images ]
##
## Bands:
## name offset scale unit nodata image_count
## 1 B01 0 1 10
## 2 B02 0 1 10
## 3 B03 0 1 10
## 4 B04 0 1 10
## 5 B05 0 1 10
## 6 B06 0 1 10
## 7 B07 0 1 10
## 8 B08 0 1 10
## 9 B09 0 1 10
## 10 B11 0 1 10
## 11 B8A 0 1 10
## 12 SCL 0 1 10Resultado obtenido: 10 imágenes con menos de 20% de nubosidad, de las 58 encontradas.
6 Definición de la vista del cubo de datos
primero se reproyecta el AOI al sistema de coordenadas oficial de Colombia (EPSG:9377, en metros). Luego, cube_view() define la “receta” con la que se construirá el cubo: resolución espacial de 10 m, agregación temporal mensual (P1M) usando la mediana de todos los píxeles disponibles en cada mes (esto ayuda a reducir el ruido de nubes residuales), y el extent espacial/temporal completo.
aoi_bbox <- st_bbox(aoi)
aoi_9377 <- st_bbox(st_transform(aoi, "EPSG:9377"))
gdalcubes_options(parallel = 4)
v <- cube_view(
srs = "EPSG:9377", dx = 10, dy = 10, dt = "P1M",
aggregation = "median", resampling = "average",
extent = list(
t0 = "2025-01-01", t1 = "2025-12-31",
left = aoi_9377["xmin"], right = aoi_9377["xmax"],
top = aoi_9377["ymax"], bottom = aoi_9377["ymin"]
)
)
v## A data cube view object
##
## Dimensions:
## low high count pixel_size
## t 2025-01-01 2025-12-31 12 P1M
## y 1535829.96541385 1536779.96541385 95 10
## x 4802007.57683429 4802877.57683429 87 10
##
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"Nota: usamos parallel = 4 (u 8, según los núcleos de tu equipo) — valores muy altos como 24 pueden hacer que el proceso se vuelva más lento en vez de más rápido, en este caso mi computador no tiene esa cantidad de núcleos y al usar parallel = 24 estaba tardando alrededor de 40 min y no llego a terminar el proceso.
7 Visualización del cubo de datos
(color natural, RGB 4-3-2): calcula la mediana anual de las bandas azul (B02), verde (B03) y roja (B04), y las combina en una imagen de color natural — así se ve la escena aproximadamente como la vería el ojo humano.
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))
(falso color, RGB 8-11-4): combina la banda roja (B04), infrarrojo de onda corta (B11) e infrarrojo cercano (B08). La vegetación sana resalta en tonos rojo/naranja intenso porque refleja mucho en el infrarrojo cercano — es la combinación clásica para distinguir vegetación de suelo/agua/construcciones.
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))
8 Descarga del cubo de datos (GeoTIFF mensuales)
crea una máscara que descarta los píxeles clasificados como sombra de nube (código SCL 3), nube de probabilidad media (8) y nube de probabilidad alta (9). Luego write_tif() materializa el cubo, generando un archivo GeoTIFF por cada mes del año, guardado en la carpeta indicada con dir y el prefijo indicado con prefix.
s2_mask <- image_mask("SCL", values = c(3, 8, 9))
gdalcubes_options(parallel = 4)
raster_cube(s2_collection, v, mask = s2_mask) |>
write_tif(dir = "./puerto_leguizamo_10", prefix = "cube_")
list.files("puerto_leguizamo_10")## [1] "cube_2025-01-01.tif" "cube_2025-02-01.tif" "cube_2025-03-01.tif"
## [4] "cube_2025-04-01.tif" "cube_2025-05-01.tif" "cube_2025-06-01.tif"
## [7] "cube_2025-07-01.tif" "cube_2025-08-01.tif" "cube_2025-09-01.tif"
## [10] "cube_2025-10-01.tif" "cube_2025-11-01.tif" "cube_2025-12-01.tif"9 Verificación de cobertura de datos por mes
recorre los 12 archivos mensuales y calcula qué porcentaje de píxeles tiene datos válidos (no NA por nubes o ausencia de imagen), para identificar qué meses son útiles para el análisis de NDVI.
meses <- sprintf("2025-%02d-01", 1:12)
for (m in meses) {
f <- sprintf("puerto_leguizamo_10/cube_%s.tif", m)
if (file.exists(f)) {
r <- rast(f)
pct_validos <- sum(!is.na(values(r[[4]]))) / ncell(r) * 100
cat(m, "-", round(pct_validos, 1), "% píxeles válidos\n")
} else {
cat(m, "- archivo no existe\n")
}
}## 2025-01-01 - 0 % píxeles válidos
## 2025-02-01 - 0 % píxeles válidos
## 2025-03-01 - 0 % píxeles válidos
## 2025-04-01 - 0 % píxeles válidos
## 2025-05-01 - 0 % píxeles válidos
## 2025-06-01 - 0 % píxeles válidos
## 2025-07-01 - 13.1 % píxeles válidos
## 2025-08-01 - 100 % píxeles válidos
## 2025-09-01 - 100 % píxeles válidos
## 2025-10-01 - 0 % píxeles válidos
## 2025-11-01 - 0 % píxeles válidos
## 2025-12-01 - 100 % píxeles válidos## [1] "2025-01-01" "2025-02-01" "2025-03-01" "2025-04-01" "2025-05-01"
## [6] "2025-06-01" "2025-07-01" "2025-08-01" "2025-09-01" "2025-10-01"
## [11] "2025-11-01" "2025-12-01"Resultado obtenido: agosto, septiembre y diciembre con 100% de cobertura; el resto de meses en 0% (o parcial en julio) por nubosidad.
10 Cálculo del NDVI para dos fechas
para cada fecha, extrae la banda roja (B04) y la infrarroja cercana (B08) del GeoTIFF correspondiente, y aplica la fórmula del NDVI: (NIR - Red) / (NIR + Red). El resultado va de -1 a 1: valores altos (>0.5) indican vegetación densa y sana; valores cercanos a 0 o negativos indican suelo desnudo, agua, o construcciones.
f_ago <- "puerto_leguizamo_10/cube_2025-08-01.tif"
r_ago <- rast(f_ago)
red_ago <- r_ago[[4]]
nir_ago <- r_ago[[8]]
ndvi_ago <- (nir_ago - red_ago) / (nir_ago + red_ago)
ndvi_ago## class : SpatRaster
## size : 95, 87, 1 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : 4802008, 4802878, 1535830, 1536780 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## varname : cube_2025-08-01
## name : B08
## min value : -0.01286
## max value : 0.922815f_dic <- "puerto_leguizamo_10/cube_2025-12-01.tif"
r_dic <- rast(f_dic)
red_dic <- r_dic[[4]]
nir_dic <- r_dic[[8]]
ndvi_dic <- (nir_dic - red_dic) / (nir_dic + red_dic)
ndvi_dic## class : SpatRaster
## size : 95, 87, 1 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : 4802008, 4802878, 1535830, 1536780 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## varname : cube_2025-12-01
## name : B08
## min value : -0.081167
## max value : 0.91205911 Visualización del NDVI
generamos un mapa estático con paleta de color Rojo-Amarillo-Verde para cada fecha, facilitando la interpretación visual: verde oscuro = vegetación sana, blanco/rojo = suelo, agua o construcciones.
tm_shape(ndvi_ago) +
tm_raster(col_palette = "RdYlGn", style = "cont", title = "NDVI") +
tm_layout(main.title = "NDVI - Puerto Leguízamo, Agosto 2025")
tm_shape(ndvi_dic) +
tm_raster(col_palette = "RdYlGn", style = "cont", title = "NDVI") +
tm_layout(main.title = "NDVI - Puerto Leguízamo, Diciembre 2025")
Interpretación (agosto vs. diciembre):
El NDVI de Puerto Leguízamo se mantiene alto (0.6-0.9) en ambas fechas, consistente con cobertura de selva húmeda tropical, donde no existe una estacionalidad fenológica marcada como en cultivos anuales. La zona urbana se distingue claramente por valores bajos de NDVI (0-0.3) en ambas imágenes, manteniendo su extensión aproximadamente constante. Entre agosto y diciembre no se observan cambios drásticos en la cobertura vegetal, lo cual es esperable en bosque primario/secundario amazónico; la ligera variación observada en la zona centro-derecha podría estar asociada a diferencias fenológicas menores o a contaminación residual de nubes no completamente enmascarada.
12 Mapa interactivo NDVI.
crea un mapa web interactivo con leaflet, superponiendo el NDVI de diciembre sobre imagen satelital de fondo, con control de capas y leyenda — útil para explorar el resultado con zoom/pan.
ndvi_palette <- colorNumeric(
palette = c("red", "yellow", "green"),
domain = c(-1, 1),
na.color = "transparent"
)
centro <- st_coordinates(st_centroid(st_transform(aoi, 4326)))
leaflet() |>
addProviderTiles("Esri.WorldImagery", group = "Esri Satellite") |>
setView(lng = centro[1], lat = centro[2], zoom = 15) |>
addRasterImage(ndvi_dic, colors = ndvi_palette, opacity = 0.6, group = "NDVI") |>
addLegend(pal = ndvi_palette, values = c(-1, 1), title = "NDVI",
position = "bottomright", opacity = 0.8) |>
addLayersControl(baseGroups = c("Esri Satellite"),
overlayGroups = c("NDVI"),
options = layersControlOptions(collapsed = FALSE))13 Conclusiones
El presente ejercicio permitió acceder a series de tiempo de imágenes Sentinel-2 mediante los paquetes rstac y gdalcubes, para el municipio de Puerto Leguízamo, Putumayo, ubicado en la Amazonía colombiana. De un total de 58 escenas disponibles para el año 2025, solo 10 cumplieron el criterio de nubosidad menor al 20%, concentradas principalmente en los meses de agosto, septiembre y diciembre — reflejo de la alta cobertura de nubes característica del bioma amazónico, que limita fuertemente la disponibilidad de imágenes ópticas utilizables en la región.
A partir de los compuestos mensuales de mediana, se calculó el índice NDVI para dos fechas representativas (agosto y diciembre de 2025), encontrando valores máximos de 0.92 y 0.91 respectivamente, típicos de cobertura de bosque húmedo tropical denso y saludable. El patrón espacial del NDVI se mantuvo consistente entre ambas fechas, permitiendo identificar con claridad el área urbana del municipio (NDVI bajo, 0-0.3) frente a la matriz boscosa circundante (NDVI alto, 0.6-0.9). No se observaron cambios drásticos en la cobertura vegetal entre agosto y diciembre, lo cual es coherente con la baja estacionalidad fenológica de los bosques amazónicos primarios y secundarios.
Este ejercicio evidencia tanto el potencial como las limitaciones del uso de Sentinel-2 en zonas de alta nubosidad persistente: si bien el filtrado por porcentaje de nubosidad y el enmascaramiento con la banda SCL permiten obtener compuestos limpios, la disponibilidad temporal efectiva de datos se reduce considerablemente, restringiendo el monitoreo continuo de la vegetación en regiones como la Amazonía colombiana.
14 Bibliografia.
Lizarazo, I., 2025. Accessing time series of Sentinel-2 imagery Ivan Lizarazo. https://rpubs.com/ials2un/sentinel2_ts
15 Entorno de trabajo
## R version 4.6.0 (2026-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8 LC_CTYPE=Spanish_Colombia.utf8
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] terra_1.9-27 ggplot2_4.0.3 leaflet_2.2.3 mapview_2.11.4
## [5] tmaptools_3.3 tmap_4.3 stars_0.7-2 sf_1.1-1
## [9] abind_1.4-8 gdalcubes_0.7.4 rstac_1.0.1
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.57 bslib_0.10.0
## [4] raster_3.6-32 htmlwidgets_1.6.4 lattice_0.22-9
## [7] leaflet.providers_3.0.0 generics_0.1.4 vctrs_0.7.3
## [10] tools_4.6.0 crosstalk_1.2.2 curl_7.1.0
## [13] stats4_4.6.0 parallel_4.6.0 tibble_3.3.1
## [16] proxy_0.4-29 pkgconfig_2.0.3 KernSmooth_2.23-26
## [19] satellite_1.0.6 data.table_1.18.4 RColorBrewer_1.1-3
## [22] S7_0.2.2 lifecycle_1.0.5 farver_2.1.2
## [25] compiler_4.6.0 codetools_0.2-20 leafsync_0.1.0
## [28] ncdf4_1.24 leaflegend_1.2.8 htmltools_0.5.9
## [31] class_7.3-23 sass_0.4.10 yaml_2.3.12
## [34] pillar_1.11.1 crayon_1.5.3 jquerylib_0.1.4
## [37] classInt_0.4-11 cachem_1.1.0 lwgeom_0.2-16
## [40] wk_0.9.5 tidyselect_1.2.1 digest_0.6.39
## [43] dplyr_1.2.1 fastmap_1.2.0 grid_4.6.0
## [46] colorspace_2.1-2 cli_3.6.6 logger_0.4.2
## [49] magrittr_2.0.5 maptiles_0.11.0 base64enc_0.1-6
## [52] XML_3.99-0.23 cols4all_0.10 leafem_0.2.5
## [55] e1071_1.7-17 withr_3.0.2 scales_1.4.0
## [58] sp_2.2-1 rmarkdown_2.31 httr_1.4.8
## [61] jpeg_0.1-11 otel_0.2.0 png_0.1-9
## [64] evaluate_1.0.5 knitr_1.51 s2_1.1.9
## [67] rlang_1.2.0 Rcpp_1.1.1-1.1 glue_1.8.1
## [70] DBI_1.3.0 rstudioapi_0.18.0 jsonlite_2.0.0
## [73] R6_2.6.1 spacesXYZ_1.6-0 units_1.0-1