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
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
s = stac("https://earth-search.aws.element84.com/v0")
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"
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 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
(union_box = st_bbox(union))
## xmin ymin xmax ymax
## -74.33083 10.89692 -74.13611 11.03113
(union_b9377 <- st_bbox(st_transform(union, "EPSG:9377")))
## xmin ymin xmax ymax
## 4854676 2762382 4875941 2777306
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"
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))
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("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"))
}