Combine multiple multi-band rasters into a single VRT

stack_rasters.R

stack_rasters <- function(rasters,
                          raster_path,
                          ...,
                          reference_layer = 1,
                          resampling_method = "bilinear") {

  check_type_and_length(
    reference_layer = integer(1)
  )

  out_dir <- dirname(raster_path)
  var_names <- unlist(lapply(rasters, \(r) terra::rast(r) |> names()))

  ref_rast <- terra::rast(rasters[reference_layer])
  ref_ext <- (terra::ext(ref_rast) |>
                as.vector())[c("xmin", "ymin", "xmax", "ymax")] |>
    setNames(NULL)
  ref_res <- terra::res(ref_rast)

  vrt_container_file <- tempfile(tmpdir = out_dir, fileext = ".vrt")
  sf::gdal_utils(
    "buildvrt",
    rasters[reference_layer],
    vrt_container_file,
    options = c(
      "-b", 1,
      "-te", ref_ext,
      "-tr", ref_res,
      "-r", resampling_method
    )
  )
  vrt_container <- readLines(vrt_container_file)
  file.remove(vrt_container_file)

  intermediate_vrt <- lapply(
    rasters,
    \(raster) {
      r_lyrs <- terra::nlyr(terra::rast(raster))
      out_files <- replicate(
        r_lyrs,
        tempfile(tmpdir = out_dir, fileext = ".vrt")
      )

      lapply(
        seq_len(r_lyrs),
        \(b) {
          sf::gdal_utils(
            "buildvrt",
            raster,
            out_files[[b]],
            options = c(
              "-b", b,
              "-te", ref_ext,
              "-tr", ref_res,
              "-r", resampling_method
            )
          )
        }
      )
      out_files
    }
  )

  intermediate_vrt <- unlist(intermediate_vrt)

  band_no <- 1
  vrt_bands <- lapply(
    intermediate_vrt,
    \(vrt) {
      vrt <- readLines(vrt)
      band_def <- grep("VRTRasterBand", vrt)
      vrt <- vrt[seq(band_def[[1]], band_def[[2]])]
      vrt[1] <- gsub("band=\"1\"", paste0("band=\"", band_no, "\""), vrt[1])
      vrt <- c(
        vrt[1],
        paste0("    <Description>", var_names[[band_no]], "</Description>"),
        vrt[2:length(vrt)]
      )

      band_no <<- band_no + 1

      vrt
    }
  )

  band_def <- grep("VRTRasterBand", vrt_container)
  writeLines(
    c(
      vrt_container[1:(band_def[[1]] - 1)],
      unlist(vrt_bands),
      vrt_container[(band_def[[2]] + 1):length(vrt_container)]
    ),
    raster_path
  )

  invisible(file.remove(intermediate_vrt))
  raster_path
}