This package offers a set of functions to work with the GEDI Level 4A data. This dataset contains GEDI Level 4A version 2.1 (L4A) predictions of the aboveground biomass density (AGBD; in Mg/ha) and estimates of the prediction standard error within each sampled geolocated laser footprint derived from parametric models that relate simulated GEDI Level 2A (L2A) waveform relative height (RH) metrics to field plot estimates of AGBD. The datasets is available for the period 2019-04-18 to 2021-11-23.
More information on Level 4A data can be found here. To download the GEDI L4A user guide plese visit this link The package follows a simple name convention: all functions names start with l4_ and are followed by a verb indicating the purpose of the function.

#install.packages("devtools")
devtools::install_github("VangiElia/GEDI4R")
#loading GEDI4R package
library(GEDI4R)

l4_download uses the gedifinder to find GEDI Level 4A data intersecting a user-defined extent and date range. The list of the resulting path is downloaded in parallel using the foreach package. Note that if the number of files to download is less than the available number of cores, the function l4_download will uses as many cores as the numebr of file, for efficency reason. At the first run of l4_download, users will need to enter their Earth Explorer Login Information (if you do not have an account, sign up at this link) for downloading the data. The function will create a netrc file in the directory specified by the argument outdir that stores the access credentials to the ORNL DAAC database, where GEDI Level 4 datasets are stored. l4_download will search for the netrc file in outdir at each run. If it does not find it, it recreates it by asking the user to enter the login credentials again. As long as this file is in outdir, the user does not need to enter this information again. If the user changes outdir, a new netrc file will be created after entering the access credentials.

#using Italy as bounding box for search GEDI data
ita <- sf::st_as_sf(raster::getData('GADM', country = 'ITA', level = 1))
#extract extent
e <- raster::extent(ita)
ul_lat <- e@ymax
lr_lat <- e@ymin
ul_lon <- e@xmin
lr_lon <- e@xmax
outdir = tempdir()
#Get the list of all files available for the study area in the period selected,
#using just_path = T
file_path <- l4_download(
  ul_lat,
  lr_lat,
  ul_lon,
  lr_lon,
  outdir = outdir,
  from = "2020-01-01",
  to = "2020-01-31",
  just_path = T
)
#download the first 4 files
file_download <- l4_download(
  ul_lat,
  lr_lat,
  ul_lon,
  lr_lon,
  ncore = parallel::detectCores()-1,
  outdir = outdir,
  from = "2020-01-01",
  to = "2020-01-31",
  just_path = F,
  subset=1:4
)

After downloading, files can be read from the original file format (h5) with l4_getmulti as data.table objects. The function can accept a list or a vector of file paths (as the output of l4_download) and read them in parallel, using the snowfall package. If the list of file path has lenght=1 the file will be read in single thread mode. The function remove by default footprints with AGBD values corrupted (agbd<0), and can be used to filter footprints based on the tree cover threshold derived for the year 2010, from Hansen et al. (2013) and encoded as a percentage per output grid cell.
The functions can list the dataset names inside the structure of h5 files. This is useful for adding other columns to the functions's default output datasets by specifying the argument add_col. See the Details section of ?l4_getmulti for the default dataset extracted from the h5 file.

outdir = tempdir()
l4_zip <- system.file("extdata",
                      c("GEDI04_A_2020036151358_O06515_02_T00198_02_002_01_V002.zip",
                        "GEDI04_A_2021150031254_O13948_03_T06447_02_002_01_V002.zip"
                      ),
                      package="GEDI4R")
#Unzipping GEDI level4A data
l4 <- lapply(l4_zip,unzip,exdir = outdir)
#list all dataset in h5 file
dataname <- l4_getmulti(l4[[1]],just_colnames = T)
head(dataname,10)
#read all footprint and merge them together.
gediL4_path <- l4_getmulti(l4,merge=T)
#select other columns to add to the default output.
#if columns are already present in the default output they will be dropped
col <-
  c("land_cover_data/leaf_off_flag",
    "agbd_pi_lower",
    "agbd_pi_upper",
    "agbd"#this will be dropped as it is already included by default in the output.
    )
#get level 4 data with the user defined column binded and with the source path of each observation
#with source=T a column with the source path for each observation will be added
gediL4 <- l4_getmulti(l4,add_col = col,source=T)

knitr::kable(head(gediL4[,c("date","tree_cover","agbd","agbd_se")]))

