rstac_PlanetaryComp_Cop30.R

library(rstac)
library(sf)
library(raytrix)
library(rayshader)

# set region and search stac!
region <- st_bbox(st_buffer(st_point(c(-151.739,-16.501)), 0.063)) %>%
  st_as_sfc() %>%
  st_set_crs(4326)

stac_region <- region%>%
  st_bbox()

s_obj <- stac("https://planetarycomputer.microsoft.com/api/stac/v1")
get_request(s_obj)

it_obj <- s_obj %>%
  stac_search(collections = "cop-dem-glo-30",
              bbox = stac_region) %>%
  get_request()

src_list <-assets_list(it_obj)$path %>%
  .[grep(".tif$", .)]

#raytrix/rayshader
set_canvas_sf(region %>% st_transform(3857))

cop30 <- topo_matrix(30, src=src_list)

texture <- cop30 %>%
  sphere_shade(texture = "imhof3") %>%
  add_water(detect_water(cop30), color='#2A6E93') %>%
  add_shadow(ray_shade(cop30, zscale = 30, sunaltitude = 25), 0.5) %>%
  add_shadow(ambient_shade(cop30), 0)

plot_map(texture)
plot_3d(texture, cop30,zscale = 20)

Thank you for sharing this code.

I am replicating your code, and when I run:

src_list <-rstac::assets_list(it_obj)$path %>% .[grep(".tif$", .)]

I get the following error:

Error: 'assets_list' is not an exported object from 'namespace:rstac'_

Could you please help me on this.

Point:

I am trying to use your package raytrix which I find much interesting.

The above dataset has a 30m resolution, any advice for a more detailed dataset for free.

Thank you,

Luis

Hi Luis, So I can help with the DEM url issues - the rstac API has updated and I tend to use functions like this now:

library(rstac)
library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.6.4, PROJ 9.1.1; sf_use_s2() is TRUE

rstac_bbox <- function(x) {
  x <- sf::st_transform(x, "EPSG:4326") |>
    sf::st_bbox()
  x[c(1:4)]
}

#' Get the tile sources for Copernicus GLO-30, ALOS World 3D-30m, NASADEM HGT
#' v001 DTM using the planetary computer STAC catalogue.
#'
#' @param aoi a spatial extent which can be provided as a numeric bounding box
#' or using a spatial object from either the {terra}, {sf} or {stars} packages.
#' @param collection Any of: "cop-dem-glo-30", "cop-dem-glo-90", "alos-dem", "nasadem".
#' @param vsicurl logical with default TRUE. If TRUE the '/vsicurl/' prefix is appended
#' to returned urls which increases read speed for gdal. see https://gdal.org/user/virtual_file_systems.html#vsicurl-http-https-ftp-files-random-access
#' for more information.
#'
#' @return A character vector of url paths.
mpc_dtm_src <- function(aoi,
                        collection = c(
                          "cop-dem-glo-30", "cop-dem-glo-90",
                          "alos-dem", "nasadem"
                        ),
                        vsicurl = TRUE) {
  checkmate::assertChoice(
    x = collection[1],
    choices = c("cop-dem-glo-30", "cop-dem-glo-90", "alos-dem", "nasadem")
  )

  stac_region <- rstac_bbox(aoi)

  s_obj <- rstac::stac("https://planetarycomputer.microsoft.com/api/stac/v1")
  rstac::get_request(s_obj)

  it_obj <- s_obj |>
    rstac::stac_search(
      collections = collection[1],
      bbox = stac_region
    ) |>
    rstac::get_request()

  src_list <- rstac::assets_url(it_obj)

  .urls <- src_list[grep(".tif$", src_list)]

  if (isTRUE(vsicurl)) {
    .urls <- sapply(.urls, function(x) paste0("/vsicurl/", x), USE.NAMES = FALSE)
  }

  return(.urls)
}


# set region and search stac!
region <- st_bbox(st_buffer(st_point(c(-151.739, -16.501)), 0.063)) %>%
  st_as_sfc() %>%
  st_set_crs(4326)


x <- mpc_dtm_src(aoi = region)
r <- terra::rast(x)
terra::plot(r)

^{Created on 2023-10-23 with reprex v2.0.2}

However, regarding {raytrix} I'm afraid much of this package is very unlikely to work - I haven't had the time to maintain it and several of the dependencies have changed considerably since I was last working on this. I'd like to get this up and running again soon but unfortunately I don't ahve the time at the moment.

I believe {rayvista} should still work but it less flexible than {raytrix} intended to be. However, many of the ideas in {raytrix} could just be extacted and used with just {rayshader} anyway.

re. higher resolutions, if you're after global data then usual 30 m is as good as it gets but if you check out https://opentopography.org/, they have an amazing collection of very high resolution data from various locations around the world.

Good luck and thanks for the interest in the package.

Hello,
thank you very much for all the tips and solutions presented.
As far as {raytrix} is concerned, in my opinion the concept is very promising.
Thanks again.