remote sensing Improving Heterogeneous Forest Height Maps by Integrating GEDI-Based Forest Height Information in a Multi-Sensor Mapping Process

: Forests are one of the key elements in ecological transition policies in Europe. Sustainable forest management is needed in order to optimise wood harvesting, while preserving carbon storage, biodiversity and other ecological functions. Forest managers and public bodies need improved and cost-effective forest monitoring tools. Research studies have been carried out to assess the use of optical and radar images for producing forest height or biomass maps. The main limitations are the quantity, quality and representativeness of the reference data for model training. The Global Ecosystem Dynamics Investigation (GEDI) mission (full waveform LiDAR on board the International Space Station) has provided an unprecedented number of forest canopy height samples from 2019. These samples could be used to improve reference datasets. This paper aims to present and validate a method for estimating forest dominant height from open access optical and radar satellite images (Sentinel-1, Sentinel-2 and ALOS-2 PALSAR-2), and then to assess the use of GEDI samples to replace ﬁeld height measurements in model calibration. Our approach combines satellite image features and dominant height measurements, or GEDI metrics, in a Support Vector Machine regression algorithm, with a feature selection process. The method is tested on mixed uneven-aged broadleaved and coniferous forests in France. Using dominant height measurements for model training, the cross-validation shows 7.3 to 11.6% relative Root Mean Square Error (RMSE) depending on the forest class. When using GEDI height metrics instead of ﬁeld measurements for model training, errors increase to 12.8–16.7% relative RMSE. This level of error remains satisfactory; the use of GEDI could allow the production of dominant height maps on large areas with better sample representativeness. Future work will focus on conﬁrming these results on new study sites, improving the ﬁltering and processing of GEDI data, and producing height maps at regional or national scale. The resulting maps will help forest managers and public bodies to optimise forest resource inventories, as well as allow scientists to integrate these cartographic data into climate models. from 1.9 to 3.1 m (7.3 to 11.6% relative RMSE). Replacing field data with GEDI RH95% data for model learning provided satisfactory results, with RMSE from 3.0 to 4.6 m (12.8 to 16.7% relative RMSE). The study also shows the valuable contribution of GEDI data by using the GEDI RH95% as a reference dataset in the machine learning approach. This contribution is twofold: on one hand, the large amount of reference data using GEDI data allows the stepwise selection of input features in the processing chain to be bypassed, reducing the computational time and complexity. On the other hand, the use of GEDI data as reference data allows more representativeness and forest height mapping to be provided in areas where little or no in situ data are available.


