get Digital Earth Africa L8/L9 GEOMAD with {vrtility}

vrtility-DEA-GEOMAD.R

# how to get at Digital Eath Africa GeoMAD datasets
library(vrtility)

mirai::daemons(10)

bbox <- gdalraster::bbox_from_wkt(
  wkt = "POINT (46.33 -15.9)",
  extend_x = 0.15,
  extend_y = 0.2
)

l8l9_query <- stac_query(
  bbox = bbox,
  stac_source = "https://explorer.digitalearth.africa/stac/",
  collection = "gm_ls8_ls9_annual",
  start_date = "2024-01-01",
  end_date = "2024-12-31"
)

# if you only want selected bands you can filter the assets.
# l8l9_query_rgb <- rstac::assets_select(
#   l8l9_query,
#   asset_names = c("SR_B2", "SR_B3", "SR_B4")
# )

# Doing this once saves you having to add it to every vrt_x command.
Sys.setenv(
  AWS_NO_SIGN_REQUEST = "YES",
  AWS_S3_ENDPOINT = "s3.af-south-1.amazonaws.com"
)

l8l9_vrt <- vrt_collect(
  l8l9_query
)

# if you just want to look at one of the tiles...
# plot(l8l9_vrt[[1]][[2]], c(3, 2, 1))

bbox_proj <- bbox_to_projected(bbox, l8l9_vrt$srs)

l8l9_vrt_warped <- vrt_warp(
  l8l9_vrt,
  t_srs = attr(bbox_proj, "wkt"),
  te = bbox_proj,
  tr = c(30, 30)
) |>
  vrt_stack()


l8l9file <- vrt_compute(
  l8l9_vrt_warped,
  outfile = fs::file_temp(ext = "tif"),
  engine = "gdalraster"
)

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