After downloading and reading GEDI data, a typical pre-processing step is to clip the data to the extent of a study area. To do so, l4_clip can be used. The function clip footprints by extent or vector boundary provided by the argument clip. Currently, it accepts a path to a shp or tif file, an object of class sf, a Raster* object, a numeric vector of coordinates or other objects from which an extent can be extracted. By specifying the argument usegeometry=TRUE foorprints will be clipped on the boundary of an sf object (or path from which an sf object can be created).
GEDI coordinates are by default in lon/lat format (EPSG 4326). The function will try to convert the extent of clip to lon/lat coordinate system to ensure compatibility during the clip. The only exception is when clip is a numeric vector or a bbox object. In these cases, the user must check that the extent is in lon/lat projection.

outdir = tempdir()
l4_zip <- system.file("extdata",
                     "GEDI04_A_2020036151358_O06515_02_T00198_02_002_01_V002.zip",
                     package="GEDI4R")
l4 <- unzip(l4_zip,exdir = outdir)
#get GEDI data
l4_data <- l4_getmulti(l4)
#clip using vector of coordinates
b_box <- c(-50,35,52,37)
clipped <- l4_clip(l4_data,clip=b_box)
#using Shapefile to clip
bound <- system.file("extdata","bound4326.shp",package="GEDI4R")
#with  extension
clipped <- l4_clip(l4_data,clip=bound,usegeometry = F)
#with  polygon boundaries
clipped2 <- l4_clip(l4_data,clip=bound,usegeometry = T)

After pre-processing, the next step is converting and exporting the data in a vector format, usually an ESRI Shapefile. This can be easily accomplished with the function l4_convert, which can also reproject the data to a user-defined coordinate reference system (specified via the EPSG code by the argument epsg). Note that in converting data to ESRI Shapefile, columns names will be abbreviated with a warning. The function is called for its side effects. It return NULL unless return_path=TRUE.

converted <- l4_convert(l4_data,epsg = 4326,filename=paste0(outdir,"/example.shp"),return_path = T,append=F)
example <- sf::read_sf(converted)
file.remove(list.files(outdir,pattern = "example",full.names = T))

Finally, the package implements two functions for plotting the data: l4_plotagb and l4_plotprofile.
The former function can plots the location of footprints, the distribution of AGBD against the elevation, or both.
The latter function returns the AGBD against the elevation profile along the GEDI track. Note that plotting elevation profiles from GEDI data (l4_plotprofile) is only advisable for single beam/track pairs. Plotting profiles from a data.table concatenating multiple GEDI files (orbits) can be misleading by forcing an overlap of data from tracks at different locations.

gediL4 <- l4_getmulti(l4)
#footprints locations and AGBD distribution against elevation
l4_plotagb(gediL4,,beam_id="all",type = "both",n=100,h=c(100,100))

#along-track AGBD profile
l4_plotprofile(gediL4,beam_id="all")

All the above pre-processing steps can be performed at once in chunks of files with the function l4_process. It runs, in order, l4_getmulti, l4_clip and l4_convert to each file chunk.
By specifying parallel=TRUE, the function will process each chunk in parallel trying to guess the best number of cores to be used in l4_getmulti and in the chunks loop, based on the maximum number of cores available and the number of files to be processed. The user can override the number of cores used in l4_getmulti with ..., by specifying the argument ncore. This will also affect the number of cores used to loop over chunks. Usually, the number of cores used by default is the best option. Modifying it can slow down the function.

outdir = tempdir()
l4_zip <- system.file("extdata",
                      c("GEDI04_A_2020036151358_O06515_02_T00198_02_002_01_V002.zip"
                      ),
                      package="GEDI4R")
#Unzipping GEDI level4A data
l4 <- unzip(l4_zip,exdir = outdir)
#create 4 copy of GEDI file to test the function
file.copy(from=l4,to=paste0(tools::file_path_sans_ext(l4),"_copy",1:4,".h5"))
#path to Shapefile for clipping the data
bound <- system.file("extdata","bound4326.shp",package="GEDI4R")
#path to GEDI files
l4_path <- list.files(outdir,pattern = "h5",full.names = T)
#proces all files in chunk each of 2 files, in sequence
l4_data <- l4_process(l4_path,nfile=2,clip=bound,epsg=4326,outdir=outdir,ext="shp",parallel=F,prefix = "ex")
file.remove(list.files(outdir,full.names = T, pattern = "ex"))
#in parallel
l4_data <- l4_process(l4_path,nfile=2,clip=bound,epsg=4326,outdir=outdir,ext="shp",parallel=T)
file.remove(list.files(outdir,full.names = T, pattern = "ex"))

The output of l4_process is the vector of file path listed in outdir. The number of resulted vector files is equal to the original number of h5 files divided by the argument nfile. The vector files thus obtained can simply be merged together, both in R or in any GIS software.

Please report any issue regarding the GEDI4R package at https://github.com/VangiElia/GEDI4R/issues