rm(list = ls(all=TRUE))

1. Introduccion

Este cuaderno ilustra cómo acceder a series temporales de imágenes Sentinel-2 utilizando los paquetes rstac y gdalcubes.

rstac ofrece funciones que simplifican el proceso de adquisición de datos espaciales desde catálogos STAC (SpatioTemporal Asset Catalog).

gdalcubes es un paquete de R y una librería en C++ cuyo objetivo es facilitar, agilizar y hacer más interactivo el procesamiento de colecciones de imágenes satelitales.

Utilizaremos el catálogo de Sentinel-2 COG disponible gratuitamente en Amazon Web Services (AWS) y su correspondiente punto de acceso STAC-API.

Adicionalmente, utilizaremos terra para la lectura de archivos ráster y tmap junto con leaflet para la visualización.

2. y 3. Instalar y llamar librerias

Lalamamos los paquetes istalados

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

4. Define your area of interest (AOI)

En este cuaderno, el área de interés corresponde al municipio de Sandoná, Nariño.

A continuación, se realiza la lectura del archivo GeoJSON cargando la capa vectorial mediante la función st_read() ya con los datos descargados de GeoJson Studio

sandona <- st_read("C:/Users/Acer/Desktop/GB2R/Proyecto5/datos/sandona.geojson")
## Reading layer `sandona' from data source 
##   `C:\Users\Acer\Desktop\GB2R\Proyecto5\datos\sandona.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.50876 ymin: 1.264058 xmax: -77.44123 ymax: 1.306711
## Geodetic CRS:  WGS 84

5. Search Sentinel-2 imagery

La función stac_search() permite a los usuarios buscar datos de Sentinel-2 a partir de un área de interés (AOI) georreferenciada en el sistema EPSG:4326 y un periodo de tiempo determinado. Probremos su funcionamiento a continuación:

s = stac("https://earth-search.aws.element84.com/v0")
items <- s |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = c(st_bbox(sandona)["xmin"],
                       st_bbox(sandona)["ymin"],
                       st_bbox(sandona)["xmax"],
                       st_bbox(sandona)["ymax"]), 
              datetime = "2025-01-01/2025-12-31",
              limit = 60) |>
  post_request() 
items
## ###Items
## - matched feature(s): 103
## - features (60 item(s) / 43 not fetched):
##   - S2C_18NTG_20251211_0_L2A
##   - S2B_17NRB_20251206_0_L2A
##   - S2B_18NTG_20251206_0_L2A
##   - S2C_17NRB_20251201_0_L2A
##   - S2C_18NTG_20251201_0_L2A
##   - S2B_17NRB_20251126_0_L2A
##   - S2B_18NTG_20251126_0_L2A
##   - S2A_17NRB_20251123_0_L2A
##   - S2A_18NTG_20251123_1_L2A
##   - S2C_17NRB_20251121_0_L2A
##   - ... with 50 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

Imprimamos la fecha y hora de la primera y última imagen obtenidas:

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

6. Create a Sentinel-2 collection

Podemos convertir la respuesta de STAC en una colección de imágenes de gdalcubes utilizando la función stac_image_collection(). Esta función recibe como entrada una lista de elementos STAC y opcionalmente permite aplicar filtros sobre los metadatos y las bandas.

A diferencia de la creación de colecciones desde archivos locales, esta operación es mucho más rápida debido a que todos los metadatos se encuentran disponibles previamente y no requiere abrir ningún archivo de imagen. A continuación, creamos una colección con imágenes que presentan menos del 20% de cobertura de nubes.

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

Comprobemos lo que hemos obtenido:

s2_collection
## Image collection object, referencing 2 images with 12 bands
## Images:
##                       name      left      top   bottom     right
## 1 S2B_18NTG_20251206_0_L2A -77.69653 1.808999 0.815529 -76.70948
## 2 S2B_18NTG_20250927_0_L2A -77.69653 1.808999 0.816433 -76.70950
##              datetime        srs
## 1 2025-12-06T15:42:44 EPSG:32618
## 2 2025-09-27T15:42:48 EPSG:32618
## 
## Bands:
##    name offset scale unit nodata image_count
## 1   B01      0     1                       2
## 2   B02      0     1                       2
## 3   B03      0     1                       2
## 4   B04      0     1                       2
## 5   B05      0     1                       2
## 6   B06      0     1                       2
## 7   B07      0     1                       2
## 8   B08      0     1                       2
## 9   B09      0     1                       2
## 10  B11      0     1                       2
## 11  B8A      0     1                       2
## 12  SCL      0     1                       2

En este caso, tenemos 2 imágenes.

7. Crear y procesar cubos de datos

Una vez que disponemos de un objeto de colección de imágenes, podemos hacer uso de gdalcubes. En particular, definimos una vista del cubo de datos (cube view) y aplicamos las operaciones deseadas. En el siguiente ejemplo, creamos una imagen compósita RGB de mediana con una resolución gruesa de 100 metros a partir de un cubo de datos mensual.

