GDAL

Generally there is not much difference between GDAL and {terra} (there is only very little overhead). However, using gdal_translate directly (sf::gdal_utils()) can be beneficial for data merging (tiles) or resolution changes (downsampling / upsampling).

Code

library("sf")
library("terra")

Project

n = 3000 * 3
r = rast(nrows = n, ncols = n, nlyrs = 3, vals = runif(n ^ 2 * 3))
input = tempfile(fileext  = ".tif")
writeRaster(r, input)
output_terra = tempfile(fileext  = ".tif")
output_gdal = tempfile(fileext  = ".tif")

system.time({
  r = rast(input)
  project(r, "EPSG:3857", method = "bilinear", gdal = TRUE,
          filename = output_terra, overwrite = TRUE)
})

##    user  system elapsed 
##   10.90    0.73   11.70

system.time({
  gdal_utils(util = "warp",
             source = input,
             destination = output_gdal,
             options = c("-t_srs", "EPSG:3857",
                         "-r", "bilinear",
                         "-co", "COMPRESS=LZW",
                         "-overwrite")
             )
})

##    user  system elapsed 
##   10.17    0.75   10.92

Crop by extent

v = ext(-90, 90, -45, 45)

system.time({
  r = rast(input)
  crop(r, v, filename = output_terra, overwrite = TRUE)
})

##    user  system elapsed 
##    4.56    0.43    4.98

ext = c("-90", "45", "90", "-45")

system.time({
  gdal_utils(util = "translate",
             source = input,
             destination = output_gdal,
             options = c("-projwin", ext,
                         "-co", "COMPRESS=LZW")
             )
})

##    user  system elapsed 
##    4.28    0.39    4.68

Downsample

system.time({
  r = rast(input)
  dest = r
  res(dest) = 0.1
  resample(r, dest, method = "cubic", filename = output_terra, overwrite = TRUE)
})

##    user  system elapsed 
##   13.50    0.72   14.25

system.time({
  gdal_utils(util = "translate",
             source = input,
             destination = output_gdal,
             options = c("-tr", 0.1, 0.1,
                         "-r", "cubic",
                         "-co", "COMPRESS=LZW")
             )
})

##    user  system elapsed 
##    6.62    0.92    7.54

Merging files

r = rast(ncols = 3000, nrows = 3000)
values(r) = rnorm(ncell(r))
filename = tempfile(fileext = "_.tif")
ff = makeTiles(r, 150, filename)
paste("Number of tiles:", length(ff))

## [1] "Number of tiles: 400"

system.time({
  vrt = vrt(ff)
  writeRaster(vrt, output_terra, overwrite = TRUE)
})

##    user  system elapsed 
##    1.54    0.93    3.03

system.time({
  vrt = tempfile(fileext = ".vrt")
  gdal_utils(util = "buildvrt",
             source = ff,
             destination = vrt)
  gdal_utils(util = "translate",
             source = vrt,
             destination = output_gdal,
             options = c("-co", "COMPRESS=LZW")
             )
})

##    user  system elapsed 
##    1.26    0.46    1.85