EVAL <- isTRUE(as.logical(Sys.getenv("R_RGEEDIM_RUN_EXAMPLES"))) &&
  requireNamespace("terra", quietly = TRUE) &&
  rgeedim::gd_is_initialized(project = Sys.getenv("GOOGLE_CLOUD_QUOTA_PROJECT", "rgeedim-demo"))

litedown::reactor(
  eval = EVAL,
  collapse = TRUE,
  fig.width = 8,
  fig.align = 'center'
)

library(rgeedim)

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

r <- gd_bbox(
  xmin = -121,
  xmax = -120.5,
  ymin = 38.5,
  ymax = 39
)

x <- gd_image_from_id('CSP/ERGo/1_0/Global/SRTM_topoDiversity')

img <- gd_download(x, filename = 'image.tif',
                   region = r, scale = 900,
                   overwrite = TRUE, silent = FALSE
                  )

library(terra)
f <- rast(img)

par(mar = c(1, 1, 1, 1))
plot(f[[1]])

# inspect object
f

library(rgeedim)
library(terra)

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

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

## hillshade example
# download 10m NED DEM in AEA
x <- "USGS/SRTMGL1_003" |>
  gd_image_from_id() |>
  gd_download(
    region = b,
    scale = 10,
    crs = "EPSG:5070",
    resampling = "bilinear",
    filename = "image.tif",
    bands = list("elevation"),
    overwrite = TRUE,
    silent = FALSE
  )
dem <- rast(x)$elevation

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

# compare elevation v.s. hillshade
plot(c(dem, hillshade = hsd))

# search and download composite from USGS 1m lidar data collection
library(rgeedim)
library(terra)

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

# wkt->SpatVector->GeoJSON
b <- 'POLYGON((-121.355 37.56,-121.355 37.555,
          -121.35 37.555,-121.35 37.56,
          -121.355 37.56))' |>
  vect(crs = "OGC:CRS84")

# create a GeoJSON-like list from a SpatVector object
# (most rgeedim functions arguments for spatial inputs do this automatically)
r <- gd_region(b)

# search collection for an area of interest
a <- "USGS/3DEP/1m" |>
  gd_collection_from_name() |>
  gd_search(region = r) 

# inspect individual image metadata in the collection
gd_properties(a)

# resampling images as part of composite; before download
x <- a |>
  gd_composite(resampling = "bilinear") |> 
  gd_download(region = r,
              crs = "EPSG:5070",
              scale = 1,
              filename = "image.tif",
              overwrite = TRUE,
              silent = FALSE) |>
  rast()

# inspect
plot(terra::terrain(x$elevation))
plot(project(b, x), add = TRUE)

# search and download individual images from daymet V4
library(rgeedim)
library(terra)

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

r <- gd_bbox(
  xmin = -121,
  xmax = -120.5,
  ymin = 38.5,
  ymax = 39
)

# search collection for spatial and date range (one week in January 2020)
res <- gd_collection_from_name('NASA/ORNL/DAYMET_V4') |> 
  gd_search(region = r,
            start_date = "2020-01-21",
            end_date = "2020-01-27")

# get table of IDs and dates
p <- gd_properties(res)
p

td <- tempdir()

# download each image as separate GeoTIFF (no compositing)
# Note: `filename` is a directory
res2 <- res |>
  gd_download(
    filename = td,
    composite = FALSE,
    dtype = 'int16',
    region = r,
    bands = list("prcp"),
    crs = "EPSG:5070",
    scale = 1000
  ) 

x2 <- rast(res2)

# filter to bands of interest (if needed)
x2 <- x2[[names(x2) == "prcp"]]

# set time for each layer
time(x2) <- p$date
panel(x2)
title(ylab = "Daily Precipitation (mm)")

unlink("image.tif")
unlink(td, recursive = TRUE)