Use rstac to query the swisstopo lidar collection and load into R. plot with terra and rayshader...

rstac_swisstopo.R

#rstac example/test
library(vapour)
library(gdalio)
library(rstac)
library(sf)
library(dplyr)
library(scico)
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))

# get random swiss region
region <- st_bbox(st_buffer(st_point(c(7.403000987649576,46.34570242792326)), 0.03)) %>%
  st_as_sfc() %>%
  st_set_crs(4326)

stac_region <- region%>%
  st_bbox()

# swiss
s_obj <- stac("https://data.geo.admin.ch/api/stac/v0.9/")

# list all collections and get id.
# colls <- rstac::collections(s_obj) %>%
#   get_request()
# idlist <- purrr::map(colls$collections, ~ .x$id)

#
it_obj <- s_obj %>%
  stac_search(collections = "ch.swisstopo.swissalti3d",
              bbox = stac_region) %>%
  get_request()

# list all COGS in collection.
cog_list <- purrr::map(it_obj$features, ~ paste0('/vsicurl/',.$assets[[1]][[2]]))  %>%
  unlist()

prj = vapour_raster_info(cog_list[1])$projection

extList <- purrr::map(it_obj$features, ~ .$bbox)

# set extent for all tiles
ext1 <- extList %>% unname() %>%
  as.data.frame() %>%
  tibble::rownames_to_column() %>%
  tidyr::pivot_longer(-rowname) %>%
  tidyr::pivot_wider(names_from=rowname, values_from=value) %>%
  summarise(xmin=min(`1`), ymin = min(`2`),
            xmax=max(`3`), ymax=max(`4`)) %>%
  unlist() %>%
  st_bbox() %>% st_as_sfc() %>% st_set_crs(4326) %>%
  st_transform(prj) %>%
  st_bbox() %>%
  .[order(c(1,3,2,4))]

# set extent for requested area.
targ_ext <- region %>%
  st_transform(prj) %>%
  st_bbox() %>%
  .[order(c(1,3,2,4))]



g <- list(extent = targ_ext, dimension = c(1202, 1602), projection = prj)

gdalio_set_default_grid(g)

im <- gdalio_terra(cog_list)
terra::plot(im, col=scico(256))

# Messing about with raytrix and rayshader if you like...
# devtools::load_all()
library(raytrix)
library(rayshader)

set_canvas(bounds = ext1, projection = prj)
get_canvas(5)

tc <- topo_matrix(5, src=cog_list)
# ov <- map_drape(5, alpha = 0.6)

pal.cols <- scico(5, palette = 'romaO', direction = 1)

tc %>%
  sphere_shade(texture = create_texture(pal.cols[2], pal.cols[5], pal.cols[1],
                                        pal.cols[4], pal.cols[3])) %>%
  # add_overlay(ov) %>%
  add_shadow(lamb_shade(tc, zscale = 2)) %>%
  plot_3d(., tc, zscale=5, windowsize = 1000)

render_snapshot(filename = 'dev_tests/exports/SwisstopoTest.png')