# kNN - Bagging NFI, GEDI, Sentinel-2 and Sentinel-1 data to produce estimates of forest volumes

Anouk Schleich, Cédric Vega, Jean-Pierre Renaud, Olivier Bouriaud, Sylvie Durrieu

► To cite this version:

Anouk Schleich, Cédric Vega, Jean-Pierre Renaud, Olivier Bouriaud, Sylvie Durrieu. kNN - Bagging NFI, GEDI, Sentinel-2 and Sentinel-1 data to produce estimates of forest volumes. Forest Ecology and Management, 2025, 595, pp.122964. 10.1016/j.foreco.2025.122964. hal-05216003

HAL Id: hal-05216003

https://hal.inrae.fr/hal-05216003v1

Submitted on 20 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# kNN - Bagging NFI, GEDI, Sentinel-2 and Sentinel-1 data to produce estimates of forest volumes

Anouk Schleich<sup>a, ∗</sup>, Cédric Vega<sup>c</sup>, Jean-Pierre Renaud<sup>d, e</sup>, Olivier Bouriaud<sup>b, e</sup>,
Sylvie Durrieu<sup>a, f</sup>

ARTICLE INFO

Keywords:

National forest inventory
Multi-source forest inventory
kNN
Bagging
GEDI
Spaceborne lidar

A B S T R A C T

This study presents a k-nearest neighbor bagging approach for predicting forest attributes, specifically growing stock volume (GSV), by integrating optical Sentinel-2, radar Sentinel-1, and spaceborne lidar GEDI datasets following a two-step procedure. Over a large region (~ 4600 km²), some selected GEDI variables were imputed to national forest inventory (NFI) plots, then models were developed to predict GSV based on the combined imputed GEDI and Sentinel variables. For GEDI imputation, three strategies were followed, each one using different auxiliary datasets: (A) solely relying on Sentinel data, (B) using Sentinel data complemented by maximum height values from both GEDI and NFI data sources, (C) using Sentinel and height values from an independent canopy height map derived from GEDI and Sentinel through Deep Learning. For each strategy, the models were built using a k-nearest neighbor approach and evaluated on an independent test dataset. The results indicated that the ability of Strategy A to predict forest volume as the percentage of explained variance (R²) was considerably limited (0.08) and the relative root mean squared error (RMSE%) was high (66%). Upon adding heights, Strategy B outperformed the previous strategy, with a R² of 0.58 and a relatively low RMSE% (43%). Compared to Strategy A, the use of estimated heights from a canopy height map improved predictions, but did not outperform Strategy B; R² was 0.34 and RMSE% 53%. Imputations yielded better models when using original GEDI footprint heights as inputs, rather than relying on height estimates extracted from the canopy height map. Therefore, small area model-assisted estimation may be preferred over wall-to-wall maps to reduce area-based local estimation errors.

1. Introduction

National forest inventories (NFIs) are designed to produce precise forest resource estimates across vast areas (i.e. national to regional). Various methods for enhancing the precision of estimates over small domains while controlling costs have been developed. Most of these methods are based on remote sensing data or thematic maps. Examples include stratification approaches, where auxiliary data are used to define homogeneous estimation strata (Haakana et al., 2019; Schleich et al., 2024), and multisource forest inventory (MFI) (Tomppo et al., 2008; Saborowski et al., 2010; Westfall et al., 2019), which relies on a model to relate field attributes to auxiliary data, allowing for the derivation of high-resolution maps of forest attributes to support small-area estimation (Guldin, 2021). Among the attributes observed by NFIs, growing stock volume (GSV) plays a key role in providing essential information for policy decisions and forest management (Gschwantner

et al., 2022). GSV, which corresponds to the amount of living stand-
ing wood per hectare, indicates the availability of wood resource, an
economically and environmentally important variable (Gschwantner
et al., 2019). Every year, the French NFI publishes GSV estimates at
regional scales using 5 years of data collection. While such information
is relevant for national forestry policy purposes, sub-regional to local
estimates of GSV are desired for forest management purposes.

Within the multisource estimation framework, moderate-resolution
imagery (10-30 m) has been used for its capacity to provide national-
scale data on a yearly basis, and its ability to ensure long-term data re-
newal while having a pixel size in the order of the NFI plots size (Coops
et al., 2023; Tomppo et al., 2008). However, optical data have limited
correlation with forest attributes in diverse and complex forest struc-
tures (Irulappa Pillai Vijayakumar et al., 2019; Saarela et al., 2018).
Three-dimensional remote sensing data from airborne laser scanning

<footnote><sup>a</sup> UMR TETIS, Univ Montpellier, AgroParisTech, CIRAD, CNRS, INRAE, 500 Rue Jean-François Breton, Montpellier, F-34196, France</footnote>

<footnote><sup>b</sup> University Stefan cel Mare of Suceava, 13 University Street, Suceava, 720229, Romania</footnote>

<footnote><sup>c</sup> ENSG, IGN, Laboratoire d'Inventaire Forestier, 14 Rue Girardet, Nancy, F-54042, France</footnote>

<footnote><sup>d</sup> Office National des Forêts, Pôle Recherche et Développement Innovation, 8 Allée de Longchamp, Villier-les-Nancy, F-54600, France</footnote>

<footnote>* Corresponding author.
E-mail address: anouk.schleich@gmail.com (A. Schleich).</footnote>
(ALS) or photogrammetry, which provide key information about forest canopy structure, are more suitable for estimating structure-related forest attributes such as volume, basal area or biomass (Lim et al., 2003). Various European countries have either completed or started their national ALS coverage to support these needs (Virtanen et al., 2017; Novo-Fernández et al., 2019; Kakoulaki et al., 2021). A key challenge associated with such aerial surveys is ensuring a consistent and timely update at the national level. Nordic European countries, in particular, have engaged regular ALS acquisitions to that aim (Appiah Mensah et al., 2023). Renewal of 3D data from aerial imagery is also a sustainable solution owing to the regular surveys carried out by most national mapping institutes for several decades (Ginzler and Hobi, 2015). The resulting long time series are also beneficial for the assessment of growth trajectories, which are valuable in monitoring forest dynamics (Véga and St-Onge, 2009). Nevertheless, the time needed to collect data over an entire country continues to be an impediment. While temporal aggregation may be adequate for nationwide estimates of forest resources, it is relatively ineffective in addressing rapid changes in forests, such as those induced by significant disturbances like storms or fires, the frequency and severity of which tend to increase owing to climate change (Forzieri et al., 2021). A potential solution lies in the use of spaceborne lidar data, which allows for the monitoring of large areas at a higher temporal frequency, thereby better meeting emerging needs for monitoring changes in forest dynamics.

In 2018, the Global Ecosystem Dynamics Investigation (GEDI) space-borne lidar system was launched to collect forest structure data (Dubayah et al., 2020). GEDI employs full waveform lidar sensors to capture data within ~ 25 m diameter footprints. Its sampling scheme, designed to capture the variability of forest conditions within an inferential framework for producing statistics at the kilometer scale (Saarela et al., 2018), shows strong potential for enhancing MFIs (Pascual et al., 2023; Bruening et al., 2023; Bullock et al., 2023). In contrast with passive optical systems, GEDI data include several vegetation structural variables such as height profiles and canopy cover, which are highly correlated with forest attributes of interest such as volume or biomass.

GEDI's L2 A and L2B products provide several variables at footprint level, such as canopy heights and coverage. The L3 product provides 1×1 km gridded mean and standard deviation of canopy height (Dubayah et al., 2021), while L4 A and L4B products give AGBD respectively at footprint and at 1×1 km grid level (Dubayah et al., 2022, 2023). The underlying models were developed using GEDI waveforms simulated from ALS data and coinciding with field inventory data, obtained from on a unique compilation of forest sites distributed worldwide (Dunanson et al., 2022). Model specification issues have been reported, highlighting the need to increase the representativeness of forest conditions and to standardize data collection protocols (Labrière et al., 2023). Meanwhile, methods have been proposed to correct estimation biases at regional scales (Bruening et al., 2023; Pascual et al., 2025).

Developing local biomass models requires additional field data and faces the challenge of limited spatial overlap between GEDI footprints and field plots. For the GEDI L4 products this was addressed through ALS-based simulations. In the era of deep learning, research on approaches leveraging deep learning models tend to become the mainstream approach for spatializing GEDI attributes wall-to-wall using remote sensing imagery. Potapov et al. (2021) created a 30 m resolution global canopy height map based on GEDI and Landsat data. Lang et al. (2023) created a similar map, at a finer resolution of 10 m, using GEDI and Sentinel-2. Recently, Schwartz et al. (2023) created a 10 m resolution canopy height map, called FORMS-H, for the French mainland territory based on GEDI, Sentinel-2, and Sentinel-1, using a deep learning approach. Height information can then be combined with field plot data to develop allometric models and produce high-resolution forest attribute maps (Silva et al., 2021; Schwartz et al., 2023) to further support inference. A current limitation of these approaches is that GEDI's vertical information is summarized by a single

height metric, such as the 99th or the 95th percentile of the relative height profile.

To fully leverage the rich structural information of GEDI data to improve biomass or volume predictions, the integrity of this information must be preserved and linked to reference field measurements. Several alternative approaches to ALS-simulated GEDI attributes have been developed to that end. Indeed, the availability of ALS data contemporaneous with field plot measurements remains a major constraint, making the simulation-based approaches not always reproducible. Zhang et al. (2022) and Bullock et al. (2023) matched GEDI footprints to field plots using a nearest neighbor algorithm within a defined search radius. However, the limited availability of field plots, combined with the irregular spatial distribution of GEDI footprints, prevents many field plots from being paired with GEDI data, resulting in an important reduction in the size of the training datasets. To address this limitation, May et al. (2024) proposed a spatial probability-based approach that considers the 25 closest GEDI footprints to each NFI plot. Nonetheless, in complex forest landscapes with structural discontinuities, the optimal footprint may lie outside the selected GEDI subset, potentially limiting model performance. To account for this, Besic et al. (2024) proposed a method to stratify GEDI and NFI structural classes, aiming to improve plot-footprint matching by considering both forest structure and spatial distance. A weighting function is used to balance the contributions of geographic space (i.e. latitude and longitude) and auxiliary data space (i.e., remote sensing imagery and thematic maps) in order to optimize the matching process, thereby relaxing the spatial distance constraint. A key strength of this approach lies in its capability to generate outputs suitable for either stratification or regression, depending on the requirements for inference.

We consider that these types of matching approaches could be highly valuable for augmenting NFI data with GEDI footprints and either support directly model-assisted small area estimation or serve as a primary step toward the propagation of NFI attribute spatially using continuous remote sensing data such as Sentinel or GEDI through computationally intensive deep learning approaches. Namely, while several studies have used machine learning or deep learning algorithms in order to propagate GEDI height on continuous data such as Sentinel images, before modeling volume or biomass (Schwartz et al., 2023; Potapov et al., 2021), we consider that it would be more efficient to predict NFI attributes at GEDI footprints before training deep learning methods to produce continuous maps.

From this statement, the main objective of this study is to propose a matching approach to propagate NFI data to GEDI footprints. We propose a flexible matching approach, taking the computational efficiency and simplicity of a k-nearest neighbor (kNN) strategy, such as in Zhang et al. (2022), with a bagging approach (Cover and Hart, 1967), to both generate a distribution of matches (May et al., 2024) and relax the distance constraints (Besic et al., 2024) through resampling. A secondary objective is to evaluate different sources of auxiliary data to support the matching.

To test our approach, we use the NFI variable GSV as the target. Therefore, the objective is to predict GSV using GEDI and Sentinel data. Sentinel data have been successfully used to predict a specific GEDI height, such as the 95th or 100th percentiles of relative heights (e.g. RH95 or RH100) measured between the ground and the maximum height (Schwartz et al., 2023; Pereira-Pires et al., 2021). However, the use of a single height variable is limited in its capability to predict structure-based NFI parameters such as GSV, as height is only a component of GSV. Propagating multiple GEDI heights and other GEDI vegetation derivatives would be more appropriate to that end. However, owing to the tenuous link between Sentinel information and forest structure, as Besic et al. (2024), we hypothesized that Sentinel data used alone might fail to efficiently bridge this gap and that it needs to be completed with structural information. We assumed that the use of forest height in addition to Sentinel information could constitute an
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Footprints were filtered to intersect with BD Forêt polygons. Moreover, multiple quality-based filters were applied, resulting in ~105 000 GEDI footprints. The NFI data consist of 675 plots.</figcaption></figure>

efficient strategy for successfully propagating NFI information at the level of GEDI footprints.

To evaluate these hypotheses, we tested three matching strategies. Strategy A only relies on Sentinel information. Strategy B adds structural information by using height information available both in NFI and GEDI datasets. Specifically, Schleich et al. (2024) found that, among various relative height (RH) candidates, GEDI RH100 and NFI maximum height were correlated best to one another and could be used as a shared bridge variable between the two datasets. Strategy C also combines Sentinel and height information, but unlike Strategy B, it uses height values from an independent canopy height map (FORMS-H) derived from GEDI and Sentinel data using Deep Learning (Schwartz et al., 2023). While Strategy B is constrained to predict GSV only at the level of GEDI footprints, preventing the realization of wall to wall predictions, Strategies A and C enable wall-to-wall predictions by leveraging spatially continuous information.

## 2. Study site and data

### 2.1. Study site

Based on bioclimatic and ecological criteria, metropolitan France is divided into 91 sylvoecoregions (Cavaignac, 2009). This geographical subdivision serves as a national reference for forest management. Our study site is located in “Central Vosges Massif”, a mountainous sylvoecoregion of North-Eastern France. Elevations range from 200 to 1400 m. The stand compositions differ considerably, depending on site conditions (altitude, slope, exposure, soil type, available water, etc.) The forest is primarily composed of European beech, silver fir, and Norway spruce trees, either in pure or mixed stands (beech-fir stands) (IGN, 2023c). The Vosges study site used to develop the proposed MFI approach is defined by the intersection of the extent of a Sentinel-2 tile with the “Central Vosges Massif” sylvoecoregion (Fig. 1(a)). It covers an area of 4634 km².

### 2.2. Data

Three main datasets were used: GEDI products, NFI field plots and Sentinel (1 & 2) images (Fig. 1(b)). We also used a forest map (IGN, 2023a) to define a forest mask, and a 1-m digital elevation model (IGN, 2022) to improve GEDI footprint locations. Furthermore, a 10 m resolution height map (FORMS-H, Schwartz et al., 2023) was used to obtain spatial information on height, a crucial structural variable, which is only available in discrete form with GEDI footprints. Finally, a 30 m resolution wood volume map (FORMS-V, Schwartz et al., 2023) was used to compare final GSV estimates. All data are summarized in Table 1.

#### 2.2.1. National forest inventory plots

