Aboveground biomass mapping by integrating ICESat-2, SENTINEL-1, SENTINEL-2, ALOS2/PALSAR2, and topographic information in Mediterranean forests

ABSTRACT The Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) provides an extraordinary opportunity to support global large-scale forest carbon mapping, but further research is needed in order to obtain wall-to-wall forest aboveground biomass (AGB) maps with this technology. The effects of vegetation structure on the performance of canopy height and AGB modeling using ICESat-2 photon-counting light detection and ranging (LiDAR) data in Mediterranean forest areas have not been previously studied in the literature. In this study, we combined recent ICESat-2 vegetation (ATL08) data, Airborne Laser Scanning (ALS)- and field-based estimates, and a multi-sensor earth observation composite for extrapolation of AGB estimates and AGB mapping. A diverse gradient of forest Mediterranean ecosystems, distributed over 19,744.15 km2 of forest area in the region of Extremadura (Spain), with different species and structural complexity forming 5 different forest types (3 Quercus spp. dominated and 2 Pinus spp. dominated forests), was used to (i) evaluate the precision of ICESat-2 canopy height estimations, (ii) develop ICESat-2-based AGB models, and (iii) generate a spatially continuous prediction of AGB by using data from the satellite missions Sentinel-1 (S1), Sentinel-2 (S2), Phased Array L-band Synthetic Aperture Radar (ALOS2/PALSAR2), and Shuttle Radar Topography Mission (SRTM). First, ALS- and ICESat-2-derived metrics that best described canopy height (p98 and rh98, respectively) were compared at the ATL08 segment level. Second, ALS-based AGB values were derived at the ATL08 segment scale. Third, ALS-based AGB estimates at the ICESat-2 segment level were used as dependent variables to fit ICESat-2-based AGB models. Fourth, a multi-sensor approach was then implemented to predict ICESat-2-derived AGB, by means of a Random Forest (RF) modeling technique, with predictors retrieved from S1, S2, ALOS2/PALSAR2, and SRTM. Finally, RF was used to generate wall-to-wall AGB maps that were compared with field-, ALS- and ICESat-2-based observations. The agreement between the ALS- and ICESat-2-derived metrics related to the canopy height distribution was higher for Pinus spp. forest than for the Quercus spp-dominated forests. The ICESat-2-based AGB models yielded model efficiency (Mef) values between 0.56 and 0.80, with a RMSE ranging from 7.76 to 17.71 Mg ha−1 and rRMSE from 19.04 to 55.21%. The multi-sensor RF models provided the following results when compared with the ICESat-2- and ALS-based AGB observations: R 2 values of 0.63 and 0.64, and RMSE values of 11.10 Mg ha−1(rRMSE = 28.15%) and 12.28 Mg ha−1 (rRMSE = 31.45%), respectively, and an approximately unbiased result (0.03 Mg ha−1 and 0.09 Mg ha−1). When applied to the field-based validation data set (4th Spanish National Forest Inventory (SNFI-4) plots = 508), the RF-derived AGB model showed a relatively lower predictive capacity (R 2 = 0.45), a higher RMSE value (25.88 Mg ha−1) and slightly biased results (−1.47 Mg ha−1), especially for larger field-derived AGB intervals. The results of this study serve to provide an initial quantitative assessment of the ICESat-2 ATL08 data for large-scale AGB estimation. The findings suggest that a multi-sensor approach may be feasible for extrapolating ICESat-2-derived AGB estimates over areas where field or ALS reference data are not available.


