Preparación del Entorno

En esta primera sección preparamos nuestro espacio de trabajo en R. Es una buena práctica limpiar el entorno para evitar conflictos de memoria.

Paquetes clave utilizados: * rstac: Para conectarnos a la base de datos de AWS y buscar imágenes. * gdalcubes: El motor principal que procesará nuestras imágenes satelitales como un “cubo de datos” espacio-temporal. * terra y sf: Para el manejo de datos espaciales y reproyecciones. * leaflet y tmap: Para la visualización estática e interactiva de nuestros resultados.

#install.packages("rstac")
#install.packages("tidyterra")
#install.packages("gdalcubes")
#install.packages("stars")
#install.packages("mapview")
#install.packages("tmap")

# clean workspace
rm(list = ls(all=TRUE))
# load libraries
library(rstac)
library(gdalcubes)
library(stars)
library(tmap)
library(tmaptools)
library(mapview)
library(leaflet)
library(sf)
library(ggplot2)
library(terra)

Cargamos el polígono que delimita nuestra zona de interés (Argelia, Cauca).

Importante: Verifica siempre que la ruta del archivo .geojson

argelia <- st_read("C:\\Users\\pc\\OneDrive\\Documentos\\GB2\\RSTUDIO\\cuaderno6\\data\\argelia2.geojson")
## Reading layer `argelia2' from data source 
##   `C:\Users\pc\OneDrive\Documentos\GB2\RSTUDIO\cuaderno6\data\argelia2.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.25845 ymin: 2.271004 xmax: -77.21494 ymax: 2.330608
## Geodetic CRS:  WGS 84

Aquí definimos los parámetros de búsqueda de nuestras imágenes satelitales con un rango de tiempo de 2025-01-01/2025-12-31 .

s = stac("https://earth-search.aws.element84.com/v0")
items <- s |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = c(st_bbox(argelia)["xmin"],
                       st_bbox(argelia)["ymin"],
                       st_bbox(argelia)["xmax"],
                       st_bbox(argelia)["ymax"]), 
              datetime = "2025-01-01/2025-12-31",
              limit = 60) |>
  post_request() 
items
## ###Items
## - matched feature(s): 52
## - features (52 item(s) / 0 not fetched):
##   - S2C_18NTH_20251211_0_L2A
##   - S2B_18NTH_20251206_0_L2A
##   - S2C_18NTH_20251201_0_L2A
##   - S2B_18NTH_20251126_0_L2A
##   - S2A_18NTH_20251123_0_L2A
##   - S2C_18NTH_20251121_0_L2A
##   - S2B_18NTH_20251116_0_L2A
##   - S2C_18NTH_20251111_0_L2A
##   - S2B_18NTH_20251106_0_L2A
##   - S2C_18NTH_20251101_0_L2A
##   - ... with 42 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

Observamos la fecha de la primera y última captura registrada:

range(sapply(items$features, function(x) {x$properties$datetime}))
## [1] "2025-06-04T15:42:53Z" "2025-12-11T15:42:41Z"

Creación de la Colección y Filtro de Nubes

Convertimos los resultados STAC a un formato compatible con gdalcubes.

Filtro de Nubosidad: Hemos establecido eo:cloud_cover < 30. Esto significa que el código descartará cualquier imagen satelital donde más del 30% de la escena esté tapada por nubes. Es un balance ideal: permite entrar suficientes imágenes, pero descarta las tormentas severas.

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"]] < 30})
s2_collection
## Image collection object, referencing 5 images with 12 bands
## Images:
##                       name     left      top   bottom     right
## 1 S2C_18NTH_20250922_0_L2A -77.6982 2.712948 1.903651 -76.71015
## 2 S2A_18NTH_20250805_0_L2A -77.6982 2.712948 1.718886 -76.70998
## 3 S2A_18NTH_20250716_0_L2A -77.6982 2.712948 1.718886 -76.70998
## 4 S2C_18NTH_20250714_0_L2A -77.6982 2.712948 1.718886 -76.70998
## 5 S2A_18NTH_20250706_0_L2A -77.6982 2.712948 1.718886 -76.70998
##              datetime        srs
## 1 2025-09-22T15:42:45 EPSG:32618
## 2 2025-08-05T15:42:49 EPSG:32618
## 3 2025-07-16T15:42:52 EPSG:32618
## 4 2025-07-14T15:42:56 EPSG:32618
## 5 2025-07-06T15:42:51 EPSG:32618
## 
## Bands:
##    name offset scale unit nodata image_count
## 1   B01      0     1                       5
## 2   B02      0     1                       5
## 3   B03      0     1                       5
## 4   B04      0     1                       5
## 5   B05      0     1                       5
## 6   B06      0     1                       5
## 7   B07      0     1                       5
## 8   B08      0     1                       5
## 9   B09      0     1                       5
## 10  B11      0     1                       5
## 11  B8A      0     1                       5
## 12  SCL      0     1                       5

