rGEDI: An R Package for NASA's Global Ecosystem Dynamics Investigation (GEDI) Data Visualizing and Processing.

Authors: Carlos Alberto Silva and Caio Hamamura.

The original rGEDI workflow provides functions for downloading, visualizing, clipping, gridding, simulating, and exporting GEDI data. The expanded package also searches, streams, filters, summarizes, models, and maps NASA GEDI data.

rGEDI searches, downloads, streams, reads, filters, clips, summarizes, simulates, models, and maps NASA Global Ecosystem Dynamics Investigation (GEDI) data. It supports footprint products GEDI01_B, GEDI02_A, GEDI02_B, and GEDI04_A, as well as raster products GEDI03 and GEDI04_B.

The examples below form one workflow. Local examples use the real GEDI subsets bundled with the package. Cloud and Google Earth Engine examples use real NASA and Earth Engine services. The scripts that regenerated the displayed assets are readme/build-local-examples.R and readme/build-modern-examples.R. A generated function reference manual is available at output/pdf/rGEDI-reference-manual.pdf.

Getting Started

1 Installation

# Development release
install.packages(
  "rGEDI",
  repos = c("https://carlos-alberto-silva.r-universe.dev",
            "https://cloud.r-project.org")
)

# Current GitHub master
# install.packages("remotes")
# remotes::install_github("carlos-alberto-silva/rGEDI")

library(rGEDI)
library(data.table)
library(sf)
library(terra)

Configure Python and Google Earth Engine

As in the ICESat2VegR workflow, rGEDI uses three Python packages through reticulate:

1. earthaccess finds and authenticates Earthdata Cloud objects.
2. h5py reads HDF5 content streamed from the cloud.
3. earthengine-api samples predictors and creates wall-to-wall maps in Google Earth Engine.

Run the configuration once. It installs Miniconda when needed and creates the package Python environment. Restart R if the installer asks you to do so.

rGEDI_configure(install = TRUE)
ee_initialize(project = "your-google-cloud-project")

Verify the active Python and required modules:

safely <- function(expr, default = NA) tryCatch(expr, error = function(e) default)
have_reticulate <- requireNamespace("reticulate", quietly = TRUE)
py_ready <- have_reticulate &&
  !inherits(try(reticulate::py_config(), silent = TRUE), "try-error")

status <- list(
  rGEDI_loaded = "package:rGEDI" %in% search(),
  python_used = if (have_reticulate)
    safely(reticulate::py_discover_config()$python) else NA,
  h5py = if (py_ready) reticulate::py_module_available("h5py") else FALSE,
  earthaccess = if (py_ready)
    reticulate::py_module_available("earthaccess") else FALSE,
  ee = if (py_ready) reticulate::py_module_available("ee") else FALSE
)
print(status)

2 Example study site

The package contains Cerrado forest stands and small matched GEDI Level 1B, Level 2A, and Level 2B granules.

stands_path <- system.file("extdata", "stands_cerrado.shp", package = "rGEDI")
study_area <- sf::st_read(stands_path, quiet = TRUE)
box <- sf::st_bbox(study_area)

xmin <- unname(box["xmin"]); xmax <- unname(box["xmax"])
ymin <- unname(box["ymin"]); ymax <- unname(box["ymax"])
daterange <- c("2019-04-18", "2019-04-19")
outdir <- file.path(tempdir(), "rGEDI")
dir.create(outdir, showWarnings = FALSE)

3 Find GEDI data and download

3.1 Find GEDI granules

gedifinder() queries NASA CMR. Use the official short names and current versions. The returned character vector contains direct HTTPS links.

products <- c(
  GEDI01_B = "002", GEDI02_A = "002", GEDI02_B = "002",
  GEDI03 = "001", GEDI04_A = "003", GEDI04_B = "002"
)

granules <- lapply(names(products), function(product) {
  gedifinder(
    product,
    ul_lat = ymax, ul_lon = xmin,
    lr_lat = ymin, lr_lon = xmax,
    version = products[[product]], daterange = daterange,
    cloud_computing = FALSE, persist = TRUE
  )
})
names(granules) <- names(products)
head(granules$GEDI02_A)