Introduction
Regional-scale forest above-ground biomass (AGB) mapping using National Forest Inventories (NFI) is important for an effective forest management planning (Ene et al. 2013). Remote sensing data have been widely used for mapping forest inventory variables such as AGB (White et al. 2016;Gherardo et al. 2020;Saarela et al. 2020). Within this domain, active remote sensing is used to describe the threedimensional (3D) structure of forests (Wulder et al. 2012;McRoberts, Andersen, and Naesset 2014). Airborne laser scanning (ALS) has increasingly been used in the last decade to accurately map important forest variables in Mediterranean forests, e.g. canopy height (Guerra-Hernández et al. 2018), forest inventory variables , aboveground carbon (Guerra-Hernández et al. 2016) and canopy fuel characteristics (Botequim et al. 2019), at a high spatial resolution.
Large-scale forest mapping using NFI data has relied on the Area-Based Approach (ABA) to build estimation models from which to predict spatially explicit forest attributes over forest landscapes (Chen et al. 2016;Huang et al. 2019a;Tang et al. 2021;Guerra-Hernández et al. 2022). Spatial alignment between NFI field data and co-registration of ALS surveys are being included in many NFI programs to improve the accuracy of the estimates and to enable the ABA method to be applied under robust conditions in terms of positioning. The region of Extremadura (Central-West of Spain) is a good showcase example. New and corrected positions for the fourth Spanish NFI (SNFI-4) samples have been collected, paving the way to improving the performance of ABA at the regional scale and optimizing model fitting with less training data (Pascual et al. 2020;Guerra-Hernández et al. 2022). However, at larger scales, a design-based NFI approach may not produce accurate estimates at subregional or local scales (Tomppo et al. 2008). On the other hand, ALS 3D mapping comes with high acquisition costs for small-scale projects, and budget constraints and the availability of technology constrain the periodic mapping of AGB worldwide with the aim of monitoring carbon stocks and fluxes. Countrywide ALS-acquisitions are increasingly common in rich countries, but they remain an unreachable goal in many less developed countries, including some countries in the Mediterranean basin.
Satellite remote sensing is an alternative approach, to supporting AGB surveys and expanding coverage over remote territories without ALS data. Current spaceborne laser scanning missions have become practical and available options Narine, Popescu, and Malambo 2020). The Global Ecosystem Dynamics Investigation (GEDI)  and Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) (Markus et al. 2017;Neumann et al. 2019), both space-based laser altimeters, provide unprecedented opportunities for observing forest structures worldwide. The GEDI mission operates within the latitudinal limits of the International Space Station (ISS) orbit (i.e. ±52° latitude). The GEDI laser instrument generates a total of 8 ground tracks that are spaced approximately 600 m apart in the cross-track direction. Each track consists of ~25 m footprints spaced every 60 m along the ground track within a ~ 4.2 km swath. GEDI waveforms are processed to provide canopy height and profile metrics (GEDI Level 2 product). GEDI Level 3 product includes gridded canopy height metrics products generated from Level 2 products and GEDI level 4A and 4B products provide footprint and gridded aboveground carbon estimates ICESat-2 uses a different approach to data collection than the GEDI mission. Notably, its revisiting orbit cycle (91 days) and inclination (92°) can even yield observations in boreal forests at high latitudes (88° N/S) (Montesano et al. 2015). ICESat-2 contains a multibeam photon-counting laser altimeter (Advanced Terrain Laser Altimeter System, ATLAS), which splits a 532 nm laser beam into 6 beams arranged in 3 beam pairs. Each beam pair comprises a weak beam and a strong beam, with a detection ratio of approximately 1:4. ATLAS generates footprints at a laser repetition rate pulse of 10 kHz, resulting in a separation of 0.7 m between shots in an along-track direction, with a footprint size of 10-12 m (Neumann et al. 2019). The mission provides several along-track products, including a dedicated land and vegetation data product or ATL08 . The ATL08 product yields ground elevation above sea level, canopy height, and other descriptive variables estimated from a fixed step length of 100 m (as a segment) along the ground track (Neuenschwander and Pitts 2019). The sampling geometry of segment size 100 m × 12 m guaranteed the availability of sufficient numbers of photons for ground and canopy height estimation (Figure 1). Accuracy assessment of terrain and canopy height estimates from the GEDI and ICESat-2 missions are essential for model calibration and the development of applications to support decision-making. Different studies have assessed the accuracy of canopy height estimates produced by GEDI and ICESat-2 (Potapov et al. 2021;Li et al. 2020;Adam et al. 2020;Lang et al. 2022; and their utility for AGB estimation Malambo 2019, 2020;Duncanson et al. 2020;Neuenschwander et al. 2020;Silva et al. 2021). In most cases, temporally matched or simulated ALS data have been used to confirm canopy heights (Potapov et al. 2021;Li et al. 2020;Adam et al. 2020;Lang et al. 2022) and to calibrate GEDI-and ICESat-2-based AGB models Narine, Popescu, and Malambo 2020;Duncanson et al. 2020;Neuenschwander et al. 2020;Silva et al. 2021). As such, the accuracy assessment of post-launch ICESat-2 data and products is fundamentally important (Liu, Cheng, and Chen 2021), because results might differ from those achieved using simulated data and then to use spaceborne laser data to calibrate local and regional ALS-based AGB models (Dorado-Roda et al. 2021). Therefore, there is still a need to characterize vegetation using ICESat-2 over a wider spectrum of forest landscapes considering vertical and horizontal structural forest complexities.
The ecosystem dynamics in Mediterranean forests are particularly interesting, and the interpretation of single-photon detection technology remains more uncertain than in homogeneous boreal forest conditions Martin, White, and Coops 2021). Testing the performance of ICESat-2 to describe the forest structure and to predict AGB could provide important insights into how the vegetation structure of different forest types influences the performance of ICESat-2-based estimates. The examination of AGB estimation with ICESat-2 can also facilitate comparisons with ALS-based estimation, especially when most studies have been developed for small study areas, where variation in forest complexity was not assessed. Thus, one of the main motivations for conducting this research was to analyze different types of forests over a large region.
This study investigates Mediterranean forests with canopy cover varying from scattered to dense, with the aim of comparing model performance for ALS- predicted and ICESat-2-derived AGB. In this study, recent ALS surveys (2018-2019) and near coincident ATL08 (version 4) data (2018-2019) were used to develop and evaluate ICESat-2-based models for different forest types. Since the time interval between satellite products and ground truth information may be a source of additional errors, we validated the ICESat-2 ATL08 product (version 4) using ALS data acquired in the same time period (2018-2019). As far as we know, studies evaluating the performance of AGB models derived from the ongoing satellite ICESat-2 mission have not been conducted in Mediterranean forests. Yet, such an understanding is critical as a first step to monitoring terrestrial carbon fluxes with space-based laser data.
ICESat-2 data are acquired in transects over the landscape, and spatially continuous AGB maps are not mission data outcomes (Narine, Popescu, and Malambo 2019). For wall-to-wall estimates, ICESat-2-based AGB must be extrapolated from the segment scale to areas without ICESat-2 coverage (Shen et al. 2018). Therefore, auxiliary data from alternative spaceborne missions are needed to achieve wall-towall coverage (Huang et al. 2019a;Li et al. 2020). Quite a few studies have recently demonstrated the feasibility of combining ICESat-2 with multispectral imagery and SAR data from other satellites to produce spatially explicit information for key forest structural variables, such as canopy height and AGB (e.g. Jiang et al. 2021;Li et al. 2020;Nandy, Srinet, and Padalia 2021;Narine, Popescu, and Malambo 2019). Thus, the main goal of this study was to evaluate the usefulness of ICESat-2 data to estimate canopy height and aboveground biomass in five different types of Mediterranean forest in Central-West Spain. The motivations and needs for the study are summarized in the following four specific objectives: 1. Evaluate the accuracy of ICESat-2-derived canopy height statistics by comparing these with ALSderived metrics collected in the same time period, 2. Apply previously published ALS-based models for Mediterranean forests in the region to analyze the performance of ICESat-2-derived statistics on canopy metrics (height and cover) as exogenous variables to predict AGB, 3. Construct a wall-to-wall map of AGB at 25 m resolution by integrating ICESat-2 with multi-source remotely sensed data (Sentinel-1 (S1), Sentinel-2 (S2), Advanced Land Observing Satellite-2 (ALOS2/ PALSAR2) and the Shuttle Radar Topography Mission (SRTM)), and 4. Compare generated AGB maps with field-, ALS-, and ICESat-2-based AGB observations.