Cambiamos el sistema de coordenadas de nuestro polígono a Origen Nacional de Colombia (EPSG:9377). Esto nos permite trabajar con medidas en metros reales.

  (argelia_box = st_bbox(argelia))
##       xmin       ymin       xmax       ymax 
## -77.258448   2.271004 -77.214943   2.330608
  (argelia_b9377 <- st_bbox(st_transform(argelia, "EPSG:9377")))
##    xmin    ymin    xmax    ymax 
## 4526265 1809651 4531130 1816265

En este paso ya hemos creado el cubo DEM con las dimensiones necesarioas para el poligono, ajustremos la resolucion dimensiones de 1 pixel 10m x 10m

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=argelia_b9377["xmin"], right=argelia_b9377["xmax"],
                                       top=argelia_b9377["ymax"], bottom=argelia_b9377["ymin"])))
## A data cube view object
## 
## Dimensions:
##                low             high count pixel_size
## t       2025-01-01       2025-12-31    12        P1M
## y 1809647.94462168 1816267.94462168   662         10
## x 4526262.25155099 4531132.25155099   487         10
## 
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"

Mapa

Este bloque descarga, procesa y exporta la información en un archivo TIF.

Máscara de nubes: La capa SCL (valores 3, 8 y 9) borra las sombras y nubes altas de las imágenes.

parallel = 4: Disminuye la carga del procesador para evitar el colapso de los workers durante la descarga desde AWS.

select_bands y reduce_time: Obligamos a R a quedarse solo con la banda Roja (B04) e Infrarroja (B08) y a fusionar los 12 meses en una sola imagen perfecta, ahorrando espacio y evitando posibles errores en la empalme de la imagen satelital.

Ploteo 1

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

Ploteo 2

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

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

if (!file.exists((file.path(tutorialDir,"argelia_10")))) {
  
  s2.mask = image_mask("SCL", values = c(3,8,9))

  gdalcubes_options(parallel = 4, ncdf_compression_level = 0)
  
  raster_cube(s2_collection, v, mask = s2.mask) |>
  
    select_bands(c("B04", "B08")) |> 
    
    reduce_time(c("median(B04)", "median(B08)")) |> 
  
    write_tif(file.path(tutorialDir,"argelia_10.tif"))
}

Se exportara la unica imagen raster al equipo y se evaluara en QGIS la calidad de la imagen y los valores dados por el provedor.

f <- "C:\\Users\\pc\\OneDrive\\Documentos\\GB2\\RSTUDIO\\cuaderno6\\argelia\\argelia_10.tif\\cube_40d4b831b972025-01-01.tif"
argelia_s2 <- terra::rast(f)

Cargamos el archivo que acabamos de exportar. Dado que filtramos el cubo anteriormente, este archivo TIF ahora solo contiene 2 capas.

red <- argelia_s2[[1]] # Band 4
nir <- argelia_s2[[2]] # Band 8

ndvi <- (nir - red) / (nir + red)

Haremos una representacion grafica para vizualizar los valores del indice calculado

  tm_shape(ndvi) +
  tm_raster(
    col_palette = "RdYlGn", 
    style = "cont",
    title = "NDVI Value"
  ) +
  tm_layout(main.title = "Sentinel-2 NDVI Dec. 2025")

Visualización Interactiva del indice calculado

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 = -77.236695, lat = 2.300806, zoom = 14) %>%
  
  addRasterImage(
    x = ndvi, 
    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