To retrieve every production granule in one orbit, supply its orbit number instead of a bounding box. CMR returns the four parts in acquisition order.

orbit_granules <- gedifinder(
  "GEDI02_A", version = "002", orbit = "O01964", return = "table"
)
orbit_granules[, c("granule_id", "orbit", "granule_part")]
stopifnot(nrow(orbit_granules) == 4L)

Set cloud_computing = TRUE for cloud-hosted links suitable for streaming. An HTTPS Earthdata Cloud URL works on Windows, macOS, and Linux. Direct s3:// reads require code running in AWS us-west-2, which is a NASA bucket policy rather than an operating-system limitation.

3.2 Download the granules

Keep credentials outside scripts. earthdata_login() can use a standard ~/.netrc, the NETRC environment variable, or a path passed at run time.

# Sys.setenv(NETRC = "/secure/path/to/.netrc")
earthdata_login()

gediDownload(granules$GEDI01_B[1], outdir)
gediDownload(granules$GEDI02_A[1], outdir)
gediDownload(granules$GEDI02_B[1], outdir)
gediDownload(granules$GEDI04_A[1], outdir)

# GEDI03 and GEDI04_B searches may return several metric GeoTIFFs.
gediDownload(granules$GEDI03[1], outdir)
gediDownload(granules$GEDI04_B[1], outdir)

gediDownload() accepts the complete URL vector and invisibly returns the downloaded paths. Orbit O01964 is about 9.2 GB, so ensure sufficient disk space before running this example.

orbit_files <- gediDownload(orbit_granules$url, outdir)

4 Read GEDI products

Read downloaded products

level1b_file <- list.files(outdir, "GEDI01_B.*\\.h5$", full.names = TRUE)[1]
level2a_file <- list.files(outdir, "GEDI02_A.*\\.h5$", full.names = TRUE)[1]
level2b_file <- list.files(outdir, "GEDI02_B.*\\.h5$", full.names = TRUE)[1]
level4a_file <- list.files(outdir, "GEDI04_A.*\\.h5$", full.names = TRUE)[1]
level3_file  <- list.files(outdir, "GEDI03.*\\.tif$", full.names = TRUE)[1]
level4b_file <- list.files(outdir, "GEDI04_B.*\\.tif$", full.names = TRUE)[1]

level1b_full <- readLevel1B(level1b_file)
level2a_full <- readLevel2A(level2a_file)
level2b_full <- readLevel2B(level2b_file)
level3  <- readLevel3(level3_file)
level4a_full <- readLevel4A(level4a_file)
level4b <- readLevel4B(level4b_file)

The bundled subsets can be opened without a network connection:

level1b_path <- unzip(system.file("extdata",
  "GEDI01_B_2019108080338_O01964_T05337_02_003_01_sub.zip",
  package = "rGEDI"), exdir = outdir)
level2a_path <- unzip(system.file("extdata",
  "GEDI02_A_2019108080338_O01964_T05337_02_001_01_sub.zip",
  package = "rGEDI"), exdir = outdir)
level2b_path <- unzip(system.file("extdata",
  "GEDI02_B_2019108080338_O01964_T05337_02_001_01_sub.zip",
  package = "rGEDI"), exdir = outdir)

level1b <- readLevel1B(level1b_path)
level2a <- readLevel2A(level2a_path)
level2b <- readLevel2B(level2b_path)

Stream GEDI HDF5 data from Earthdata Cloud

cloud_urls <- gedifinder(
  "GEDI04_A",
  ul_lat = ymax, ul_lon = xmin, lr_lat = ymin, lr_lon = xmax,
  version = "003", daterange = daterange,
  cloud_computing = TRUE, persist = TRUE
)

# readLevel1B(), readLevel2A(), and readLevel2B() accept the same URL type.
level4a_cloud <- readLevel4A(cloud_urls[1])
level4a_cloud$beams

Level 3 and Level 4B are raster collections. Download the selected metric GeoTIFF first and then use readLevel3() or readLevel4B(). This avoids the Earthdata redirect restrictions encountered by GDAL virtual-file reads.