Study area, NFI, and ALS data sets
This study was carried out in the region of Extremadura (Central-West Spain) (Figure 1), covering an area of about 27,300 km 2 of diverse forested landscapes. Elevation values vary from 116 m in the Guadiana valley to 2,405 m in the Calvitero peak located in the north-eastern part of the region, presenting a mean altitude of 425 m. The climate is a Mediterranean semiarid, with a mean annual precipitation varying from <400 mm in the center part of the region to >1000 mm in the northern and eastern mountainous areas. The Spanish Forest Map (Scale: 1:25,000) (SFM25) (MAPA 2018) and the 4th Spanish National Forest Inventory (SNFI-4) datasets were used as field training data for the AGB model, representing a wide spectrum of forest structural complexities in the region (MAGRAMA 2017). The most recent SFM25 version was based on a refined classification of the main species present and forest stage of development using aerial photointerpretation. We selected five dominant Mediterranean forest ecosystems in the region: (i) Dehesas, (ii) Encinares, (iii), Pinaster, (iv) Alcornocales, and (v) Pinea (for more details about the description of the forest types, see Dorado-Roda et al. 2021).
The set of SNFI-4 training plots used to model AGB consisted of 508 concentric circular plots with 25 m of radius and accurately georeferenced plots (Table 1). These field inventory data were collected in 2017. Tree count and tree-level measurements of diameter at breast height and total height were extrapolated to area-level estimates of tree density (N, stems ha −1 ), basal area (G, m 2 ha −1 ) and AGB (Mg ha −1 ). The allometric equations of AGB used in the SNFI for each species were applied at tree level (Montero, Ruiz-Peinado, and Munoz 2005;Ricardo, Montero González, and Del Rio 2012;Ruiz-Peinado, Del Rio, and Montero 2011). For further details about the SNFI-4 field data and processing, see (Álvarez-González et al. 2014) and Dorado-Roda et al. (2021).
The ALS datasets used were covered the regions of Extremadura: EXT-N (Extremadura north, October 2018 and March 2019) and EXT-S (Extremadura south, October 2018 and July 2019). These datasets belong to the 2nd round of nationwide ALS measurements collected by the PNOA project (for more details about ALS data acquisition, see 2021).

Implemented workflow
The complete workflow is shown in Figure 2. ALS point cloud data were processed using LAStools software (Isenburg et al. 2021) following the methodological steps described by Pascual et al. (2020). As a result, 23 ALS-derived metrics ( Table 2) computed for the area of each SNFI4 plot were used as potential explanatory variables to build ALS-based AGB estimation models (the AGB estimation models were fitted using fieldbased AGB as endogenous variable and ALS-derived metrics as exogenous variables) (see details in Dorado-Roda et al. (2021)). Using an area-based approach (ABA), these models were developed by integrating two Spanish countrywide datasets: the ALS PNOA project (the Spanish National Plan for Aerial Ortophotography and LiDAR) and the fourth Spanish National Forest Inventory (SNFI4). Using the extent and location of ATL08 segments, a range of metrics were first derived from the ALS data for inclusion as explanatory variables in previously published models to produce ALS-based estimates of AGB. Lascanopy was used to compute the p98 (98th height percentile), and the previously selected ALS-derived metrics included in the ALS-based AGB estimation models at the ICESat-2 segment level ( Figure 2). Second, the ALS-and ICESat-2-derived metrics that best describe canopy height (p98 and rh98, respectively) were compared at Table 1. Summarized field data from SNFI-4 for the five dominant forest types in the region of Extremadura at plot level: aboveground biomass (AGB, Mg ha −1 ), stand basal area (G, m 2 ha −1 ), and tree density (N, trees ha −1 ).  Figure 2. Flowchart based on Lastools software, scripts in R and Python programming languages and GIS geospatial sequences used to (i) process ALS data, (ii) extract ALS-derived metrics at the level of the ICESat-2 segments, (iii) generate ICESat-2-based AGB model estimates for each forest ecosystems, and (iv) map AGB using a multi-sensor extrapolation approach. the ICESat-2 segment level. Third, the ALS-based AGB estimates at the ICESat-2 segment level were used as a dependent variable to fit ICESat-2-based AGB models, which were analyzed in terms of performance. Fourth, the multi-sensor approach was tested to predict AGB, by means of a Random Forest (RF) modeling technique, with predictors retrieved from S1, S2, ALOS2/ PALSAR2, and SRTM. Finally, RF was used to generate wall-to-wall AGB maps that were compared with field-, ALS-and ICESat-2-based observations.

