Simulating GEDI footprints from discrete-return airborn LiDAR data.
Author
Xiuyu Cao
Published
September 29, 2026
Modified
September 29, 2026
1 Intro
Large-footprint waveform lidar (e.g., NASA’s GEDI) can measure terrestrial vegeation structure by recording the full shape of a laser pulse’s return within a broad footprint. However, space lidar missions often provides sparse spatio-temporal coverage of the Earth’s surface. Simulating large-footprint waveform lidar from discrete-return airborne lidar can help address this gap, supporting the cal/val of these space missions.
This blog records the workflow of simulating GEDI-like footprints from discrete-return airborne lidar data. Blair and Hofton (1999) showed that discrete-return lidar can be used to simulate large-footprint waveforms. This blog refers to the workflow described by Hancock et al. (2019), which builds on Blair and Hofton’s work.
2 Study area and data
This blog focuses on the Wind River Experimental Forest (WREF) in Washington state, which is one of the NEON sites.
2.1 Discrete-return lidar
The discrete return airborne lidar used in this blog is from NEON. The NEON Airborne Observation Platform (AOP) discrete return lidar point cloud is provided in a compressed American Society for Photogrammetry and Remote Sensing (ASPRS) LAS format (LAZ), referenced to a UTM map projection and ITRF00 datum horizontally and NAVD88 (Geoid12A realization) vertically. It provides the X, Y, Z coordinates and intensity for each laser return point. AOP discrete LiDAR is collected at several points per square meter, and each point can have multiple returns. Source: NEON.
Get NEON lidar info
neon_folder ="../data/raw/neon/neon_wref2021/LAZ_ground"files =sorted({str(p) for pat in ("*.las", "*.laz") for p in Path(neon_folder).glob(pat)})info = lt.get_info(files[0])pprint(info)
{'crs': <Projected CRS: EPSG:32610>
Name: WGS 84 / UTM zone 10N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).
- bounds: (-126.0, 0.0, -120.0, 84.0)
Coordinate Operation:
- name: UTM zone 10N
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
,
'max_x': np.float64(571592.527),
'max_y': np.float64(5076984.737),
'max_z': np.float64(1648.8700000000001),
'min_x': np.float64(570364.084),
'min_y': np.float64(5076249.277),
'min_z': np.float64(533.995),
'point_count': 6883298,
'time': datetime.date(2021, 9, 10)}
I downloaded real spaceborne waveform lidar (GEDI L1B and L2A) from NASA Earthdata using the NEON lidar coverage as the bounding box, with the temporal frame similar to when the NEON data was collected.
Bounding box:
45.5, -122.2
46.0, -121.5
Temporal range: 2021-07-01 - 2021-10-01
Get GEDI shots
gedi_l1b_folder ="../data/raw/gedi/neon_wref/GEDI01_B_003-20260907_162920/"gedi_l2a_folder ="../data/raw/gedi/neon_wref/GEDI02_A_003-20260907_175459/"# Get and quality filter all shotsshots = gt.get_filtered_shots_batch(gedi_l1b_folder, get_shots_only=True)# Filter the shots to the boundaryshots_gpd = gt.get_shot_geometry(shots)coverage2filter = coverage_combined_gdf.to_crs(shots_gpd.crs)shots_gpd_within = gpd.sjoin(shots_gpd, coverage2filter, how='inner', predicate='within')valid_ids = shots_gpd_within["shot_number"]shots_filtered = [s for s in shots if s["shot_number"] in valid_ids.values]# Prepare the simulation function inputshots_input = gt.get_simulation_input(shots_filtered, crs_to ='32610')
The shots_filtered variable contains the complete shot data, including the waveform, shot id, as well as geolocations. While the shots_input is modified for the later simulation function input, and thus only include the beam number and geolocation information. They are one-by-one matched to each other.
The simulation workflow described here follows Hancock et al. (2019). The simulation package used in this section is developed by John Armston.
3.1 Simulating the clean waveform
Given a footprint geometry, the simulator pulls the corresponding discrete-return lidar points, and then simulates the waveform by convolving the weight of each point, \(I_{w,i}\) with the system pulse shape, \(p(z - z_i)\).
The laser point intensity distribution can be modeled as a Gaussian, weighting the contribution of each ALS point by its distance from the footprint center.
\(x_i, y_i\) = horizontal coordinates of that point
\(x_0, y_0\) = horizontal coordinates of the footprint center
\(\sigma_f\) = width of the footprint (see Table 1)
\(I_i\) = relative weighting to account for partial hits
“count”: \(I_i = 1\) –> all points be weighted equally ignoring partial hits (e.g., Blair and Hofton, 1999).
“frac”: \(I_i = 1/nHits\) –> points be weighted by the number of hits each beam records, assuming that each hit along a laser beam intersects a surfaces of equal area (e.g., Armston et al., 2013)
“int” = proportional to surface area –> assumes that the return laser intensity recorded by ALS systems is proportional to the surface area intersected (e.g., Hancock et al., 2017). This assumption is valid for full-waveform lidar but is often not the case for discrete-return systems over diffuse targets.
This weight for each ALS point, \(I_{w,i}\), can be further adjusted based on the local pulse density. For a given simulated footprint, the ALS pulse density will be variable due to varying scan angles and flight-line overlap. This can be corrected by weighting the contribution of each ALS point by the inverse of the pulse density at that area.
Vertical binning is required for the simulation. Real lidar instruments don’t output continuous, infinite-precision waveforms. Instead, they output digitized signals at a fixed sampling resolution. Therefore, the simulated waveform is binned into a fixed vertical resolution. The binning can be performed before or after the convolution. Convolving before binning is more accurate but more computationally expensive, while binning before convolution is much faster.
The \(\text{binned}[z]\) represents the total summed weight of all ALS points whose height \(z_i\) falls within height bin \(z\). ___
The emitted laser pulse has a finite temporal duration, while the detector has a finite response time. The convolution of these two effects determines the lidar system pulse, which can be approximated by a Gaussian:
Lidar’s SNR can be given in terms of a link margin. Link margin expresses the ratio of the signal threshold to the noise threshold in decibels.
Noise threshold (\(t_n\)): a threshold set to give a certain probability of background noise being above it (false positive).
Signal threshold (\(t_s\)): a threshold set to give a certain probability of a real ground-return peak being below it (false negative).
Both thresholds are found as specific locations on the Gaussian distribution. For sections of pure background noise, the relevant distribution is centered on the mean noise level; for a real ground return, it’s centered on the return’s true signal intensity (e.g., the ground-return amplitude, \(\mu_g\)). Assuming that photon shot noise is constant with varying return intensity, both distributions share the same width \(\sigma_n\). To convert a target probability (say, 5% or 10%) into an actual intensity value, we first find the corresponding z-score, then scale that z-score by \(\sigma_n\) and offset it from the relevant mean. This is how \(t_n\) and \(t_s\) are computed. After that, the link margin can be calculated as:
\(t_s\) = signal threshold –> set to give a 10% probability of a false negative. Can be calculated as an offset from the ground-return amplitude, \(\mu_g\).
\(t_n\) = noise threshold –> set to give a 5% probability of a false positive within a 30 m window.
When \(t_s = t_n\) (i.e., \(linkM = 0\)), the signal is barely distinguishable from the noise, and this is the edge case used to define beam sensitivity (\(b_s\)). \(b_s\) is the canopy cover that we would expect to be able to detect the ground through certain chance of false negatives and false positives. It can be denoted as the fraction of energy contained within a Gaussian with the ground-return amplitude, \(\mu_g\):
where \(\sigma_{eff}\) is the ground return’s effective width
\(\bar{n}\) is the mean noise level.
When the ground is sloped, the ground return’s effective width is increased. This is because different parts of the footprint are at different elevations. So instead of all the ground signal arriving at nearly the same height, it gets spread out over a range of heights. The ground return’s effective width can be modeled as:
\(\sigma_{eff}\) = effective width of the ground return (combines instrument blur and slope-induced spreading)
\(\sigma_p\) = system pulse width (the instrument’s own blur, present even on flat ground)
\(\sigma_f\) = footprint width
\(\theta\) = ground slope
if given ground return fraction information (i.e., canopy cover), \(\mu_g\) can be calculated with the effective width of the ground return, \(\sigma_{eff}\), and \(\sigma_n\) can be calculated according to the target link margin.
3.3 Example
Here I use the real GEDI shot locations retrieved in Section 2.2 to simulate the waveforms from the NEON discrete-return lidar. Note that the NEON lidar data and the GEDI data use different vertical reference systems. For GEDI, the height data is relative to the WGS-84 ellipsoid, while the NEON lidar data is referenced to NAVD88 (Geoid12A realization) vertically (GEDI L1B; NEON APOP). Left uncorrected, this produces a geoid undulation (about −30 m across most of CONUS) plus a smaller horizontal and vertical offset of order 1–2 m from the frame difference, which varies with observation date.
Example input footprint information
n =3print(f"Example - First {n} shot information:")for i inrange(n):print(shots_input[i])
# Global style for publication figuresset_plot_style(theme="ggplot", font_scale=1.2, serif=True)test = matched[50]# Plotfig, axes = plt.subplots(1, 2, figsize=(12, 6))ax1, ax2 = axesgt.plot_waveform(test[0]['rx_raw'], test[0]['tx_raw'], (test[0]["elevation_bin0"], test[0]["elevation_lastbin"]), simulation=False, title ='Real GEDI Shot Waveform', ax=ax1)gt.plot_waveform(test[1].rxwaveform, test[1].txwaveform, test[1].z, simulation=True, min_signal =0.0001, title ='Simulated GEDI Waveform', ax=ax2)
The relative height and shape are similar, but there is an offeset in the elevation. It is due to the different vertical coordinate reference system. Therefore, for precise comparison of the elevations, it is needed to convert the simulated elevations and footprint coordinates to the GEDI system in two steps: (1) adding the GEOID12B undulation to convert orthometric to ellipsoidal heights, and (2) transforming from epoch of NEON data to GEDIV003’s ITRF2020 at each shot’s acquisition date using NOAA’s HTDP.
