Copernicus Dataspace Ecosystem (CDSE) Access with vrtility

s2-cdse-vrtility.R

library(vrtility)

# only use 4 daemons - this is limit for concurrent requests: https://documentation.dataspace.copernicus.eu/Quotas.html
mirai::daemons(4) 

Sys.setenv(
  AWS_VIRTUAL_HOSTING = "FALSE",
  AWS_ACCESS_KEY_ID = Sys.getenv("CDSE_ACCESS_KEY"),
  AWS_SECRET_ACCESS_KEY = Sys.getenv("CDSE_SECRET_KEY"),
  AWS_S3_ENDPOINT = "eodata.dataspace.copernicus.eu"
)

bbox <- gdalraster::bbox_from_wkt(
  wkt = "POINT (95.415 27.78)",
  extend_x = 0.3,
  extend_y = 0.2
)
bbx_proj <- bbox_to_projected(bbox)

# run the STAC query
s2copdse <- stac_query(
  bbox = bbox,
  stac_source = "https://stac.dataspace.copernicus.eu/v1",
  collection = "sentinel-2-l2a",
  start_date = "2025-06-01",
  end_date = "2025-06-30"
) |>
  rstac::assets_select(
    asset_names = c("B02_10m", "B03_10m", "B04_10m", "B08_10m", "SCL_20m")
  ) |>
  stac_cloud_filter(max_cloud_cover = 30)

# Download the data and process
s2_median <- vrt_collect(s2copdse) |>
  vrt_move_band(1, 5) |>
  vrt_set_nodata(0, band_idx = 1:4) |>
  vrt_set_scale(scale_value = 0.0001, offset_value = -0.1) |>
  vrt_set_maskfun(
    mask_band = "SCL_20m",
    mask_values = c(0, 1, 2, 3, 8, 9, 10, 11),
    drop_mask_band = TRUE
  ) |>
  vrt_warp(
    t_srs = attr(bbx_proj, "wkt"),
    te = bbx_proj,
    tr = c(10, 10)
  ) |>
  vrt_stack() |>
  vrt_set_py_pixelfun(median_numpy()) |>
  vrt_compute(
    outfile = fs::file_temp(ext = "tif"),
    engine = "gdalraster"
  )

plot_raster_src(s2_median, bands = 4)

plot_raster_src(s2_median, bands = c(3, 2, 1))