country_dem.md

terra_out <- 
  function(x) {
    loadNamespace("terra")
    if (isS4(x) && inherits(x, "SpatRaster")) {
      x <- try(list(extent = c(terra::xmin(x), terra::xmax(x), terra::ymin(x), terra::ymax(x)),
                    dimension = dim(x)[2:1],
                    projection = terra::crs(x), 
                    lonlat = terra::is.lonlat(x), terra = TRUE), silent = TRUE)
      if (inherits(x, "try-error")) stop("cannot use terra grid")
    }
    x  
  }

raster_out <- function(x) {
  loadNamespace("raster")
  if (isS4(x) && inherits(x, "BasicRaster")) {
    ## we have a {raster}
    x <- list(extent = c(x@extent@xmin, x@extent@xmax, x@extent@ymin, x@extent@ymax),
              dimension = c(x@ncols, x@nrows),
              projection = x@crs@projargs, 
              lonlat = raster::couldBeLonLat(x), terra = FALSE)
  }
  x
}
#' Get digital elevation data
#'
#' Read elevation data for any region. 
#' 
#' Currently using GEBCO 2019 as a background, and SRTM 30m resolution for 
#' regions that fit approximately within the size of an SRTM tile (these are 1 degree wide). 
#' 
#' Note that data is streamed into memory, so don't make the dimensions of the 'x' target raster too big. 
#' @param x a terra rast object, any extent any crs and any dims (not too big! - 2000x2000, not 10000x10000)
#' @param ... 
#' @param source a GDAL raster source, to override the inbuild GEBCO + SRTM (in future we might patch in local source)
#' @param threshold a size in degrees above which no SRTM data is queried (about )
#'
#' @return a terra rast object with elevation data
#' @export
#'
#' @examples
#' get_elev(terra::rast())
#' get_elev(raster::raster())
#' 
#' get_elev(raster::raster(raster::extent(80, 120, -60, -40), res = 0.25, crs = "OGC:CRS84"))
#' 
#' get_elev(raster::raster(raster::extent(c(-1, 1, -1, 1) * 25e3), nrows = 1024, ncols = 1024, crs = "+proj=laea +lat_0=44.6371 +lon_0=-63.5923"))
get_elev <- function(x, resample = "bilinear", ..., source = NULL, threshold = 0.5) {
  xraster <- x
  if (inherits(x, "SpatRaster")) {
    x <- terra_out(x)
  } else {
    x <- raster_out(x)
  }
  wh <- c(diff(x$extent[1:2]), diff(x$extent[3:4]))
  no_srtm <- FALSE
  if (x$lonlat) {
    if (max(wh) > threshold) {
      no_srtm <- TRUE
    } 
  } else {
    ## assuming metres ...
    if (max(wh) > (threshold * 111111.12)) {
      no_srtm <- TRUE
    }
  }
  
  if (is.null(source)) {
    rso <- c("/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2019_GEOTIFF/GEBCO_2019.tif", 
             "/vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt")
    if (no_srtm) rso <- rso[1L]
    if (!no_srtm) print("SRTM in use, in addition to GEBCO")
  } else {
    rso <- source
  }
  vals <- vapour::vapour_warp_raster_dbl(rso, extent = x$extent, dimension = x$dimension, projection = x$projection, resample = resample, ...)
  if (x$terra) {
    terra::setValues(xraster, vals)
  } else {
    raster::setValues(xraster, vals)
  }
}

library(terra)

sfx <- rnaturalearth::ne_countries(returnclass = "sf")
op <- par(mfrow = n2mfrow(nrow(sfx)), mar = rep(0, 4))
dv <- dev.size("px")
mfrow <- par("mfrow")
for (i in seq_len(nrow(sfx))) {
country <- vect(sfx[i, ])
  template <- terra::rast(country, nrows = dv[1]/mfrow[2], ncols = dv[2]/mfrow[1])

  dem1 <- get_elev(template)
  image(mask(dem1, country), asp = 1, axes = F, xlab = "", ylab = "", col = hcl.colors(64))
}

projected version

library(terra)

sfx <- rnaturalearth::ne_countries(returnclass = "sf")
op <- par(mfrow = n2mfrow(nrow(sfx)), mar = rep(0, 4))
dv <- dev.size("px")
mfrow <- par("mfrow")
for (i in seq_len(nrow(sfx))) {
  country <- vect(sfx[i, ])
  centre <- geom(terra::centroids(country))
  pcountry <- project(country, sprintf("+proj=laea +lon_0=%f +lat_0=%f", centre[1, "x"], centre[1, "y"]))
  template <- terra::rast(pcountry, nrows = dv[1]/mfrow[2] * 2, ncols = dv[2]/mfrow[1] * 2)

  dem1 <- get_elev(template)
  image(mask(dem1, pcountry), asp = 1, axes = F, xlab = "", ylab = "", col = hcl.colors(64))
}