The French NFI is a two-phase continuous inventory, based on a 1 km sampling grid (Bouriaud et al., 2023). Each year, one-tenth of the grid cells is covered. The first phase sample is drawn from the yearly grid fraction and constitutes photo-interpreting aerial photographs to estimate forest cover (~100,000 points/year). The second phase sample is drawn from the first phase one classified as forest and comprises field plots measurements (~7500 points/year) for estimating forest resource attributes. Tree measurements are made in concentric plots of 6, 9 and 15 m radii according to their diameter at breast height (DBH; [7.5;22.5[, [22.5;37.5[, [37.5; +∞[(in cm), respectively). Along with the species identification and vitality information (dead or alive), three main variables are measured: DBH, total height (H), and stem height with up to a 7 cm stem diameter (Hdec, used to determine the solid stem volume). Species, vitality, and DBH information are collected for all trees, and both H and Hdec are measured for a single representative tree within each species and DBH classes.

The GSV is subsequently calculated using species-specific volume equations and imputation methods involving the tree DBH, H, and Hdec. Plot-level attribute density values are then computed by aggregating individual tree data and using tree inclusion probability. Official statistics are consolidated using five-year moving averages (IGN, 2023b). For this study, plots from 2018 to 2022 were used. Over the area of interest, 675 plots were inventoried and the mean GSV of the Sylvoecoregion was estimated to be 295 ± 20 m³ ha⁻¹. For comparison, the mean GSV over France was estimated to be 173 ± 3 m³ ha⁻¹ over the same period.

#### 2.2.2. GEDI data

GEDI is a spaceborne lidar system mounted on the International Space Station. It is an active sensor, emitting 1064 nm laser pulses and measuring the backscattered signal over ~ 25 m diameter footprints on the ground (Dubayah et al., 2020). Launched by NASA in 2018, during its initial phase, the GEDI mission collected data from April 2019 to March 2023. The instrument was subsequently moved, temporarily pausing data acquisition. GEDI resumed operations and began collecting data again in April 2024 (LP DAAC, 2023). This study uses data collected during leaf-on seasons, each spanning from May to September, of the first phase of the mission.

Various filters were applied to ensure footprint data quality. Footprints acquired using full-power lasers with good degrade and quality flags were retained. Additional filters based on GEDI Level 2 A (L2 A) data proposed by Schleich et al. (2024) were also used. Moreover, GEDI footprints were filtered using a forest mask based on BD Forêt (IGN, 2023a) to identify forest footprints and discard non-forest ones. An additional height criterion was introduced to exclude footprints with a maximum height (RH100) below 7 m or above 60 m. This allows for
<table><caption>**Table 1**
Datasets used in this study.</caption><thead><tr><th>Data</th><th>Description</th><th>Product</th><th>Source</th><th>Variables</th><th>Date</th><th>Filters or preprocessing</th><th>N</th><th>Utility in study</th></tr></thead><tbody><tr><td>RGE Alti</td><td>1-m resolution Digital Elevation Model</td><td>RGE ALTI 1M</td><td>https://geoservices.ign.fr/rgealti</td><td>elev_lowestmode, lat_lowestmode, lon_lowestmode, delta_time</td><td></td><td>quality_flag = 1<br>degrade_flag = 0<br>beam = full power (extent study site)</td><td>268 109 footprints</td><td>to apply<br>GeoGEDI</td></tr><tr><td rowspan="2">GEDI</td><td rowspan="2">full waveform spaceborne lidar 25 m circular footprints</td><td>L2A product v2</td><td>https://lpaac.usgs.gov<br>https://search.earthdata.nasa.gov/search</td><td>Relative heights :<br>RH100, RH90, RH80, RH70, RH60, RH50, RH40, RH30, RH20, RH10</td><td>Apr 2019-Mar 2023<br>Summer: Mai - Sep</td><td>quality_flag = 1<br>degrade_flag = 0<br>beam = full power<br>Additional filters (Schleich et al., 2024)<br>RH100 &lt;= 60<br>RH100 &gt;= 7<br>in BD Forêt</td><td>104 831 footprints</td><td>kNN</td></tr><tr><td>L2B product v2</td><td></td><td>canopy cover cumulative canopy covers pai vertical pai profile pavd profile fhd_normal</td><td></td><td>use footprints filtered in L2A</td><td>104 831 footprints</td><td>kNN</td></tr><tr><td>NFI</td><td>data from french forest inventory</td><td></td><td>https://inventaire-forestier.ign.fr</td><td>growing stock volume (GSV) (i.e. variable named vwac)</td><td>2018-2022</td><td></td><td>675 plots</td><td>kNN</td></tr><tr><td rowspan="2">Sentinel-2 (S2)</td><td rowspan="2">optical multispectral spaceborne images spectral bands in the visible, near infrared (NIR), and short-wave infrared (SWIR)</td><td>L2A product</td><td>https://earthengine.google.com</td><td>bands: B3, B4, B8</td><td>12 Aug 2022</td><td>calculate lai, fCover and fAPAR with R package prosail</td><td>lai, fCover, fAPAR</td><td>kNN</td></tr><tr><td>Theia L3A product</td><td>https://theia.cnes.fr/</td><td>bands: B2, B3, B4, B5, B6, B7, B8, B8A, B11, B12</td><td>Aug 2022<br>Aug 2017<br>Jun 2022</td><td>calculate vegetation indices and time series of indices</td><td>10 bands, vegetation indices, time series of indices</td><td>kNN</td></tr><tr><td>Sentinel-1 (S1)</td><td>Synthetic Aperture Radar (SAR)</td><td>Ground Range Detected (GRD)</td><td>https://earthengine.google.com</td><td>Ascending VV, Ascending VH, Descending VV, Descending VH</td><td>August 2022</td><td>orbits with VV and VH polarizations</td><td>total of 37 orbits: 17 ascending and 20 descending; resumed to one monthly mean by polarization and by direction</td><td>kNN</td></tr><tr><td>BD Forêt</td><td>database of polygons providing information about french forest stands</td><td>BD Forêt v2</td><td>https://geoservices.ign.fr/bdforet</td><td></td><td></td><td></td><td></td><td>to filter GEDI footprints</td></tr><tr><td>FORMS-H</td><td>National Canopy Height Map</td><td></td><td>Schwartz et al. (2023)</td><td></td><td></td><td></td><td></td><td>kNN<br>Strategy C</td></tr><tr><td>FORMS-V</td><td>National Wood Volume Map</td><td></td><td>Schwartz et al. (2023)</td><td></td><td></td><td></td><td></td><td>final comparison</td></tr></tbody></table>

exclusion of the footprints located in areas where logging has occurred and alleviates unrealistic height estimates. After applying all these filters, 104 831 footprints were retained.

Two GEDI products at the 25-m footprint level were used: L2 A data providing relative height (RH) metrics and Level 2B (L2B) data, including estimates of canopy cover (cover), plant area index (pai), plant area volume density (pavd), and foliage height diversity index (fhd_normal). Pavd and both cumulative cover and pai from height z to the ground were available at a vertical step size of $dz = 5$ m. Pai is also cumulative, therefore pai_z2 contains all pai between the ground and a height of 10 m. Pavd, however, constitutes individual layers, pavd_z2 contains pavd from 5 m to 10 m heights.

RH values were transformed to vegetation RH, representing the relative height metrics at 1% interval for the vegetation component of the backscattered signal. To that aim, for each footprint, a simplified waveform was reconstructed from RH values and a Gaussian function was adjusted to the ground return. Then, the ground component of the signal was removed by subtracting the Gaussian ground peak from the waveform. Finally, the resulting vegetation waveform was transformed back into Relative Height metrics, i.e. heights above the ground at which a certain percentile of the energy returned by the vegetation is reached. These refined vegetation RH values will further be referred to as RHv. Additional variables were derived from the RHv values to provide supplementary information related to the shape of the RHv profile and thus to vertical structure. G_cor represents the correlation between the RHv profile and a straight line connecting RHv_0 (i.e. the ground level) to RHv_100. G_Hr_chgt indicates the first height where the RHv

profile intersects the RHv_0-RHv_100 line. This height is expressed as a percentage of the maximum height (G_RHv_100). G Coeff_var is the coefficient of variation of the RHvs calculated as the ratio of the standard deviation to the mean.

Moreover, cover_z, pai_z, and pavd_z variables were used to estimate additional variables. Depending on already available data, the vertical z step values were used to estimate the variables from different perspectives: cumulative (already given for cover and pai) and by Z step values (already given for pavd), from top to bottom, and from bottom to top. Variables from the top to the bottom were organized in alphabetical order, i.e. the pavd from the highest 5 m layer detected in the footprint, were called pavd_za. Layer wise cover (and pai) variables were named cover_z3z4 when including the cover between 15 and 20 m.

Two different geolocations were used for GEDI: the one provided by GEDI data version 2, and the one obtained after correcting positions available in version 2 using the GeoGEDI algorithm (Schleich et al., 2023). GeoGEDI corrects GEDI footprint positions using a digital elevation model (DEM). In this study, a publicly available 1-m lidar DEM obtained through the French mapping agency (IGN, 2022), was used.

### 2.2.3. Sentinel-2

The Copernicus Sentinel-2 (S2) satellites acquire optical multispectral imagery measuring the Earth's reflected radiance across 13 spectral bands ranging from the visible to the shortwave infrared (SWIR) (443 - 2190 nm) (ESA European Space Agency), a).

This study uses Sentinel-2 L2 A and Sentinel-2 Theia L3 A products. The surface reflectance L2 A product (Copernicus Sentinel-2, 2021),
downloaded from Google Earth Engine (Gorelick et al., 2017), offers atmospherically corrected data and is available at 110 × 110 km² ortho-image tiles. Our study site corresponds to tile ULU32. The Theia L3 A product is a monthly (45 days) cloud-free synthesis of the L2 A product using the Weighted Average Synthesis Processor (WASP) (Hagolle et al., 2018). For this study, we downloaded the L3 A product for August 2022 as the reference image. L3 A images from August 2017 and June 2022 were also downloaded to generate indices related to temporal changes.

For these three L3 A images, vegetation indices were calculated using the 10 m resolution bands in the visible (B2 Blue, B3 Green, B4 Red) and NIR (B8). The chosen indices have proven to be effective for vegetation monitoring (Rouse et al., 1973; Gao, 1996; Louhaichi et al., 2001; Hunt et al., 2011; Eng et al., 2019). They are all defined in the Index DataBase (Henrich et al., 2012).

$$ndvi = (B8 - B4)/(B8 + B4) \qquad (1)$$

$$msavi = 0.5 * (2 * B8 + 1 - \sqrt{(2 * B8 + 1)^2 - 8 * (B8 - B4)}) \qquad (2)$$

$$ndwi = (B3 - B8)/(B3 + B8) \qquad (3)$$

$$gli = (2 * B3 - B4 - B2)/(2 * B3 + B4 + B2) \qquad (4)$$

Time series indicators, potentially indicative of stand type and condition, were calculated by comparing August 2017 to the reference image for assessing annual change, and June 2022 to the reference image for seasonal change. Images from August 2017 and June 2022 were used exclusively for these respective purposes.

In addition, using the R package prosail (Feret and de Boissieu, 2023), we calculated through an hybrid inversion process the leaf area index (lai), the fraction of vegetation cover (fCover), i.e. the gap fraction for nadir direction, and the Fraction of Absorbed Photo-synthetically Active Radiation (fAPAR) (Weiss et al., 2016). To work with a radiative transfer model, monthly syntheses are not appropriate; instead, original Level-2 A data in its native format is required to run the existing code for calculating these biophysical variables. Therefore, we used a single cloud-free L2 A image acquired on August 12, 2022, i.e. within the time period used to compute the August 2022 L3 A product.

### 2.2.4. Sentinel-1

Sentinel-1 (S1), a synthetic aperture radar (SAR) satellite mission developed by ESA, uses C-band microwave pulses (around 4 to 8 GHz) directed towards the Earth's surface (ESA (European Space Agency), b).

Using Google Earth Engine (GEE), Level 1 ground range detected interferometric wide-swath (L1 GRD IW) products, with both VV and VH polarizations, acquired during the same period as the 45-day period of Theia L3 A of S2 for August 2022 were selected, resulting in a total of 37 orbits, split between 17 ascending and 20 descending orbits (e.g. S1 AIW_GRDH_1SDV_20220729T054249_20220729T054314_044311_0549CE_5665). Each point in the study site was observed from 7 to 13 available orbit images (descending and ascending combined). The L1 GRD IW product consists of projected to the ground range SAR data, at 10 m resolution. Backscatter coefficient sigma0 was used.

To mitigate speckle, mean values for each polarization (VV and VH) and orbit direction (ascending and descending) were retained, i.e. mean values of 3 to 7 images. The four resulting images were corrected from the influence of topography on the backscattered values using an angular-based radiometric slope correction algorithm from Vollrath et al. (2020a) ready to run with GEE and available on github (Vollrath et al., 2020b). We also calculated the average of the corrected ascending and descending mean images. This results in six images: VVasc, VHasc, VVdesc, VHdesc, VVascdesc, and VHascdesc.

Moreover, ratio and radar vegetation index (rvi) were calculated after normalizing the VV and VH values by minimum and maximum.

rvi is an alternative to ndvi with optical images (Kim and Zyl, 2009; Charbonneau et al., 2005; Nasirzadehdizaji et al., 2019).

$$Ratio = VV/VH \qquad (5)$$

$$rvi = (4 * VH)/(VV + VH) \qquad (6)$$

## 3. Methods

### 3.1. Data preparation

#### 3.1.1. NFI dataset

In contrast with Strategy A, which relies solely on Sentinel information for imputing GEDI values, Strategy B (see Section 3.2) takes advantage of the maximum heights to supplement the Sentinel data. In GEDI products, RHv_100 represents the maximum height of a footprint. The default NFI data provide dominant height and not maximum height that could be related to RHv_100. To address this disparity, we augmented the NFI dataset by imputing a maximum height to each NFI plot. This imputation was performed using MissForest on the forest database, as presented in Schleich et al. (2024). The subsequently added maximum height from the NFI data is denoted as Hmax. In addition to Hmax, two other variables from the NFI dataset will be used. The first variable, labeled vwac in the official French NFI database, corresponds to the GSV. The second variable, esspre_fr, defines whether the NFI plot is dominated by coniferous or deciduous trees. The NFI dataset is composed of 415 coniferous plots with a mean and standard deviation GSV of 360 ± 207 m³ ha⁻¹, and of 260 deciduous plots with a mean and standard deviation GSV of 227 ± 155 m.

#### 3.1.2. GEDI footprint georeferencing

GEDI footprints have a low geolocation accuracy (~ 10 m for version 2 (Beck et al., 2021)), which may cause problems when coupling GEDI data with other georeferenced datasets. Therefore, we applied the GeoGEDI algorithm outlined by Schleich et al. (2023) to improve the geolocation of GEDI footprints. As suggested by Schleich et al. (2023), we opted for the GeoGEDI algorithm with beam-pair implementation. To identify the optimal location, a search window of 30 m with a step size of 1 m was used around each footprint position.

Mean, median, and standard deviation of the distances between original GEDI v2 positions and corrected GeoGEDI positions were 9.20, 7.29, and 7.19 m, respectively.

#### 3.1.3. Auxiliary variables

The following variables were considered for modeling:

- Sentinel-2:

- L3 A data August 2022

- Bands: B2, B3, B4, B5, B6, B7, B8, B8 A, B11, B12

- Indices: ndvi, msavi, ndwi, gli

- Monthly differences of indices with June 2022

- Yearly differences of indices with August 2017

- L2 A data August 12th 2022

- Biophysical variables: lai, fCover, fAPAR

- Sentinel-1 Mean of August 2022

- VVasc, VHasc, VVdesc, VVasc, VVascdesc and VHascdesc

- Indices: Ratio, Rvi for asc, desc and ascdesc

- FORMS-H
<figure><img src="image_4.png" /><figcaption>Fig. 2. Step 1: GEDI variables are imputed over NFI plots using auxiliary data. Colors are only a figurative way of suggesting the diversity in auxiliary data, for instance the different values of S2 bands and indices.</figcaption></figure>

### 3.1.4. Data framework for modeling

GEDI and NFI data were intersected with Sentinel and FORMS-H data using spatial intersection, with X and Y coordinates in Lambert 93 as the coordinate reference system (EPSG 2154). To account for the different spatial resolutions of the data, zonal statistics were used to summarize Sentinel and FORMS-H data at both GEDI footprint and NFI plot centers. Zonal statistics consisted in the mean and standard deviations of all the variables described in Section 3.1.3. The statistics were computed using radii of 15 m and 50 m. The first radius was set to match the NFI plot radius, and the scale at which dendrometric measurements are made in the field. The second radius was set to consider the spatial context around the plot.

To account for the accuracy of GEDI positions, two versions of GEDI were used: one using initial GEDI v2 positions and another using GeoGEDI-corrected positions.

For clarity, variable names were formed using a prefix characterizing the data source (S1 for Sentinel-1 and S2 for Sentinel-2), followed by the variable name (e.g. B4 for S2 band 4), followed by a suffix made of the aggregation radius (i.e. 15 or 50) and the statistic (mean or sd) when appropriate (for example, S2_B4_15_mean is the mean of Sentinel-2 band 4 computed using a 15 m radius). For the monthly and yearly differences of indices, prefixes MD and YD were added, respectively. Notably, in one modeling strategy (see Section 3.2.1 hereafter), height information was obtained from both GEDI (RHv_100) and NFI (maximum height) as a bridging variable. The variable is denoted Hmax. Similarly, when GEDI data are used as auxiliary data (see Section 3.2.2), a similar naming convention is used with a prefix G.

### 3.2. kNN bagging approach

Modeling is based on random patch kNN bagging regression. It consists in an ensemble of kNN imputations, constructed from a random subset on both instances and features. While kNN has been widely used in MFIs owing to their capacity to predict multiple attributes with a single model, random patched kNN bagging was found to be an efficient approach toward competitive predictive performance with high dimensional data (Louppe and Geurts, 2012; Gomes et al., 2021). In this study, the kNN models are trained using an Euclidean distance, a neighborhood (k) of 1, and an inverse weighted distance. The algorithm is applied n times (n = 1000). For each ensemble, the random patch is formed using sampling with replacement for the instances (i.e. the lines) and sampling without replacement (set to 3 to 6) for the features (i.e. the columns) (Ho, 1998). From the resulting ensemble of predictions, an aggregated predictor is generated for each point (Breiman, 1996).

The modeling framework constitutes two steps described below. The first step involves spatializing GEDI data using Sentinel (Fig. 2). The second step involves spatializing NFI attributes using both Sentinel data and GEDI imputed values (Fig. 3). In each step, the predictions were evaluated using an independent test dataset.

#### 3.2.1. Step I: Imputing GEDI variables using sentinel data

For all strategies, a test dataset of 500 GEDI footprints (approximately 0.5% of the total dataset of (N = 104 831)) was randomly extracted for validation purposes. Although this represents a relatively small proportion, the fixed test sample size was chosen to balance computational efficiency with stable and reliable error estimation. The remaining GEDI footprints (≈99.5%) were used as the training dataset for the kNN approach.

Before modeling, dimension-reduction is conducted in the common variable space to create the feature vector, similar to the method mentioned by Sagar et al. (2022). Auxiliary variables (i.e. S1, S2, and possibly auxiliary height depending on the strategy) correlated by more than 0.33 to one of the GEDI variables and correlated by less than 0.85 between each other were selected. The latter criterion was aimed at limiting correlation between the selected variables, which otherwise may pose a risk of overfitting (Moser et al., 2017).

kNN bagging was performed under R using the package yaImpute (Crookston and Finley, 2007). From the ensemble predictions, for each predicted variable, the median value is calculated and regarded as the predicted value. The median was selected here rather than the mean due to the risk of outliers caused by issues in lidar signal analysis. Analyzing the test dataset comprising 500 GEDI footprints, theoretically, the imputed median values should align with the observed GEDI values. Nevertheless, a recognized phenomenon in regression analysis is the tendency for values to regress towards the mean, implying that for small values, imputations tend to overestimate the observed values, and for large values, imputations tend to underestimate the observed values. This introduces an additional source of correlated imputation errors, as suggested by Ehlers et al. (2018). To address this issue and enhance the established relationship, we applied a linear regression for each RHv variable on the 500 test footprints. In the context of GSV, this classical calibration procedure, reducing correlated errors, has also been used by Lindgren et al. (2021). First, using the median imputations and observed values, for each variable g, linear regression parameters $A_g$ and $B_g$ are estimated for the GEDI test dataset according to Eq. (7) :

$$ \hat{y}_{i,g} = A_g + B_g y_{i,g} + \varepsilon \quad (7) $$

where:

<table><tr><td>ŷ<sub>i,g</sub></td><td>median imputation of variable g</td></tr><tr><td>A<sub>g</sub></td><td>linear regression intercept of variable g</td></tr><tr><td>B<sub>g</sub></td><td>linear regression slope of variable g</td></tr><tr><td>y<sub>g</sub></td><td>observed value of variable g</td></tr><tr><td>ε</td><td>the error term supposed to be normally distributed.</td></tr></table>

Subsequently, using the linear regression parameters $A_g$, and $B_g$, the systematic errors are removed by classical calibration (Osborne, 1991) so that calibrated imputations ($ŷ_{i,g,c}$) are obtained as presented in Eq. (8):

$$ \hat{y}_{i,g,c} = \frac{(\hat{y}_{i,g} - A_g)}{B_g} \quad (8) $$

where:
<figure><img src="image_5.png" /><figcaption>Fig. 3. Step 2: predicting GSV based on auxiliary data and imputed GEDI data from step 1. Colors are only a figurative way of suggesting the variability in auxiliary data, for instance the different values of S2 bands and indices.</figcaption></figure>

<table><tr><td>ŷ<sub>l,g,c</sub></td><td>calibrated median imputations of variable g</td></tr><tr><td>ŷ<sub>l,g</sub></td><td>median imputation of variable g</td></tr><tr><td>A<sub>g</sub></td><td>linear regression intercept estimated in Eq. (7)</td></tr><tr><td>B<sub>g</sub></td><td>linear regression slope estimated in Eq. (7).</td></tr></table>

The optimal GEDI footprint is determined by selecting the footprint minimizing residuals of its RHv profile to the calibrated median profile of the 1000 candidates. The different RHv percentiles were scaled, to ensure equal weight.

Step I therefore allows for the matching of an existing GEDI footprint to each NFI plot. Therefore, all GEDI variables, and not only the variables used for the selection of the best imputation, are imputed.

Before moving to Step II, we filtered out poorly imputed GEDI variables. In the assessment of observed and imputed values on the test dataset, variables with a correlation below 0.6 were considered to be poorly imputed and subsequently excluded from consideration in step II.

### 3.2.2. Step II: Predicting GSV at the level of GEDI footprints

Similar to Step I, the enlarged common variables space is used to develop a model for predicting GSV from GEDI, S1, and S2 features, and from FORMS-H for Strategy C. The model is developed using the GEDI-Sentinel-NFI dataset (Fig. 3). To evaluate the method, a 20% subset of NFI plots was used as test dataset. The remaining GEDI-Sentinel-NFI plots serve as training dataset. Variable reduction is the same as in Step I. Among all the variables, i.e. well imputed GEDI variables from step I (i.e. correlation >= 0.6), S1 and S2 and, for Strategy C, FORMS-H metrics, variables correlated by more than 0.33 to NFI GSV and correlated by less than 0.85 between each other were selected. The kNN algorithm is run 1000 times. Similar to step I, the number of features to be used was randomized to retain 3 to 6 features, if there were more than 3. If there were fewer than 6 features, the number of features to use was randomized between 3 and the number of features, and if there were exactly 3 features, for each ensemble, the number of features to use was randomized to be either 2 or 3.

GSV estimates at NFI plot level are considered accurate. Given the improved accuracy (compared to GEDI signal analysis) and the considerably smaller NFI plot dataset, compared to the GEDI dataset, the mean value of the predictions of the 1000 kNNs was used to predict the GSV. This kNN bagging approach results in GSV estimates for all plots in the test data and in a model that can be further applied to predict GSV at the level of each GEDI footprint (Strategy B) or possibly at each cell of a 30 m grid, provided that all auxiliary data are available wall-to-wall (Strategies A and C).

### 3.3. Analysis of results

The different strategies were compared to assess if Strategy A (using only S1 and S2 data) allows for the improvement of GSV estimates, and how the other two strategies performed. For this comparison, we performed the following analysis in particular:

* The variables selected for the first step depending on the strategies and their ability to impute GEDI variables. Their ability to impute GEDI variables is assessed through the analysis of residuals and correlations between observed and imputed values for the test dataset.

* The variables chosen for the second step and their ability to predict GSV at the footprint level. To analyze the ability to predict GSV, using the test dataset, Mean Absolute Error (MAE), Root Mean Square Error (RMSE), MAE%, RMSE% as the ratio of the RMSE with the mean, R², Pearson's correlation coefficient and boxplots of errors are calculated between observed and predicted GSV values. The distribution of standard deviation errors of the three strategies will also be compared. The equations of these evaluation metrics are provided in Appendix A.

* We compared the GSV predictions generated by our three strategies with those from FORMS-V (Schwartz et al., 2023). The GSV map derived by Schwartz et al. (2023) was derived from their height map (FORMS-H), using two distinct allometric equations — one for deciduous and another for coniferous plots.

* To enhance the comparability of our strategies with FORMS-V, given that the map was created without georeferencing improvement and relied on two allometric equations, we implemented the strategy that performed best using GEDI v2 coordinates and developed separate models for deciduous and coniferous plots. For both v2 and improved footprint coordinates, we will present results using the overall kNN including all data, as well as results from the same overall kNN model categorized by forest stand types. Furthermore, we will present outcomes from the individually conducted coniferous kNN and deciduous kNN. Results will be presented according to stand types and aggregated for direct comparison with the overall kNN regarding prediction accuracy on the entire test dataset.

## 4. Results

This Section is divided into three main subsections. First, results from Step I and II are presented in Sections 4.1 and 4.2. In addition, alternative setups using stand-specific kNN and offering insights on results without geolocation improvement are presented in Section 4.3

### 4.1. Step I: Imputing GEDI variables

#### 4.1.1. Auxiliary variables selected for kNN in step I

The size of the feature space used to impute GEDI variables has been significantly reduced using the variable reduction strategy based on correlations 3.2. In Table 2 the selected variables are listed in decreasing order of their correlation with GEDI variables.

In Strategy B, Hmax (i.e. RHv_100) is used as the auxiliary height. Consequently, the maximum correlation between the auxiliary height and the GEDI variables is 1. When using FORMS-H data, i.e. in Strategy C, the strongest correlation is observed with RHv_90, which is 0.71. All other variables are consistent across strategies. Only six variables were retained among the more than 150 initial S1 and S2 variables. Evidently, no S1 variable was selected, all of them exhibited a weak correlation with GEDI variables (−0.33 < r < + 0.33), the strongest correlation being equal to −0.22 for S1_vvAscDesc_15_sd.

The S2 variable most correlated to GEDI data is S2.greenness_15_mean, followed by non-correlated variables S2_fAPAR_15_mean, S2_ndwi_15_mean, S2_YD_ndvi_15_sd,
<caption>**Table 2**</caption>

Variables selected for step I. The first column contains the selected variables, the second column contains the GEDI variable to which the auxiliary variable is the most correlated, and the third column indicates the correlation between the two. Variables start with “S1” or “S2” if they originate from Sentinel-1 or Sentinel-2 data, followed by the variable name, followed by 15 or 50 depending if the 15 m or 50 m radius was used for the extraction of the variable. Variables end with “.mean” or “.sd” depending on the zonal extraction of mean or standard deviation. Auxiliary height for Strategy C is FORMS-H_15_mean. Auxiliary height for Strategy B is Hmax from GEDI and NFI.

<table><thead><tr><th>Auxiliary variable</th><th>Most correlated GEDI variable</th><th>Corr</th></tr></thead><tbody><tr><td>Hmax for Strategy B</td><td>RHv_100</td><td>1.00</td></tr><tr><td>FORMS-H_15_mean Strategy C</td><td>RHv_90</td><td>0.71</td></tr><tr><td>S2.greenness_15_mean</td><td>cover_zf</td><td>0.40</td></tr><tr><td>S2.fAPAR_15_mean</td><td>pavd_z3</td><td>0.40</td></tr><tr><td>S2.ndwi_15_mean</td><td>cover_zg</td><td>-0.35</td></tr><tr><td>S2.YD_ndvi_15_sd</td><td>cover_z2</td><td>-0.35</td></tr><tr><td>S2.ndvi_50_mean</td><td>cover_z2</td><td>0.34</td></tr><tr><td>S2.B4_15_sd</td><td>cover_z2</td><td>-0.34</td></tr></tbody></table>

S2_ndvi_50_mean and S2_B4_15_sd. Not surprisingly, S2 variables are more strongly correlated with GEDI cover variables than to the height ones. The most correlated GEDI variables to S2 data are cover_zf, pavd_z3, cover_zg and cover_z2.

Notably, correlated variables were excluded (abs(Cor) > 0.85). Therefore, while other variables (lai_mean, fCover_mean, msavi_mean and ndvi_sd) had a correlation exceeding 0.33 with GEDI variables, they were not selected as their correlation to a previously chosen variable exceeded 0.85.

### 4.1.2. Evaluation of the kNN regression model used to impute GEDI variables

Fig. 4 shows results obtained on the 500 GEDI test footprints dataset, comparing imputed to observed values for G_RHv_100.

For Strategy A, we observe that median imputations of G_RHv_100 are restricted to heights ranging from 15 to 36 m, whereas the observed G_RHv values range from 9 to 52 m. When a calibration is applied, the range is stretched (Fig. 4(d)) included with negative predicted values, while values should ideally follow a 1:1 trend. The best imputation then chooses the closest existing GEDI footprint, thereby minimizing the difference to the calibrated median value based on all RHv deciles. For Strategy A, the best imputations for G_RHv_100 have a correlation of only 0.10 to observed G_RHv_100 values. Residual dispersion is vast and increased from G_RHv_10 to G_RHv_100, as seen in Fig. 4(j). This increase subsequently translates to an increase in both the mean and range of observed heights from G_RHv_10 to G_RHv_100 combined, with a very poor performance of imputations over this range, regardless of the observed value. From G_RHv_40 and onwards, imputations underestimate observed heights.

For Strategy B, G_RHv_100 imputations are highly correlated with the observed values (r = 0.99). Median height imputations are however slightly overestimated for high values and underestimated for low values. The calibration allows for the correction of this slight bias. For this strategy, unlike Strategies A and C, residual distributions tend to narrow from G_RHv_10 to G_RHv_100. This pattern underlines the more optimized predictions for higher deciles compared to lower ones, highly influenced by the use of RHv_100 as an auxiliary variable.

As observed in the third column of Fig. 4, imputations with Strategy C outperform that of Strategy A. The correlation between best imputations and observed G_RHv_100 values is 0.592, which is similar to the correlation based on median imputations (0.594). Globally, the residuals of the imputed values are considerably lower than those for Strategy A. The overall median difference is low (<1 m) for all G_RHvs (Fig. 4(k)). Residual distributions only slightly narrowed from G_RHv_100 to G_RHv_10, suggesting an influence of the observed RHv decile distribution on the imputed value distributions.

All GEDI variables were imputed in the GEDI-Sentinel dataset, not only G_RHv_100 and other height deciles. However, we restricted further processing to the ones with a correlation coefficient with observed

values that are greater than 0.6 (Table 3) to avoid using wrongly predicted GEDI variables in the subsequent step.

In summary, Strategy A did not allow for the retention of any GEDI variable for step II, while Strategy B allowed for the retention of 18 and Strategy C the retention of 4 GEDI variables.

## 4.2. Step II: Predicting GSV

### 4.2.1. GEDI-sentinel variables for step II

The imputed GEDI variables kept at step I, if any, as well as those of Sentinel (and FORMS-H for Strategy C) were further used to predict GSV. To reduce the number of variables used in the kNN for GSV predictions, and similarly to step I, a variable reduction step based on their correlation with GSV was introduced. The retained variables are displayed in Table 4.

Common to all strategies, three S2 features computed with a 15 m buffer were selected: two bands (B6 and B11, i.e. vegetation red edge and SWIR, respectively) and fCover. Similar to step I, variables that were highly correlated with the retained ones were discarded. As no GEDI variables were “correctly” imputed in step I for Strategy A, only S2 data were used in this case.

For Strategy B, 5 among the 18 correctly imputed GEDI variables were chosen. The most correlated variable was G_RHv_60 (r = 0.69), followed closely by other G_RHv variables. However, the latter were discarded based on their correlation with G_RHv_60. Then G_pai_z5, GCover_z3, G_fhd_normal, and GCover_z5z6 were also selected.

For Strategy C, the same variables as in Strategy A and the auxiliary height FORMS-H were selected. Notably, FORMS-H height was used because of its correlation (r = 0.59) with GSV. Imputed G_RHv_60 was right behind (r = 0.58), but not selected because of its strong correlation with FORMS-H height.

### 4.2.2. Evaluation of the kNN regression model used to predict GSV

Results from kNN bagging using previously selected variables 4.2.1, and applied on the NFI-Sentinel-GEDI dataset, are illustrated for the three strategies in Figs. 5 and B.8(a). The test dataset was also intersected with the volume map FORMS-V. Therefore, we also compared our results to an existing GSV map.

Comparing the observed with the predicted GSV, a gradual improvement is observed from Strategy A to C to B. All strategies have difficulty predicting GSV above ~600 m³ ha⁻¹. Predictions tend to saturate at high GSV values.

In terms of errors, Strategy B has a small MAE of 101.23 m compared to 160.73 m and 126.51 m for Strategies A and C, respectively. The RMSE% improves by 23.23 percentage points from Strategy A to B, and by 10.13 percentage points from Strategy C to B. R² and correlation coefficients also improve from A to C to B. The median value, as observed in Fig. B.8(a), is closer to 0 for Strategy B, and the interquartile range also reduces from A to C to B.

When comparing results with estimates from FORMS-V, the MAE for FORMS-V is 125.38 m. Interestingly, results obtained with our Strategy C yield similar outcomes to estimates from FORMS-V.

All results presented in Sections 4.1 and 4.2 refer to data used with corrected GEDI positions.

## 4.3. Comparing several setups for the optimal strategy

Since Strategy B outperformed the other two strategies in Section 4.2, only its results are presented for these more specific tests. This Section presents results obtained by implementing Strategy B while using stand-type-specific models, with and without GeoGEDI correction. For both GEDI footprint location versions (GEDI v2 and GEDI v2 corrected by GeoGEDI), an overarching model (overall kNN) and individual models for each stand type (Coniferous/Deciduous kNN) emerged. For each setup, error metrics are computed using NFI plots of the test dataset (N = 135), both aggregated and separately for
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** G_RHv_100 imputation vs observed data for the 500 test GEDI footprints. From the 1000 iterations of the bagging process, the median imputation is first selected (first line), then the median predictions are compared to observed values to calibrate the median imputations to account for the regression toward the mean trend (second line). The best imputation is further defined by selecting the real existing GEDI footprint with variables minimizing the distance to the calibrated median imputations of 10 Rhv decile variables (third line). The last line shows the distribution of the residuals between the best imputations and the real Rhv values. *Correlation*, *R*², *RMSE* and *RMSE%* are reported in the upper left corner of each plot.</figcaption></figure>
<table><caption>**Table 3**
GEDI variables with a correlation coefficient greater than 0.6.</caption><thead><tr><th>Strategy A</th><th>Strategy B</th><th>Strategy C</th></tr></thead><tbody><tr><td>/</td><td><ul><li>all the G_RHv deciles except G_RHv_10</li><li>G_cover_z3, G_cover_z4, G_cover_z5</li><li>G_fhd_normal, G_pai_z4, G_pai_z5, G_pavd_z6, G_cover_z5z6, G_pai_z5z6</li></ul></td><td><ul><li>G_RHv_90, G_RHv_80, G_RHv_70, G_RHv_60</li></ul></td></tr></tbody></table>

<table><caption>**Table 4**
Selected variables for step II. Variables were chosen based on their correlation with GSV.</caption><thead><tr><th colspan="2">Strategy A</th><th colspan="2">Strategy B</th><th colspan="2">Strategy C</th></tr><tr><th>Variable</th><th>Corr</th><th>Variable</th><th>Corr</th><th>Variable</th><th>Corr</th></tr></thead><tbody><tr><td>S2_B6_15_mean</td><td>-0.42</td><td>G_RHv_60</td><td>0.69</td><td>FORMS-H_15_mean</td><td>0.59</td></tr><tr><td>S2_fCover_15_mean</td><td>-0.37</td><td>G_pai_z5</td><td>0.57</td><td>S2_B6_15_mean</td><td>-0.42</td></tr><tr><td>S2_B11_15_mean</td><td>-0.34</td><td>G_cover_z3</td><td>0.54</td><td>S2_fCover_15_mean</td><td>-0.37</td></tr><tr><td></td><td></td><td>G_fhd_normal</td><td>0.53</td><td>S2_B11_15_mean</td><td>-0.34</td></tr><tr><td></td><td></td><td>G_cover_z5z6</td><td>0.44</td><td></td><td></td></tr><tr><td></td><td></td><td>S2_B6_15_mean</td><td>-0.42</td><td></td><td></td></tr><tr><td></td><td></td><td>S2_fCover_15_mean</td><td>-0.37</td><td></td><td></td></tr><tr><td></td><td></td><td>S2_B11_15_mean</td><td>-0.34</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** GSV predictions diagnosis based on test NFI plots (N = 135 plots) according to the estimation strategy. Correlation, R², RMSE, RMSE%, MAE and MAE% are reported on each plot. Errors were calculated as observed GSV - predicted GSV.</figcaption></figure>

each stand type. For the test dataset the mean GSV of coniferous NFI plots was 356.15 m³ ha⁻¹ and the mean GSV of deciduous NFI plots was 234.76 m³ ha⁻¹.

Results of step I and step II without geolocation correction are presented in Appendix B. For step I, the Sentinel features, extracted at uncorrected GEDI positions, are less correlated with GEDI variables than the ones extracted at corrected GEDI positions, i.e. the greenness indicator correlation coefficient is 0.40 for corrected positions, and 0.37 for uncorrected positions. Consequently, applying the same correlation thresholds compared to the corrected positions, lead to a selection of less Sentinel features for the kNN, i.e. 3 instead of 6.

For step II, with and without GeoGEDI correction, no imputed GEDI variables were retained to predict GSV for Strategies A and C, therefore only Sentinel variables (and FORMS-H) were used. For the test data

set these variables were extracted at NFI plot positions, and therefore results do not change between GeoGEDI corrected and uncorrected test NFI datasets. For Strategy B, results on test NFI datasets change, as GEDI variables have been imputed differently. The most correlated variable to GSV for GeoGEDI corrected footprints was G_RHv_60 (r = 0.69) whereas for uncorrected datasets, the highest correlation was for G_RHv_100 (r = 0.68). Five GEDI variables were retained with GeoGEDI correction and four without GeoGEDI correction. The G_fhd_normal had the same correlation without and with GeoGEDI correction (r = 0.54 and 0.53), while the G_cover_z5z6 variable was better predicted with the GeoGEDI correction (r = 0.44 vs 0.37). Final results showed that GeoGEDI corrected datasets performed better to predict GSV than uncorrected data (Table 5 and Fig. B.7). RMSE% with GeoGEDI was 42.60%, and without GeoGEDI it was 47.75%.
<caption>**Table 5**</caption>

Strategy B with and without GeoGEDI correction, for combined (All), coniferous and deciduous test datasets, using overall, and stand-specific models. The “With overall kNN” corresponds to the method described in the methodology. All NFI plots were used to run the kNN, and at the end we split NFI predictions by their dominant stand type. The “With Conif/Decid” kNN columns correspond to two kNNs run distinctly for step II, and aggregating them for “All”. The last column is for comparison with estimations from FORMS-V.

<table><thead><tr><th rowspan="2">Error</th><th colspan="3">With overall kNN</th><th colspan="3">With Conif/Decid kNN</th><th colspan="3">FORMS-V</th></tr><tr><th>All</th><th>Conif</th><th>Decid</th><th>All</th><th>Conif</th><th>Decid</th><th>All</th><th>Conif</th><th>Decid</th></tr></thead><tbody><tr><td colspan="10"><strong>with GeoGEDI correction</strong></td></tr><tr><td>MAE</td><td>101.23</td><td>106.27</td><td>91.82</td><td>109.58</td><td>119.36</td><td>91.30</td><td>125.38</td><td>127.82</td><td>120.82</td></tr><tr><td>MAE%</td><td>32.26</td><td>29.84</td><td>39.11</td><td>34.92</td><td>33.51</td><td>38.89</td><td>39.95</td><td>35.89</td><td>51.47</td></tr><tr><td>RMSE</td><td>133.70</td><td>144.18</td><td>111.52</td><td>142.86</td><td>155.25</td><td>116.21</td><td>166.07</td><td>172.44</td><td>153.44</td></tr><tr><td>RMSE%</td><td>42.60</td><td>40.48</td><td>47.50</td><td>45.52</td><td>43.59</td><td>49.50</td><td>52.91</td><td>48.42</td><td>65.36</td></tr><tr><td>R<sup>2</sup></td><td>0.58</td><td>0.60</td><td>0.34</td><td>0.52</td><td>0.51</td><td>0.30</td><td>0.38</td><td>0.42</td><td>0.16</td></tr><tr><td>Correlation</td><td>0.77</td><td>0.78</td><td>0.60</td><td>0.72</td><td>0.72</td><td>0.56</td><td>0.62</td><td>0.65</td><td>0.42</td></tr><tr><td colspan="10"><strong>without GeoGEDI correction</strong></td></tr><tr><td>MAE</td><td>114.52</td><td>130.10</td><td>85.41</td><td>110.67</td><td>129.34</td><td>75.79</td><td>125.38</td><td>127.82</td><td>120.82</td></tr><tr><td>MAE%</td><td>36.49</td><td>36.53</td><td>36.38</td><td>35.26</td><td>36.31</td><td>32.28</td><td>39.95</td><td>35.89</td><td>51.47</td></tr><tr><td>RMSE</td><td>149.87</td><td>168.94</td><td>105.31</td><td>149.72</td><td>170.33</td><td>100.44</td><td>166.07</td><td>172.44</td><td>153.44</td></tr><tr><td>RMSE%</td><td>47.75</td><td>47.44</td><td>44.86</td><td>47.70</td><td>47.82</td><td>42.78</td><td>52.91</td><td>48.42</td><td>65.36</td></tr><tr><td>R<sup>2</sup></td><td>0.47</td><td>0.44</td><td>0.41</td><td>0.47</td><td>0.41</td><td>0.46</td><td>0.38</td><td>0.42</td><td>0.16</td></tr><tr><td>Correlation</td><td>0.69</td><td>0.67</td><td>0.65</td><td>0.69</td><td>0.65</td><td>0.69</td><td>0.62</td><td>0.65</td><td>0.42</td></tr></tbody></table>

With Strategy B of our approach, improving georeferencing had a notable impact on predictions, particularly for coniferous trees. For this stand type, we observed a substantial increase in correlation of 0.11 and 0.07 and a decrease in RMSE% of 6.96 and 4.23, with overall and coniferous-specific kNN, respectively. However, this georeferencing improvement led to a slight degradation in predictions for deciduous trees, with a reduction in correlation of 0.05 to 0.13 and an increase in RMSE% from 2.64 to 6.72, with overall and deciduous-specific kNN, respectively. Surprisingly, using a single model (i.e. overall kNN) amplified the positive impact of georeferencing improvement, enhancing prediction accuracy for coniferous trees while mitigating the degradation of predictions for deciduous trees. Consequently, overall results indicate more significant improvement with a single model compared to a model for each stand type, showing an increase in correlation of 0.08 instead of 0.03 and a decrease in RMSE% of 5.15 percentage points instead of 2.18 between non corrected and GeoGEDI corrected datasets for overall kNN and stand-specific kNN, respectively.

Compared to FORMS-V, our overall-kNN model without GeoGEDI correction exhibited similar performance for coniferous trees but demonstrated significantly improved results for deciduous trees (Fig. B.8). Similar trends were observed when comparing our stand-specific kNN model, without improved georeferencing, to FORMS-V. Notably, there was an even more substantial improvement for deciduous trees without GeoGEDI correction than seen with the overall kNN model. Consequently, our two models without improved georeferencing slightly superior predictions compared to the FORMS-V approach, showing a correlation increase of +0.07 and an RMSE% reduction of -5.2 for both models. In contrast, the overall-kNN model with improved georeferencing significantly enhanced predictions, with a correlation increase of +0.15 and an RMSE% reduction of -10.31.

# 5. Discussion

## 5.1. GEDI to NFI matching strategy

An effective matching strategy between GEDI structural information and NFI plots is essential for building robust calibration models that fully exploit both GEDI and NFI attributes. In this study, the matching was based solely on auxiliary data, which introduces the risk of associating mismatched pairs. This risk is particularly relevant over large spatial extents, where stands growing under different ecological conditions may display similar auxiliary data values, potentially leading to inaccurate linkages. Introducing spatial distance dependencies could help address this issue, either by setting a maximum search distance (Zhang et al., 2022; Bullock et al., 2023) or by fixing the number of neighboring footprints to consider (May et al., 2024). A key

limitation of using a fixed search distance is that it may significantly reduce the number of available matches, especially when the average minimum distance between NFI plots and high-quality GEDI footprints (i.e., those retained after filtering) exceeds the specified threshold. This can be mitigated by instead fixing the number of neighboring footprints. For instance, May et al. (2024) selected the 25 nearest GEDI footprints for each NFI plot and developed a geostatistical approach to generate probabilistic GEDI profiles at NFI plot locations, enabling robust model training for forest attribute prediction. However, a potential limitation of using fixed neighborhood — whether defined by spatial distance or footprint count — is the risk of including footprints located beyond the spatial autocorrelation range. This concern is especially relevant in structurally heterogeneous forests (Hunka et al., 2023). For example, in the study by Zhang et al. (2022) in France, spatial correlation was lost at approximately 500 m.

An alternative approach may involve balancing geographic distance with distance in the auxiliary data space, as proposed in Besic et al. (2024). In this study, the optimal weight between the two distances was found to be approximately equal (50/50). This strategy maximized the correlation between GEDI attributes and auxiliary variables, while maintaining spatial proximity between GEDI footprints and NFI plots, along with likely similar environmental conditions.

In Pascual et al. (2025), the authors compared three matching strategies: one based on GEDI-like waveforms simulated from ALS data; another using spatial nearest-neighbor matching — which led to a reduction in NFI-GEDI footprint pairs by over 50% with a 100 m distance threshold, and to just one-eighth with a 25 m threshold; and a third implementing the imputation-based method proposed by May et al. (2024). The best performance in biomass model calibration, linking field-plot biomass to GEDI attributes, was achieved using ALS-simulated GEDI waveforms, with normalized RMSE values approximately 10% lower than those of the other two strategies and no observed bias. Nonetheless, the two other strategies also yielded satisfactory results. The nearest-neighbor approach, despite the reduced number of training pairs, showed no bias and achieved NRMSE values of 82% and 84% for enhanced and standard NFI plot geolocation accuracy, respectively. Additionally, the findings of Besic et al. (2024) — in the same study area as ours — indicate that relying solely on auxiliary variable space for matching leads to suboptimal performance, though it still outperforms strategies based solely on geographic distance. The potential of our approach is likely comparable to the matching methods proposed by Besic et al. (2024) and May et al. (2024), but with a reduced computational cost.

Furthermore, we observed that improving the georeferencing accuracy of GEDI footprints led to enhanced predictions. Schwartz et al. (2022) also reported notable improvements of their height predictions
when transitioning from GEDI version 1 to GEDI version 2, largely due to reduced georeferencing bias. A similar result was found by Pascual et al. (2025) regarding the enhanced geolocation of NFI plots. Georeferencing issues are a well-known challenge in the field of forest attribute modeling with remote sensing data (Bouvier et al., 2019; Mcroberts et al., 2018; Saarela et al., 2016), and various methods have been proposed to address this issue, often relying on high-resolution ALS data (Fadili et al., 2019). By ensuring more accurate geolocation of both the GEDI and NFI data sources, the integration with other auxiliary data sources is expected to be more robust, leading to more reliable and coherent predictions. However, achieving such improvements is either costly (e.g. enhancing NFI plot geolocation) or reliant on auxiliary datasets that may not always be available, thus making them challenging to generalize.

## 5.2. Advantages and limitations of sentinel and GEDI

Results from Strategy A indicate that Sentinel data alone is insufficient for accurately imputing GEDI variables on our study site. We also observed that, whatever the strategy, features coming from S1 data were never selected. This result contrasts to various studies. In a similar context of mountainous forest in China, Guo et al. (2023) reported a significant contribution of S1 metrics using Recursive Feature Elimination - Support Vector Machine variable selection. While modeling canopy height using IceSat-2, Li et al. (2020) also indicated that backscattering coefficients from S1 could positively contribute to the prediction. The results of Kacic et al. (2021) over a tropical forest in Paraguay are more balanced, limiting the interest of S1 data to the detection of permanent water bodies and high biomass level. Shendryk (2022) in Australian forests reported that S1 performed less than S2, but could contribute to improving results when combined. This result was confirmed by Ge et al. (2022), using S1 times-series and S2 data to map forest heights. Considering the importance on S1 data in the literature cited, the reason for the non-selection in our study may be explained by the selection approach and the overall low correlation of S1 data with the attribute of interest, that is more difficult to model than height.

GSV predictions in Strategy B performed better than FORMS-V. A parallel study by Sagar et al. (2022) on a similar study site in the Vosges region using S2 and ALS data reported an RMSE% of 41.7%, comparable to our RMSE% of 42.60% with Strategy B. The observed saturation effect of S2 data at approximately 600 m³ ha⁻¹ appeared to be reduced compared to the other tested strategies, but persisted.

Incorporating FORMS-H largely improved imputations (Strategy C) compared to Strategy A, although a regression toward the mean persisted. This was partly mitigated by a calibration process (Lindgren et al., 2021). However, FORMS-H did not efficiently replace GEDI height values and its use led to reduced accuracy in predictions. Notably, FORMS-H has been produced from Sentinel 1 and 2 data using a U-Net deep learning model trained with GEDI RH95 data to map canopy height. Results showed a trend of underestimating height above 25 m (Schwartz et al., 2023). Therefore, the choice of the auxiliary data used to obtain information on forest structure, thereby complementing Sentinel information, appears to be decisive. In our approach, a height product can be used in addition to Sentinel data as a bridge variable, but this product has to be consistent enough with observed forest height and avoid saturation issues to fully play its role.

The comparison between deciduous and coniferous stands yielded more optimized predictions for coniferous data. Models had more difficulty in predicting GSV for deciduous plots. This trend can also be observed with FORMS-V (Schwartz et al., 2023). Schwartz et al. (2023) used allometric equations based on NFI data to create a GSV map from their height map, having one equation for deciduous and another for coniferous forest stands. Over the French territory and at the sylvoecoregional scale, the authors reported an R² of 0.63 with an MAE of 30 m³ ha⁻¹. However, aggregating results at the sylvoecoregion

level masked the high local heterogeneity in plot level prediction accuracy. For our study site, i.e. within an area covering most of the “Central Vosges Massif” sylvoecoregion, FORMS-V showed an R² of 0.38 with an MAE of 125.38 m³ ha⁻¹ compared to NFI test plots. With Strategies B and C, results were significantly improved compared to FORMS-V. Our study site had challenging conditions, mountainous and mixed stand. Other study sites, might have performed better. Previous research by Pereira-Pires et al. (2021) assessed the use of different S2 features to estimate GEDI RH100 and demonstrated strongly varying correlations between S2 and GEDI RH100, depending on the vegetation type.

Interestingly, Strategy C did not select the highest RHv and the overall structure variables, favoring variables reflecting middle layers of the canopy. These layers have more vegetation density than the top layers, and therefore might be better in relation with GSV.

Despite promising results for Strategy B, the method exhibited some instability. Stand-specific kNNs were expected to better fit to their respective datasets, but they did not improve estimates compared to the overall kNN method. When combining results from the two different kNNs, the overall correlation with observed GSV values was 0.72, compared to 0.77 with the overall kNN including both forest stand types. Moreover, predictions are resulting from the aggregation of 1000 kNN runs, and as illustrated in Appendix C, individually they differ considerably. The mean standard deviation within the 1000 kNN imputations for Strategies A and C is 144.65 m, compared to 126.62 m for Strategy B, as shown in Fig. C.9. Fig. C.10 provides examples of the 1000 predictions of GSVs by NFI plot. It shows a considerable variation in prediction quality and modes. While few were adequately predicted by a numerous kNN runs (C.10(a),C.10(b),C.10(c)), some had two modes (C.10(i),C.10(j)), some had very variable predictions resulting in mean estimates (C.10(l),C.10(k),C.10(e)), and others were simply inaccurately predicted (C.10(f),C.10(d),C.10(g), Fig. C.10(h)).

## 5.3. Consequences for MFI estimations

Strategy B allows for the creation of a densified GSV point dataset, with GSV estimated at each GEDI footprint. As Strategy B outperformed the other strategies, the results suggest that small area estimation (SAE) approaches, e.g. at the municipality or forest stand level, should be preferred to high resolution maps with pixel-level GSV predictions. SAE creates local predictions using model-assisted approaches, allowing for the correction of bias and improvement of variance. Alternatively, spatialized wall-to-wall maps obtained through model-based MFI approaches present a risk of important local errors, at least requiring accurate quantification of uncertainty at pixel level. Ståhl et al. (2024) also suggests that design-based inference might be a more conceptually appropriate framework, less exposed to biases.

To achieve accurate predictions of GSV, both Sentinel-2 optical data and high-quality height information are necessary. The quality of height data is a critical factor influencing the overall performance of GSV predictions. Improving the height map used (here, FORMS-H) would allow for direct application of the kNN bagging approach to create a wall-to-wall map. Advanced deep-learning approaches, such as vision transformers (Aleissaee et al., 2023), should be able to learn even better spatial features. Such features have proved useful to mitigating the saturation effect (Lang et al., 2023) to improve wall-to-wall height products. For example, Ge et al. (2022) used an improved semi-supervised deep learning approach to map forest height with S1 and S2 data and the RMSE% of tree heights was estimated to be 24.1% at pixel-level. The integration of multivariate techniques also holds potential for significant improvements.

However, as shown in studies with stands exceeding 20–30 m in height (e.g. Potapov et al. (2021), Lang et al. (2023), Schwartz et al. (2023), Sothe et al. (2022)) and, Gupta and Sharma (2022) revealed that lidar-based information can minimize the saturation issue from
optical and SAR C sensors, but cannot solve it entirely. This satura-
tion is indeed inherent to the signal itself. Fully spatializing forest
structure features, such as dominant or maximum height, with only
high-resolution S2 and S1 data even calibrated using dense datasets
such as GEDI data continues to be extremely challenging and might
not be the optimal solution for obtaining a bridge auxiliary structure
variable that is accurate enough.

With the data sources we used and the proposed double kNN
bagging approach, mapping GSV at the pixel level with its associated
precision could be counterproductive. Based on standard deviations
of plot level GSV predictions (see Appendix C), 95% confidence in-
tervals could be assessed for each pixel. These confidence intervals
would range between 4 and 22 m³ ha⁻¹, with a mean around 9 m³ ha⁻¹
with Strategy C. However, despite these reasonable precisions, model
predictions remain highly inaccurate and local accuracy cannot be as-
sessed. Recently, Ståhl et al. (2024) discussed emerging concerns about
model-based wall-to-wall maps and recommended that great caution be
exercised in the use of such products. They demonstrated the existence
of systematic over- or under-prediction that might call into question
the usefulness of the remote-sensing based predictions. Indeed, the use
of maps with uncontrolled errors or results from scenario modeling
based on such maps may lead to erroneous policy decisions. The risk
of negative consequences of the use of remote-sensing based maps is
higher when the goodness-of-fit of estimated models is intermediate or
poor (Ståhl et al., 2024), which is the case for our GSV models. This
is why, even though Strategy C is likely to provide a map with more
optimized predictions than those available in FORMS-V, we will not
retain this strategy.

Moreover, in our kNN bagging approach, variables representing the
middle layers of canopy were chosen over top canopy features. While
top of canopy features may capture the maximum height, which is
correlated with GSV, the GSV is also influenced by vegetation density
present in the middle layers of the canopy. This choice emphasizes the
significance of the middle layers in GSV estimation, as they include
crucial information related to vegetation density. This highlights an
advantage of GEDI data over the utilization of simple top of canopy
height maps derived from photogrammetry.

Last but not least, a key advantage of the proposed approach lies
in its potential to be readily extended to other NFI stand attributes.
KNN methods are robust, well-established methods that enable the
prediction of multiple plot-level attributes, albeit with varying levels
of accuracy (Tomppo et al., 2008; Chirici et al., 2016). In contrast, ad-
vanced machine learning approaches show a more limited application;
it is worth noting that among the studies reviewed by Hamedianfar
et al. (2022), multi-output models were employed in only one out of
ten cases for continuous forest variable estimation, and in only one out
of eleven cases for segmentation and classification tasks.

In this study, we focused on GSV as an illustration; however, the
method could also simultaneously estimate other forest variables, with
better results expected for those related to structure. This contrasts
with most recent studies using GEDI data, where biomass is central
due to its significance in carbon accounting and climate change-related
policies (Hunka et al., 2025; May et al., 2024; Pascual et al., 2025). In
our study, however, we did not focus on biomass for two main reasons.
First, standardized biomass data at the NFI plot level is currently
unavailable, as the agency responsible for NFI is in the process of
transitioning from field measurements to biomass estimation. Second,
our focus was on addressing resource management at the regional or
sub-regional level, where GSV is a key variable. While this makes direct
comparison challenging (our RMSE = 42.6% using overall kNN with
GeoGEDI correction compared for example to 71% or higher at the plot
level in Pascual et al. (2025)), the findings of Ortiz-Reyes et al. (2024)
suggest that some loss of accuracy in AGB estimates is to be expected.
However, their reported RMSE values (24.6% for volume and 38% for
biomass) are context-specific and may not be directly applicable to our
study.

Another advantage of KNN methods is that no assumption on the
relationship between forest and lidar variables is to be assumed. This
contrasts with linear regression models, that have been widely used
in the field of multi-source forest inventory to link forest attributes to
lidar-based variables and are still in use in recent studies (e.g. Pascual
et al. (2025), May et al. (2024)). Several studies have highlighted the
superiority of nonlinear models over linear ones in effectively capturing
the spatial complexity and inherent heterogeneity of forest ecosystems
(e.g., Moser et al. (2017), Lu and Jiang (2024)).

5.4. Possible improvements in terms of modeling and auxiliary data sources

For step I, where the goal is to impute multiple GEDI variables,
owing to the large dimension of the data, variable selection is required.
The efficiency of the selection is constrained by the large number of
predictors in both X and Y dimensions and the relatively low cor-
relation of Sentinel data with the vertical information contained in
GEDI. While for this study, simple selection using correlation was
found to perform better than more complex methods, additional work
would be required to optimize variable selection and associated model
performance with such data. Toward that goal, various approaches
have been tested in step II, where the goal is to predict a single variable
(GSV). These included stepwise regression, partial least squares, or
lasso regressions. However, these methods resulted in the selection of
numerous variables, which is not adequate for kNN models, as a large
number of predictors tends to homogenize distances. Consequently,
simple filtering based on correlation was retained.

This approach is more of a first step of variable reduction rather
than a final variable selection. It simply relies on correlations. The
thresholds of 0.85 and 0.33 are questionable, and changing these
thresholds affects final results. Future improvements may involve the
use of this correlation method with less restrictive thresholds as a
preliminary variable reduction step, followed by a more advanced
variable selection.

The correlation-based method may not capture the potential bene-
fits of less correlated variables. For example, Strategy C's retention of
FORMS-H could be revisited, with potential exploration of an imputed
RHv variable for improved results.

With Strategy C, adding a supplementary common variable of aux-
iliary height has opened up the possibility to predict GSV at any
point of the area of interest. However, the addition of a height map
raises concerns about its dependence to Sentinel data and its potential
reliability. Moreover, data could be misgeoreferenced, and errors might
be transferred to the GSV estimations. Furthermore, FORMS-H shows
important errors in mountainous areas (Schwartz et al., 2023), which
is the case for our study site. Additionally, Schwartz et al. (2023)
attempted to map RH95, considered as a measure of dominant height,
introducing a slight departure from Strategy B, where the maximum
height is considered. Improving existing height maps could enhance the
predictions. However, even with a more optimized height map, results
are not expected to be better than those obtained with Strategy B; to
improve steps I and II models, other levels, rather than the production
of a more optimized height map, should initially be the focus. For ex-
ample, introducing additional features could present a possible source
of improvement. Calculating variables over an extended time span or
S1 time series (Ge et al., 2022), could offer insights into temporal
dynamics. Extracting textural features from optical images, as described
by Couteron et al. (2015), or exploring other GEDI waveform derived
indicators may enhance predictions. Furthermore, using a radiative
transfer model to simulate GEDI data at the level of NFI plots could
provide insight into the link between forest characteristics and GEDI
measurements.
## Conclusion

In this study, we applied a kNN bagging approach using Sentinel and GEDI datasets to predict forest GSV. The approach involved imputing GEDI variables using Sentinel auxiliary variables, to subsequently predict GSV. Three strategies were tested to impute GEDI data using different sets of auxiliary data: (A) using only Sentinel data, (B) using Sentinel and a maximum height variable available on GEDI footprints and on NFI plots, and (C) using Sentinel and a national height map. The outcomes emphasized the inadequacy of relying only on Sentinel data for imputing GEDI variables. The use of maximum height variables largely improved estimations. However, the height map used in Strategy C seems to contain uncertainties that tend to reduce its ability to support GEDI imputations. GSV is influenced by more than just height; therefore, the inclusion of various structure-related variables, as available in GEDI data, contribute to improving predictive performance.

Strategy B outperformed all others tested in this study, and appeared to be the most promising strategy. However, it does not allow for pixel-wise mapping of GSV. Even so, the densified GSV points estimated at each GEDI position hold potential for small area estimation, thereby offering insights into local assessments.

In summary, our study highlights the complexities involved in predicting NFI volumes using a kNN bagging approach with Sentinel and GEDI datasets. The need for careful consideration of auxiliary data, geolocation refinement, and methodological stability is emphasized, thereby providing valuable insights for future research and applications in forest volume estimation.

## CRediT authorship contribution statement

**Anouk Schleich:** Writing – original draft, Software, Methodology, Formal analysis, Conceptualization. **Cédric Vega:** Writing – review & editing, Supervision, Software, Methodology, Formal analysis, Conceptualization. **Jean-Pierre Renaud:** Writing – review & editing, Methodology, Formal analysis, Conceptualization. **Olivier Bouriaud:** Writing – review & editing, Methodology, Formal analysis, Conceptualization.
**Sylvie Durrieu:** Writing – review & editing, Supervision, Methodology, Formal analysis, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This study was carried out in the framework of the SLIM project, supported by the TOSCA Continental Surface program of the Centre National d'Etudes Spatiales (CNES) (order N° 4500066524). It also benefited from the cofounding by INRAE and IGN for the PhD thesis of Anouk Schleich. Olivier Bouriaud acknowledges funding from EU Horizon Europe Research and Innovation Programme under Grant Funded by the European Union, project INFORMA Grant agreement No. 101060309. We also thank Jean-Baptiste Féret for the quick introduction to the R libraries preprocS2 and prosail.

## Appendix A. Equations of evaluation metrics

$$ \text{MAE} = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i| \qquad (A.1) $$

$$ \text{MAE}\% = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{|y_i - \hat{y}_i|}{y_i} \right) \times 100 \qquad (A.2) $$

$$ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} \qquad (A.3) $$

$$ \text{RMSE}\% = \frac{\text{RMSE}}{\bar{y}} \times 100 \qquad (A.4) $$

$$ R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \qquad (A.5) $$

$$ r = \frac{\sum_{i=1}^{n} (y_i - \bar{y})^2 \sqrt{\sum_{i=1}^{n} (\hat{y}_i - \bar{y})^2}}{\sqrt{\sum_{i=1}^{n} (y_i - \bar{y})^2} \sqrt{\sum_{i=1}^{n} (\hat{y}_i - \bar{y})^2}} \qquad (A.6) $$

## Appendix B. Results using datasets without GeoGEDI correction

### B.1. Step I

Other S2 variables show correlations above 0.33 with one or more GEDI variables (i.e. lai_mean, fCover_mean, ndvi_mean, msavi_mean) but are not selected because they are correlated by more than 0.85 to a previously chosen variable. Fig. B.6 shows results on the 500 GEDI footprint test dataset, comparing imputed to observed values for G\_RHv\_100 (see Table B.6).

### Table B.6

Variables selected for step I. The first column contains the selected variables, the second column contains the GEDI variable to which the auxiliary variable is the most correlated, and the third column indicates the correlation between the two. Variables start with “S1” or “S2” if they come from Sentinel-1 or Sentinel-2 data, followed by the variable name, followed by 15 or 50 depending if the 15 m or 50 m radius was used for the extraction of the variable. Variables end with “mean” or “sd” depending on the zonal extraction of mean or standard deviation. Auxiliary height for Strategy C is FORMS-H_15_mean. Auxiliary height for Strategy B is RHv_100 from GEDI.

<table><thead><tr><td>Auxiliary variable</td><td>Most correlated GEDI variable</td><td>Corr</td></tr></thead><tbody><tr><td>Hmax for Strategy B</td><td>RHv_100</td><td>1.00</td></tr><tr><td>FORMS-H_15_mean for Strategy C</td><td>RHv_90</td><td>0.70</td></tr><tr><td>S2_fAPAR_15_mean</td><td>G_pavd_z3</td><td>0.37</td></tr><tr><td>S2.greenness_15_mean</td><td>G_cover_zf</td><td>0.37</td></tr><tr><td>S2_ndvi_50_mean</td><td>G_cover_z2</td><td>0.33</td></tr></tbody></table>

Subsequently, for further processing only the GEDI variables which were considered well imputed (correlation between observed and best imputations on the test dataset > 0.6), were kept.

For Strategy A: /

For Strategy B: G\_RHv\_100, G\_RHv\_90, G\_RHv\_80, G\_RHv\_70, G\_RHv\_60, G\_RHv\_50, G\_RHv\_40, G\_cover\_z4, G\_cover\_z5, G\_fhd\_normal, G\_pai\_z5, G\_pavd\_z6, and G\_cover\_z5z6.

For Strategy C: /

### B.2. Step II

This results in graphics in Fig. B.7 (see Table B.7).

For errors of Strategy B using uncorrected GEDI positions, refer to Table 5.

## Appendix C. Step II: the 1000 predictions from the 1000 kNN runs

This Section shows results of the 1000 kNN runs, before aggregating the single predictions into mean values. Fig. C.9 shows the distribution of standard deviations within the 1000 predictions for all NFI plots.

Fig. C.10 shows examples of the 1000 kNN predicted GSVs for some NFI plots.

While few were adequately predicted and by numerous kNN runs (C.10(a),C.10(b),C.10(c)), some exhibited two modes (C.10(i),C.10(j)), some yielded very variable predictions resulting in mean estimates (C.10(l),C.10(k),C.10(e)), and others were simply erroneously predicted (C.10(f),C.10(d),C.10(g),C.10(h)). C.10(k) and C.10(l) have very similar-looking prediction histograms, resulting in similar aggregated final predictions. However, their observed GSV values (in blue) are quite different.
<figure><img src="image_8.png" /><figcaption>**Fig. B.6.** G_RHv_100 imputations for 500 test GEDI footprints - Without GeoGEDI correction. Correlation, *R*², RMSE and RMSE% are reported in the upper left corner of each plot.</figcaption></figure>
<table><caption>Table B.7</caption><caption>Selected variables for step II - Without GeoGEDI correction.</caption><thead><tr><th colspan="2">Strategy A</th><th colspan="2">Strategy B</th><th colspan="2">Strategy C</th></tr><tr><th>Variable</th><th>Corr</th><th>Variable</th><th>Corr</th><th>Variable</th><th>Corr</th></tr></thead><tbody><tr><td>S2_B6_15_mean</td><td>-0.42</td><td>G_RHv_100</td><td>0.68</td><td>FORMS-H_15_mean</td><td>0.59</td></tr><tr><td>S2_fCover_15_mean</td><td>-0.37</td><td>G_fhd_normal</td><td>0.54</td><td>S2_B6_15_mean</td><td>-0.42</td></tr><tr><td>S2_B11_15_mean</td><td>-0.34</td><td>G Cover_z5</td><td>0.53</td><td>S2_fCover_15_mean</td><td>-0.37</td></tr><tr><td></td><td></td><td>S2_B6_15_mean</td><td>-0.42</td><td>S2_B11_15_mean</td><td>-0.34</td></tr><tr><td></td><td></td><td>S2_fCover_15_mean</td><td>-0.37</td><td></td><td></td></tr><tr><td></td><td></td><td>G Cover_z5z6</td><td>0.37</td><td></td><td></td></tr><tr><td></td><td></td><td>S2_B11_15_mean</td><td>-0.34</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. B.7. GSV predictions for test NFI plots with Strategy B.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. B.8.** Boxplot of errors using overall and stand-specific models, compared to FORMS-V estimates. The boxes constitute medians, 1st quartiles (Q1) and 3rd quartiles (Q3) errors and 95% confidence intervals are represented by horizontal lines.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. C.9.** Distribution of standard deviations within the 1000 kNN predictions. The red line presents the mean.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. C.10.** A thousand predictions by NFI plot. The blue line represents the observed GSV value, the red line the mean of the 1000 predicted GSVs used in our methodology, and the green line represents the median of the 1000 predicted GSVs.</figcaption></figure>

## Data availability

The data that has been used is confidential.

## References

Aleissaee, A.A., Kumar, A., Anwer, R.M., Khan, S., Cholakkal, H., Xia, G.S., Khan, F.S., 2023. Transformers in remote sensing: A survey. Remote. Sens. 15 (7), 1860. http://dx.doi.org/10.3390/rs15071860.

Appiah Mensah, A., Jonzén, J., Nyström, K., Wallerman, J., Nilsson, M., 2023. Mapping site index in coniferous forests using bi-temporal airborne laser scanning data and field data from the Swedish national forest inventory. Forest Ecol. Manag. 547, 121395. http://dx.doi.org/10.1016/j.foreco.2023.121395.

Beck, J., Wirt, B., Luthcke, S., Hofton, M., Armston, J., 2021. Global Ecosystem Dynamics Investigation (GEDI) level 1B user guide. Document version 2.0. SGS Earth Resources Observation and Science (EROS) Center: NASA's Land Pro- cesses Distributed Active Archive Center (LP DAAC), URL https://lpdaac.usgs.gov/ documents/987/GEDI01B_User_Guide_V2.pdf.

Besic, N., Durrieu, S., Schleich, A., Vega, C., 2024. Using structural class pairing to address the spatial mismatch between GEDI measurements and NFI plots. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 17, 12854–12867. http://dx.doi.org/10. 1109/JSTARS.2024.3425431.

Bouriaud, O., Morneau, F., Bontemps, J.D., 2023. Square-grid sampling support to reconcile systematicity and adaptivity in periodic spatial surveys of natural populations. J. Veg. Sci. 34 (3), http://dx.doi.org/10.1111/jvs.13195.

Bouvier, M., Durrieu, S., Fournier, R., Saint-Geours, N., Guyon, D., Grau, E., De Boissieu, F., 2019. Influence of sampling design parameters on biomass predictions derived from airborne LiDAR data. Can. J. Remote Sens. 45, 1–23. http://dx.doi.org/10.1080/07038992.2019.1669013.

Breiman, L., 1996. Bagging predictors. Mach. Learn. 24 (2), 123–140. http://dx.doi.org/10.1007/BF00058655.

Bruening, J., May, P., Armston, J., Dubayah, R., 2023. Precise and unbiased biomass estimation from GEDI data and the US forest inventory. Front. For. Glob. Chang. 6, http://dx.doi.org/10.3389/ffgc.2023.1149153.

Bullock, E.L., Healey, S.P., Yang, Z., Acosta, R., Villalba, H., Insfrán, K.P., Melo, J., Wilson, S., Duncanson, L.I., Næsset, E., Armston, J., Saarela, S., Ståhl, G., Patter- son, P.L., Dubayah, R., 2023. Estimating aboveground biomass density using hybrid statistical inference with GEDI lidar data and paraguay's national forest inventory. Environ. Res. Lett. http://dx.doi.org/10.1088/1748-9326/acdf03.

Cavaignac, S., 2009. Les sylvoécorégions (SER) de France métropolitaine : Étude de définition. IGN URL https://inventaire-forestier.ign.fr/IMG/pdf/Part1_rapport-ser. pdf.

Charbonneau, F., Trudel, M., Fernandes, R., 2005. Ruse of dual polarization and multi- incidence SAR for soil permeability mapping. In: In Proceedings of the 2005 Advanced Synthetic Aperture Radar (ASAR) Workshop. St-Hubert, QC, Canada.

Chirici, G., Mura, M., McInnerney, D., Py, N., Tomppo, E., Waser, L., Travaglini, D., Mcroberts, R., 2016. A meta-analysis and review of the literature on the k-nearest neighbors technique for forestry applications that use remotely sensed data. Remote Sens. Environ. 176, 282–294. http://dx.doi.org/10.1016/j.rse.2016.02.001.
Coops, N.C., Tompalski, P., Goodbody, T.R.H., Achim, A., Mulverhill, C., 2023. Framework for near real-time forest inventory using multi source remote sensing data. For.: An Int. J. For. Res. 96 (1), 1-19. http://dx.doi.org/10.1093/forestry/cpac015.

Copernicus Sentinel-2, 2021. (Processed by ESA), MSI level-2A BOA reflectance product. Collection 1. http://dx.doi.org/10.5270/S2-znk9xsj

Couteron, P., Barbier, N., Deblauwe, V., Pelissier, R., Ploton, P., 2015. Texture analysis of very high spatial resolution optical images as a way to monitor vegetation and forest biomass in the tropics. In: Murthy, M., Wesselman, S., Gilani, H. (Eds.), Multi-Scale Forest Biomass Assessment and Monitoring in the Hindu Kush Himalayan Region: A Geospatial Perspective. In: ICIMOD special publication, International Centre for Integrated Mountain Development (ICIMOD), pp. 157-164, URL https://hal.umontpellier.fr/hal-02303735.

Cover, T., Hart, P., 1967. Nearest neighbor pattern classification. IEEE Trans. Inform. Theory 13 (1), 21-27. http://dx.doi.org/10.1109/TIT.1967.1053964.

Crookston, N.L., Finley, A.O., 2007. yalmpute: An R package for kNN imputation. J. Stat. Softw. 23 (10), http://dx.doi.org/10.18637/jss.v023.i10.

Dubayah, R., Armston, J., Healey, S., Yang, Z., Patterson, P., Saarela, S., Stahl, G., Duncanson, L., Kellner, J., Bruening, J., Pascual, A., 2023. GEDI L4B gridded aboveground biomass density, version 2.1. http://dx.doi.org/10.3334/ORNLDAAC/2299, ORNL Distributed Active Archive Center URL https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=2299.

Dubayah, R., Armston, J., Kellner, J., Duncanson, L., Healey, S., Patterson, P., Hancock, S., Tang, H., Bruening, J., Hofton, M., Blair, J., Luthcke, S., 2022. GEDI L4A footprint level aboveground biomass density, version 2.1. http://dx.doi.org/10.3334/ORNLDAAC/2056, ORNL Distributed Active Archive Center, URL https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=2056.

Dubayah, R., Blair, J.B., Goetz, S., Fatoyinbo, L., Hansen, M., Healey, S., Hofton, M., Hurtt, G., Kellner, J., Luthcke, S., Armston, J., Tang, H., Duncanson, L., Hancock, S., Jantz, P., Marselis, S., Patterson, P.L., Qi, W., Silva, C., 2020. The global ecosystem dynamics investigation: High-resolution laser ranging of the earth's forests and topography. Sci. Remote. Sens. 1, 100002. http://dx.doi.org/10.1016/j.srs.2020.100002.

Dubayah, R., Luthcke, S., Sabaka, T., Nicholas, J., Preaux, S., Hofton, M., 2021. GEDI L3 gridded land surface metrics, version 2. http://dx.doi.org/10.3334/ORNLDAAC/1952, ORNL Distributed Active Archive Center, URL https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1952.

Duncanson, L., Kellner, J.R., Armston, J., Dubayah, R., Minor, D.M., Hancock, S., Healey, S.P., Patterson, P.L., Saarela, S., Marselis, S., Silva, C.E., Bruening, J., Goetz, S.J., Tang, H., Hofton, M., Blair, B., Luthcke, S., Fatoyinbo, L., Abernethy, K., Alonso, A., Andersen, H.E., Aplin, P., Baker, T.R., Barbier, N., Bastin, J.F., Biber, P., Boeckx, P., Bogaert, J., Boshetti, L., Moucher, P.B., Boyd, D.S., Burslem, D.F.R.P., Calvo-Rodriguez, S., Chave, J., Chazdon, R.L., Clark, D.B., Clark, D.A., Cohen, W.B., Coomes, D.A., Corona, P., Cushman, K.C., Cutler, M.E.J., Dalling, J.W., Dalponte, M., Dash, J., de-Miguel, S., Deng, S., Ellis, P.W., Erasmus, B., Fekety, P.A., Fernandez-Landa, A., Ferraz, A., Fischer, R., Fisher, A.G., García-Abril, A., Gobakken, T., Hacker, J.M., Heurich, M., Hill, R.A., Hopkinson, C., Huang, H., Hubbell, S.P., Hudak, A.T., Huth, A., Imbach, B., Jeffery, K.J., Katoh, M., Kearsley, E., Kenfack, D., Kljun, N., Knapp, N., Král, K., Krůček, M., Labrière, N., Lewis, S.L., Longo, M., Lucas, R.M., Main, R., Manzanera, J.A., Martínez, R.V., Mathieu, R., Memiaghe, H., Meyer, V., Mendoza, A.M., Montanez, P., Morsdorf, F., Næsset, E., Naidoo, L., Nilus, R., O'Brien, M., Orwig, D.A., Papathanassiou, K., Parker, G., Philipson, C., Phillips, O.L., Pisek, J., Poulsen, J.R., Pretzsch, H., Rüdiger, C., Saatchi, S., Sanchez-Azofeifa, A., Sanchez-Lopez, N., Scholes, R., Silva, C.A., Simard, M., Skidmore, A., Stereiczczak, K., Tanase, M., Torresan, C., Valbuena, R., Verbeeck, H., Vrska, T., Wessels, K., White, J.C., White, L.J.T., Zahabu, E., Zgraggen, C., 2022. Aboveground biomass density models for NASA's global ecosystem dynamics investigation (GEDI) lidar mission. Remote Sens. Environ. 270, 112845. http://dx.doi.org/10.1016/j.rse.2021.112845.

Ehlers, S., Saarela, S., Lindgren, N., Lindberg, A., Nyström, M., Persson, H., Olsson, H., Ståhl, G., 2018. Assessing error correlations in remote sensing-based estimates of forest attributes for improved composite estimation. Remote. Sens. 10, 667. http://dx.doi.org/10.3390/rs10050667.

Eng, L., Ismail, R., Hashim, W., Baharum, A., 2019. The use of VARI, GLI, and Vigreen formulas in detecting vegetation in aerial images. Int. J. Technol. 10, 1385. http://dx.doi.org/10.14716/ijtech.v10i17.3275.

ESA (European Space Agency), URL https://sentinels.copernicus.eu/web/sentinel-missions/sentinel-2. (Accessed 14 November 2023).

ESA (European Space Agency), URL https://sentinels.copernicus.eu/web/sentinel-technical-guides/sentinel-1-sar. (Accessed 14 November 2023).

Fadili, M., Renaud, J.P., Bock, J., Vega, C., 2019. RegisTree: A registration algorithm to enhance forest inventory plot georeferencing. Ann. For. Sci. 76, http://dx.doi.org/10.1007/s13595-019-0814-2.

Feret, J.B., de Boissieu, F., 2023. prosail: PROSAIL leaf and canopy radiative transfer model and inversion routines. R package version 2.2.2, URL https://gitlab.com/jbferet/prosail.

Forzieri, G., Girardello, M., Ceccherini, G., Spinoni, J., Feyen, L., Hartmann, H., Beck, P.S.A., Camps-Valls, G., Chirici, G., Mauri, A., Cescatti, A., 2021. Emergent vulnerability to climate-driven disturbances in European forests. Nat. Commun. 12 (1), 1081. http://dx.doi.org/10.1038/s41467-021-21399-7.

Gao, B.c., 1996. NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sens. Environ. 58 (3), 257-266. http://dx.doi.org/10.1016/S0034-4257(96)00067-3.

Ge, S., Gu, H., Su, W., Praks, J., Antropov, O., 2022. Improved semisupervised unet deep learning model for forest height mapping with satellite SAR and optical data. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 15, 5776-5787. http://dx.doi.org/10.1109/JSTARS.2022.3188201.

Ginzler, C., Hobi, M.L., 2015. Countrywide stereo-image matching for updating digital surface models in the framework of the Swiss National Forest Inventory. Remote. Sens. 7 (4), 4343-4370. http://dx.doi.org/10.3390/rs70404343.

Gomes, H.M., Read, J., Bifet, A., Durrant, R.J., 2021. Learning from Evolving Data Streams through Ensembles of Random Patches. Open Access Te Herenga Waka-Victoria University of Wellington, http://dx.doi.org/10.25455/wgtn.20104412.

Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., Moore, R., 2017. Google earth engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. http://dx.doi.org/10.1016/j.rse.2017.06.031.

Gschwantner, T., Alberdi, I., Balázs, A., Bauwens, S., Bender, S., Borota, D., Bosela, M., Bouriaud, O., Cañellas, I., Donis, J., Freudenschuß, A., Hervé, J.C., Hladnik, D., Jan-sons, J., Kolozs, L., Korhonen, K.T., Kucera, M., Kulbokas, G., Kuliešis, A., Lanz, A., Jeuneve, P., Lind, T., Marin, G., Morneau, F., Nagy, D., Nord-Larsen, T., Nunes, L., Pantić, D., Paulo, J.A., Pikula, T., Redmond, J., Rego, F.C., Riedel, T., Saint-André, L., Šebeň, V., Sims, A., Skudnik, M., Solti, G., Tomter, S.M., Twomey, M., Westerlund, B., Zell, J., 2019. Harmonisation of stem volume estimates in European national forest inventories. Ann. For. Sci. 76 (1), 1-23. http://dx.doi.org/10.1007/s13595-019-0800-8.

Gschwantner, T., Alberdi, I., Bauwens, S., Bender, S., Borota, D., Bosela, M., Bouriaud, O., Breidenbach, J., Donis, J., Fischer, C., Gasparini, P., Heffernan, L., Hervé, J.-C., Kolozs, L., Korhonen, K.T., Koutsias, N., Kovácsévics, P., Kučera, M., Kulbokas, G., Kuliešis, A., Lanz, A., Jeuneve, P., Lind, T., Marin, G., Morneau, F., Nord-Larsen, T., Nunes, L., Pantić, D., Redmond, J., Rego, F.C., Riedel, T., Šebeň, V., Sims, A., Skudnik, M., Tomter, S.M., 2022. Growing stock monitoring by European national forest inventories: Historical origins, current methods and harmonisation. Forest Ecol. Manag. 505, 119868. http://dx.doi.org/10.1016/j.foreco.2021.119868.

Guldin, R.W., 2021. A systematic review of small domain estimation research in forestry during the twenty-first century from outside the United States. Front. For. Glob. Chang. 4, http://dx.doi.org/10.3389/ffgc.2021.695929.

Guo, Q., Du, S., Jiang, J., Guo, W., Zhao, H., Yan, X., Zhao, Y., Xiao, W., 2023. Combining GEDI and sentinel data to estimate forest canopy mean height and aboveground biomass. Ecol. Informatics 78, 102348. http://dx.doi.org/10.1016/j.ecoinf.2023.102348.

Gupta, R., Sharma, L., 2022. Mixed tropical forests canopy height mapping from spaceborne LiDAR GEDI and multisensor imagery using machine learning models. Remote. Sens. Appl.: Soc. Environ. 27, 100817. http://dx.doi.org/10.1016/j.rsase.2022.100817.

Haakana, H., Heikikinen, J., Katila, M., Kangas, A., 2019. Efficiency of post-stratification for a large-scale forest inventory—case Finnish NFI. Ann. For. Sci. 76, http://dx.doi.org/10.1007/s13595-018-0795-6.

Hagolle, O., Morin, D., Kadiri, M., 2018. Detailed Processing Model for the Weighted Average Synthesis Processor (WASP) for Sentinel-2. Zenodo, http://dx.doi.org/10.5281/zenodo.1401360.

Hamedianfar, A., Mohamedou, C., Kangas, A., Vauhkonen, J., 2022. Deep learning for forest inventory and planning: a critical review on the remote sensing approaches so far and prospects for further applications. For.: An Int. J. For. Res. 95 (4), 451-465. http://dx.doi.org/10.1093/forestry/cpac002, arXiv:https://academic.oup.com/forestry/article-pdf/95/4/451/45293980/cpac002.pdf.

Henrich, V., Krauss, G., Göte, C., Sandow, C., 2012. IDB - www.indexdatabase.de, entwicklung einer datenbank für fernerkundungsindizes. AK Fernerkundung, Bochum URL https://www.indexdatabase.de/info/credits.php.

Ho, T.K., 1998. Nearest neighbors in random subspaces. In: Amin, A., Dori, D., Pudil, P., Freeman, H. (Eds.), Advances in Pattern Recognition. In: Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, pp. 640-648. http://dx.doi.org/10.1007/BFb0033288.

Hunka, N., May, P., Babcock, C., de la Rosa, J.A.A., de los Ángeles Soriano-Luna, M., Saucedo, R.M., Armston, J., Santoro, M., Suarez, D.R., Herold, M., Málaga, N., Healey, S.P., Kennedy, R.E., Hudak, A.T., Duncanson, L., 2025. A geostatistical approach to enhancing national forest biomass assessments with earth observation to aid climate policy needs. Remote Sens. Environ. 318, 114557. http://dx.doi.org/10.1016/j.rse.2024.114557.

Hunka, N., Santoro, M., Armston, J., Dubayah, R., McRoberts, R.E., Næsset, E., Quegan, S., Urbazaev, M., Pascual, A., May, P.B., Minor, D., Leitold, V., Basak, P., Liang, M., Melo, J., Herold, M., Málaga, N., Wilson, S., Durán Montesinos, P.,
Arana, A., Ernesto De La Cruz Paiva, R., Ferrand, J., Keoka, S., Guerra-Hernández, J., Duncanson, L., 2023. On the NASA gedi and ESA cci biomass maps: Aligning for uptake in the UNFCCC global stocktake. Environ. Res. Lett. 18 (12), 124042. http://dx.doi.org/10.1088/1748-9326/ad0b60.

Hunt, Jr., E.R., Daughtry, C.S.T., Eitel, J.U.H., Long, D.S., 2011. Remote sensing leaf chlorophyll content using a visible band index. Agron. J. 103, 1090-1099. http://dx.doi.org/10.2134/agronj2010.0395.

IGN, 2022. RGE ALTI. https://geoservices.ign.fr/rgealti. (Accessed 6 March 2022).

IGN, 2023a. BD Forêt version 2. URL https://geoservices.ign.fr/bdforet. (Accessed 21 June 2023).

IGN, 2023b. **Méthodologie : Résultats d'inventaire forestier** – pour bien comprendre les résultats publiés. URL https://inventaire-forestier.ign.fr/IMG/pdf/methodologie-2023.pdf. (Accessed 11 December 2023).

IGN, 2023c. **Sylvéoécorégion.** URL https://inventaire-forestier.ign.fr/spip.php?article773. (Accessed 21 June 2023).

Irulappa Pillai Vijayakumar, D.B., Renaud, J.P., Morneau, F., Mcroberts, R., Vega, C., 2019. Increasing precision for french forest inventory estimates using the K-NN technique with optical and photogrammetric data and model-assisted estimators. Remote. Sens. 11, 991. http://dx.doi.org/10.3390/rs11080991.

Kacic, P., Hirner, A., Canova, E., 2021. Fusing sentinel-1 and 2 to model GEDI-derived vegetation structure characteristics in GEE for the paraguayan chaco. Remote. Sens. 13, 5105. http://dx.doi.org/10.3390/rs12425105.

Kakoulaki, G., Martinez, F.A.M., Florio, P., 2021. Non-commercial light detection and ranging (LiDAR) data in Europe. JRC Publ. Repos. http://dx.doi.org/10.2760/212427, https://publications.jrc.ec.europa.eu/repository/handle/JRC126223.

Kim, Y., Zyl, J., 2009. A time-series approach to estimate soil moisture using polarimetric radar data. IEEE Trans. Geosci. Remote Sens. 47, 2519-2527. http://dx.doi.org/10.1109/TGRS.2009.2014944.

Labrière, N., Davies, S.J., Disney, M.I., Duncanson, L.I., Herold, M., Lewis, S.L., Phillips, O.L., Quegan, S., Saatchi, S.S., Schepaschenko, D.G., Scipal, K., Sist, P., Chave, J., 2023. Toward a forest biomass reference measurement system for remote sensing applications. Global Change Biol. 29 (3), 827-840. http://dx.doi.org/10.1111/gcb.16497.

Lang, N., Jetz, W., Schindler, K., Wegner, J.D., 2023. A high-resolution canopy height model of the Earth. Nat. Ecol. Evol. 7 (11), 1778-1789. http://dx.doi.org/10.1038/s41559-023-02206-6.

Li, W., Niu, Z., Shang, R., Qin, Y., Wang, L., Chen, H., 2020. High-resolution mapping of forest canopy height using machine learning by coupling ICESat-2 LiDAR with Sentinel-1, Sentinel-2 and Landsat-8 data. Int. J. Appl. Earth Obs. Geoinf. 92, 102163. http://dx.doi.org/10.1016/j.jag.2020.102163.

Lim, K., Treitz, P., Wulder, M., St-Onge, B., Flood, M., 2003. LIDAR remote sensing of forest structure. Prog. Phys. Geogr. 27, 88-106. http://dx.doi.org/10.1191/0309133303pp360ra.

Lindgren, N., Olsson, H., Nyström, K., Nyström, M., Ståhl, G., 2021. Data assimilation of growing stock volume using a sequence of remote sensing data from different sensors. Can. J. Remote Sens. 48 (2), 127-143. http://dx.doi.org/10.1080/07038992.2021.1988542.

Louhaichi, M., Borman, M., Johnson, D., 2001. Spatially located platform and aerial photography for documentation of grazing impacts on wheat. Geocart Int. 16, http://dx.doi.org/10.1080/1010604075842184.

Louppe, G., Geurts, P., 2012. Ensembles on random patches. In: Flach, P.A., De Bie, T., Cristianini, N. (Eds.), Machine Learning and Knowledge Discovery in Databases. In: Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, pp. 346-361. http://dx.doi.org/10.1007/978-3-642-33460-3_28.

LP DAAC, 2023. NASA announces pause in GEDI mission. URL https://www.earthdata.nasa.gov/news/nasa-announces-pause-gedi-mission?. (Accessed 27 March 2023).

Lu, D., Jiang, X., 2024. A brief overview and perspective of using airborne Lidar data for forest biomass estimation. Int. J. Image Data Fusion 15 (1), 1-24. http://dx.doi.org/10.1080/19479832.2024.2309615.

May, P.B., Dubayah, R.O., Bruening, J.M., Gaines, G.C., 2024. Connecting spaceborne Lidar with NFI networks: A method for improved estimation of forest structure and biomass. Int. J. Appl. Earth Obs. Geoinf. 129, 103797. http://dx.doi.org/10.1016/j.jag.2024.103797.

Mcroberts, R., Chen, Q., Walters, B., Kaisershot, D., 2018. The effects of global positioning system receiver accuracy on airborne laser scanning-assisted estimates of aboveground biomass. Remote Sens. Environ. 207, 42-49. http://dx.doi.org/10.1016/j.rse.2017.09.036.

Moser, P., Vibrans, A.C., McRoberts, R.E., Naesset, E., Gobakken, T., Chirici, G., Mura, M., Marchetti, M., 2017. Methods for variable selection in LiDAR-assisted forest inventories. For.: An Int. J. For. Res. 90 (1), 112-124. http://dx.doi.org/10.1093/forestry/cpw041.

Nasizadehdihazaji, R., Balik Sanli, F., Abdikan, S., Cakir, Z., Sekertekin, A., Ustuner, M., 2019. Sensitivity analysis of multi-temporal Sentinel-1 SAR parameters to crop height and canopy coverage. Appl. Sci. 9 (4), 655. http://dx.doi.org/10.3390/app9040655.

Novo-Fernández, A., Barrio-Anta, M., Recondo, C., Cámara-Obregón, A., López-Sánchez, C.A., 2019. Integration of national forest inventory and nationwide airborne laser scanning data to improve forest yield predictions in North-Western Spain. Remote. Sens. 11 (14), 1693. http://dx.doi.org/10.3390/rs11141693.

Ortí-Reyes, A.D., Barrera-Ortega, D., Velasco-Bautista, E., Romero-Sánchez, M.E., Díaz, A.C., 2024. Predicting forest parameters through generalized linear mixed models using GEDI metrics in a temperate forest in Oaxaca, Mexico. Int. J. Remote Sens. 45 (22), 8037-8060. http://dx.doi.org/10.1080/01431161.2024.2396565.

Osborne, C., 1991. Statistical calibration: A review. Rev. de Stat. (Int. Stat. Rev.) 59 (3), 309-336. http://dx.doi.org/10.2307/1403690, arXiv:1403690.

Pascual, A., Guerra-Hernández, J., Armston, J., Minor, D.M., Duncanson, L.I., May, P.B., Kellner, J.R., Dubayah, R., 2023. Assessing the performance of NASA's GEDI L4A footprint aboveground biomass density models using national forest inventory and airborne laser scanning data in Mediterranean forest ecosystems. Forest Ecol. Manag. 538, 120975. http://dx.doi.org/10.1016/j.foreco.2023.120975.

Pascual, A., May, P.B., Cárdenas-Martínez, A., Guerra-Hernández, J., Hunka, N., Bruening, J.M., Healey, S.P., Armston, J.D., Dubayah, R.O., 2025. Calibration of GEDI footprint aboveground biomass models in Mediterranean forests with NFI plots: A comparison of approaches. J. Environ. Manag. 375, 124313. http://dx.doi.org/10.1016/j.jenvman.2025.124313.

Pereira-Pires, J.E., Mora, A., Aubard, V., Silva, J.M.N., Fonseca, J.M., 2021. Assessment of Sentinel-2 spectral features to estimate forest height with the new GEDI data. In: Camarinha-Matos, L.M., Ferreira, P., Brito, G. (Eds.), Technological Innovation for Applied AI Systems. In: IFIP Advances in Information and Communication Technology, Springer International Publishing, Cham, pp. 123-131. http://dx.doi.org/10.1007/978-3-030-78288-7_12.

Potapov, P., Li, X., Hernandez-Serna, A., Tyukavina, A., Hansen, M.C., Kommareddy, A., Pickens, A., Turubanova, S., Tang, H., Silva, C.E., Armston, J., Dubayah, R., Blair, J.B., Hofton, M., 2021. Mapping global forest canopy height through integration of GEDI and Landsat data. Remote Sens. Environ. http://dx.doi.org/10.1016/j.rse.2020.112165.

Rouse, J., Haas, R.H., Schell, J.A., Deering, D., 1973. Monitoring vegetation systems in the great plains with ERTS. In: Third ERTS Sympos NASA.

Saarela, S., Holm, S., Healey, S.P., Andersen, H.-E., Petersen, W., Patterson, P.L., Naesset, E., Gregoire, T.G., Ståhl, G., 2018. Generalized hierarchical model-based estimation for aboveground biomass assessment using gedi and Landsat data. Remote. Sens. 10 (11), 1832. http://dx.doi.org/10.3390/rs10111832.

Saarela, S., Schnell, S., Tuominen, S., Balazs, A., Hyypäpp, J., Grafström, A., Ståhl, G., 2016. Effects of positional errors in model-assisted and model-based estimation of growing stock volume. Remote Sens. Environ. 172, 101-108. http://dx.doi.org/10.1016/j.rse.2015.11.002.

Saborowski, J., Marx, A., Nagel, J., Böckmann, T., 2010. Double sampling for stratification in periodic inventories—Infinite population approach. Forest Ecol. Manag. 260 (10), 1866-1895. http://dx.doi.org/10.1016/j.foreco.2010.08.035.

Sagar, A., Vega, C., Bouriaud, O., Piedallu, C., Renaud, J.P., 2022. Multisource forest inventories: A model-based approach using k-NN to reconcile forest attributes statistics and map products. ISPRS J. Photogramm. Remote Sens. 192, 175-188. http://dx.doi.org/10.1016/j.isprsjprs.2022.08.016.

Schleich, A., Bouriaud, O., Vega, C., Durrieu, S., 2024. Showcasing the Effectiveness of GEDI Footprints as First-Phase Sample for Forest Inventories Based on Double Sampling for Post-Stratification. Institute of Electrical and Electronics Engineers (IEEE), http://dx.doi.org/10.36227/techxiv.173092115.59963750/v1.

Schleich, A., Durrieu, S., Soma, M., Vega, C., 2023. Improving GEDI footprint geolocation using a high resolution digital elevation model. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 16, 7718-7732. http://dx.doi.org/10.1109/JSTARS.2023.3298991.

Schwartz, M., Ciais, P., Ottlé, C., De Truchis, A., Vega, C., Fayad, I., Brandt, M., Fensholt, R., Baghdadhi, N., Morneau, F., Morin, D., Guyon, D., Dayau, S., Wigneron, J.P., 2022. High-resolution canopy height map in the lands forest (France) based on GEDI, Sentinel-1, and Sentinel-2 data with a deep learning approach. http://dx.doi.org/10.48550/arXiv.2212.10265, arXiv arXiv:arXiv:2212.10265.

Schwartz, M., Ciais, P., Truchis, A., Chave, J., Ottlé, C., Vega, C., Wigneron, J.P., Nicolas, M., Jouaber, S., Liu, S., Brandt, M., Fayad, I., 2023. FORMS: forest multiple source height, wood volume, and biomass maps in France at 10 to 30 m resolution based on Sentinel-1, Sentinel-2, and GEDI data with a deep learning approach. Earth Syst. Sci. Data http://dx.doi.org/10.5194/essd-2023-196.

Shendryk, Y., 2022. Fusing GEDI with earth observation data for large area above-ground biomass mapping. Int. J. Appl. Earth Obs. Geoinf. 115, 103108. http://dx.doi.org/10.1016/j.jag.2022.103108.

Silva, C., Duncanson, L., Hancock, S., Neuenschwander, A., Thomas, N., Hofton, M., Fatoyinbo, L., Simard, M., Marshak, C., Armston, J., Lutchke, S., Dubayah, R., 2021. Fusing simulated GEDI, ICESat-2 and NISAR data for regional aboveground biomass mapping. Remote Sens. Environ. 253, 112234. http://dx.doi.org/10.1016/j.rse.2020.112234.
Sothe, C., Gonsamo, A., Lourenço, R.B., Kurz, W.A., Snider, J., 2022. Spatially continuous mapping of forest canopy height in Canada by combining GEDI and ICESat-2 with PALSAR and Sentinel. Remote. Sens. 14 (20), 5158. http://dx.doi.org/10.3390/rs14205158.

Ståhl, G., Gobakken, T., Saarela, S., Persson, H.J., Ekström, M., Healey, S.P., Yang, Z., Holmgren, J., Lindberg, E., Nyström, K., Papucci, E., Ulvdal, P., Ørka, H.O., Næset, E., Hou, Z., Olsson, H., McRoberts, R.E., 2024. Why ecosystem characteristics predicted from remotely sensed data are unbiased and biased at the same time – and how this affects applications. For. Ecosyst. 11, 100164. http://dx.doi.org/10.1016/j.fecs.2023.100164.

Tomppo, E., Olsson, H., Ståhl, G., Nilsson, M., Hagner, O., Katila, M., 2008. Combining national forest inventory field plots and remote sensing data for forest databases. Remote Sens. Environ. 112 (5), 1982–1999. http://dx.doi.org/10.1016/j.rse.2007.03.032.

Véga, C., St-Onge, B., 2009. Mapping site index and age by linking a time series of canopy height models with growth curves. Forest Ecol. Manag. 257 (3), 951–959. http://dx.doi.org/10.1016/j.foreco.2008.10.029.

Virtanen, J.P., Kukko, A., Kaartinen, H., Jaakkola, A., Turppa, T., Hyyppä, H., Hyyppä, J., 2017. Nationwide point cloud—The future topographic core data. ISPRS Int. J. Geo-Inf. 6 (8), 243. http://dx.doi.org/10.3390/ijgi6080243.

Vollrath, A., Mullissa, A., Reiche, J., 2020a. Angular-based radiometric slope correction for Sentinel-1 on google earth engine. Remote. Sens. 12 (11), 1867. http://dx.doi.org/10.3390/rs12111867.

Vollrath, A., Mullissa, A., Reiche, J., 2020b. Angular-based radiometric slope correction of Sentinel-1 on google earth engine. Google Colab python script URL https://github.com/ESA-PhiLab/radiometric-slope-correction.

Weiss, M., Baret, F., Jay, S., 2016. S2ToolBox level 2 products:LAI, FAPAR, FCOVER. URL https://step.esa.int/docs/extra/ATBD_S2ToolBox_V2.0.pdf. (Accessed 25 October 2023).

Westfall, J., Lister, A., Scott, C., Weber, T., 2019. Double sampling for post-stratification in forest inventory. Eur. J. For. Res. 138, http://dx.doi.org/10.1007/s10342-019-01171-9.

Zhang, S., Vega, C., Deleuze, C., Durrieu, S., Barbillon, P., Bouriaud, O., Renaud, J.P., 2022. Modelling forest volume with small area estimation of forest inventory using GEDI footprints as auxiliary information. Int. J. Appl. Earth Obs. Geoinf. 114, http://dx.doi.org/10.1016/j.jag.2022.103072.