5 Extract, quality-filter, clip, rasterize, summarize, and simulate footprint-level products

5.1 Extract the Reference Ground Track and plot the GIF animation

getGEDITrack() standardizes coordinates, beam names, and acquisition order from a Level 1B, 2A, 2B, or 4A object, an extracted table, or a list of open granules. NASA partitions one GEDI orbit into four consecutive production granules (_01 through _04). gedifinder(orbit = ...) retrieves every part; every thins the display without changing acquisition order.

stopifnot(length(orbit_files) == 4L)

orbit_h5 <- lapply(sort(orbit_files), readLevel2A)
rgt <- getGEDITrack(orbit_h5, every = 200, segment_gaps = FALSE)
head(rgt)

plot_gedi_orbit_animation(
  orbit_h5,
  output_file = file.path(outdir, "gedi-orbit-animation.gif"),
  title = "GEDI aboard the International Space Station",
  duration = 8, every = 200, launch = FALSE
)

# Use .html for the interactive 3D globe.
plot_gedi_orbit_animation(
  orbit_h5, output_file = file.path(outdir, "gedi-orbit-animation.html"),
  every = 200,
  track_speed = 2, earth_rotation_speed = 2,
  launch = interactive()
)

The displayed example was generated from the four complete Level 2A V002 granules for orbit O01964. Its 22,350 sampled footprints span all eight GEDI beams, four granule parts, 81.6 minutes, and latitudes from -50.346 to 51.825 degrees. The reproducible cloud-reading and animation script is readme/build-complete-orbit.R.

The GIF follows the ICESat2VegR presentation while showing the correct GEDI platform: GEDI is mounted on the International Space Station. GEDI operates at 1064 nm in the near infrared; the red beam and track visualize that invisible laser pulse and connect the ISS payload to the accumulating reference ground track. Every displayed ISS position follows the time-ordered geolocation from the most complete beam in the open HDF5 granule. Following the ICESat2VegR animation, the textured Earth rotates independently while the NASA ISS image, GEDI payload label, red laser, and red orbit track move together in 3D. The HTML controls adjust both track playback and Earth rotation speed.

5.2 Extract the data

Get GEDI pulse geolocation (GEDI Level 1B)

level1b_geo <- getLevel1BGeo(level1b)
head(level1b_geo)

The original rGEDI map is retained below. It shows the Level 1B footprints over high-resolution imagery.

Get GEDI full waveform (GEDI Level 1B)

shot <- level1b_geo$shot_number[1]
waveform <- getLevel1BWF(level1b, shot_number = shot)
plot(waveform)

Get GEDI elevation and height metrics (GEDI Level 2A)

level2a_metrics <- getLevel2AM(level2a)
level2a_good <- level2a_metrics[
  quality_flag == 1 & degrade_flag == 0 & sensitivity >= 0.9
]
level2a_good[, .(beam, shot_number, elev_lowestmode, rh50, rh90, rh98, rh100)]

Plot waveform with RH metrics

plotWFMetrics(level1b, level2a, shot_number = shot,
              rh = c(25, 50, 75, 90, 98))

Get GEDI vegetation biophysical variables (GEDI Level 2B)

level2b_vpm <- getLevel2BVPM(level2b)
level2b_good <- level2b_vpm[l2b_quality_flag == 1 & sensitivity >= 0.9]
level2b_good[, .(beam, shot_number, rh100, pai, fhd_normal, cover)]

Get and plot PAI and PAVD profiles (GEDI Level 2B)

pai_profile  <- getLevel2BPAIProfile(level2b)
pavd_profile <- getLevel2BPAVDProfile(level2b)
profile_beam <- unique(pai_profile$beam)[1]
plotPAIProfile(pai_profile, beam = profile_beam)
plotPAVDProfile(pavd_profile, beam = profile_beam)

Get GEDI aboveground biomass at footprint level (GEDI Level 4A)

