dsn <- "https://raw.githubusercontent.com/OSGeo/gdal/refs/heads/master/frmts/wms/frmt_wms_bluemarble_s3_tms.xml"
library(terra)
r <- rast(dsn)
## first I did
m <- do.call(cbind, maps::map(plot = F)[1:2])
plot(project(m, to = "+proj=isea", from = "EPSG:4326"), pch = ".", asp = 1)
## which gives an idea of the overall extent par("usr") which looks a lot like c(-2, 2, -1, 1) *1e7
template <- rast(ext(c(-1.925, 1.925, -1, 1) *1e7), ncols = 1024, nrows = 512, crs = "+proj=isea")
pp <- project(r, template, by_util = TRUE)
plotRGB(pp)
with gdalwarp would do something like (set one of ts to "0" so GDAL works out the right aspect ratio, or use -tr)
not sure about setting nodata with cli
export dsn=https://raw.githubusercontent.com/OSGeo/gdal/refs/heads/master/frmts/wms/frmt_wms_bluemarble_s3_tms.xml
gdalwarp $dsn out.tif -te -1.925e7 -1e7 1.925e7 1e7 -ts 1024 0 -of COG -t_srs "+proj=isea"
(note that ext in terra is xmin,xmax,ymin,ymax but -te is xmin,ymin,xmax,ymax)