Hillshade Example

library(rgeedim)
library(terra)

project_id <- Sys.getenv("GOOGLE_CLOUD_QUOTA_PROJECT", "rgeedim-demo")
gd_initialize(project = project_id)

## hillshade example
b <- gd_bbox(
  xmin = -120.296,
  xmax = -120.227,
  ymin = 37.9824,
  ymax = 38.0071
)

# download 1m DEM in AEA
# need to set resampling method in composite step
x <- "USGS/3DEP/1m" |>
  gd_collection_from_name() |>
  gd_search(region = b) |>
  gd_composite(
    resampling = "bilinear"
  ) |> 
  gd_download(
    region = b,
    bands = list("elevation"),
    crs = "EPSG:5070",
    scale = 10,
    filename = "dem.tif"
  ) |>
  rast()

# calculate slope, aspect, and hillshade with terra
slp <- terrain(x, "slope", unit = "radians")
asp <- terrain(x, "aspect", unit = "radians")
hs <- shade(slp, asp)

# compare elevation v.s. hillshade
plot(c(x, hillshade = hs))