level4a_footprints <- getLevel4A(
  level4a_cloud,
  cols = c("beam", "shot_number", "delta_time", "lat_lowestmode",
           "lon_lowestmode", "agbd", "agbd_se", "sensitivity",
           "l4a_quality_flag_rel3", "degrade_include_flag",
           "elev_highestreturn_outlier_flag"),
  quality = TRUE
)
head(level4a_footprints)

5.3 Clip

Use the generic clip() function for every GEDI product. A bounding box uses the rGEDI/terra coordinate order c(xmin, xmax, ymin, ymax); polygons can be sf, sfc, or SpatVector objects.

bbox <- c(xmin, xmax, ymin, ymax)

# Open HDF5 products
level1b_clip <- clip(level1b, bbox,
  output = file.path(outdir, "level1b-clip.h5"))
level2a_clip <- clip(level2a, bbox,
  output = file.path(outdir, "level2a-clip.h5"))
level2b_clip <- clip(level2b, study_area,
  output = file.path(outdir, "level2b-clip.h5"))
level4a_clip <- clip(level4a, bbox,
  output = file.path(outdir, "level4a-clip.h5"))

# Extracted footprint tables and rasters
level2a_bbox <- clip(level2a_metrics, bbox)
level4a_geom <- clip(level4a_footprints, study_area)
level3_aoi <- clip(level3, vect(study_area))
level4b_bbox <- clip(level4b, bbox)

The output type follows the input type: open HDF5 objects return open clipped HDF5 objects, extracted tables return data.table objects, and Level 3/4B rasters return SpatRaster objects. The product-specific functions shown below remain available when an explicit function name is useful in a script.

Clip open GEDI HDF5 objects

Level 1B, 2A, 2B, and 4A clippers write valid subset HDF5 files and return open GEDI objects. Passing an extracted table instead returns a data.table.

level1b_clip <- clipLevel1B(level1b, xmin, xmax, ymin, ymax,
                            output = file.path(outdir, "level1b-clip.h5"))
level2a_clip <- clipLevel2A(level2a, xmin, xmax, ymin, ymax,
                            output = file.path(outdir, "level2a-clip.h5"))
level2b_clip <- clipLevel2B(level2b, xmin, xmax, ymin, ymax,
                            output = file.path(outdir, "level2b-clip.h5"))
level4a_clip <- clipLevel4A(level4a, xmin, xmax, ymin, ymax,
                            output = file.path(outdir, "level4a-clip.h5"))
level4a_table_clip <- clipLevel4A(level4a_footprints, xmin, xmax, ymin, ymax)

Clip the HDF5 products by geometry:

level1b_geom <- clipLevel1BGeometry(level1b, study_area,
  output = file.path(outdir, "level1b-geometry"))
level2a_geom <- clipLevel2AGeometry(level2a, study_area,
  output = file.path(outdir, "level2a-geometry"))
level2b_geom <- clipLevel2BGeometry(level2b, study_area,
  output = file.path(outdir, "level2b-geometry"))
level4a_geom <- clipLevel4AGeometry(level4a, study_area,
  output = file.path(outdir, "level4a-geometry.h5"))
level4a_table_geom <- clipLevel4AGeometry(level4a_footprints, study_area)

Clip extracted data.table objects

level1b_geo_bbox <- clipLevel1BGeo(level1b_geo, xmin, xmax, ymin, ymax)
level2a_bbox <- clipLevel2AM(level2a_metrics, xmin, xmax, ymin, ymax)
level2b_bbox <- clipLevel2BVPM(level2b_vpm, xmin, xmax, ymin, ymax)

level1b_geo_geom <- clipLevel1BGeoGeometry(level1b_geo, study_area)
level2a_geom_dt <- clipLevel2AMGeometry(level2a_metrics, study_area)
level2b_geom_dt <- clipLevel2BVPMGeometry(level2b_vpm, study_area)

plot(st_geometry(study_area))
points(level2a_geom_dt$lon_lowestmode, level2a_geom_dt$lat_lowestmode,
       pch = 16, col = "#762A83")

5.4 Compute descriptive statistics

metric_set <- function(x) c(
  n = sum(is.finite(x)), mean = mean(x, na.rm = TRUE),
  sd = sd(x, na.rm = TRUE), min = min(x, na.rm = TRUE),
  max = max(x, na.rm = TRUE)
)

