# 1. CARGAR LAS LIBRERIAS.
library(rstac)
## Warning: package 'rstac' was built under R version 4.5.3
library(gdalcubes)
## Warning: package 'gdalcubes' was built under R version 4.5.3
library(stars)
## Warning: package 'stars' was built under R version 4.5.3
## Cargando paquete requerido: abind
## Cargando paquete requerido: sf
## Warning: package 'sf' was built under R version 4.5.3
## 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.5.3
library(tmaptools)
## Warning: package 'tmaptools' was built under R version 4.5.3
library(mapview)
## Warning: package 'mapview' was built under R version 4.5.3
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.5.3
library(sf)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
# 2. DEFINICIÓN DEL ÁREA DE INTERÉS (AOI).
#Cargo mi zona de estudio (la Ciénaga). Lo tengo en un archivo geojson que preparé antes.
union <- st_read("C:\\Users\\suare\\OneDrive\\Documentos\\GB2\\proyecto5\\cienaga.geojson")
## Reading layer `cienaga' from data source
## `C:\Users\suare\OneDrive\Documentos\GB2\proyecto5\cienaga.geojson'
## using driver `GeoJSON'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.33083 ymin: 10.89692 xmax: -74.13611 ymax: 11.03113
## Geodetic CRS: WGS 84
# 3. BÚSQUEDA EN EL CATÁLOGO STAC
#Me conecto al catálogo de Earth Search para buscar imágenes Sentinel-2.
s = stac("https://earth-search.aws.element84.com/v0")
# Hago la consulta pasándole la caja delimitadora (bounding box) de la Ciénaga,durante el 2025.
items <- s |>
stac_search(collections = "sentinel-s2-l2a-cogs",
bbox = c(st_bbox(union)["xmin"],
st_bbox(union)["ymin"],
st_bbox(union)["xmax"],
st_bbox(union)["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_18PWT_20251213_0_L2A
## - S2C_18PWT_20251208_0_L2A
## - S2B_18PWT_20251203_0_L2A
## - S2A_18PWT_20251130_0_L2A
## - S2C_18PWT_20251128_0_L2A
## - S2B_18PWT_20251123_0_L2A
## - S2A_18PWT_20251120_0_L2A
## - S2C_18PWT_20251118_0_L2A
## - S2B_18PWT_20251113_0_L2A
## - S2A_18PWT_20251110_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, type
range(sapply(items$features, function(x) {x$properties$datetime}))
## [1] "2025-06-01T15:30:23Z" "2025-12-13T15:30:00Z"
# 4. FILTRADO DE BANDAS Y NUBES.
# Selecciono solo las bandas que necesito.
assets = c("B01","B02","B03","B04","B05","B06", "B07","B08","B8A","B09","B11","SCL")
#Armo la colección, pero filtro para quedarme únicamente con imágenes que tengan menos del 20% de nubes.
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 12 images with 12 bands
## Images:
## name left top bottom right
## 1 S2C_18PWT_20251208_0_L2A -74.75914 11.75966 10.76541 -73.99249
## 2 S2B_18PWT_20251113_0_L2A -74.74231 11.75964 10.76541 -73.99249
## 3 S2C_18PWT_20251108_0_L2A -74.75529 11.75966 10.76541 -73.99249
## 4 S2C_18PWT_20250929_0_L2A -74.75612 11.75966 10.76541 -73.99249
## 5 S2A_18PWT_20250911_0_L2A -74.74295 11.75964 10.76541 -73.99249
## 6 S2C_18PWT_20250909_0_L2A -74.74587 11.75964 10.76541 -73.99249
## datetime srs
## 1 2025-12-08T15:30:11 EPSG:32618
## 2 2025-11-13T15:29:59 EPSG:32618
## 3 2025-11-08T15:30:14 EPSG:32618
## 4 2025-09-29T15:30:18 EPSG:32618
## 5 2025-09-11T15:30:15 EPSG:32618
## 6 2025-09-09T15:30:18 EPSG:32618
## [ omitted 6 images ]
##
## Bands:
## name offset scale unit nodata image_count
## 1 B01 0 1 12
## 2 B02 0 1 12
## 3 B03 0 1 12
## 4 B04 0 1 12
## 5 B05 0 1 12
## 6 B06 0 1 12
## 7 B07 0 1 12
## 8 B08 0 1 12
## 9 B09 0 1 12
## 10 B11 0 1 12
## 11 B8A 0 1 12
## 12 SCL 0 1 12
# 5. ESTRUCTURACIÓN DEL CUBO DE DATOS
(union_box = st_bbox(union))
## xmin ymin xmax ymax
## -74.33083 10.89692 -74.13611 11.03113
# Reproyecto el polígono al Origen Nacional (EPSG:9377)
(union_b9377 <- st_bbox(st_transform(union, "EPSG:9377")))
## xmin ymin xmax ymax
## 4854676 2762382 4875941 2777306
#Defino las reglas de mi cubo de datos.
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=union_b9377["xmin"], right=union_b9377["xmax"],
top=union_b9377["ymax"], bottom=union_b9377["ymin"])))
## A data cube view object
##
## Dimensions:
## low high count pixel_size
## t 2025-01-01 2025-12-31 12 P1M
## y 2762379.05425008 2777309.05425008 1493 10
## x 4854673.6878443 4875943.6878443 2127 10
##
## SRS: "EPSG:9377"
## Temporal aggregation method: "median"
## Spatial resampling method: "average"
# 6. VISUALIZACIÓN: COLOR VERDADERO
#Ploteo la composición RGB normal (Rojo=B04, Verde=B03, Azul=B02) para ver la Ciénaga tal cual ve el ojo humano.
gdalcubes_options(threads = 2)
## Warning in gdalcubes_options(threads = 2): 'threads' option is deprecated;
## please use 'parallel' instead
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))
#Al armar la imagen en color verdadero (mezclando las bandas 4, 3 y 2), el resultado es básicamente lo que veríamos a simple vista si estuviéramos sobrevolando la zona. Observando el mapa de la Ciénaga, lo primero que noto es que los cuerpos de agua se ven súper oscuros, entre un azul muy profundo y negro, lo cual tiene sentido por la profundidad y la cantidad de sedimentos que maneja el lugar. Por otro lado, toda la zona de manglares y vegetación es muy fácil de identificar porque resalta en diferentes tonos de verde. Finalmente, esos parches más grises, blanquitos o color arena corresponden a los suelos desnudos, los playones o las áreas urbanas que están alrededor.
# 7. VISUALIZACIÓN: FALSO COLOR.
#Ahora hago una composición de falso color usando infrarrojo (B08, B11, B04).Con esto diferencio muchísimo mejor los cuerpos de agua de los manglares y vegetación.
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))
#Cuando se cambia a la composición en falso color (usando las bandas B08, B11 y B04), la imagen dio un giro total. Como estoy usando las bandas infrarrojas, la vegetación más densa y saludable, como los bosques de manglar, ahora 'brilla' en tonos anaranjados. Esto pasa porque las hojas reflejan muchísimo esa luz infrarroja que el ojo huemano no ve. Lo que más me gusta de esta vista es lo fácil que es separar la tierra del agua; como el agua no refleja casi nada de infrarrojo , la ciénaga y los humedales quedan casi negros o azules oscuros, logrando un contraste súper marcado que en el color verdadero era mucho más difícil de notar.
# 8. LIMPIEZA FINAL Y EXPORTACIÓN
tutorialDir = path.expand("C:\\Users\\suare\\OneDrive\\Documentos\\GB2\\proyecto561")
if (!file.exists((file.path(tutorialDir,"cienaga_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,"cienaga_10.tif"))
}
#¿Qué aprendí hoy?
#1. Sobre el color verdadero (RGB 4-3-2): Me di cuenta de que al mezclar las bandas roja, verde y azul vemos el mapa tal cual como si estuviéramos ahí. Es súper útil y fácil de entender, pero también noté que a veces tiene sus límites; por ejemplo, si hay un humedal muy oscuro, se confunde fácilmente con un cuerpo de agua normal porque los tonos se ven casi iguales.
#2. Sobre el falso color: Aprendí que usar las bandas infrarrojas (como la B08 y B11) es increíble porque nos deja ver cosas que el ojo humano normalmente se perdería. Me pareció muy práctico cómo el infrarrojo hace que la vegetación sana "grite" en color naranja, y cómo nos ayuda a delimitar los bordes del agua con mucha exactitud porque absorbe la luz y queda oscura.
#3.Sobre la resolución de Sentinel-2: Algo nuevo para mí fue entender que Sentinel-2 no toma toda la imagen con el mismo nivel de detalle. Maneja distintas resoluciones: tiene bandas muy nítidas a 10 metros (que son geniales para el color y el infrarrojo cercano), otras a 20 metros que sirven más para ver detalles agrícolas o de humedad, y unas menos detalladas a 60 metros que se usan más que todo para corregir cosas de la atmósfera, como las nubes.