## random city on earth
pt <- maps::world.cities |> dplyr::sample_n(1) |> dplyr::select(long, lat, "name")
library(gdalraster)
prj <- sprintf("+proj=laea +lon_0=%i +lat_0=%i +datum=WGS84", as.integer(pt[1]), as.integer(pt[2]))
xy <- transform_xy(pt[,1:2, drop = F], srs_to = prj, srs_from = "EPSG:4326")
bb <- c(-1, -1, 1, 1) * 5000
warp(sds::cop30(), dstwarp <- tempfile(fileext = ".vrt", tmpdir = "/vsimem"),
prj, cl_arg = c("-ts", 1024, 0, "-te", bb[1], bb[2], bb[3], bb[4], "-r", "bilinear"))
args <- list(input = ds <- new(GDALRaster, dstwarp),
output_format = "MEM",
output = "", zfactor = "10")
alg <- gdal_run("raster hillshade", args)
(ds_hillshade <- alg$output())
alg$release()
warp(sds::wms_arcgis_mapserver_ESRI.WorldImagery_tms(), imwarp <- tempfile(fileext = ".vrt", tmpdir = "/vsimem"),
prj, cl_arg = c("-ts", 1024, 0, "-te", bb[1], bb[2], bb[3], bb[4], "-r", "bilinear"))
im <- new(GDALRaster, imwarp, read_only = FALSE)
im$addBand(im$getDataTypeName(1L), NULL)
translate(im, tf <- tempfile(fileext = ".tif"))
im2 <- new(GDALRaster, tf, read_only = FALSE)
im2$fillRaster(4L, .65*255, 0)
im2$close()
#plot_raster(im2)
par(mfrow = c(1, 2), mar = rep(0, 4))
plot_raster(ds_hillshade, main = pt$name)
ximage::ximage(vapour::gdal_raster_image(tf, bands = 1:3), add = T)
plot_raster(ds_hillshade, main = pt$name)
ximage::ximage(vapour::gdal_raster_image(tf), add = T)
using a mix of ways about this