rh98_stats <- polyStatsLevel2AM(
  level2a_geom_dt, func = metric_set(rh98), id = NULL
)
cover_stats <- polyStatsLevel2BVPM(
  level2b_geom_dt, func = metric_set(cover), id = NULL
)
agbd_stats <- polyStatsLevel4A(
  level4a_footprints, study_area, metric = "agbd", fun = metric_set
)

5.5 Compute grids with descriptive statistics

rh98_grid <- gridStatsLevel2AM(
  level2a_metrics, func = mean(rh98, na.rm = TRUE), res = 0.002
)
cover_grid <- gridStatsLevel2BVPM(
  level2b_vpm, func = mean(cover, na.rm = TRUE), res = 0.002
)
agbd_grid <- gridStatsLevel4A(
  level4a_footprints, metric = "agbd", fun = mean,
  res = 0.002, na.rm = TRUE
)

5.6 Convert to SpatVector

to_vect() detects the standard coordinate fields for Level 2A, Level 2B, and Level 4A tables.

level2a_vect <- to_vect(level2a_metrics)
level2b_vect <- to_vect(level2b_vpm)
level4a_vect <- to_vect(level4a_footprints)

6 Predicting and rasterizing local GEDI-derived forest attributes with machine learning

6.1 Create a model for GEDI data

This example reads the real Level 2A table generated from the bundled granule.

model_data <- fread("readme/gedi-level2a-example.csv")
predictors <- c("rh50", "rh75", "rh90", "rh100", "elev_lowestmode")

selection <- varSel(
  model_data[, ..predictors], model_data$rh98,
  method = "rfe", threshold = 0, seed = 42, ntree = 200
)
plot(selection, which = "importance", main = "RFE predictor importance")
selected_predictors <- selection$selvars

fit <- fit_model(
  model_data[, ..selected_predictors], model_data$rh98,
  method = "randomForest", test = "kfold", k = 5,
  seed = 42, ntree = 300
)
fit$stats_test

6.2 Predict field or fuel data

For a field response such as fuel load, join measured plots to GEDI footprints by spatial proximity or shot ID, then pass the measured column as y. The following call predicts the real rh98 response in this reproducible example.

predicted <- predictGEDI(fit, model_data, name = "predicted_rh98")

# For a measured fuel-load table, the corresponding call is:
# fuel_fit <- fit_model(fuel_training[, ..predictor_names],
#                       fuel_training$fuel_load,
#                       method = "randomForest", test = "kfold", k = 5)

6.3 Rasterize predicted data

rh98_prediction <- rasterizeGEDI(
  predicted, metric = "predicted_rh98", res = 0.002,
  lon = "lon_lowestmode", lat = "lat_lowestmode", fun = mean
)
plot(rh98_prediction)

7 Extract, clip, and summarize GEDI-derived grid products

level3 <- readLevel3(level3_file)
level4b <- readLevel4B(level4b_file)

7.1 Visualization

plotLevel3(level3, main = "GEDI Level 3 metric")
plotLevel4B(level4b, main = "GEDI Level 4B AGBD")

7.2 Clip using a bounding box or geometry

level3_bbox <- clipLevel3(level3, c(xmin, xmax, ymin, ymax))
level4b_bbox <- clipLevel4B(level4b, c(xmin, xmax, ymin, ymax))
level3_aoi <- clipLevel3(level3, vect(study_area))
level4b_aoi <- clipLevel4B(level4b, vect(study_area))
level3_values <- extractLevel3(level3, vect(study_area))
level4b_values <- extractLevel4B(level4b, vect(study_area))

7.3 Compute Level 4B descriptive statistics within an AOI

level4b_stats <- polyStatsLevel4B(
  level4b, vect(study_area), fun = mean, na.rm = TRUE
)
level4b_stats

8 Simulating GEDI full-waveform data from ALS point clouds

