get_elev.md

Get elevation data from Copernicus 30m or GEBCO 2023. Run get_elev(rast) with any raster to get elevation data from online.

The function chooses which source based on the distance across the centre cell of the query raster, if that is less than 'threshold' it goes for the 30m SRTM source "cop30".

Requires in-dev terra, see install instructs in code below

  library(terra)
#> terra 1.7.41
get_elev <- function(x, method = "bilinear", maxcell = 25e6, silent = TRUE, threshold = 50) {
  if (terra::ncell(x) > maxcell) {
    stop("number of cells in x exceeds 'maxcell'")
  }
  
  sources <- c(gebco =  "/vsicurl/https://gebco2023.s3.valeria.science/gebco_2023_land_cog.tif", 
               cop30 =  "/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt")
  
  dm <- dim(x)
  ## gebco is about 500m
  centre <- terra::as.polygons(terra::ext(x, cells =  terra::cellFromRowCol(x, dm[1] %/% 2, dm[2] %/% 2)), crs = terra::crs(x))
  ds <- sqrt(terra::expanse(centre))
    
  if (ds < threshold) {
    src <- sources[["cop30"]]
  } else {
    src <- sources[["gebco"]]
  }
  if (!silent) {
    message(
      sprintf("choosing source %s: %s\n based on distance across centre cell %f (metres)", names(src), src, ds)
      
    )
  }
  terra::project(terra::rast(src), x, by_util = TRUE, method = method)
}


## an example from elevatr
ex <- c( 1802081, 1807610, 622060, 633833.8 )
crs <- "+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" 
plot(r <- get_elev(rast(ext(ex), crs = crs, res = 5)))

## get any random location
cit <- maps::world.cities[sample(nrow(maps::world.cities), 1L), ]
ex <- c(-1, 1, -1, 1) * 5e3
crs <- sprintf("+proj=laea +lon_0=%f +lat_0=%f +datum=WGS84", cit$long, cit$lat)
plot(r <- get_elev(rast(ext(ex), crs = crs, ncols = 1024, nrows = 1024), silent = FALSE ))
#> 
title(sprintf("%s (%s)", cit$name, cit$country.etc))

r
#> class       : SpatRaster 
#> dimensions  : 1024, 1024, 1  (nrow, ncol, nlyr)
#> resolution  : 9.765625, 9.765625  (x, y)
#> extent      : -5000, 5000, -5000, 5000  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=50.53 +lon_0=16.23 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> source(s)   : memory
#> name        : COP30_hh 
#> min value   : 367.9446 
#> max value   : 724.8431

## anywhere on earth
get_elev(rast())
#> class       : SpatRaster 
#> dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source(s)   : memory
#> name        : gebco_2023_land_cog 
#> min value   :               -8104 
#> max value   :                5528

^{Created on 2023-07-08 with reprex v2.0.2}