ICESat-2 data acquisition and processing
The boundary of the Extremadura region was used to select all ICESat-2 granules in the study area. The ATL08 V004 data acquired by the ICESat-2 satellite between October 2018 and September 2019 were downloaded and processed using the NSIDC DAAC Data Access and Service API (https://nsidc.org/data/ ATL08/versions/4) ). The ATL08 data product is delivered in HDF5 file format. The pysl4landICESat-2.py (https://github.com/remo tesensinginfo/pysl4land) tool was used to generate ATL08 *.gpkg file from HDF5 files for only strong beams. A total of 95 files were created in HDF5 format. ICESat-2 segments that were completely included in SFM25 polygons were first selected for the five forest types. A final filter was used to exclude ICESat-2 segments where forest height estimates were less than 2 m and where the 98 th height percentile (rh98) values were higher than the maximum tree height measured in the SNFI4 plots within Extremadura. After spatial and data quality filtering (night_flag = 1 and no scattering or msw = 0), we obtained 10,009 sample segments containing valid measurements (Table 3): 7,712 for Dehesas, 1,055 for Encinares, 436 for Alcornocales, 620 for Pinaster, and 186 for the Pinea forest ecosystem.
The ATL08 algorithm has been updated various times, and version 4 represents an improvement in canopy height estimation at high vegetation densities in order to classify photons as ground, top of the canopy, canopy, and noise . ATL08 contains three canopy height variables, specifically, h_max_canopy (maximum canopy height), h_canopy (98% relative height), and cano-py_h_metrics (i.e. height metrics from the cumulative distribution of relative heights, calculated at 5% intervals for the range 10-95%) ( Table 4). Although the h_max_canopy is equivalent to the rh100 metric , its use can lead to errors because solar background noise may not have been completely removed. Therefore, h_canopy (rh98, the 98% relative canopy height) can be used as the height at the top of the canopy. The ATL08 algorithm can be affected by solar background noise. Thus, the solar elevation angle, seasonal variation, and beam intensity may impact the signal radiometry which will affect the number of photons in the segment (n_seg_ph) classified as ground (n_te_photons) and canopy (n_ca_photons). In addition to canopy height metrics, two canopies cover metrics (Narine, Malambo, and Popescu 2022) were computed from the ATL08 product, considering a combination of the variables n_ca_photons (number of canopy photons), n_toc_photons (number of top canopy photons), and n_te_photons (number of terrain photons identified in segment) (Table 4).  To assess the strength of the relationship between the ICESat-2 and ALS-derived canopy height metrics (i.e. between rh98 and p98), Pearson's correlation coefficient (r) (Eq. 1), overall root mean square error (RMSE, Eq. 2), relative root mean square error (rRMSE, Eq. 3), Bias (Eq. 4), and rBias % (Eq. 5) were used.
Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi P n i¼1 y i À � y i ð Þ 2 q (1)

RMSE ¼
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where n is the number of ICESat-2 segments, x i is the p98 (m) computed at the level (100 × 12-m) of the ICESat-2 i-segment, y i is the rh98 (m) computed from ICESat-2 ATL08 product at the level of the ICESat-2 i-segment, and � x is the mean of p98 observed at the levels of the ICESat-2 segment.

ICESat-2-based AGB models
ALS-based AGB estimation models specifically developed for the five forest types are described in detail by Dorado-Roda et al. (2021) (see Section 2.5 for methodology and Table A1 (appendix A) to check the performance of ALS-based AGB prediction models and evaluation of the plot-level accuracy). The ALSbased AGB models were applied to each of the 10,009 ICESat-2 segments for which the complete set of ALSderived metrics (Table 2) was computed. After predicting the ALS-based AGB (using the previously mentioned ALS-based AGB models), this was used as an endogenous variable to fit the ICESat-2-based AGB models, while the ICESat-2-derived metrics were used as exogenous variables in this process. The steps were similar to those described in Dorado-Roda et al. (2021) (see Section 2.6. GEDI-derived AGB models), but with ICESat-2-derived metrics instead of GEDI-derived metrics, and at the ICESat-2 segment level instead of GEDI footprints. In this case, we proposed restricting the models to two explanatory variables in order for them to be parsimonious, and thus Table 4. Set of number of photon statistics derived from ATL08 data for each ICESat-2 segment ranging within the study area. only two exogenous variables computed from ATL08 (one based on height, and one based on a cover metric) were retained in the final models to estimate AGB. We also applied a 3-fold cross-validation procedure to each potential regression model, using the functions available in R (R Core Team 2020). We then calculated the performance of each forest-type model at the segment level. We used the model efficiency (Mef) statistic (Eq 7), which yields a simple index of relative performance, where Mef values close to 1 indicate "a perfect" model fit, and negative values reveal poor model performance (Vanclay and Peter Skovsgaard 1997). Finally, we computed the overall root mean square error (RMSE, Eq. 8), relative root mean square error (rRMSE, Eq. 9), Bias (Eq. 10), and rBias (Eq. 11) to establish the accuracy of using ICESat-2-based models.
RMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where n is the number of segments, y i is the ALSbased AGB estimates at the ICESat-2 i-segment level, ŷ i is the value of AGB estimated from the species specific ICESat-2-based AGB model, � y is the mean observed value of the ALS-based AGB estimates at ICESat-2 segment level, and p is the number of model parameters.

Predictor variables for wall-to-wall AGB mapping
Wall-to-wall statistics from three different satellite remote sensing products were retrieved from (i) multispectral data from S2, (ii) SAR data from S1 and ALOS2/PALSAR2, and (iii) topographical data from SRTM. S1 and S2 images from October 2018 to October 2019 were processed using the Google Earth Engine (GEE) platform (Gorelick et al. 2017). The S2 data corresponds to the Level-2A product which was already atmospherically corrected using the sen2cor processor (Louis et al. 2016). From the 13 S2 spectral bands, only bands at 10 m and 20 m spatial resolution were used: Blue, Green, Red, Red edge 1, 2, and 3, NIR1, NIR2, SWIR1, and SWIR2. Each S2 image was pre-processed to remove clouds and poor-signal pixels using the QA60 band, which is a quality flag developed to identify and mask pixels affected by clouds and cirrus clouds. Based on the cloud-masked S2 imagery, six vegetation indices were calculated, namely the normalized difference vegetation index (NDVI), the normalized difference red edge index 1, 2, and 3 (NDRE1, NDRE2, and NDRE3), the normalized difference infrared index (NDII), and the normalized difference water index (NDWI) ( Table A2, appendix A). S1 is a C-band radar system providing images in dual polarizations (HH+HV, VV+VH) (Torres et al. 2012). In this study, the VV (transmitter-vertical and receiver-vertical) and VH (transmitter-vertical, receiver-horizontal) polarizations were extracted from S1 Level-1 Interferometric Wide Swath (IW) Ground Range Detected (GRD) at a spatial resolution of 10 m. In GEE, the S1 imagecollection is processed using the SNAP (Sentinel Application Platform) Toolbox to produce a calibrated and orthocorrected product. The pre-processing steps include thermal noise removal, radiometric calibration that converts the intensity into normalized backscatter coefficient (σ°) in decibels (dB), and terrain correction using the Digital Elevation Model (DEM) data from the SRTM (Rosenqvist et al. 2007). After the preprocessing steps for both S1 and S2 data sets, a median composite image was computed for each band by temporally aggregating high-quality observations for each of the four seasons: autumn, winter, spring, and summer (Table A2, appendix A).
Two annual HH and HV mosaics (2018 and 2019) were obtained from the Phased Array L-band Synthetic Aperture Radar sensor on the Advanced Land Observing Satellite-2 (ALOS2/PALSAR2) of the Japanese Aerospace Exploration Agency (JAXA). These data, processed in GEE, were orthorectified and slope-corrected using the 90 m SRTM DEM, and they were converted to backscatter gamma-naught (γ0) values in decibels (dB) using the equation provided by (Rosenqvist et al. 2007). For the S1 and ALOS2/ PALSAR2 polarizations, three textural metrics were computed using the Gray-level Co-occurrence Matrix (GLCM) algorithm (Haralick 1979), specifically Mean, Variance, and Homogeneity. These textural metrics were calculated with a window size of 3 × 3 pixels using the R package glcm (v.1.6.1) (Chen et al. 2016). As topographical predictor variables, elevation, slope, and aspect metrics were derived from the SRTM DEM data set at a 30-m resolution using GEE. All the abovementioned datasets were resampled to a 25-m resolution.

Wall-to-wall forest biomass estimation
The AGB estimates derived from ICESat-2 segmentlevel data (Table 3) were used as dependent variables and were related to multispectral, SAR, and SRTM topographical parameters as predictors, to construct a spatially continuous forest AGB map for the entire area. For this purpose, the Random Forest (RF) algorithm (Breiman 2001) was used to establish the relationship between ICESat-2-derived AGB estimates and up to 91 candidate predictor variables from the multisource catalog (Table A2, appendix A). The good performance of RF in predicting forest attributes with remote sensing data as explanatory variables has been widely demonstrated (e.g. Gherardo et al. 2020;H. Huang, Liu, and Wang 2019;2019a).
The RF algorithm has two important tuning parameters, mtry and ntree. The first refers to the number of predictor variables randomly selected at each split, and p p, where p the number of predictor variables, was used to determine the value of mtry. The second parameter is the final number of independent trees to be grown. In this study, the ntree was set to 1000 trees to produce stabilized variable importance estimates (Liaw and Wiener 2002). For RF modeling, the ICESat-2 segment dataset was divided into training (80%) and test (20%) sets using createDataPartition function in the caret R package (Kuhn 2015). This function creates balanced splits in the data by ensuring that random sampling occurs within each forest type while also preserving the overall class distribution over the data set (Kuhn and Johnson 2013). The training dataset was used for RF model building, while the test dataset was used to evaluate model performance by examination of R 2 , RMSE, rRMSE, Bias, and rBias. The spatially continuous AGB maps for the five ecosystems were assessed using the ALS-based AGB estimates (section 2.2) as reference data for computing R 2 and RMSE.

Accuracy of ICESat-2-vs ALS-derived metrics
For the five forest ecosystems considered, the p98 -h_canopy (rh98) comparison produced r Pearson values between 0.70 and 0.83 for Dehesas, Encinares, and Alcornocales and a value of 0.93 for Pinaster and Pinea forests (Table 5, Figure 3). For Dehesas, Encinares, Alcornocales, Pinaster, and Pinea, the RMSE values for the p98 -h_canopy (rh98) comparisons were 0.95, 1.64, 1.71, 2.24, and 1.24 m, respectively. The respective rRMSE values were 12.13, 25.62, 22.31, 17.99, and 11.94%, and bias values were −0.26, −0.14, −0.46, −0.56, and −0.30 m. Finally, the mean differences between h_canopy and p98 metrics for CC ALS -based classification of the values (see Table 2) for all the forest types are shown in Figures 3b, d, f, h and j. For all forest types, negative differences in average bias were found for the rh98 and p98 comparison, which shows that ICESat-2 underestimates canopy height (rh98) relative to the estimate derived from ALS p98.

Performance of ICESat-2-based AGB models
The performance of the ICESat-2-based AGB models for the five forest types considered is shown in Table 6. Scatterplots of ALS-observed against ICESat-2 segment-level AGB estimates are presented in Figure 4 for the best-performing forest-  Table A1 and Appendix A. The final models included two variables (one related to the height distribution and another to canopy cover), which were highly statistically significant (P < 0.001) explanatory variables. The positive and negative mean values for bias (Mg ha −1 ) and rBias (%) indicate that ICESat-2 overestimated (Encinares, Alcorncoles, Pinaster, and Pinea) and underestimated (Dehesa) the AGB relative to the ALS-based observations, respectively. Regression models for the five forest types produced Mef values ranging from 0.56 to 0.80.
The rRMSE values were slightly lower for Dehesas (19.05%) Pinaster (34.09%) and Pinea (37.54%) than for Encinares (44.29%) and Alcornocales (55.29%). The results obtained for the best models (histograms Figure 4b, d, f, h, j) indicated that the modelderived AGB estimates were slightly negatively biased at lower and higher intervals, which correspond to low and high canopy closures, respectively. The models, using mean canopy height (hmean ICS2 ) and CC1 ICS2 and CC2 ICS2 as predictors, were the most accurate and least unbiased models for all forest types, except the Encinares model, which included quadratic mean relative height (hquad ICS2 ). The models for pure homogeneous Pinaster and Pinea coniferous forest performed

Wall-to-wall forest biomass modeling
The RF model using 91 candidate predictors showed an R 2 value of 64% and a RMSE value of 10.98 Mg/ha. From that model, the importance of each independent variable in predicting AGB values was assessed by computing the percentage increase in mean squared error (%IncMSE), where higher values of % IncMSE indicate greater importance of the independent variable. The variables that made a greater contribution to AGB estimation among the 91 independent variables (dashed vertical line in Figure  A1, appendix A) are shown in Figure 5. The importance ranking showed that the summer mean composite of the near-infrared band (B8a)   Evaluation of the relevance of the Sentinel-1 C-band polarizations (VV and VH) showed that these contributed little (<6.5%) to forest AGB prediction relative to the L-band polarizations. Regarding the topographical variables, slope was the most important, being ranked 18 th among the 91 variables tested ( Figure 5). Considering the results from the relative importance assessment, a final RF model was built using only the top 24 most important variables, resulting in a substantial reduction in computation time, while suffering only a minor decrease in accuracy. The predictive accuracy of the RF model for estimating ICESat-2-derived AGB is plotted in Figure 6 for three different independent validation data sets: i) ICESat-2-derived AGB test set; ii) ALS-derived AGB; and iii) field-based AGB estimates. The estimated AGB values of the RF model showed a reasonable accuracy in predicting AGB from both ICESat-2 and ALSderived AGB data sets, with R 2 values of 0.63 and 0.64 and RMSE values of 11.10 Mg/ha and 12.28 Mg/ ha, respectively (Figure 6a and 6c). When applied to the field-based validation data set, the RF model showed a relatively low predictive capacity (R 2 = 0.43) and a higher RMSE value (25.88 Mg/ha). The observed mean bias of the RF predictions was 0.03 Mg/ha and 0.09 Mg/ha when compared with the AGB observations from the ICESat-2 and ALSbased independent data sets, respectively. The calculation of the mean bias for RF predictions and the field-based AGB estimates showed that the RF model produced a negative bias (−1.47 Mg/ha), suggesting greater underestimation of the RF predictions. Overall, the histograms in Figure 6 (b, d, and f) reveal that the RF model tended to overestimate forest AGB in the range 0 to 45 Mg/ha, while the higher (>100 Mg/ha) levels of biomass are underestimated by the RF model. To assess the effects of spatial autocorrelation on the predictive performance of the wall-to-wall mapping model, a spatial 5-fold cross-validation was also performed. To do so, the blockCV R package (Roozbeh et al. 2018) was used to compute the variogram and build the spatial blocks (more details are provided in the Appendix). Testing the RF model in a spatial 5-fold crossvalidation resulted in an R 2 value of 0.58 and an RMSE of 11.53 Mg/ha, which are similar to those obtained from the ICESat-2-derived AGB test set (R 2 = 0,63; RMSE = 11.10 Mg/ha). These results suggest that the existing spatial autocorrelation had an effect, yet small (0.048%), on the predictive performance of the RF model.

