1. Introducción

En este cuaderno veremos cómo acceder a imágenes Sentinel-2. Con la ayuda de la librería preprocS2 podemos acceder a los Cátálogos de Activos Espaciales y Temporales (STAC) cuya colección “sentinel-2-l2a es del catálogo de Microsoft Planetary Computers, utilizamos también paquetes como rstac, bidiovMapR y spinR:

2. Instalación de paquetes y configuración

Instalamos las librerías previamente mencionadas, las cargamos y limpiamos la memoria:

rm(list = ls(all=TRUE))
library(rstac)
## Warning: package 'rstac' was built under R version 4.6.1
library(gdalcubes)
## Warning: package 'gdalcubes' was built under R version 4.6.1
library(stars)
## Warning: package 'stars' was built under R version 4.6.1
## Cargando paquete requerido: abind
## Cargando paquete requerido: sf
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(tmap)
## Warning: package 'tmap' was built under R version 4.6.1
library(tmaptools)
## Warning: package 'tmaptools' was built under R version 4.6.1
library(mapview)
## Warning: package 'mapview' was built under R version 4.6.1
library(leaflet)
library(sf)
library(ggplot2)
library(preprocS2)
## Warning: replacing previous import 'kernlab::predict' by 'stats::predict' when
## loading 'prosail'
library(biodivMapR)
## 
## Adjuntando el paquete: 'biodivMapR'
## The following objects are masked from 'package:preprocS2':
## 
##     get_si_tiles_from_raster, mosaic_tiles, split_line
library(spinR)
## 
## Adjuntando el paquete: 'spinR'
## The following object is masked from 'package:biodivMapR':
## 
##     IQR_outliers

3. Definimos el área de interes, AOI (area of interest)

list.files()
## [1] "aoi_bbox.GPKG"     "invg.nb.html"      "rsconnect"        
## [4] "sentinel2.html"    "sentinel2.nb.html" "sentinel2.Rmd"    
## [7] "sentinel2_files"   "velez"             "velez.geojson"

Leemos el archivo GeoJSON:

velez <- st_read("velez.geojson")
## Reading layer `velez' from data source 
##   `C:\Users\isabe\OneDrive\Documentos\GB2\Proyecto5\velez.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -73.69009 ymin: 5.992513 xmax: -73.65366 ymax: 6.026015
## Geodetic CRS:  WGS 84

5. Buscamos las imágenes de Sentinel- 2

La función stac_search le permite a los usuarios buscar datos de Sentinel-2 avpartir de una AOI, una georreferencia EPSG:4326 y un intervalo de tiempo:

s = stac("https://earth-search.aws.element84.com/v0")
items <- s |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = c(st_bbox(velez)["xmin"],
                       st_bbox(velez)["ymin"],
                       st_bbox(velez)["xmax"],
                       st_bbox(velez)["ymax"]), 
              datetime = "2025-07-31/2025-08-31",
              limit = 60) |>
  post_request() 
items
## ###Items
## - matched feature(s): 10
## - features (10 item(s) / 0 not fetched):
##   - S2C_18NXM_20250830_0_L2A
##   - S2B_18NXM_20250825_0_L2A
##   - S2A_18NXM_20250822_1_L2A
##   - S2C_18NXM_20250820_0_L2A
##   - S2B_18NXM_20250815_1_L2A
##   - S2B_18NXM_20250815_0_L2A
##   - S2A_18NXM_20250812_0_L2A
##   - S2B_18NXM_20250805_0_L2A
##   - S2A_18NXM_20250802_1_L2A
##   - S2C_18NXM_20250731_0_L2A
## - 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

Ahora, vamos a ver la fecha y hora de la súltimas imágenes:

range(sapply(items$features, function(x) {x$properties$datetime}))
## [1] "2025-07-31T15:31:48Z" "2025-08-30T15:31:45Z"

5. Creando una colección de Sentinel-2