(sandona_box = st_bbox(sandona))
##       xmin       ymin       xmax       ymax 
## -77.508757   1.264058 -77.441229   1.306711
(sandona_b9377 <- st_bbox(st_transform(sandona, "EPSG:9377")))
##    xmin    ymin    xmax    ymax 
## 4498097 1698146 4505622 1702861
gdalcubes_options(parallel = 8)
(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=sandona_b9377["xmin"], right=sandona_b9377["xmax"],
                                       top=sandona_b9377["ymax"], bottom=sandona_b9377["ymin"])))
## A data cube view object
## 
## Dimensions:
##                low             high count pixel_size
## t       2025-01-01       2025-12-31    12        P1M
## y 1698143.67743554 1702863.67743554   472         10
## x 4498094.35133148 4505624.35133148   753         10
## 
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"

8.Visualizing the data cube

Primero con color natural (i.e. 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 con falso color (i.e. 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. Download the data

Descargamos los arhivos para usarlos y verificamos lo que se descargó:

# Ruta absoluta hacia tu carpeta "datos"
tutorialDir = "C:/Users/Acer/Desktop/GB2R/Proyecto5/datos"

# Definimos la subcarpeta específica para Sandoná dentro de datos
outputFolder = file.path(tutorialDir, "sandona_10")

if (!file.exists(outputFolder))
{
  # Máscara para nubes (8, 9) y sombras (3)
  s2.mask = image_mask("SCL", values = c(3, 8, 9))
  
  # Hilos de procesamiento seguros y barra de progreso activa
  gdalcubes_options(parallel = 4, progress = TRUE)
  
  # Generamos y exportamos las imágenes
  raster_cube(s2_collection, v, mask = s2.mask) |> 
    write_tif(outputFolder)
}

# Verificar archivos creados dentro de la carpeta
list.files(outputFolder)
##  [1] "cube_23e4127a56152025-01-01.tif" "cube_23e4127a56152025-02-01.tif"
##  [3] "cube_23e4127a56152025-03-01.tif" "cube_23e4127a56152025-04-01.tif"
##  [5] "cube_23e4127a56152025-05-01.tif" "cube_23e4127a56152025-06-01.tif"
##  [7] "cube_23e4127a56152025-07-01.tif" "cube_23e4127a56152025-08-01.tif"
##  [9] "cube_23e4127a56152025-09-01.tif" "cube_23e4127a56152025-10-01.tif"
## [11] "cube_23e4127a56152025-11-01.tif" "cube_23e4127a56152025-12-01.tif"

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

En primer lugar, leemos un archivo GeoTIFF exportado previamente.

# First, read one GTiff file:
f <- file.path(outputFolder, "cube_23e4127a56152025-12-01.tif")
sandona_s2 <- terra::rast(f)

A continuación, calculamos el Índice de Vegetación de Diferencia Normalizada (NDVI) utilizando la banda del Rojo (B4) y del Infrarrojo Cercano (B8).

# Then, compute the NDVI index:

# 1. Define the bands
red <- sandona_s2[[4]] # Band 4
nir <- sandona_s2[[8]] # Band 8

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

Comprobemos las propiedades y estadísticas del ráster NDVI resultante:

ndvi
## class       : SpatRaster
## size        : 472, 753, 1  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : 4498094, 4505624, 1698144, 1702864  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s)   : memory
## varname     : cube_23e4127a56152025-12-01
## name        :       B08
## min value   : -0.279392
## max value   :  0.940999

Ahora es momento de la visualización. Comencemos con un mapa estático sencillo:

# Map NDVI with a continuous Red-Yellow-Green palette
tm_shape(ndvi) +
  tm_raster(
    col_palette = "RdYlGn",
    style = "cont",
    title = "NDVI Value"
  ) +
  tm_layout(main.title = "Sentinel-2 NDVI Dec. 2025")

Interpretación del NDVI para Sandoná (Diciembre 2025)

De acuerdo con el mapa obtenido y los rangos teóricos del NDVI (-1 a 1):

  • Zona Urbana y Construcciones: La cabecera municipal de Sandoná (visible en la zona central) presenta valores bajos/cercanos a 0 o ligeramente positivos (0.1 - 0.2), correspondientes a infraestructura y suelo desnudo.
  • Cobertura Vegetal: Las áreas circundantes presentan valores predominantemente altos (mayores a 0.5 y alcanzando hasta 0.8 - 0.94), lo que indica la presencia de vegetación densa y saludable (zonas agrícolas y vegetación nativa de la región).

Otra visualización puede ayudar a interpretar los valores de NDVI:

library(tmap)
library(terra)

# 1. Establecer el modo interactivo de tmap
tmap_mode("view")

# 2. Construir el mapa interactivo
mapa_ndvi <- tm_shape(ndvi) +
  tm_raster(
    style = "cont",
    palette = c("red", "yellow", "green"),
    title = "NDVI Index",
    alpha = 0.6
  ) +
  tm_basemap("Esri.WorldImagery")

# 3. Mostrar el mapa explícitamente en el cuaderno
mapa_ndvi