Here get raw tiles and calculate median across time. We're only plotting in a faceted index space, but it's lined up relatively.

library(rustycogs)  ## hypertidy/rustycogs  - we need dev-version of async_tiff which is in the Cargo.toml (#PR 265)
library(ximage)
js <- jsonlite::fromJSON(sds::stacit("55GDN", c("2025-01-01", format(Sys.Date()))))  ## hypertidy/sds relies on recent MGRS grid code rather than extent
hrf <- unlist(lapply(js$features$assets[1:19], function(x) x[["href"]]))
all(grepl("tif$",hrf))
                     
## index all the tifs (like virtualizarr, we get the IFD tables)
system.time(refs <- tiff_refs(hrf, region = "", anon = TRUE, concurrency = 120))   

##calculate median after reading all the data for a given tile across time steps
calc_med <- function(x0, tile = c(0L, 0L), concurrency = 180L) {
  tile00 <- x0 |> dplyr::filter(tile_col == tile[1L], tile_row == tile[2L]) 
  tilesize <- prod(unlist(tile00[1L, c("tile_w", "tile_h")]))
 dat <- array(unlist(rustycogs::tiff_read_tiles(tile00, region = "", anon = TRUE, concurrency = concurrency)$data), 
               dim = c(tilesize, dim(tile00)[1L]))
  dat[(dat == 0) | (dat > 5000)] <- NA
  matrixStats::rowMedians(dat, na.rm = TRUE)
}
ifd_1 <- 2L
x1 <- refs |>  dplyr::filter(ifd == ifd_1, path == path[1])

## clean up for ximage plot down below
setna <- function(x) {x[is.na(x)] <- 0; x}

## loop calc and plot at nominated IFD level                      
par(mfrow = n2mfrow(nrow(x1)), mar = rep(0, 4) ,bg = "black") 
tile = c(0L, 0L);  concurrency = 180L; i <- 1; j <- 1
for (j in seq_len(max(x1$tile_row) + 1)) {
  for (i in seq_len(max(x1$tile_col) +  1)) {
    tile <- c(i - 1, j - 1)
  red <- calc_med(dplyr::filter(refs, ifd == ifd_1, grepl("B04\\.tif$", path)), tile = tile)
  grn <- calc_med(dplyr::filter(refs, ifd == ifd_1, grepl("B03\\.tif$", path)), tile = tile)
  blu <- calc_med(dplyr::filter(refs, ifd == ifd_1, grepl("B02\\.tif$", path)), tile = tile)
  ximage(scales::rescale(aperm(array(setna(c(red, grn, blu)), c(x1$tile_w[1], x1$tile_h[1], 3L)), c(2, 1, 3))), asp = 1)
  box(col = "white")
}}))


dim(refs)
#[1] 1094613      23
#length(hrf)
#[1] 4199

dplyr::filter(refs, ifd == ifd_1, grepl("B04\\.tif$", path)) |> dplyr::distinct(path) |>  dplyr::tally()
#    n
#221