library(gdalraster)
library(ximage) ## remotes::install_github("hypertidy/ximage")
dsn <- "/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/USGS_Seamless_DEM_1.vrt"
## we want a grid of a local area
ex <- c(-2, 2, -1, 1) * 85000
args <- c("-te", ex[1], ex[3], ex[2], ex[4], "-ts", 1024, 0)
crs <- "+proj=laea +lon_0=-112.6 +lat_0=36.3"
tf <- tempfile(fileext = ".vrt", tmpdir = "/vsimem")
gdalraster::warp(dsn, tf, crs, cl_arg = args)
ds <- new(GDALRaster, tf)
## gdal raster has a read_ds()
dat1 <- read_ds(ds)
## now mars
dsn <- "/vsicurl/https://planetarymaps.usgs.gov/mosaic/Mars/HRSC_MOLA_Blend/Mars_HRSC_MOLA_BlendDEM_Global_200mp_v2.tif"
## we want a grid of a local area
ex <- c(-2, 2, -1, 1) * 1400000
args <- c("-te", ex[1], ex[3], ex[2], ex[4], "-ts", 1024, 0)
crs <- "+proj=laea +lon_0=-52.2 +lat_0=-8.9 +R=3396190"
tf <- tempfile(fileext = ".vrt", tmpdir = "/vsimem")
gdalraster::warp(dsn, tf, crs, cl_arg = args)
ds <- new(GDALRaster, tf)
## gdal raster has a read_ds()
dat <- read_ds(ds)
par(mfrow = c(2, 1), mar = rep(0, 4))
ximage::ximage(dat1, asp = 1)
ximage::ximage(dat, asp = 1)
bb <- attr(dat1, "gis")$bbox + c(-1e6, 1e6, -1e6, 1e6)
rect(bb[1], bb[2], bb[3], bb[4], border = "red")
I really like this gist. Thanks for posting it. The picture is awesome, and it's instructive for ximage. I like how it reads through an in-memory VRT that is a virtual warp, whose source is then a virtual mosaic on a remote file system. I don't know how the coordinate systems work here but very cool that you can line everything up. The
gisattributes fromgdalraster::read_ds()were recently added in the development version in case anyone tries this code.