gediWFSimulator() accepts LAS/LAZ files, lidR::LAS objects, SpatVector points, or an X/Y/Z table. Its focused implementation follows Steven Hancock's gediRat core: Gaussian footprint weighting, local ALS-density correction, vertical binning, Gaussian pulse convolution, LAS class 2 ground separation, and unit-integral normalization. Its supported waveform path is tested bin by bin against a C reference extraction from upstream commit 49ad4f26b03083b61ee0bd424e8a7e4a8f2c301d, with agreement within 5e-7. Vectorized bin accumulation keeps this path portable and fast without the native library stack that caused CRAN problems. gediWFMetrics() supplies the focused RH, canopy-cover, energy, ground-height, and foliage-height-diversity metrics used here; it does not reproduce every optional gediMetric mode.

set.seed(11)
als <- data.frame(
  X = c(runif(1500, -12, 12), runif(3500, -12, 12)),
  Y = c(runif(1500, -12, 12), runif(3500, -12, 12)),
  Z = c(rnorm(1500, 0, 0.25), pmax(0, rnorm(3500, 17, 6))),
  Classification = c(rep(2L, 1500), rep(5L, 3500))
)

8.1 Extract metrics without adding waveform noise

sim_clean <- gediWFSimulator(
  als, output = file.path(outdir, "sim-clean.h5"),
  coords = c(0, 0), noise = 0, seed = 11
)
metrics_clean <- gediWFMetrics(sim_clean)
metrics_clean[, .(cover, rh50, rh90, rh100, waveEnergy)]

8.2 Extract metrics after adding waveform noise

sim_noisy <- gediWFSimulator(
  als, output = file.path(outdir, "sim-noisy.h5"),
  coords = c(0, 0), noise = 0.03, seed = 11
)
metrics_noisy <- gediWFMetrics(sim_noisy)
metrics_noisy[, .(cover, rh50, rh90, rh100, waveEnergy)]

9 Upscaling GEDI products using AlphaEarth Embeddings

9.1 Introduction

This workflow models real quality-filtered GEDI Level 4A aboveground biomass density (agbd) with the 64-band AlphaEarth annual embedding plus terrain predictors in Google Earth Engine (GEE). It follows the ICESat2VegR sequence, adapted to GEDI footprints: discover Level 4A data, apply product quality flags, create a spatially balanced footprint sample, extract AlphaEarth and terrain values, select predictors, validate a Random Forest model, train the equivalent GEE regressor, predict the AOI, visualize it, and export GeoTIFF.

AlphaEarth bands A00–A63 are general-purpose satellite embeddings. The terrain variables (elevation, slope, and aspect) provide explicit topographic context. Longitude and latitude can also be included, but their importance should be interpreted carefully because they may encode geographic location rather than transferable ecological relationships.

9.2 Code availability

A clean end-to-end script is installed with the package and can also be opened directly from the repository:

GEDI AlphaEarth upscaling workflow: inst/scripts/upscaling_alphaearth_workflow.R

The script covers GEDI discovery and cloud reading, Level 4A quality filtering, spatial sampling, ancillary extraction, RFE, independent validation, wall-to-wall prediction, visualization, and Drive export. The separate readme/build-modern-examples.R script regenerates the figures and data artifacts displayed on this page.

9.3 Install and load required packages

repos <- c(
  rgedi = "https://carlos-alberto-silva.r-universe.dev",
  CRAN = "https://cloud.r-project.org"
)
need <- c(
  "rGEDI", "reticulate", "sf", "terra", "data.table",
  "randomForest", "leaflet"
)
missing <- need[!vapply(need, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing)) install.packages(missing, repos = repos, dependencies = TRUE)

suppressPackageStartupMessages({
  library(rGEDI)
  library(data.table)
  library(sf)
  library(terra)
})

9.4 Read AOI and define the date range

aoi <- sf::st_make_valid(sf::st_transform(study_area, 4326))
aoi_box <- sf::st_bbox(aoi)
study_extent <- unname(aoi_box[c("xmin", "xmax", "ymin", "ymax")])
daterange <- c("2019-04-18", "2019-04-19")
start_year <- 2019
end_year <- 2019

9.5 Package configuration