Forest biomass mapping
The map of the predicted forest AGB exhibited general consistency with the ALS-based AGB map across the entire study area (Figure 7). The predicted RF AGB values were highly correlated with the ALS-AGB map, with an r value of 0.71. However, in accordance with the results shown in the histograms in Figure 6, underestimation by the RF model of higher biomass levels (>100 Mg/ha) (Figure 8, local 3) and overestimation of moderate biomass levels (30-45 Mg/ha) (Figure 8, local 2)

Discussion
Our study evaluated ICESat-2ʹs performance using robust AGB models under temporally coincident Figure 6. Relationship between AGB RF-predictions (y-axis) and three different independent data sets (x-axis): a) ICESat-2-derived AGB test set, c) ALS-derived AGB, and e) field-based AGB estimates. Right column represents the histograms of the reference (gray) and the predicted AGB (yellow).
ALS and ICESat-2 datasets, and over five forest types with different vegetation structures. Despite recent literature evaluating the accuracy of real ICESat-2 data (Liu, Cheng, and Chen 2021;Nandy, Srinet, and Padalia 2021) or ALS-simulated ICESat-2 data Narine, Popescu, and Malambo 2020;Neuenschwander et al. 2020;Silva et al. 2021), to the best of our knowledge, our study is the first one evaluating the performance of on-orbit ICESat-2 ATL08 (Version 4) products for characterizing AGB by comparing with spatially and temporally coincident ALS coverage across vast areas of diverse Mediterranean forests. Our study area covers a variety of forest types and thus provides a basis to explore the impacts of forest structure on canopy height and forest AGB biomass estimations. Additionally, this work aimed to map AGB by integrating multi-sensor earth observation data to extrapolate the AGB derived at the ICESat-2 segment level to other areas without ICESat-2 or ALS coverage.

