rGEDI: An R Package for NASA's Global Ecosystem Dynamics Investigation (GEDI) Data Visualizing and Processing.
Authors: Carlos Alberto Silva, Caio Hamamura, Ruben Valbuena, Steven Hancock, Adrian Cardil, Eben N. Broadbent, Danilo R. A. de Almeida, Celso H. L. Silva Junior, and Carine Klauberg.
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
rGEDI function reference manual is also
available for review.
# 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)
As in the ICESat2VegR workflow, rGEDI uses three Python packages through
reticulate:
earthaccess finds and
authenticates Earthdata Cloud objects.h5py reads HDF5 content streamed from the
cloud.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)
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)
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.
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)
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)
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.
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.
level1b_geo <- getLevel1BGeo(level1b)
head(level1b_geo)
The original rGEDI map is retained below. It shows the Level 1B footprints over high-resolution imagery.
shot <- level1b_geo$shot_number[1]
waveform <- getLevel1BWF(level1b, shot_number = shot)
plot(waveform)
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)]
plotWFMetrics(level1b, level2a, shot_number = shot,
rh = c(25, 50, 75, 90, 98))
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)]
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)
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)
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.
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)
data.table objectslevel1b_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")
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
)
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
)
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)
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
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)
rh98_prediction <- rasterizeGEDI(
predicted, metric = "predicted_rh98", res = 0.002,
lon = "lon_lowestmode", lat = "lat_lowestmode", fun = mean
)
plot(rh98_prediction)
level3 <- readLevel3(level3_file)
level4b <- readLevel4B(level4b_file)
plotLevel3(level3, main = "GEDI Level 3 metric")
plotLevel4B(level4b, main = "GEDI Level 4B AGBD")
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))
level4b_stats <- polyStatsLevel4B(
level4b, vect(study_area), fun = mean, na.rm = TRUE
)
level4b_stats
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. Vectorized bin accumulation keeps this path
portable and fast without the native library stack that caused CRAN problems.
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))
)
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)]
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)]
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.
A clean end-to-end script is installed with the package and can also be opened directly from the repository:
Download the GEDI AlphaEarth upscaling workflow
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.
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)
})
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
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)
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])
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
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))
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)
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")
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
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
)
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)
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.
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.
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
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.