Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active April 19, 2025 21:33
Show Gist options
  • Save mdsumner/b7f35d3af366332806e347dc68b14b20 to your computer and use it in GitHub Desktop.
Save mdsumner/b7f35d3af366332806e347dc68b14b20 to your computer and use it in GitHub Desktop.
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)

@mdsumner
Copy link
Author

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment