Acceso a series temporales de imágenes de Sentinel-2
Andres Felipe Ramirez Ramirez
1. Introducción
Este cuaderno ilustra cómo acceder a series temporales de imágenes Sentinel-2 utilizando los paquetes rstac y datacubes.
rstac ofrece funciones que simplifican la adquisición de datos espaciales de STAC (Catálogo de Activos Espaciotemporales). Más información aquí.
gdalcubes es un paquete de R y una biblioteca de C++ que facilita, agiliza y hace más interactivo el procesamiento de colecciones de imágenes satelitales. Más información aquí.
Utilizaremos el catálogo Sentinel-2 COG, disponible gratuitamente en Amazon Web Services (AWS), y el endpoint correspondiente de la API de STAC en https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a.
Además, utilizaremos terra para leer archivos ráster y tmap para la visualización.
2. Instalación de paquetes
En primer lugar, necesitamos instalar el paquete rstac. Utilizamos la consola para ejecutar el siguiente código, línea por línea:
#Utilizar la consola para correr este codigo
#install.packages("rstac")
#install.packages("tidyterra")
#install.packages("gdalcubes")
#install.packages("stars")
#install.packages("mapview")
#install.packages("tmap")3. Setup
# limpiamos el espacio de trabajo
rm(list = ls(all=TRUE))
# Cargar las librerias
library(rstac)
library(gdalcubes)
library(stars)
library(tmap)
library(tmaptools)
library(mapview)
library(leaflet)
library(sf)
library(ggplot2)4. Definir area de interes
Sigue estas instrucciones para crear un archivo GeoJSON con tu AOI:
• Ve a GeoJSON Studio.
• Busca el municipio de tu departamento con mayor producción agrícola.
• Cambia el fondo a Vista Satélite.
• Haz clic en Dibujar para crear un polígono con tu área de interés.
• Guarda el AOI como GeoJSON en tu directorio de datos.
En este cuaderno, el AOI es *Chinchiná en el departamento de Caldas.
Luego, lee el archivo GeoJSON.
# Utlizamos la función list.files() para comprobar el directorio y el nombre de los archivos
Chinchina <- st_read("Chinchina.geojson")## Reading layer `Chinchina' from data source
## `/Users/felipe/Documents/GB2/Proyecto5/Chinchina/Chinchina.geojson'
## using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -75.65196 ymin: 4.942417 xmax: -75.54635 ymax: 5.023619
## Geodetic CRS: WGS 845. Buscar imágenes de Sentinel-2
La función stac_search() permite a los usuarios buscar datos de Sentinel-2 a partir de un área de interés (AOI) con georreferencia EPSG:4326 y un período de tiempo. Probemos:
s = stac("https://earth-search.aws.element84.com/v0")
items <- s |>
stac_search(collections = "sentinel-s2-l2a-cogs",
bbox = c(st_bbox(Chinchina)["xmin"],
st_bbox(Chinchina)["ymin"],
st_bbox(Chinchina)["xmax"],
st_bbox(Chinchina)["ymax"]),
datetime = "2025-01-01/2025-12-31",
limit = 60) |>
post_request()
items## ###Items
## - matched feature(s): 59
## - features (59 item(s) / 0 not fetched):
## - S2B_18NVL_20251213_0_L2A
## - S2A_18NVL_20251210_0_L2A
## - S2C_18NVL_20251208_0_L2A
## - S2B_18NVL_20251203_0_L2A
## - S2A_18NVL_20251130_0_L2A
## - S2C_18NVL_20251128_0_L2A
## - S2B_18NVL_20251123_1_L2A
## - S2A_18NVL_20251120_0_L2A
## - S2C_18NVL_20251118_0_L2A
## - S2B_18NVL_20251113_0_L2A
## - ... with 49 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, typeImprimir la fecha y hora de la primera y la última imagen:
range(sapply(items$features, function(x) {x$properties$datetime}))## [1] "2025-06-01T15:32:08Z" "2025-12-13T15:31:46Z"6. Crear una colección de Sentinel-2
Podemos convertir la respuesta STAC en colecciones de imágenes
gdalcubes usando stac_image_collection(). Esta función
espera una lista de características STAC como entrada y, opcionalmente,
puede aplicar filtros a los metadatos y las bandas. Cabe destacar que
esta operación es mucho más rápida que la creación de colecciones de
imágenes a partir de archivos locales, ya que todos los metadatos están
disponibles y no es necesario abrir ningún archivo de imagen. A
continuación, creamos una colección para imágenes con menos del 50 % de
nubosidad.
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"]] < 25})Veamos qué tenemos:
s2_collection## Image collection object, referencing 4 images with 12 bands
## Images:
## name left top bottom right
## 1 S2C_18NVL_20251208_0_L2A -75.90301 5.42821 4.434363 -74.91191
## 2 S2A_18NVL_20250921_0_L2A -75.90301 5.42821 4.434363 -74.91191
## 3 S2C_18NVL_20250909_0_L2A -75.90301 5.42821 4.434363 -74.91191
## 4 S2C_18NVL_20250731_0_L2A -75.90301 5.42821 4.434363 -74.91191
## datetime srs
## 1 2025-12-08T15:31:56 EPSG:32618
## 2 2025-09-21T15:32:03 EPSG:32618
## 3 2025-09-09T15:32:03 EPSG:32618
## 4 2025-07-31T15:32:09 EPSG:32618
##
## Bands:
## name offset scale unit nodata image_count
## 1 B01 0 1 4
## 2 B02 0 1 4
## 3 B03 0 1 4
## 4 B04 0 1 4
## 5 B05 0 1 4
## 6 B06 0 1 4
## 7 B07 0 1 4
## 8 B08 0 1 4
## 9 B09 0 1 4
## 10 B11 0 1 4
## 11 B8A 0 1 4
## 12 SCL 0 1 4La colección contiene las imágenes que cumplen el criterio de nubosidad establecido.
7. Creación y procesamiento de cubos de datos
Con un objeto de colección de imágenes, podemos usar gdalcubes. En concreto, definimos una vista de cubo de datos y, posiblemente, algunas operaciones adicionales. En el siguiente ejemplo, creamos una imagen RGB compuesta de mediana simple con baja resolución (100 m) a partir de un cubo de datos mensual.
(Chinchina_box = st_bbox(Chinchina))## xmin ymin xmax ymax
## -75.651962 4.942417 -75.546346 5.023619(Chinchina_b9377 <- st_bbox(st_transform(Chinchina, "EPSG:9377")))## xmin ymin xmax ymax
## 4706043 2104713 4717723 2113699gdalcubes_options(parallel = 8)
(v = cube_view(srs="EPSG:9377", dx=10, dy=10, dt="P1M",
aggregation="median", resampling = "average",
extent=list(t0 = "2025-09-01", t1 = "2025-09-30",
left=Chinchina_b9377["xmin"], right=Chinchina_b9377["xmax"],
top=Chinchina_b9377["ymax"], bottom=Chinchina_b9377["ymin"])))## A data cube view object
##
## Dimensions:
## low high count pixel_size
## t 2025-09-01 2025-09-30 1 P1M
## y 2104710.73603885 2113700.73603885 899 10
## x 4706037.91083781 4717727.91083781 1169 10
##
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"8. Visualización del cubo de datos
Primero, el color natural (es decir, RGB 432):
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 (es decir, 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 un archivo netcdf o tif.
tutorialDir = path.expand("./Chinchina/")
if (!file.exists((file.path(tutorialDir,"Chinchina_10"))))
{
s2.mask = image_mask("SCL", values = c(3,8,9))
gdalcubes_options(parallel = 24, ncdf_compression_level = 0)
raster_cube(s2_collection, v, mask = s2.mask) |>
write_tif(file.path(tutorialDir,"Chinchina_10.tif"))
}Compruebe los nombres de los archivos descargados:
list.files("Chinchina/Chinchina_10.tif")## [1] "cube_16536628bbd972025-09-01.tif" "cube_176b320af62e92025-09-01.tif"
## [3] "cube_17a53631761522025-09-01.tif" "cube_17c4cd5ad4332025-09-01.tif"10. Calcula un índice de vegetación y visualízalo.
Primero, vamos a leer un archivo GTiff:
f <- "Chinchina/Chinchina_10.tif/cube_16536628bbd972025-09-01.tif"
Chinchina_s2 <- terra::rast(f)A continuación, calcule el índice NDVI:
# 1. Definimos las bandas
red <- Chinchina_s2[[4]] # Band 4
nir <- Chinchina_s2[[8]] # Band 8
# 2. Calculamos NDVI
ndvi <- (nir - red) / (nir + red)¿Qué tenemos?
ndvi## class : SpatRaster
## size : 899, 1169, 1 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : 4706038, 4717728, 2104711, 2113701 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## varname : cube_16536628bbd972025-09-01
## name : B08
## min value : -0.836578
## max value : 0.976065Ahora, es momento de visualizar. Comencemos con un gráfico sencillo:
# Mapa NDVI con una paleta continua de rojo, amarillo y verde
tm_shape(ndvi) +
tm_raster(
col_palette = "RdYlGn",
style = "cont",
title = "NDVI Value"
) +
tm_layout(main.title = "Sentinel-2 NDVI Sep. 2025")
Los resultados del cálculo del NDVI oscilan entre -1 y 1. Los valores negativos corresponden a zonas con agua, estructuras artificiales, rocas, nubes y nieve; el suelo desnudo suele estar entre 0,1 y 0,2; y las plantas siempre tendrán valores positivos entre 0,2 y 1.
Una cubierta vegetal sana y densa debería superar los 0,5, mientras que la vegetación escasa probablemente se sitúe entre 0,2 y 0,5. Sin embargo, se trata solo de una regla general y siempre debe tenerse en cuenta la estación del año, el tipo de planta y las particularidades regionales para comprender con exactitud el significado de los valores del NDVI.
Otra visualización puede ayudar a interpretar los valores NDVI:
# 2. Defina la paleta de colores solicitada: Rojo -> Amarillo -> Verde
# Ajustamos los valores del dominio para que coincidan con los límites teóricos del NDVI (-1 a 1).
ndvi_palette <- colorNumeric(
palette = c("red", "yellow", "green"),
domain = c(-1, 1),
na.color = "transparent" # Mantiene invisibles los píxeles con datos faltantes
)
# 3. Crea el mapa interactivo de Leaflet.
map <- leaflet() %>%
# Añadir la capa de imágenes satelitales de Esri de fondo
addProviderTiles("Esri.WorldImagery", group = "Esri Satellite") %>%
# Definir la vista del mapa: longitud, latitud y nivel de zoom específico
setView(lng = -76.0, lat = 4.50, zoom = 14) %>%
# Agregar la superposición de capa ráster NDVI personalizada
addRasterImage(
x = ndvi,
colors = ndvi_palette,
opacity = 0.6, # Ajusta la opacidad para que la imagen satelital subyacente sea parcialmente visible.
group = "NDVI Layer" # Identificador de grupo para el panel de control de capas
) %>%
# Agregar una leyenda de color visual estática al mapa
addLegend(
pal = ndvi_palette,
values = c(-1, 1),
title = "NDVI Index",
position = "bottomright",
opacity = 0.8
) %>%
# Agrega interruptores estructurales para que puedas activar o desactivar la capa del mapa NDVI.
addLayersControl(
baseGroups = c("Esri Satellite"),
overlayGroups = c("NDVI Layer"),
options = layersControlOptions(collapsed = FALSE)
)
# Display the map
mapTenga en cuenta que la fecha de la imagen satelital que se muestra en el fondo no coincide necesariamente con la fecha del NDVI.
11. Conclución
Se accedió exitosamente a imágenes Sentinel-2 mediante la API STAC utilizando los paquetes rstac y gdalcubes. Posteriormente, se construyó un cubo de datos para el área de interés correspondiente al municipio de Chinchiná (Caldas) y se generó una imagen compuesta del período de septiembre de 2025.
A partir de las bandas roja (B04) e infrarrojo cercano (B08) se calculó el índice de vegetación NDVI. Los resultados muestran predominio de valores positivos, lo que indica la presencia de cobertura vegetal en gran parte del área de estudio. Las zonas con valores más altos corresponden a vegetación más vigorosa, mientras que los valores bajos o negativos se asocian con áreas urbanas, suelo desnudo, cuerpos de agua o presencia de nubes.