Earth Engine Export Example • rgeedim

library(rgeedim)

## rgeedim v0.4.0 -- using geedim 2.0.0 w/ earthengine-api 1.7.10

library(terra)

## terra 1.8.97

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

## Using Application Default Credentials (ADC)

ee <- earthengine()

# use the earth engine API and reticulate to access
# World Database of Protected Areas <https://developers.google.com/earth-engine/datasets/catalog/WCMC_WDPA_current_polygons>
# filter to a specific area (yellowstone)
b <- ee$FeatureCollection("WCMC/WDPA/current/polygons")$filter(
  ee$Filter$And(
    ee$Filter$eq("NAME", "Yellowstone"),
    ee$Filter$eq("DESIG_ENG", "National Park")
  )
)
g <- gd_bbox(b)

# create a terra SpatVector from the GeoJSON-like list from gd_region()
yellowstone <- gd_region(b) |>
  gd_region_to_vect() |>
  project("EPSG:5070")

# inspect
plot(yellowstone)

# create an Earth Engine asset you can use again in the future
# 250m resolution for Yellowstone extent
"USGS/SRTMGL1_003" |>
  gd_image_from_id() |>
  gd_export(
    region = g,
    scale = 250,
    bands = list("elevation"),
    resampling = "bilinear",
    type = "asset",
    filename = "Yellowstone",
    folder = project_id,
    overwrite = TRUE,
    crs = "EPSG:5070"
  )

## <Task NMPBXBIINRYL7DSRFB6HX4VT Type.EXPORT_IMAGE: USGS-SRTMGL1_003 (State.UNSUBMITTED)>

# create a local GeoTIFF from the above created asset
x_path <- sprintf("projects/%s/assets/Yellowstone", project_id) |>
  gd_image_from_id() |>
  gd_download(
    bands = list("elevation"),
    resampling = "bilinear",
    region = g,
    scale = 250,
    crs = "EPSG:5070"
  )

## Error in py_call_impl(callable, call_args$unnamed, call_args$named) : 
##   ee.ee_exception.EEException: Image.load: Image asset 'projects/rgeedim-demo/assets/Yellowstone' not found (does not exist or caller does not have access).
## Run `reticulate::py_last_error()` for details.

if (!inherits(x_path, "try-error")) {
  x <- rast(x_path)
  plot(x)

  # calculate hillshade locally
  slp <- terrain(x, unit = "radians")
  asp <- terrain(x, "aspect", unit = "radians")
  shd <- shade(slp, asp)

  # plot with boundary
  plot(shd)
  plot(yellowstone, add = TRUE)
  
  # clean up
  unlink(sources(x))
}