Introduction
Forests cover one third of Earth's land surface and play an essential role in stabilising climate, regulating the water cycle, and providing habitat to a multitude of life forms [1,2]. Besides that, forests are of great economic benefit by providing wood and timber, used as fuel or as raw material for numerous products [1,3]. However, tropical forests are threatened by deforestation [1,2], and the effects of climate change are more and more important for every region (tropical, boreal and temperate) with an increase in forest disturbances (fires, diseases, storms) [1,4] and long-term changes in the geographical distribution of forest species [1,5,6].
In France, as in western Europe, forests are one of the key elements in ecological transition policies [7][8][9]. Sustainable forest management is needed in order to optimise wood harvesting, while preserving carbon storage, biodiversity and other ecological functions of the forests [1, 7,9]. Management techniques include timber extraction, planting and replanting of different species, and preventing and reducing the impacts of fires, diseases and storms.
For sustainable forest management and climate change mitigation measures, quantitative and dynamic information on 3D forest structure and biomass are needed by forest managers (local scale) and public bodies (local and regional scale). Furthermore, scientists need the information to understand the impacts of climate on forest health and growth in order to guide public policies for adaptation and mitigation measures.
Most countries in temperate regions conduct national forest inventories (NFI) on a five-year time scale, based on forest measurement plots sampled throughout the territory. The inventory database includes the following measurements (which are referred to in this paper as forest parameters): diameter at breast height (DBH), number of trees per hectare (tree number density or tree density), basal area (BA) and tree height. Known species-specific allometry equations allow the estimation of timber volume (i.e., wood resources) or above-ground biomass (AGB) (i.e., carbon storage) from these measurements. The inventory data are used to derive statistics of forest resources over large administrative regions. However, the forest inventory is not designed to produce geospatial distribution or high-resolution maps that are required by forest managers to carry out precise estimates at a local scale. For that purpose, forest managers usually conduct forest stands measurements themselves, that in most cases lack standardisation (due to estimation methods, spatial sampling and surveyor bias), require labor, and are less regularly updated than the NFIs performed by national bodies.
To produce forest parameter maps and to overcome the limitations of survey approaches, research studies have been carried out to assess the use of remote sensing data. The use of satellite images has been evaluated as a way to build operational and costeffective methods for mapping forest parameters over large areas. The focus is put on the existing open access images that are more and more diversified (radar, optical, high spatial and temporal resolution) and on the use of machine learning approaches [10][11][12][13][14]. Several satellite programs currently provide global and open access high resolution images. In the optical domain, Sentinel-2 and Landsat-8 images are acquired, respectively, every 5 and 16 days, with a spatial resolution of 10-20 m and 30 m for visible, near infrared and short-wave infrared bands (vis, nir, swir). Optical reflectance can be related to the green leaf area of the canopy and its horizontal distribution [15][16][17]. Studies have highlighted relationships between reflectance and forest parameters [18][19][20], although the signal quickly saturates [21][22][23]. Moreover, spatial texture indicators from optical sensors have been used with successful results [18,22,[24][25][26], and a few studies have shown that the temporal dimension can also improve forest parameter estimations [27]. In the radar domain, Sentinel-1 carries a synthetic aperture radar (SAR), that provides data in C-band (wavelength = 5.6 cm) every 6 days at~20 m spatial resolution, and annual global mosaics of L-band ALOS PALSAR (wavelength = 22.9 cm) imagery at~25 m spatial resolution are also available. The SAR signal penetrates into the forest canopies, with deeper penetration at longer wavelengths. The backscatter energy is a function of the geometric and dielectric properties of the volume of scatterers in a forest, which interacts with the incoming waves. Indirectly, SAR backscatter could be related to the water content of the scattering volume, which in turn is related to biomass, as highlighted in large scale studies [28][29][30][31][32][33]. However, for application at the local scale, relationships between SAR backscatter and forest biomass could be disturbed by other types of parameters (forest structure, soil and vegetation moisture, understory rebound, etc.), which need to be accounted for. A growing number Remote Sens. 2022, 14, 2079 3 of 25 of studies have shown that the combination of optical data, SAR data and textural data improves forest parameter estimations [10,14,[34][35][36][37][38]. Due to the diversity of remote sensing data and forests, most approaches use machine learning methods. In order to map large areas of forests with a diversity of tree species and management practices, these methods require statistical sampling to have a training dataset representative of forests [39][40][41]. However, field campaigns are very expensive and cannot be multiplied to represent forests with diverse and heterogeneous conditions, as found in Europe. For example, in France there are 11 large ecological regions (called GRECO) that represent the different topological and climate conditions and imply different silvicultural practices and forest growth. Because of the diversity of forest management practices (monospecific, mixed species with different compositions and age classes, different planting, thinning practices, etc.), studies often lack a sufficient volume and representativeness of data. Results rely on the forest region and forest type where training samples are measured. The main limitations of these methods are therefore the quality, quantity and representativeness of the reference field data.
Another approach to complement or replace ground survey is the use of light detection and ranging (LiDAR) methods. Airborne LiDAR (ALS) data provide highly accurate canopy height estimates, in addition to providing information on tree shape and understory. ALS data are the most widely used data operationally to assist in forest inventories [42,43]. Several studies used LiDAR data to estimate AGB, forest height and other forest parameters [42,[44][45][46][47][48][49], and showed improved results compared to photogrammetry, radargrammetry, and satellite optical and SAR images. Although ALS data are very accurate and useful to complement or replace ground survey for forest inventories, these data remain expensive. ALS campaigns are therefore often either spatially or temporally (repetitiveness) limited. Spaceborne LiDAR has also been used to improve the estimation of forest height. The Geoscience Laser Altimeter System (GLAS, ICESat satellite mission) has been used to improve the estimation of forest height [50][51][52][53][54][55][56][57]. GLAS has provided transects of footprints between 2003 and 2009. The footprints have a nearly circular shape of~60 m diameter and a footprint spacing of~170 m along the track. More recently, the Global Ecosystem Dynamics Investigation mission (GEDI, full waveform LiDAR on board the ISS) has provided an unprecedented number of forest canopy height samples from 2019 and its mission is scheduled until beginning of 2023. The size of GEDI footprints (~25 m diameter) is smaller than ICESat, providing improved measurements over forest areas with relief higher than 10%, and is closer to the spatial resolution of satellite images from the Sentinel and ALOS PALSAR programs. GEDI forest height metrics have already been studied to estimate forest parameters [58][59][60][61][62].
The downside of the GLAS and GEDI is that the data are footprints (15-60 m radius) spaced along transects; thus, they are not spatially continuous (i.e., not a map with data on every forest pixel). On the other hand, the large amount of LiDAR footprints could be used as forest height reference data. Studies have shown that the combination of satellite images with spaceborne LiDAR data could be used to improve the production of forest height or AGB maps at national or continental scale. GLAS data have been used in combination with optical images such as MODIS [63,64] or Landsat [65][66][67][68][69][70]. The newly available GEDI data have already been used to map forest structure [71][72][73]. A global forest canopy height map at 30 m spatial resolution was produced [72] and is freely available [74] (hereafter called GLAD-2019). GEDI data offer an opportunity to improve representativeness in methods using multi-sensor satellite images to estimate and map forest height.
We built a processing chain [37,75] for the estimation and mapping of forest structure parameters using field measurements and freely available spaceborne images (Sentinel-1, Sentinel-2 and ALOS-2 PALSAR-2 annual mosaics) in supervised learning techniques. The chain was first tested on a case study with single species and even-aged maritime pine (Pinus pinaster) plantation forests, for the estimation of AGB, BA, DBH, mean height (Hmean), dominant height (Hdom) and density. In this study, we assess the method on mixed broadleaved and coniferous forests in France, and we test the possibility of integrating GEDI height metrics as reference data for forest dominant height estimation Remote Sens. 2022, 14, 2079 4 of 25 (i.e., replacing field measurements by GEDI height metrics in model learning). In the pine forest plantations, Hdom is roughly similar to Hmean since all the trees are of the same age and species. However, in uneven-aged forests and with several species, the height of trees is often more difficult to estimate in the field (lack of a direct view of the top of the canopy due to more vertical complexity), but also the dominant or average height no longer mean the same things for forest management. The dominant height has more significance in these uneven-aged forests since it takes into account the largest trees, which are those which have the most value for logging. From the perspective of satellite data, it is also harder to build models for estimating the forest parameters of uneven-aged and mixed forests. Satellite remote sensing of forests is partly based on allometric relationships, but these relationships change from one species to another and from one forest management activity to another. In addition, in pixel-based mapping procedures, the allometric equations need to account for the mixture of species and age classes in a pixel. Working with machine learning, the most direct way to improve Hdom estimates on these highly diverse forests is to increase the number and representativeness of reference samples. Since field measurements are time-consuming and not always representative, the GEDI data could constitute a large reference dataset that can help to overcome the lack of sufficient reference data to cover the European temperate forest diversity.
In this study, we assess the use of GEDI height metrics as training data for the estimation of forest dominant height (Hdom) with multi-sensor satellite imagery. We first conduct the estimation of forest parameters (BA, DBH, Density, Hdom) using field measurements for training and validation (Section 3.2, Annex A). Secondly, we conduct the estimation of Hdom using GEDI RH95% for training and validation (Section 3.3.2). Finally, we conduct the estimation of Hdom using GEDI RH95% for training and field measurement for validation (Section 3.3.2). Results are compared to the existing canopy height map, GLAD-2019 (Section 3.1).

Study Sites and In Situ Data
The study was conducted over two test sites in northern France (Figure 1). Both sites contain public forests managed by the Office National des Forêts (ONF). The Orléans forest contains coniferous high stands and broadleaved uneven-aged high stands. Both can contain understory and coppice, but coniferous stands are generally more regular. Predominant species are oaks (Quercus robur and Quercus petraea) and Scots pines (Pinus sylvestris), with about ten other secondary species. The climate is temperate with~640 mm annual rain and a mean annual temperature of 11.2 • C. The soil is poor and acid, and the topography is flat. The Saint-Gobain Forest (hereafter called StGobain) contains mainly broadleaved uneven-aged high stands. Predominant species are oaks (Quercus robur and Quercus petraea), beech (Fagus sylvatica) and hornbeam (Carpinus betulus), with about 15 other secondary broadleaved species. The climate is also temperate with~700 mm annual rain and a mean annual temperature of 10.4 • C. The soil is more nutrient-rich than in Orléans, so the forests are denser and more productive. The topography is generally flat, but some slopes can reach 20 • locally. Data collection was conducted by the ONF in early 2015 in StGobain and during the 2016-2017 winter in the Orléans forests. Tree count and diameter at breast height (DBH) measurements were conducted on circular plots from 15 to 17 m radius; only trees with a DBH larger than 17.5 cm were considered. On part of the plots, the height of the largest trees was measured to calculate the dominant height (height of the 100 largest trees per hectare). Additionally, airborne LiDAR (ALS) data were acquired on each site. LiDAR parameters for Orléans and StGobain are, respectively, 37 and 18 points/m 2 scanning density; 6 and 15 cm vertical accuracy; and 20 and 30 cm absolute positional accuracy. ALS data enabled a canopy height model (CHM) to be obtained at 1 m spatial resolution. The dominant height of all the samples was estimated by the ONF at plot level with Computree software (v5.0) [76,77] using a crown segmentation approach on the CHM, and then, the density; 6 and 15 cm vertical accuracy; and 20 and 30 cm absolute positional accuracy. ALS data enabled a canopy height model (CHM) to be obtained at 1 m spatial resolution. The dominant height of all the samples was estimated by the ONF at plot level with Computree software (v5.0) [76,77] using a crown segmentation approach on the CHM, and then, the heights of the 6 highest trees in the plot were extracted to obtain the dominant height. Results have been validated and corrected by the ONF with the samples where in situ height measurements were made. In the end, at plot level, the values of the basal area (BA), mean diameter at breast height (DBH), tree density (Density) and dominant height (Hdom) are available. We selected two classes for each site in order to work on consistent sets of samples (i.e., similar allometric relationships within each class). The classes are based on the dominant species that can be found in the forest cover cartographic database. At the Orléans test site, the "oak" class contains 75 plots with at least 75% of the basal area being oaks; secondary species are mainly hornbeams and pines. The "Scots pine" class contains 28 plots with at least 75% of the basal area being Scots pines; secondary species are mainly oaks and other pines. At the StGobain test site, the "oak" class contains 21 plots (still 75% basal area threshold); the secondary species is manly hornbeam. The "beech" class contains 25 plots (50% basal area threshold); secondary species are mainly oaks, hornbeams and ash (Fraxinus excelsior). Field inventory values are summarised in Table 1. The forest structure statistics for each class are broadly similar except for the dominant height which is higher at the StGobain site. We selected two classes for each site in order to work on consistent sets of samples (i.e., similar allometric relationships within each class). The classes are based on the dominant species that can be found in the forest cover cartographic database. At the Orléans test site, the "oak" class contains 75 plots with at least 75% of the basal area being oaks; secondary species are mainly hornbeams and pines. The "Scots pine" class contains 28 plots with at least 75% of the basal area being Scots pines; secondary species are mainly oaks and other pines. At the StGobain test site, the "oak" class contains 21 plots (still 75% basal area threshold); the secondary species is manly hornbeam. The "beech" class contains 25 plots (50% basal area threshold); secondary species are mainly oaks, hornbeams and ash (Fraxinus excelsior). Field inventory values are summarised in Table 1. The forest structure statistics for each class are broadly similar except for the dominant height which is higher at the StGobain site. Sentinel-2 images: Sentinel-2 data consist of optical time series acquired in the visible, near infrared (nir) and short-wave infrared (swir), with a spatial resolution of 10 to 20 m. The raw images are provided by the European Space Agency (ESA). The Theia Land Data Center [78] produces the level 2A images (top of canopy reflectance: TOC) with orthorectification, atmospheric corrections and cloud mask, using MAJA processing [79,80].
where B2, B3, B4 and B8 are, respectively, blue, green, red and nir spectral bands at 10 m spatial resolution; B8A and B11 are, respectively, nir and swir at 20 m spatial resolution. NDWI is then resampled to 10 m spatial resolution. These indices from winter and summer acquisitions were used in the following processing. Sentinel-1: C-band SAR time series resampled at 10 m spatial resolution are provided by the ESA and downloaded from the Sentinel Product Exploitation Platform [81] (PEPS). We used the Level-1 Ground Range Detected (GRD) products on the descending orbit (DES, morning acquisition) with~33 • incidence angle (IA) on the test sites. The 2016 Sentinel-1 time series were calibrated, orthorectified and corrected from slope effect, and then, we applied multi-image filtering [82,83] on the 2016 time series, leading to an equivalent number of looks (ENL) of~92 for each acquisition. Then, average images for VH and VV polarisations were calculated using all the dates included in an interval of 2 months, centred on the Sentinel-2 winter and summer dates. VH and VV winter and summer averages were used in the following processing.
ALOS-2 PALSAR-2: L-band SAR annual mosaics at 25 m spatial resolution are produced by the Japan Aerospace Exploration Agency (JAXA) through large scale mosaicking [84] that includes orthorectification, slope correction and radiometric calibration between neighbouring strips. Five annual mosaics were sourced from JAXA for the years 2015 through 2019. We converted digital numbers (DN) to gamma naught values. Then, we applied multi-image filtering [82,83] using the 2 polarisations and the 5 years, leading to an ENL of 135 for each image. In the following, we used resulting HV and HH images from 2016.
Textural indices: spatial texture metrics have been generated from both Sentinel-2 10 m spatial resolution indices (NDVI, BI) and Sentinel-1 polarisations (VH, VV), using a Grey Level Co-occurrence Matrix (GLCM) with the Orfeo Toolbox library (OTB, [85]). Three second order Haralick metrics were calculated on GLCM: homogeneity (also called Inverse Difference Moment), contrast (also called Inertia), and Haralick correlation. Based on a previous sensitivity analysis [75], the following parameters were used to obtain the textural metrics for further processing: offset set to 1 pixel; 3 pixels radius; 100 grey levels (number of bins); minimum and maximum calculated separately between Sentinel-2 indices and dates, and between Sentinel-1 polarisations; and orientations 0 • , 45 • , 90 • and 135 • were computed then averaged.
Final set of features: after the processing steps explained above, we obtained a total of 36 features: 16 Sentinel-1 features (VH and VV for winter and summer dates and derived textural metrics), 18 Sentinel-2 features (BI, NDVI and NDWI for winter and summer dates, BI and NDVI textural metrics), and 2 ALOS-2 PALSAR-2 metrics (HV and HH polarisations).

GEDI Data
The Global Ecosystem Dynamics Investigation (GEDI) LiDAR is a laser instrument on board the International Space Station. It has three lasers; two of the lasers are full power and the other one is split into two beams (called coverage beams), producing a total of four beams and eight ground tracks (four power and four cover tracks). Footprints average 25 m in diameter and they are separated by 60 m along track and 600 across track.
The GEDI L2A data from May 2019 to August 2020 were downloaded [86] and preprocessed as presented in [71]. On this data, the first settings group is available to derive relative height metrics (RH): Smoothing Width (Noise) = 6.5σ; Smoothing Width (Signal) = 6.5σ; Waveform Signal Start Threshold = 3σ; Waveform Signal End Threshold = 6σ. σ represents the standard deviation of the background noise level. Relative height metrics (RH n ) were available with n from 0 to 100% (with a step of 1%). In order to select the best RH n , we compared different GEDI RH percentiles to the airborne LiDAR (ALS) data available on the two test sites. Based on the results (see Section 3.2) and previous studies, [59,72,87] we chose to use RH95 as the GEDI dominant height estimates. We also tested the relevance of using all the beams or only the full power beams.
After filtering the data as performed in [71] (heights not consistent with tree height), we considered different options to reduce GEDI height metric errors. We kept only waveforms where at least 2 modes were detected to ensure a top of canopy and soil return signal. Finally, we tested to remove the noisiest waveforms, keeping only the 25, 50 or 75% waveforms with the best signal-noise ratio.

Processing of Ancillary Products
The forest species on the GEDI samples were unknown; therefore, we used a combination of two datasets for tree species classification: the OSO land cover map (website: https://www.theia-land.fr/en/product/land-cover-map/, accessed on 16 March 2022) and BDForêtV2 (website: https://geoservices.ign.fr, accessed on 16 March 2022).
OSO: annual land cover maps are available for France on the Theia platform and produced with IOTA2 software (website: https://framagit.org/iota2-project/iota2, accessed on 16 March 2022); the main stages are described in [88]. The nomenclature includes two classes of forests: broadleaved and coniferous. This nomenclature is not sufficient for our study, but this map offers a high spatial accuracy (10 × 10 m pixels) and up-to-date land cover information (annual product). The 2018 OSO map is used.
BDForetV2: the French national institute of geographic and forest information (IGN) has produced forest species maps using airborne optical images from 2004 to 2016, depend- ing on localisation. Over the two sites, Orléans and StGobain, forest maps were produced, respectively, with 2010 and 2016 images. The nomenclature provides knowledge of the location of the pure species (>75% of the surface is occupied by a given species), which is enough to assign tree species classes to GEDI footprints. However, the spatial accuracy, the updating and the minimum mapping unit are somewhat limited for our use.
Forest maps combination: for each GEDI footprint, we assigned a species (or group of species) class when the OSO and BDForetV2 agreed (e.g., OSO indicates broadleaved and BDForetV2 indicates oak). We kept only footprints with pure species.
In the end, we had 826 GEDI oak dominant footprints and 1701 Scots pine dominant footprints in Orléans, and 1195 oak dominant footprints and 1090 beech dominant footprints in StGobain.

Estimation of Forest Parameters
Our processing chain is based on a machine learning approach; the full methodology is described in detail in [37,75]. An example with our number of features and a forward approach is schematised in Figure 2. The key points are:

•
ONF Hdom measurements or GEDI RH95% are used as reference data, and the remote sensing features from Section 2.3 are used as predictive variables.

•
Reference data and remote sensing features are fed to a Support Vector Machine (SVR) algorithm for regression; the SVR has 3 key parameters that need to be optimised: cost parameter, gamma and epsilon.

•
Validation is made with a leave-one-out (LOO) cross-validation method when using the ONF measurement data because the number of samples is very low for these data (21 to 75 samples); when the reference data is GEDI RH95, we use a stratified K-fold cross-validation method.

•
We used a stepwise approach for feature selection. It is an iterative process that starts with 1 feature or N features and iteratively adds or removes features, respectively, for forward and backward selection. Within the 1 to N feature combinations, the one with the lowest RMSE is considered the best predictive set we can obtain on the dataset. BDForetV2: the French national institute of geographic and forest information (IGN) has produced forest species maps using airborne optical images from 2004 to 2016, depending on localisation. Over the two sites, Orléans and StGobain, forest maps were produced, respectively, with 2010 and 2016 images. The nomenclature provides knowledge of the location of the pure species (>75% of the surface is occupied by a given species), which is enough to assign tree species classes to GEDI footprints. However, the spatial accuracy, the updating and the minimum mapping unit are somewhat limited for our use.
Forest maps combination: for each GEDI footprint, we assigned a species (or group of species) class when the OSO and BDForetV2 agreed (e.g., OSO indicates broadleaved and BDForetV2 indicates oak). We kept only footprints with pure species.
In the end, we had 826 GEDI oak dominant footprints and 1701 Scots pine dominant footprints in Orléans, and 1195 oak dominant footprints and 1090 beech dominant footprints in StGobain.

Estimation of Forest Parameters
Our processing chain is based on a machine learning approach; the full methodology is described in detail in [37,75]. An example with our number of features and a forward approach is schematised in Figure 2. The key points are:

•
ONF Hdom measurements or GEDI RH95% are used as reference data, and the remote sensing features from Section 2.3 are used as predictive variables.

•
Reference data and remote sensing features are fed to a Support Vector Machine (SVR) algorithm for regression; the SVR has 3 key parameters that need to be optimised: cost parameter, gamma and epsilon.

•
Validation is made with a leave-one-out (LOO) cross-validation method when using the ONF measurement data because the number of samples is very low for these data (21 to 75 samples); when the reference data is GEDI RH95, we use a stratified K-fold cross-validation method.

•
We used a stepwise approach for feature selection. It is an iterative process that starts with 1 feature or N features and iteratively adds or removes features, respectively, for forward and backward selection. Within the 1 to N feature combinations, the one with the lowest RMSE is considered the best predictive set we can obtain on the dataset. In this study, three ways to estimate forest parameters from satellite imagery are tested. First, is the estimation of BA, DBH, Density and Hdom using ONF field In this study, three ways to estimate forest parameters from satellite imagery are tested. First, is the estimation of BA, DBH, Density and Hdom using ONF field measurement for learning and validation (LOO cross-validation). Second, is the estimation of Hdom using GEDI RH95 for learning and validation (K-fold cross-validation). Finally, the estimation of Hdom using GEDI RH95 for learning and ONF field measurements for validation is conducted.
The statistics used to evaluate performances are Pearson's correlation (r, r 2 and p0), RMSE, relative RMSE, Mean Absolute Error (MAE) and relative MAE.
where N is the number of samples. Relative scores (in percentages) are these statistics divided by the mean of reference data values. RMSE gives more weight to outliers compared to MAE.

GLAD-2019 Canopy Height Map Evaluation
A global map of forest canopy height at 30 m spatial resolution was produced by [72], based on GEDI RH95% and Landsat-8 2019 times series. We evaluated the accuracy of this map over our test sites in temperate forests (France). Figure 3 shows the scatterplot between the ONF field measurements (x-axis) based on dominant trees and the GLAD-2019 forest canopy height map (y-axis) over the four forest classes. The scatterplots show poor agreement between the two. The dynamic range of the GLAD-2019 values are restricted compared with the field data in the StGobain forest, and multiple outliers are found in the Orléans test site. Several possible causes can explain the poor agreement for these mixed forests (to silvicultural interventions at Orléans or to atypical forest stands, etc.). More importantly, the GLAD-2019 height estimates were not able to retrieve the forest height dynamic range of the reference data, as indicated by the low correlation coefficients. For the Orléans oak class, reference values range from 14 to 30 m, while the GLAD-2019 height values range from~19 to~25 m. The dynamics of the GLAD-2019 values are even more restricted for other classes compared to reference value dynamics. This is due both to the type of forests that we observed which are mixed (different ages, diameters and species within the plots) and spatially complex, and to the satellite image data used for height map production: Landsat-8 time series are optical data; several studies showed that they have a significant but limited relationship with forest parameters compared to radar and textural metrics [19,22,27]. measurement for learning and validation (LOO cross-validation). Second, is the estimation of Hdom using GEDI RH95 for learning and validation (K-fold cross-validation). Finally, the estimation of Hdom using GEDI RH95 for learning and ONF field measurements for validation is conducted.
The statistics used to evaluate performances are Pearson's correlation (r, r 2 and p0), RMSE, relative RMSE, Mean Absolute Error (MAE) and relative MAE.
where N is the number of samples. Relative scores (in percentages) are these statistics divided by the mean of reference data values. RMSE gives more weight to outliers compared to MAE.

GLAD-2019 Canopy Height Map Evaluation
A global map of forest canopy height at 30 m spatial resolution was produced by [72], based on GEDI RH95% and Landsat-8 2019 times series. We evaluated the accuracy of this map over our test sites in temperate forests (France). Figure 3 shows the scatterplot between the ONF field measurements (x-axis) based on dominant trees and the GLAD-2019 forest canopy height map (y-axis) over the four forest classes. The scatterplots show poor agreement between the two. The dynamic range of the GLAD-2019 values are restricted compared with the field data in the StGobain forest, and multiple outliers are found in the Orléans test site. Several possible causes can explain the poor agreement for these mixed forests (to silvicultural interventions at Orléans or to atypical forest stands, etc.). More importantly, the GLAD-2019 height estimates were not able to retrieve the forest height dynamic range of the reference data, as indicated by the low correlation coefficients. For the Orléans oak class, reference values range from 14 to 30 m, while the GLAD-2019 height values range from ~19 to ~25 m. The dynamics of the GLAD-2019 values are even more restricted for other classes compared to reference value dynamics. This is due both to the type of forests that we observed which are mixed (different ages, diameters and species within the plots) and spatially complex, and to the satellite image data used for height map production: Landsat-8 time series are optical data; several studies showed that they have a significant but limited relationship with forest parameters compared to radar and textural metrics [19,22,27].

Dominant Height Estimation Using Field Campaign as Reference Data
The processing chain described in Section 2.3.2 and [37] was used to estimate the dominant height (Hdom), basal area (BA), mean diameter at breast height (DBH) and tree density (Density) on the ONF dataset for the forest classes in the Orléans and StGobain test sites. In this section, we focus on the Hdom results. BA, DBH and Density estimations are presented in Appendix A ( Figures A1 and A2, Table A1) for informational purposes: the overall results on the broadleaved and Scots pines forests are good, with RMSE from 11% for dominant height to 30% for tree Density. The results by class show slightly less accurate predictions on the Orleans oaks for all the forest parameter estimates; the other classes have better and fairly similar results.
Using only in situ data (ONF field measurements for model training and validation with a LOO cross-validation), a forward stepwise feature selection process was applied, since the number of reference samples was low compared to the number of input features in the learning model. Figure 4 presents the evolution of the dominant height validation statistics (r 2 ) during the forward feature selection process. The r 2 first increases for all the four forest classes until it reaches a maximum (i.e., best combination) of between 5 and 15 features. The features selected are different for each class. At least two textural indices from Sentinel-1 and two others from Sentinel-2 are selected in the optimal combination, even if they differ between classes (sensor indices, textural indices or acquisition date). Spectral indices and L-or C-band polarisations are used even more variably. After reaching the maximum, precision decreases and shows poor estimation quality with all features in the model. This may be due to the low number of samples compared to the maximum number of features; there is over-fitting with 36 features and therefore the cross-validation results show large errors. This result shows the importance of feature selection when using small numbers of samples; forward stepwise selection allows very good results to be obtained where the total feature set causes over-fitting. However, the selection process can also cause excessive method specialisation and therefore difficulties in generalising the results. Figure 5 presents the validation scatterplots of the four classes for dominant height estimation with the ONF samples as references for training and testing using a LOO crossvalidation. Good results are obtained for StGobain and for Scots pine in Orléans with 0.74-0.78 r 2 and low errors (1.9-3.1 m RMSE, 7.3-11.4% relative RMSE). For Orléans oak forest, the relative RMSE is still satisfactory (2.6 RMSE, 11.6% relative RMSE), but the r 2 (0.48) is not as good as for the other classes. There might be a saturation effect; dominant height beyond 25 m is underestimated.

Dominant Height Estimation Using Field Campaign as Reference Data
The processing chain described in Section 2.3.2 and [37] was used to estimate the dominant height (Hdom), basal area (BA), mean diameter at breast height (DBH) and tree density (Density) on the ONF dataset for the forest classes in the Orléans and StGobain test sites. In this section, we focus on the Hdom results. BA, DBH and Density estimations are presented in Appendix A ( Figures A1 and A2, Table A1) for informational purposes: the overall results on the broadleaved and Scots pines forests are good, with RMSE from 11% for dominant height to 30% for tree Density. The results by class show slightly less accurate predictions on the Orleans oaks for all the forest parameter estimates; the other classes have better and fairly similar results.
Using only in situ data (ONF field measurements for model training and validation with a LOO cross-validation), a forward stepwise feature selection process was applied, since the number of reference samples was low compared to the number of input features in the learning model. Figure 4 presents the evolution of the dominant height validation statistics (r 2 ) during the forward feature selection process. The r 2 first increases for all the four forest classes until it reaches a maximum (i.e., best combination) of between 5 and 15 features. The features selected are different for each class. At least two textural indices from Sentinel-1 and two others from Sentinel-2 are selected in the optimal combination, even if they differ between classes (sensor indices, textural indices or acquisition date). Spectral indices and L-or C-band polarisations are used even more variably. After reaching the maximum, precision decreases and shows poor estimation quality with all features in the model. This may be due to the low number of samples compared to the maximum number of features; there is over-fitting with 36 features and therefore the cross-validation results show large errors. This result shows the importance of feature selection when using small numbers of samples; forward stepwise selection allows very good results to be obtained where the total feature set causes over-fitting. However, the selection process can also cause excessive method specialisation and therefore difficulties in generalising the results. Figure 5 presents the validation scatterplots of the four classes for dominant height estimation with the ONF samples as references for training and testing using a LOO crossvalidation. Good results are obtained for StGobain and for Scots pine in Orléans with 0.74-0.78 r 2 and low errors (1.9-3.1 m RMSE, 7.3-11.4% relative RMSE). For Orléans oak forest, the relative RMSE is still satisfactory (2.6 RMSE, 11.6% relative RMSE), but the r 2 (0.48) is not as good as for the other classes. There might be a saturation effect; dominant height beyond 25 m is underestimated.

Dominant Height Estimation Using GEDI Height Metrics as Reference Data
As it has been discussed, the processing chain is able to estimate forest dominant height from a set of satellite features. But the main limitation in the performance results is the small number of samples for model learning. We are thus interested in the use of GEDI data to increase the number of reference samples for the height predictor. In this section, we first present the results of our analysis of the GEDI height metrics compared to ALS data. Secondly, we present the application of our processing chain to estimate dominant height using GEDI as reference data and the satellite image features.

GEDI RH95% Comparison with ALS Data
GEDI RH percentiles were compared to tree height metrics extracted from the ALS canopy height model (CHM) at 1 m spatial resolution. CHM metrics were calculated on the pixels inside the footprints of the GEDI sample (circles 25 m in diameter). Among different CHM metrics, we selected the 95% percentile (H95) as a proxy for dominant height [89,90], and compared it to GEDI RH metrics. The analysis of GEDI RHn with n from 50 to 100% (steps (50,75,85,90,95,99,100)) revealed that RH95% was the best metric overall (according to r 2 and RMSEs of the four classes studied). Figure 6 shows the comparison between CHM-H95 (x-axis) and GEDI RH95% (yaxis) for the four classes. The best result was obtained on the only coniferous class (Scots pine), with r 2 equal to 0.65 and 3.7 m RMSE. Broadleaved classes have lower coefficients of determination of around 0.48. In France, coniferous forests are generally managed as uniform high stands and thus easier to characterise. We can observe outliers overestimated in the low range and underestimated in the high range. This can be caused partly by growth or silvicultural interventions in the period between the ALS and GEDI acquisitions, or the acquisition date (GEDI winter dates can underestimate height due to leaf loss), or by not enough filtering of GEDI footprints. For the purposes of the study and to use the GEDI samples in the best conditions, outliers were removed before the application of our processing chain for forest parameter estimation with satellite image features. The threshold was set at 5 m away from the linear regression line between CHM-H95 and GEDI RH95% for each class. In the end, we had 479 GEDI oak dominant footprints and 958 Scots pine dominant footprints in Orléans, and 808 oak dominant footprints and 723 beech dominant footprints in StGobain. Apart from outliers, there is generally a bias or saturation for high values, more visible in StGobain because the heights are higher.

Dominant Height Estimation Using GEDI Height Metrics as Reference Data
As it has been discussed, the processing chain is able to estimate forest dominant height from a set of satellite features. But the main limitation in the performance results is the small number of samples for model learning. We are thus interested in the use of GEDI data to increase the number of reference samples for the height predictor. In this section, we first present the results of our analysis of the GEDI height metrics compared to ALS data. Secondly, we present the application of our processing chain to estimate dominant height using GEDI as reference data and the satellite image features.

GEDI RH95% Comparison with ALS Data
GEDI RH percentiles were compared to tree height metrics extracted from the ALS canopy height model (CHM) at 1 m spatial resolution. CHM metrics were calculated on the pixels inside the footprints of the GEDI sample (circles 25 m in diameter). Among different CHM metrics, we selected the 95% percentile (H95) as a proxy for dominant height [89,90], and compared it to GEDI RH metrics. The analysis of GEDI RHn with n from 50 to 100% (steps (50,75,85,90,95,99,100)) revealed that RH95% was the best metric overall (according to r 2 and RMSEs of the four classes studied). Figure 6 shows the comparison between CHM-H95 (x-axis) and GEDI RH95% (yaxis) for the four classes. The best result was obtained on the only coniferous class (Scots pine), with r 2 equal to 0.65 and 3.7 m RMSE. Broadleaved classes have lower coefficients of determination of around 0.48. In France, coniferous forests are generally managed as uniform high stands and thus easier to characterise. We can observe outliers overestimated in the low range and underestimated in the high range. This can be caused partly by growth or silvicultural interventions in the period between the ALS and GEDI acquisitions, or the acquisition date (GEDI winter dates can underestimate height due to leaf loss), or by not enough filtering of GEDI footprints. For the purposes of the study and to use the GEDI samples in the best conditions, outliers were removed before the application of our processing chain for forest parameter estimation with satellite image features. The threshold was set at 5 m away from the linear regression line between CHM-H95 and GEDI RH95% for each class. In the end, we had 479 GEDI oak dominant footprints and 958 Scots pine dominant footprints in Orléans, and 808 oak dominant footprints and 723 beech dominant footprints in StGobain. Apart from outliers, there is generally a bias or saturation for high values, more visible in StGobain because the heights are higher.

Using GEDI RH95% Instead of Field Data
In this section, we present the results of the processing chain presented in Section 2.3.2 to estimate the dominant height using GEDI RH95% as reference samples in Orléans and StGobain test sites (GEDI RH95% for model training and validation using a K-fold cross-validation). Figure 7 shows the evolution of r 2 on the four forest classes with the feature selection process. The r 2 is stable from 36 features (maximum) to ~15 features. Having a large number of samples makes it possible to have "unnecessary" features in the model without affecting the results. Among these ~15 features, the variability between classes is lower than when working with the small ONF datasets (all feature types are always present in the 15 features). However, differences remain when looking at precisely which features are selected for each class. When this threshold of 15 features is exceeded, the precision starts to fall faster and faster when there are a very few numbers of features in the model. The models no longer have enough information to allow a good estimation of the dominant height. The feature selection step might not be necessary when working with GEDI samples. Moreover, directly using all the features without specific selection for each class is better in order to be more easily generalisable. Indeed, working by class implies having cartographic information on these classes to assign them to the GEDI footprints, then to produce maps with the prediction models. This cartographic information is not always accurate, it is therefore preferable to have identical feature sets between classes to avoid having artefacts when there are species errors. Figures 8-11 present the scatterplots of the estimation of Hdom using one model per forest class (algorithm parametrisation and model learning are made separately for each forest class). The graphs on the left side show GEDI RH95% cross-validation results, meaning the ability of satellite images to retrieve GEDI height metrics. The graphs on the right side show the application of GEDI RH95% based models on ONF measurement plots, meaning an independent validation of the dominant height predictions using GEDI data and satellite image features. GEDI RH95% K-fold cross-validation results are roughly similar for the four classes, just slightly worse on Orléans oak. Statistics are good, showing r 2 from 0.54 to 0.65, and RMSEs from 2.7 to 3.4 m (~12.9 to 16.7%). Despite the results being scattered, there is no significant bias and the satellite image features are able to retrieve the range of GEDI RH95% value dynamics. Learning on GEDI data and testing on the ONF measurement plots also provide good overall results (right graph on Figures 8-11). Validation statistics are however more variable, partly because of the lower numbers of samples and, less importantly, the larger variability in the ONF value dynamics between

Using GEDI RH95% Instead of Field Data
In this section, we present the results of the processing chain presented in Section 2.3.2 to estimate the dominant height using GEDI RH95% as reference samples in Orléans and StGobain test sites (GEDI RH95% for model training and validation using a K-fold crossvalidation). Figure 7 shows the evolution of r 2 on the four forest classes with the feature selection process. The r 2 is stable from 36 features (maximum) to~15 features. Having a large number of samples makes it possible to have "unnecessary" features in the model without affecting the results. Among these~15 features, the variability between classes is lower than when working with the small ONF datasets (all feature types are always present in the 15 features). However, differences remain when looking at precisely which features are selected for each class. When this threshold of 15 features is exceeded, the precision starts to fall faster and faster when there are a very few numbers of features in the model. The models no longer have enough information to allow a good estimation of the dominant height. The feature selection step might not be necessary when working with GEDI samples. Moreover, directly using all the features without specific selection for each class is better in order to be more easily generalisable. Indeed, working by class implies having cartographic information on these classes to assign them to the GEDI footprints, then to produce maps with the prediction models. This cartographic information is not always accurate, it is therefore preferable to have identical feature sets between classes to avoid having artefacts when there are species errors. Figures 8-11 present the scatterplots of the estimation of Hdom using one model per forest class (algorithm parametrisation and model learning are made separately for each forest class). The graphs on the left side show GEDI RH95% cross-validation results, meaning the ability of satellite images to retrieve GEDI height metrics. The graphs on the right side show the application of GEDI RH95% based models on ONF measurement plots, meaning an independent validation of the dominant height predictions using GEDI data and satellite image features. GEDI RH95% K-fold cross-validation results are roughly similar for the four classes, just slightly worse on Orléans oak. Statistics are good, showing r 2 from 0.54 to 0.65, and RMSEs from 2.7 to 3.4 m (~12.9 to 16.7%). Despite the results being scattered, there is no significant bias and the satellite image features are able to retrieve the range of GEDI RH95% value dynamics. Learning on GEDI data and testing on the ONF measurement plots also provide good overall results (right graph on Figures 8-11). Validation statistics are however more variable, partly because of the lower numbers of samples and, less importantly, the larger variability in the ONF value dynamics between classes. The Orléans oak class has the worst r 2 (0.23), but the error and bias are still satisfactory. The best result is on the StGobain oak class, with 0.64 r 2 and 12.8% relative RMSE.
classes. The Orléans oak class has the worst r 2 (0.23), but the error and bias are still satisfactory. The best result is on the StGobain oak class, with 0.64 r 2 and 12.8% relative RMSE.   classes. The Orléans oak class has the worst r 2 (0.23), but the error and bias are still satisfactory. The best result is on the StGobain oak class, with 0.64 r 2 and 12.8% relative RMSE.

Discussion
In a previous study, we built and validated a processing chain on forest coniferous plantations. In the study presented in this paper, our processing chain was able to accurately estimate forest parameters on more complex temperate semi-natural forests (uneven ages, mixed species, spatial variability). All the results obtained with the ONF data, GEDI data and the comparison with GLAD-2019 height map are summarised in Table 2.

Discussion
In a previous study, we built and validated a processing chain on forest coniferous plantations. In the study presented in this paper, our processing chain was able to accurately estimate forest parameters on more complex temperate semi-natural forests (uneven ages, mixed species, spatial variability). All the results obtained with the ONF data, GEDI data and the comparison with GLAD-2019 height map are summarised in Table 2.
The best results we obtained are relative RMSEs from 7.3% to 11.6% for dominant height estimation using the ONF measurements as a reference dataset in a LOO crossvalidation. The error rates obtained in this study with satellite imagery are similar to the results seen in the literature with photogrammetry and radargrammetry (6.6% to 14% relative RMSE depending on test sites), and slightly higher than the results seen with ALS data (4.6% to 7.8% relative RMSE) [43,47,49,91]. Thereby, the combination of optical and SAR sensors with textural indices at 10 m spatial resolution, in a forward stepwise feature selection process, allowed accurate models for different coniferous and broadleaved forests to be built. Satellite images used in this study are freely available within long-term programs, acquired globally and regularly at high spatial resolution (10-25 m). It is very promising to be able to use these data to estimate forest parameters that can be used for forest management or carbon modelling. However, there are limitations to generalising the method. Field data that meet the needs of mapping by satellite imagery (precise geolocation, the distribution of values, the size of the sampled surfaces, etc.) are often of limited number and concentrated in small areas. The accuracy and the ability to generalise the results are therefore limited by field measurement availability. In this study, with small measurement datasets, we could produce accurate estimations by using an automated feature selection process (forward stepwise selection) and showed that this selection was important to improve the results. However, a selection of features, whether automatic or manual (operator analysis), also increases the risk of being unable to obtain high accuracy on other types of forest or when dealing with large areas (regional or country scale), with forests not represented in the reference data. GEDI data availability therefore presents a huge opportunity for methods based on machine learning. The data acquisition will be carried out only until 2023, but the coverage and the large number of samples have already the potential to be integrated into learning approaches for the mapping of local forests worldwide. We have shown that the heights provided by GEDI are in good agreement with ALS heights. When GEDI RH95% is compared with ALS-CHM 95% ( Figure 6, Section 3.3.1), r 2~0 .48 and RMSE ranges from 5.1 to 5.6 m for the broadleaved classes, and r 2 is 0.65 and RMSE = 3.7 m for the coniferous class (Orléans Scots pines). These statistics are improved when removing outliers that may be due to time shift, silvicultural interventions and acquisition conditions (r 2 = 0.78 and RMSE = 4.5 m for the Orléans oak class, r 2 = 0.84-0.86 and RMSE = 2.8-2.9 m for the other classes). The lower accuracy on the Orléans oak class can be explained by a larger presence of coppice forests and open forests in this class, while other classes are more dense forests.
A global map of forest canopy height has been produced by [72] using GEDI RH95% and Landsat-8 time series. Yet, the validation on our test sites showed very poor results, making it necessary to produce new estimates of forest height that are accurate enough at local scale to allow forest management and monitoring of these temperate complex forests. We have shown that our processing chain using more diversified satellite images and textural indices could produce good estimates of the dominant height with GEDI data as reference samples for model learning. Moreover, the quantity of GEDI samples allows the use of all input features together without a feature selection process, which will facilitate the generalisation of the method. The GEDI RH95% cross-validation showed no saturation; r 2 from 0.54 to 0.65 and RMSE from 2.7 m to 3.4 m (12.9 to 16.7% relative RMSE). The validation of the GEDI's predictions on the measurement plots showed roughly equivalent errors. This is a very good result given that GEDI RH95% comparison to ALS data had RMSE from 3.7 to 5.6 m ( Figure 6). Furthermore, the comparison of the map produced by our method and the CHM-H95 at 10 m spatial resolution (Table 2, bottom column) shows that the prediction values are similar or even better than the GEDI RH95% values. Bias and saturation trends observed between GEDI RH95% and CHM-H95 on footprints ( Figure 6) are still observed between the prediction map and CHM-H95 at the pixel level. The validation results obtained in this study for height estimation combining GEDI and remote sensing images (r 2 from 0.54 to 0.65) are slightly better than the results seen in the literature combining GLAS and Landsat (r 2~0 .46) [67,69], and similar to [72] which uses GEDI and Landsat-8 (r 2~0 .61), although the GLAD-2019 map shows poor results on our validation data (r 2 from 0.01 to 0.21). More investigation on GEDI height metrics and their relationships with reference data (measurements or ALS data) could make it possible to apply a correction of the GEDI data on these temperate forests and to use more accurate values to train the models. Alongside this, the filtering of GEDI samples could be improved (outliers, noise, acquisition dates). Thereby, GEDI data could replace or complete field measurements in machine learning methods. The accuracy of the predictions slightly decreases using GEDI as training data instead of field measurements (relative RMSE 7.3-11.4% with plot measurements against 12.9-16.7% with GEDI on the same validation plots). Nonetheless, the quantity and spatial distribution of GEDI data can compensate for the lower accuracy than that of field measurements in large areas where several forest types are present.
Once the Hdom predictor model is built, the next step is to apply it to the whole satellite image to generate forest height maps. Figure 12 shows Hdom maps over Orléans (top) and StGobain (bottom). The graphs on the left side are extracted from ALS data, the graphs in the centre are produced with GEDI and satellite imagery (methodology presented in this paper), and the graphs on the right side are the GLAD-2019 maps [72]. The forest mask presented in Section 2.3.1 has been used to select only pixels from forest classes. We see that the forest height map in StGobain shows a more mature and dense forest, with higher trees and not many gaps, whereas the Orleans forest height map shows a more open forest, with higher diversity and height levels. Predictions from models trained on GEDI and satellite images allow the retrieval of general patterns of dominant height in forest stands, although the dynamic range and spatial precision are not as good as those of the map from the CHM. Looking at the GLAD-2019 maps, we see, as in Figure 3 (Section 3.1), that the dynamic range is very low compared to the ALS data and even compared to the maps produced in this paper. This shows that it is possible to obtain better results from GEDI and satellite images when adapting the method to these managed temperate forests and with a larger diversity of satellite sensors. Beyond that, several issues must be taken into account to produce accurate maps over larger territories: standardise the satellite image data (temporal shifts between orbits, different cloud frequencies among regions, etc.) and split the territory between eco-climatic regions, improve the quality of land cover maps (forest localisation and species), and design a nomenclature that allows the best precision of the models by class.  [67,69], and similar to [72] which uses GEDI and Landsat-8 (r 2~0 .61), although the GLAD-2019 map shows poor results on our validation data (r 2 from 0.01 to 0.21). More investigation on GEDI height metrics and their relationships with reference data (measurements or ALS data) could make it possible to apply a correction of the GEDI data on these temperate forests and to use more accurate values to train the models. Alongside this, the filtering of GEDI samples could be improved (outliers, noise, acquisition dates). Thereby, GEDI data could replace or complete field measurements in machine learning methods. The accuracy of the predictions slightly decreases using GEDI as training data instead of field measurements (relative RMSE 7.3-11.4% with plot measurements against 12.9-16.7% with GEDI on the same validation plots). Nonetheless, the quantity and spatial distribution of GEDI data can compensate for the lower accuracy than that of field measurements in large areas where several forest types are present.
Once the Hdom predictor model is built, the next step is to apply it to the whole satellite image to generate forest height maps. Figure 12 shows Hdom maps over Orléans (top) and StGobain (bottom). The graphs on the left side are extracted from ALS data, the graphs in the centre are produced with GEDI and satellite imagery (methodology presented in this paper), and the graphs on the right side are the GLAD-2019 maps [72]. The forest mask presented in Section 2.3.1 has been used to select only pixels from forest classes. We see that the forest height map in StGobain shows a more mature and dense forest, with higher trees and not many gaps, whereas the Orleans forest height map shows a more open forest, with higher diversity and height levels. Predictions from models trained on GEDI and satellite images allow the retrieval of general patterns of dominant height in forest stands, although the dynamic range and spatial precision are not as good as those of the map from the CHM. Looking at the GLAD-2019 maps, we see, as in Figure 3 (Section 3.1), that the dynamic range is very low compared to the ALS data and even compared to the maps produced in this paper. This shows that it is possible to obtain better results from GEDI and satellite images when adapting the method to these managed temperate forests and with a larger diversity of satellite sensors. Beyond that, several issues must be taken into account to produce accurate maps over larger territories: standardise the satellite image data (temporal shifts between orbits, different cloud frequencies among regions, etc.) and split the territory between eco-climatic regions, improve the quality of land cover maps (forest localisation and species), and design a nomenclature that allows the best precision of the models by class.
(a) The fact that GEDI data are a one-time operation (acquired only over a period of 3 years) is a problem for continuing medium-or long-term forest monitoring. However, the Multi-footprint Observation LiDAR and Imager system (MOLI) is planned by JAXA to provide continuity in the 2022-2027 period.
The amount of data and the large spatial coverage will help to better understand forest parameters estimation using satellite images. Moreover, models learned from a large representativeness of forests should be able to be applied over the following years without new learning processes and thus allow annual map production with our processing chain. Future missions providing global scale open access images, such as NISAR (L-band time series at 10 m spatial resolution) or Sentinel-HR (optical images at very high spatial resolution), will help to provide accurate forest information maps and make sustainable forest monitoring possible with spaceborne images.

Conclusions
This study confirms the interest of machine learning techniques for forest parameters estimation in European temperate forests. These forests are highly variable (heterogeneous stands, uneven age, mixed species, management, environment), making it difficult to predict forest parameters on a large scale by simple laws relating remote sensing observables to a specific forest parameter or physical based models. There is a need for more complex systems such as machine learning methods. Our processing chain was validated on three semi-natural managed broadleaved forests and one pine forest in temperate areas. We demonstrated that the combination of Sentinel-1 backscatters, Sentinel-2 indices and textural metrics at 10 m spatial resolution from both sensors, together with ALOS-2 PALSAR-2 backscatter, was able to provide accurate estimations of forest dominant height with RMSE from 1.9 to 3.1 m (7.3 to 11.6% relative RMSE). Replacing field data with GEDI RH95% data for model learning provided satisfactory results, with RMSE from 3.0 to 4.6 m (12.8 to 16.7% relative RMSE). The study also shows the valuable contribution of GEDI data by using the GEDI RH95% as a reference dataset in the machine learning approach. This contribution is twofold: on one hand, the large amount of reference data using GEDI data allows the stepwise selection of input features in the processing chain to be bypassed, reducing the computational time and complexity. On the other hand, the use of GEDI data as reference data allows more representativeness and forest height mapping to be provided in areas where little or no in situ data are available. The fact that GEDI data are a one-time operation (acquired only over a period of 3 years) is a problem for continuing medium-or long-term forest monitoring. However, the Multi-footprint Observation LiDAR and Imager system (MOLI) is planned by JAXA to provide continuity in the 2022-2027 period.
The amount of data and the large spatial coverage will help to better understand forest parameters estimation using satellite images. Moreover, models learned from a large representativeness of forests should be able to be applied over the following years without new learning processes and thus allow annual map production with our processing chain. Future missions providing global scale open access images, such as NISAR (L-band time series at 10 m spatial resolution) or Sentinel-HR (optical images at very high spatial resolution), will help to provide accurate forest information maps and make sustainable forest monitoring possible with spaceborne images.

Conclusions
This study confirms the interest of machine learning techniques for forest parameters estimation in European temperate forests. These forests are highly variable (heterogeneous stands, uneven age, mixed species, management, environment), making it difficult to predict forest parameters on a large scale by simple laws relating remote sensing observables to a specific forest parameter or physical based models. There is a need for more complex systems such as machine learning methods. Our processing chain was validated on three semi-natural managed broadleaved forests and one pine forest in temperate areas. We demonstrated that the combination of Sentinel-1 backscatters, Sentinel-2 indices and textural metrics at 10 m spatial resolution from both sensors, together with ALOS-2 PALSAR-2 backscatter, was able to provide accurate estimations of forest dominant height with RMSE from 1.9 to 3.1 m (7.3 to 11.6% relative RMSE). Replacing field data with GEDI RH95% data for model learning provided satisfactory results, with RMSE from 3.0 to 4.6 m (12.8 to 16.7% relative RMSE). The study also shows the valuable contribution of GEDI data by using the GEDI RH95% as a reference dataset in the machine learning approach. This contribution is twofold: on one hand, the large amount of reference data using GEDI data allows the stepwise selection of input features in the processing chain to be bypassed, reducing the computational time and complexity. On the other hand, the use of GEDI data as reference data allows more representativeness and forest height mapping to be provided in areas where little or no in situ data are available.
Building on the results shown in this study, the next step is to apply our processing chain to larger areas such as regions or at country scale and to produce new and accurate