get_elev_right_scale.md

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.9.0dev-cb4d30f56d, PROJ 9.3.1; sf_use_s2() is
#> TRUE
library(raster)
#> Loading required package: sp


## threshold is roughly side length in metres of the roi, i.e. sqrt(area)
dem_source <- function(threshold = 1) {
  ## threshold should be about the width of the region, in metres
  switches <- tibble::tibble(thresh = c(0, 5e3, 5e4), 
                             dsn = c("/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt", 
                                     "/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP90/COP90_hh.vrt", 
                                     "/vsicurl/https://gebco2023.s3.valeria.science/gebco_2023_land_cog.tif"
                                     ))
  idx <- findInterval(threshold, switches$thresh)
  if (is.na(idx)) idx <- nrow(switches)
  switches$dsn[idx]
}
## x is sf
## dsn is string DSN
get_elev0 <- function(x, dm = 1024L, dsn = NULL) {

  ## what area do we have (will dictate resolution of source data)
  A = as.numeric(units::set_units(sf::st_area(sf::st_as_sfc(sf::st_bbox(x))), "m^2"))
  ## pick a side length dm, use the res trick to get the aspect ratio right
   if (is.null(dsn)) dsn <- dem_source(sqrt(A))
  ## extra homework to align to a nice clean grid maybe ...
  dem_r = terra::rast(terra::ext(sf::st_bbox(x)[c(1, 3, 2, 4)]), ncols = dm, 
                      crs = sf::st_crs(x)$wkt)
  terra::res(dem_r) = rep(terra::res(dem_r)[1L], each = 2L)
  ## now hit up that VRT
  print(dsn)
  terra::project(terra::rast(dsn), dem_r, by_util = TRUE, method = "bilinear")
}

get_elev0(slopes::cyclestreets_route)
#> [1] "/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt"
#> class       : SpatRaster 
#> dimensions  : 498, 1024, 1  (nrow, ncol, nlyr)
#> resolution  : 3.249023e-05, 3.249023e-05  (x, y)
#> extent      : -1.55598, -1.52271, 53.80275, 53.81893  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        : COP30_hh 
#> min value   : 27.93542 
#> max value   : 98.10830


trace_dsn = "/vsizip//vsicurl/https://ndownloader.figshare.com/files/14331185/3DGRT_AXIS_EPSG25830_v2.zip/3DGRT_AXIS_EPSG25830_v2.shp"

trace = sf::read_sf(trace_dsn)
get_elev0(trace)
#> [1] "/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt"
#> class       : SpatRaster 
#> dimensions  : 1302, 1024, 1  (nrow, ncol, nlyr)
#> resolution  : 3.003385, 3.003385  (x, y)
#> extent      : 443496.7, 446572.2, 4165588, 4169498  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830) 
#> source(s)   : memory
#> name        : COP30_hh 
#> min value   : 588.7089 
#> max value   : 944.1055

## it should scale out ok, because the cog has overviews (COP does have overview within a 1-degree tile)
get_elev0(sf::st_as_sfc(sf::st_bbox(c(xmin = 0, ymin = 0, xmax = 50, ymax = 50), crs = "EPSG:4326")))
#> [1] "/vsicurl/https://gebco2023.s3.valeria.science/gebco_2023_land_cog.tif"
#> class       : SpatRaster 
#> dimensions  : 1024, 1024, 1  (nrow, ncol, nlyr)
#> resolution  : 0.04882812, 0.04882812  (x, y)
#> extent      : 0, 50, 0, 50  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        : gebco_2023_land_cog 
#> min value   :               -5128 
#> max value   :                5134

^{Created on 2024-04-17 with reprex v2.0.2}