Configure the Python environment once. Authentication remains in the standard user-level stores managed by Earthdata and the Earth Engine API; it is never placed in this README or the workflow script.

rGEDI_configure(install = TRUE)

# Reads NETRC or ~/.netrc. Do not put usernames or passwords in scripts.
earthdata_login()

# Set this in the R session or the user's environment, outside the repository.
# Sys.setenv(EE_PROJECT = "your-google-cloud-project")
ee_project <- Sys.getenv("EE_PROJECT", unset = "")
if (!nzchar(ee_project)) stop("Set EE_PROJECT to your Google Cloud project ID.")
ee <- ee_initialize(project = ee_project)

9.6 Build the AlphaEarth predictor stack and extract footprint values

Start from the quality-filtered Level 4A footprints extracted in Section 5.2. The spacedSampling() step reduces spatial clustering before predictors are sampled from the 10 m AlphaEarth embedding at a 30 m working scale.

level4a_aoi <- clipLevel4AGeometry(level4a_footprints, aoi)
level4a_aoi <- level4a_aoi[is.finite(agbd) & agbd >= 0]

set.seed(42)
footprint_sample <- sampleGEDI(
  level4a_aoi,
  spacedSampling(size = min(500L, nrow(level4a_aoi)), radius = 25)
)
footprint_sample$sample_id <- seq_len(nrow(footprint_sample))

stack <- ee_build_AlphaEarth_embedding_terrain_stack(
  aoi, start_year, end_year, add_lonlat = TRUE
)
bands <- reticulate::py_to_r(stack$bandNames()$getInfo())
predictor_names <- c(
  grep("^A[0-9]{2}$", bands, value = TRUE),
  intersect(
    c("elevation", "slope", "aspect", "longitude", "latitude"),
    bands
  )
)

training <- extractEE(
  stack$select(as.list(predictor_names)),
  footprint_sample, scale = 30, chunk_size = 250
)

# Restore coordinates if Earth Engine returned only feature properties.
coordinates <- as.data.table(footprint_sample)[, .(
  sample_id, lon_lowestmode, lat_lowestmode
)]
if (!all(c("lon_lowestmode", "lat_lowestmode") %in% names(training))) {
  training <- merge(training, coordinates, by = "sample_id",
                    all.x = TRUE, sort = FALSE)
}

complete <- training[
  complete.cases(training[, c("agbd", predictor_names), with = FALSE])
]
fwrite(complete, file.path(outdir, "gedi-alphaearth-training.csv"))
head(complete[, c("agbd", head(predictor_names, 6)), with = FALSE])

9.7 Visualize the predictor stack as false-color RGB

rgb <- stack$select(c("A00", "A20", "A40"))
rgb_map <- map_view(
  list(`AlphaEarth RGB` = rgb),
  vis = list(`AlphaEarth RGB` = list(
    bands = c("A00", "A20", "A40"), min = -0.06, max = 0.12
  )), aoi = aoi
)
rgb_map

9.8 Variable selection with RFE

RFE repeatedly removes the least useful predictor and chooses the smallest subset within one standard error of the minimum out-of-bag error. Green bars in the importance plot identify the retained variables.

selection <- varSel(
  complete[, ..predictor_names], complete$agbd,
  method = "rfe", threshold = 0, seed = 42,
  ntree = 200
)
best_predictors <- selection$selvars
print(best_predictors)

par(mfrow = c(1, 2), mar = c(4.2, 7, 3, 1))
plot(selection, which = "importance")
plot(selection, which = "rfe")
par(mfrow = c(1, 1))

9.9 Train/test split and fit a Random Forest model

Use a reproducible 70%/30% train/test split to estimate performance on GEDI footprints that were not used to fit each validation model. The final model stored in fit$model is refitted with all complete observations for subsequent prediction.

fit <- fit_model(
  complete[, ..best_predictors], complete$agbd,
  method = "randomForest", test = "split", test_size = 0.30,
  seed = 42, ntree = 500
)
fit$stats_train
fit$stats_test

