rgeedim Demo

This vignette shows how to extract a Google Earth Engine asset by name for an arbitrary extent and visualize it in R.

First, we load {rgeedim}.

library(rgeedim)

If you have the necessary Python dependencies installed (geedim, earthengine-api), you will see the versions printed out when the package is loaded.

If this is your first time using any Google Earth Engine tools, authenticate with gd_authenticate().

You can pass arguments to use several different authorization methods. Perhaps the easiest to use is auth_mode="notebook" in that does not rely on an existing GOOGLE_APPLICATION_CREDENTIALS file nor an installation of the gcloud CLI tools. However, the other options are better for non-interactive use.

gd_authenticate(auth_mode = "notebook")

You only need to authenticate periodically, depending on the method you used.

In each session we need to call gd_initialize(). This is a wrapper function around geedim.Initialize() that must be run before using the Python Google Earth Engine API.

gd_initialize()

Note that with auth_mode="gcloud" you need to specify the project via project= argument, in your default configuration file or via system environment variable GOOGLE_CLOUD_QUOTA_PROJECT. The authorization tokens generated for auth_mode="notebook" are always associated with a specific project.

Perhaps the simplest way to specify the target extent is using the xmin/xmax/ymin/ymax arguments to gd_bbox(). This function returns a Python object equivalent to GeoJSON, which is interchangeably represented as a simple list object in R using {reticulate}.

Determine Target Region

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

As is standard for GeoJSON, coordinates of the bounding box are expressed in WGS84 decimal degrees ("OGC:CRS84"). Note that longitude, latitude (X, Y) coordinate pair order is implied.

Access Images by ID

We can find IDs of assets of interest using the Google Earth Engine data catalog: https://developers.google.com/earth-engine/datasets/catalog

To obtain an R object reference to the asset we pass the "id" to gd_image_from_id(). For example here we use Global SRTM Topographic Diversity:

Download SRTM Topographic Diversity product

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

gd_image_from_id() will return geedim.mask.MaskedImage and gd_collection_from_name() will return geedim.collection.MaskedCollection objects.

Now we pass the image result to gd_download(). We can specify output filename and target area as region arguments. See gd_bbox() for examples of making a region argument from bounding coordinates or a {terra} SpatExtent object.

Other options that can be passed to the BaseImage.download() method include scale which allows warping of the result to a target resolution. Try modifying this example to use scale=90 (~native SRTM resolution):

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

gd_download() (invisibly) returns the filename on successful download, which helps to “pipe” into functions that might read the result.

So we can use the {terra} rast() function to read the GeoTIFF gd_download() result.

library(terra)
f <- rast(img)
par(mar = c(1, 1, 1, 1))
plot(f[[1]])

# inspect object
f

Create a local Hillshade from DEM

This example demonstrates the download of a section of the USGS NED seamless 10m grid. This DEM is processed locally with {terra} to calculate some terrain derivatives (slope, aspect) and a hillshade.

library(rgeedim)
library(terra)

gd_initialize()

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/NED" |>
  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))

Subsets of the "USGS/NED" image result in multi-band GeoTIFF with "elevation" and "FILL_MASK" bands. In the contiguous US we know the DEM is continuous so the FILL_MASK is not that useful. With geedim >1.7 we retrieve only the "elevation" band by specifying argument bands = list("elevation"). This cuts the raw image size that we need to download in half.

Working with Image Collections

Many Google Earth Engine assets of interest are “collections.” These assets contain data for several bands, dates, processes, resolutions, etc. There are standard ways that these collections can be aggregated; also, custom processing can be achieved through the Earth Engine API or by downloading data to your computer and working offline.

Filters on date and data quality are possible with gd_search(); many processes on collections cannot proceed without first “searching” it. We can quickly inspect search results with gd_properties()

A key step in many processes involving image collections is the use of gd_composite() to re-sample (or otherwise “combine”) component images of interest from the collection into a single layer. These operations are often performed on the server-side prior to download. A composite image can also be made automatically through gd_download() (which defaults to composite=TRUE). Setting composite=FALSE will allow individual images to be downloaded as separate GeoTIFF without aggregation to a single layer.

Composite LiDAR DEM

This example demonstrates how to access 1-meter LiDAR data from the USGS 3D Elevation Program (3DEP). LiDAR data are not available everywhere, and are generally available as collections of tiles for specific dates. We will use gd_search() to narrow down the options.

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

gd_initialize()

# 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)

The small sample extent covers only one tile, but in larger examples LiDAR data sources co-occur or span larger extents with variable coverage.

Downloading individual images in a collection

This example downloads each component image of the DAYMET climate collection “daily precipitation” layers into a temporary folder as separate GeoTIFFs. These individual images can be further aggregated or processed outside of Google Earth Engine.

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

gd_initialize()

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)
gd_collection_from_name('NASA/ORNL/DAYMET_V4') |> 
  gd_search(region = r,
            start_date = "2020-01-21",
            end_date = "2020-01-27") -> res

# get table of IDs and dates
p <- gd_properties(res)
td <- file.path(tempdir(), "DAYMET_V4")

# create a new collection using gd_collection_from_list()
# download each image as separate GeoTIFF (no compositing)
# Note: `filename` is a directory
gd_collection_from_list(p$id) |>
  gd_download(
    filename = td,
    composite = FALSE,
    dtype = 'int16',
    region = r,
    bands = list("prcp"),
    crs = "EPSG:5070",
    scale = 1000
  ) |> 
  rast() -> x2

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

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