Convertimos la respuesta STAC en una colección de imágenes de gdalcubes con la función stac_image_collection(),se le pueden aplicar filtros a los metadatos y las bandas. Esto es más rápido que la creación de colecciones a partir de archivos locales, pondremos a demás una nubosidad del 20%:

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 1 images with 12 bands
## Images:
##                       name      left      top   bottom     right
## 1 S2C_18NXM_20250820_0_L2A -74.09746 6.332246 5.338704 -73.10435
##              datetime        srs
## 1 2025-08-20T15:31:46 EPSG:32618
## 
## Bands:
##    name offset scale unit nodata image_count
## 1   B01      0     1                       1
## 2   B02      0     1                       1
## 3   B03      0     1                       1
## 4   B04      0     1                       1
## 5   B05      0     1                       1
## 6   B06      0     1                       1
## 7   B07      0     1                       1
## 8   B08      0     1                       1
## 9   B09      0     1                       1
## 10  B11      0     1                       1
## 11  B8A      0     1                       1
## 12  SCL      0     1                       1

6. Creando y procesando cubos de datos

Con el objeto de colección de imágenes podemos utilizar gdalcubes, vamos a crear una imagen RGB con resolución de 100m usando la técina de la mediana compuesta a partir de un cubo de datos:

(velez_box = st_bbox(velez))
##       xmin       ymin       xmax       ymax 
## -73.690090   5.992513 -73.653661   6.026015
(velez_b9377 <- st_bbox(st_transform(velez, "EPSG:9377")))
##    xmin    ymin    xmax    ymax 
## 4923656 2220205 4927690 2223911
gdalcubes_options(parallel = 8)
(v = cube_view(srs="EPSG:9377", dx=10, dy=10, dt="P1M", 
                           aggregation="median", resampling = "average",
                           extent=list(t0 = "2025-07-31", t1 = "2025-08-31",
                                       left=velez_b9377["xmin"], right=velez_b9377["xmax"],
                                       top=velez_b9377["ymax"], bottom=velez_b9377["ymin"])))
## A data cube view object
## 
## Dimensions:
##                low             high count pixel_size
## t       2025-07-01       2025-08-31     2        P1M
## y 2220202.99798432 2223912.99798432   371         10
## x 4923653.25192526 4927693.25192526   404         10
## 
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"

7. Visualizando el cubo de datos

Primero con una composición de color natural:

cube<- raster_cube(s2_collection,v)
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:

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. Descargando los datos

tutorialDir = path.expand("./velez/")

if (!file.exists((file.path(tutorialDir,"velez_10"))))

{

  s2.mask = image_mask("SCL", values = c(3,8,9))

  gdalcubes_options(parallel = 2, ncdf_compression_level = 0)

  raster_cube(s2_collection, v, mask = s2.mask) |>

    write_tif(file.path(tutorialDir,"velez_10.tif"))

}

Revisamos los nombres de los archivos descargados:

list.files("velez/velez_10")
## [1] "cube_4e483b7148ec2025-07-01.tif" "cube_4e483b7148ec2025-08-01.tif"

9. Calculando el indice de vegetación y visualizándolo

Primero leemos uno de los archivos GTiff:

f <- "velez/velez_10/cube_4e483b7148ec2025-08-01.tif"
velez_s2 <- terra::rast(f)

Luego, calculamos el NDVI:

red <- velez_s2[[4]]
nir <- velez_s2[[8]] 
ndvi <- (nir - red) / (nir + red)
ndvi
## class       : SpatRaster
## size        : 3672, 3679, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 599811.6, 710181.6, 590165.9, 700325.9  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s)   : memory
## varname     : cube_4e483b7148ec2025-08-01
## name        :       B08
## min value   : -0.602484
## max value   :  0.996377

Ahora lo visualizamos:

tm_shape(ndvi) +
  tm_raster(
    col_palette = "RdYlGn", 
    style = "cont",
    title = "NDVI Value"
  ) +
  tm_layout(main.title = "Sentinel-2 NDVI Jul. 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 = )`
## SpatRaster object downsampled to 3160 by 3166 cells.
## 
## [`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.

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

ndvi_raster <- raster::raster(ndvi)
ndvi_pequeno <- terra::aggregate(ndvi, fact = 2, fun = mean)
## |---------|---------|---------|---------|=========================================                                          
ndvi_palette <- colorNumeric(
  palette = c("red", "yellow", "green"), 
  domain = c(-1, 1),
  na.color = "transparent"
)
options(leaflet.maxBytes =  20000)
map <- leaflet() %>%
  addProviderTiles("Esri.WorldImagery", group = "Esri Satellite") %>%
  setView(lng = -73.9, lat = 6.0, zoom = 14) %>%
  addRasterImage(
    x = ndvi_pequeno, 
    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)
  )
map