ok <- is.finite(fit$validation)
plot(
  fit$response[ok], fit$validation[ok],
  pch = 21, bg = "#1fa187", col = "#173f5f",
  xlab = "Observed GEDI AGBD (Mg/ha)",
  ylab = "Holdout prediction (Mg/ha)"
)
abline(0, 1, col = "#d1495b", lwd = 2)

9.10 Create a wall-to-wall aboveground biomass map in GEE

The local holdout model above provides accuracy diagnostics. For scalable wall-to-wall prediction, train an equivalent regression forest inside Earth Engine using the same response and selected predictors.

ee_columns <- c(
  "agbd", "lon_lowestmode", "lat_lowestmode", best_predictors
)
training_vect <- to_vect(
  complete[, ..ee_columns],
  lon = "lon_lowestmode", lat = "lat_lowestmode"
)
training_ee <- vect_as_ee(training_vect)
training_ee <- training_ee$filter(
  ee$Filter$notNull(as.list(c("agbd", best_predictors)))
)
forest <- build_ee_forest(
  training_ee, response = "agbd", predictors = best_predictors,
  trees = 500, seed = 42
)
agbd_map <- map_create(forest, stack$select(as.list(best_predictors)),
                       aoi = aoi, name = "agbd")

9.11 Visualize the aboveground biomass map

agbd_view <- map_view(
  list(`Predicted AGBD` = agbd_map),
  vis = list(`Predicted AGBD` = list(
    min = 0, max = 200,
    palette = c("#f7fcf5", "#74c476", "#00441b")
  )), aoi = aoi
)
agbd_view

9.12 Export the map to GeoTIFF via Google Drive

Use a Drive task for large exports. start = TRUE starts it immediately and ee_check_task_status() reports its state.

drive_task <- ee_image_to_drive(
  agbd_map,
  description = "rGEDI_AGBD_2019", folder = "EE_Exports",
  file_name_prefix = "rGEDI_AGBD_2019", region = aoi,
  scale = 30, start = TRUE
)
ee_check_task_status(drive_task, quiet = FALSE)

# Small images can be downloaded directly:
map_download(
  agbd_map, file.path(outdir, "rGEDI_AGBD_2019.tif"),
  region = aoi, scale = 30, overwrite = TRUE
)

10 Close the files

Close every local or streamed HDF5 object after use.

close(level1b)
close(level2a)
close(level2b)
close(level1b_full)
close(level2a_full)
close(level2b_full)
close(level4a_full)
close(level4a_cloud)
invisible(lapply(orbit_h5, close))
close(level1b_clip)
close(level2a_clip)
close(level2b_clip)
lapply(level1b_geom, close)
lapply(level2a_geom, close)
lapply(level2b_geom, close)
close(sim_clean)
close(sim_noisy)

References

• Dubayah, R. et al. (2020). The Global Ecosystem Dynamics Investigation: High-resolution laser ranging of the Earth's forests and topography. Science of Remote Sensing, 1, 100002.
• Hancock, S. et al. (2019). The GEDI simulator: A large-footprint waveform lidar simulator for calibration and validation of spaceborne missions. Earth and Space Science, 6, 294–310.
• Steven Hancock's reference implementation: https://bitbucket.org/StevenHancock/gedisimulator/src/master/
• GEDI project and product guides: https://www.earthdata.nasa.gov/data/projects/gedi

Acknowledgements

GEDI data are provided by NASA's Land Processes Distributed Active Archive Center. The package includes simulator work originally developed by Caio Hamamura and contributors to the GEDI waveform-simulation ecosystem.

Reporting issues

Report reproducible problems at https://github.com/carlos-alberto-silva/rGEDI/issues. Include the rGEDI version, operating system, product/version, a minimal example, and the complete error message. Never include Earthdata or Google credentials.

Citing rGEDI

citation("rGEDI")

Silva, C. A. et al. (2020). rGEDI: NASA's Global Ecosystem Dynamics Investigation (GEDI) data visualization and processing. Remote Sensing, 12(19), 3208. https://doi.org/10.3390/rs12193208

Disclaimer

rGEDI is an independent open-source project and is not an official NASA software product. Users remain responsible for checking product quality flags, release notes, model assumptions, and the suitability of outputs for their application.