Accuracy of ICESat-2-derived canopy height
Our results are consistent with those reported in previously published pre-launch studies that analyzed simulated ICESat-2 data (e.g. Amy and Magruder 2016), which found simulated top-ofcanopy heights from ICESat-2, underestimating true top-of-canopy returns for all types analyzed with negative bias ranging from 0.28 m (1.39 m RMSE) to 1.25 m (2.63 m RMSE). However, it is important to evaluate the accuracies of post-launch ICESat-2 data across other forest types beyond boreal forests Martin, White, and Coops 2021). Our results on canopy structure metrics show that the ATL08 relative height metrics (rh98) are reasonably accurate, with RMSE values stabilizing at 1.5-2 m in Mediterranean forest (RMSE varied from 0.95 m in Dehesas to 2.24 m in Pinaster), and substantially better than the values obtained in terms of performance from p98 -rh98 at the GEDI footprint level (level 2A using Version 0001) for the same study area (Dorado-Roda et al. 2021). In terms of bias, the p98 -rh98 relationship was slightly better in Dehesas than the other forest types. The results highlight that the accuracy of ICESat-2 canopy heights is influenced by the vegetation structure. ICESat-2 segment estimates were better for regular-even-aged coniferous forests of P. pinaster and P. pinea species (r = 0.93) than in sparse canopy Quercus-dominated forests with values of r ranging from 0.70 (irregular-uneven-aged-multi-layered Encinares and Alcornocales forests) to 0.83 (sparse homogeneous Dehesas forests).
The comparison between ALS-and ICESat-2-derived metrics showed that the results for coniferous-dominated forest were better in terms of r than those obtained by Martin, White, and Coops (2021) who compared ICESat-2 ATL03 top of canopy photon height and 90 th height percentile ALS in boreal forest (r = 0.84, RMSE = 2.5 m, rRMSE = 19.2%, bias = −2 m). The performance between metrics was better in our study in terms of RMSE, rRMSE, and bias, when compared to a recent study over 40 sites covering the mainland USA using nighttime acquisitions and strong beam data (same as our study) and ALS data (1-4 point m −2 ) collected by the National Ecological Observatory Network (NEON) in 2019 (Liu, Cheng, and Chen 2021) (RMSE = 3.93 m, rRMSE = 24.8%, and negative bias = −0.87, n = 7543 segments). Only the rRMSE for Encinares forest produced less accurate values than those reported in Martin, White, and Coops (2021). The comparison between ALS and ICESat-2-derived metrics in this study was slightly better in terms of bias and slightly Figure 8. Maps comparing the observed ALS-derived AGB and the RF predicted AGB for three different locations. The third row depicts the difference between ALS-derived AGB with the RF predicted AGB. White pixels represent non-forest areas.
worse in terms of rRMSE than the results reported in Neuenschwander et al. (2020). Neuenschwander et al. (2020) compared ICESat-2-derived canopy height (h_canopy = rh98) with the same metric from the ALS data in boreal forests in Finland (strong beam/ night/summer acquisitions, underestimate with a negative bias = −0.56 m, bias% = 3.18%; RMSE % = 13.75%). Our results confirmed negative differences on average bias, between rh98 and p98 (Figure 3), indicating that ICESat-2 underestimates canopy height (rh98) at the segment level when compared with ALS p98 computed at the same extent. This negative bias compared to Neuenschwander et al. (2020), and Liu, Cheng, and Chen (2021) (see figure 11 b) is likely related to the forest vertical structure (2-20 m interval in our study). Overall, our bias and RMSE values were better than those obtained by Neuenschwander et al. (2020) and Liu, Cheng, and Chen (2021), respectively. Mediterranean vegetation has, in general, different structural complexities when compared to boreal forests and the range of type of forest analyzed by Liu, Cheng, and Chen (2021) in the USA. This may be the main reason for the relatively low errors compared with these recently published studies Liu, Cheng, and Chen 2021;Lonesome and Popescu 2021).
The higher precision on the estimation of canopy cover with ICESat-2 was achieved in the 40-90% range of canopy cover. The highest errors were computed in dense canopy closure conditions (>90%) consistent with the findings of Liu, Cheng, and Chen (2021). This confirms that, at low canopy cover (<20%) conditions, both ICESat-2 photon-counting and GEDI full-waveform (FW) LiDAR sensors are more likely to record returns reflected from the terrain rather than canopies, which impedes precise estimation of canopy height (Dubayah et al. 2020a(Dubayah et al. , 2020b. On the contrary, for high canopy cover conditions (CC > 90%), the terrain-reflected signal received by both sensors is weaker than the canopy signal, leading to errors in canopy height measurements . Therefore, rh metrics from both missions may be biased especially in low and high canopy cover conditions. Negative bias values exist between canopy height residuals for ICESat-2 when canopy cover exceeds 90%, indicating that the underestimation of canopy may be due to the overestimation of terrain elevation (Liu, Cheng, and Chen 2021). In summary, our results on canopy height estimation confirmed that (i) ICESat-2 was more accurate in forests with canopy cover ranging between 40% and 85% and (ii) canopy height was underestimated on average by 0.14 to 0.56 m, slightly better than recent studies in boreal forest that showed that canopy height was underestimated on average by 0.5 m -0.6 m (Amy and Magruder 2019;Neuenschwander et al. 2020;Martin, White, and Coops 2021) and by −0.87 covering the US forest (Liu, Cheng, and Chen 2021).

Performance of ICESat-2 AGB-derived models
ICESat-2-derived AGB models using mean and quadratic canopy heights and CC1 ICS2 and CC2 ICS2 as predictors represent a satisfactory description of vegetation structure compared to ALS-based AGB estimates. In terms of RMSE and rRMSE for the five forest ecosystems, the precision of the ICESat-2-derived AGB models were similar or better in the case of Pinaster (rRMSE = 36.67%) and Dehesas (rRMSE = 28.79%) than those values reported by . In , the simulated ICESat-2-derived 10 th and 90 th height percentiles and canopy cover statistics resulted, under linear regression modeling, an R 2 value of 0.74, and an rRMSE of 32% with a RMSE of 25.50 Mg ha −1 , using 85 segments as training data under a nighttime acquisition scenario. In our study, ICESat-2-based AGB estimates for Encinares, Alcornocales, and Pinea yielded rRMSE values slightly worse than those reported by . The rRMSE values achieved were also similar or better than those reported by Silva et al. (2021) (rRMSE of 54% for Sonoma County (US), using GEDI and ICESat-2 fused AGB and GEDI's AGB models from Duncanson et al. (2020)). In general, the range of values in terms of Mef or adjR 2 from ICESat-2-based models (Mef = 0.56-0.80) at the 100 × 12 m segment level were better than the values obtained by GEDI-based models (Mef = 0.31-0.46) at 25-m footprint level using version 001 in the same Mediterranean formations. The results confirmed the more complex and uneven-aged-multilayered forests as Encinares and Alcornocales (rRMSE = 44.29% and 55.21%) were, the lower the accuracies in modeling AGB for Mediterranean forests with GEDI and ICESat-2.
Our ICESat-2-derived AGB models were less unbiased in terms of bias and rbias than those values obtained by Duncanson et al. (2020) (bias = −26.3 Mg ha −1 and bias% = −18.7%) using US-wide GEDI-based AGB models that rely only on rh metrics and simulated AGB estimates at the footprint level. In general, the ICESat-2-based models were less biased at lower and higher AGB intervals than GEDI-based models, although ICESat-2 models underestimated at lower and higher canopy cover, similar to findings reported with simulated or real GEDI data (Dorado-Roda et al. 2021;Duncanson et al. 2020). Our models achieved slightly better values of Mef (0.56 to 0.80) (similar to adj. R 2 ) than those reported by Silva et al. (2021) Silva et al. (2021) used ALSsimulated GEDI and ICESat-2 data and the ALS point cloud density was higher than in our study.
The outcomes of our ICESat-2 exercise are better than those reported by (Dorado-Roda et al. 2021) at the 25 m GEDI footprint level over the same region, using the same ALS data for benchmarking. There are important reasons for the differences in the performance between AGB models based on ICESat-2 or GEDI. The impacts of the GEDI (Version 0001) geolocation errors in the study of (Dorado-Roda et al. 2021) were higher, as already discussed, than those from ICESat-2 (geolocation horizontal accuracy of ICESat-2 is <5 m and it gets close to 2-3 m , especially in scattered tree ecosystem as Dehesas (Dorado-Roda et al. 2021).
The metrics derived from ATL08 products at the segment level, such as mean canopy height (hmean ICS2 = rh98), quadratic mean relative height (hquad ICS2 ), and canopy cover metrics (CC1 ICS2 and CC1 ICS2 ), were found to be significant predictors of AGB. The exhaustive search step to identify the best metric included hmean ICS2 , which is more stable to changes in both the vertical and horizontal canopy structures (Lefsky et al. 2002;Ni-Meister et al. 2010;Asner et al. 2012;Guerra-Hernández et al. 2016). The use of canopy height metrics alone may ignore some information in profiles with more vegetation structural heterogeneities, such as in natural Mediterranean forest. The canopy cover metric (CC1 ICS2 and CC1 ICS2 ) also improved the fit in all the ICESat-2-derived AGB models, consistent with a previous study using ICESat-2 simulated data . Our results also demonstrate that a second metric related to canopy cover from photon-count LiDAR is important for improving ICESat-2 AGB-derived models (Table 6).

Wall-to-wall forest biomass mapping
The approach of combining wall-to-wall multispectral imagery, SAR, and topographic data to extrapolate estimates of forest AGB values from ICESat-2 tracks to obtain wall-to-wall coverage showed a reasonable agreement when benchmarking with ALS-based AGB predictions. The moderate agreement between RF predictions and ALS-based model estimates of AGB (R 2 = 0.64) reflects the challenging study area, and a mosaic of structurally complex Mediterranean forest ecosystems which imposes limitations in acquiring vegetation structure information from space-borne lidar systems (Glenn et al. 2016;Gwenzi et al. 2016). The performance of our RF model in predicting AGB values agrees with previous studies that utilize similar methodological approaches. For instance, , who mapped AGB by combining ICESat-2 and Landsat-8 OLI data, found an RF model explaining 58% of the variations in AGB. In a study carried out by Huang et al 2019b in China, an AGB map was produced with an R 2 value of 0.64 by integrating ICESat GLAS sensor data with Landsat and PALSAR-derived predictor variables.
To the best of our knowledge, the present study is the first to predict AGB for areas outside of ICESat-2 coverage by combining S-2, S-1, ALOS2/PALSAR2, and topographical data. The variable importance ranking indicated that both spectral information from Sentinel-2 and L-band SAR data from ALOS2/PALSAR2 played an important role in predicting AGB. From the multispectral source, it was found that near-infrared and shortwave region (B8a and NDWI) were key in predicting AGB, which is consistent with other studies focused on mapping AGB and tree canopy height (e.g. (Campbell et al. 2021;Nandy, Srinet, and Padalia 2021;Zhu et al. 2020). With regard to the L-band derived variables, this study also corroborates with findings from Huang et al. (2019a), who showed that texture features of HV and HH polarizations contribute more to AGB predictions than the original backscatter PALSAR data. This may suggest that SAR texture measurements have a higher capability to discriminate spatial information, as well as to reduce the noise in the SAR data (e.g. (Sarker et al. 2012;Laurin et al. 2017). When accounting for the combined variable importance of the top five most important variables, it was found that SAR L-band outperformed the Sentinel-2 derived variables, which demonstrates the usefulness of ALOS-2/PALSAR-2 L-band in predicting AGB. In terms of the Sentinel-1 data, this study showed that SAR C-band polarizations had a limited effect on the RF model accuracy. Previous studies assessing the utility of L-and C-band SAR data for AGB estimation emphasize that L-band data generally produce better results than C-band due to its higher penetration ability through the canopy (e.g. Huang et al. 2018;Naidoo et al. 2015). Therefore, the overall R 2 value of 0.63 and the outcomes from RF-based variable importance might justify the synergetic use of ICESat-2, Sentinel-2, and ALOS2/PALSAR2 L-band SAR satellite data to derive AGB over large areas lacking complete coverage of ALS data and temporally co-registered NFI data to build AGB estimation models.
The RF-approach to generate wall-to-wall maps of AGB may lead to improved temporal resolution as weekly/monthly data become available. As a result, changes from growth or silvicultural activities, for instance, can be traceable. Results from variable importance in the RF modeling exercise highlight the potentialities and capabilities of the upcoming NASA-ISRO Synthetic Aperture Radar (NISAR) satellite mission in 2023 (NISAR 2022) that will deliver denser L-band time-series data at a higher spatial resolution (12 m) with the potential to further boost the approach presented (Khati, Lavalle, and Singh 2021). Nevertheless, and although this methodological approach can be applied in other ecological contexts, the direct application of our RF model to other regions with different ecological characteristics should be approached with caution, particularly for areas where field and ALS reference data are not available.
The accuracy of satellite-based AGB estimations is influenced by many factors, such as the characteristics of the satellite systems used, the methodological approach implemented, the edaphoclimatic and topographic characteristics of the study area, and the human and equipment errors when measuring vegetation structural parameters from the field (e.g. Arnan et al. 2022). In this study, with the exception of the expected uncertainty associated with tree measurements in the field and the allometric equations used for AGB estimations, a potential source of uncertainty may be related to the sampling nature of the ICESat-2 mission. The sparse spatial distribution and the gaps between ICESat-2 segments over the landscape may have impacted the accuracy of the RF model used for wall-to-wall AGB mapping (e.g. Huang et al. 2019). Therefore, detailed analysis should be carried out to identify the processes of error propagation in the extrapolation processes.

Conclusions
Real ICESat-2 data were used to assess the use of canopy forest height for estimating AGB in five different Mediterranean forest types representing a wide spectrum of forest structural complexity in the region. The workflow comprises the use of ICESat-2 data and several satellite earth observation (EO) sensors (S1, S2, ALOS2/PALSAR, and STRM), in addition to two Spanish countrywide data sets (ALS PNOA project and SNFI4). The Spanish countrywide data sets were used in the training and validation phases of the ICESat-2-based AGB model development, while the EO data sets were used to scale up the AGB estimates and to generate spatially explicit AGB maps for the study area. This initial assessment of the capacity of ICESat-2 to estimate AGB indicated feasibility of the approach. The results serve as a basis for further extrapolation efforts and are of particular interest for areas or countries for which field and ALS reference data are not available. In this respect, it is important to highlight the strong correlations between the ALS-and the ICESat-2-derived metrics that best describe the top of the canopy height (p98 and rh98, respectively), together with the high percentage of the AGB variability explained by the ICESat-2-based models, since ALSbased AGB (and other key forest variables) observations are currently the most reliable ground truth reference data available for model fitting and extrapolation of estimates. These features are particularly relevant to scientists investigating methods to better understand AGB dynamics and associated uncertainties using ongoing wall-to-wall satellite image and SAR missions, especially considering the high acquisition costs of field reference data. Moreover, the high correlation between ALS-and the ICESat-2-derived metrics opens the door to further examination of the utility of previously developed species-specific or global AGB models derived in former ALS surveys and expensive field surveys. This is important in relation to country-or region-wide sampling designs, especially in countries that cannot afford to acquire new ALS or field data. Finally, the multi-sensor EO composite was analyzed and validated, producing promising results and revealing an approach that could be considered for extrapolation of AGB estimates and mapping. Nevertheless, further research should be conducted to optimize the design of field data sampling and ALS acquisition in order to reduce costs and to maintain the robustness of the estimates.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This research was supported by the project "Extensión del cuarto inventario forestal nacional mediante técnicas LiDAR para la gestión sostenible de los montes de Extremadura" from the Extremadura Forest Service (FEADER nº 1952SE1FR435) and the FUEL-SAT project "Integration of multi-source satellite data for wildland fuel mapping: the role of remote sensing for an effective wildfire fuel management" from Fundação para a Ciência e a Tecnologia (PCIF/GRF/0116/ 2019).

Data Availability
The data that support the findings of this study are available from the corresponding author, Guerra-Hernandez., J., upon reasonable request.