# A fused canopy height map of Italy (2004–2024) from spaceborne and airborne LiDAR, and Landsat via deep learning and Bayesian averaging

Yang Su, Nikola Besic, Xianglin Zhang, Yidi Xu, Saverio Francini, Giovanni d'Amico, Gherardo Chirici, Martin Schwartz, Ibrahim Fayad, Sarah Brood, et al.

► To cite this version:

Yang Su, Nikola Besic, Xianglin Zhang, Yidi Xu, Saverio Francini, et al. A fused canopy height map of Italy (2004–2024) from spaceborne and airborne LiDAR, and Landsat via deep learning and Bayesian averaging. 2025. hal-05347031

HAL Id: hal-05347031

https://hal.science/hal-05347031v1

Preprint submitted on 4 Nov 2025

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

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

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

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

1 A fused canopy height map of Italy (2004–2024) from spaceborne and airborne LiDAR, and Landsat

2 via deep learning and Bayesian averaging

3

4 Authors

5 Yang Su ¹,²,³, Nikola Besic ⁴, Xianglin Zhang ³, Yidi Xu ², Saverio Francini ⁵, Giovanni D'Amico ⁶, Gherardo

6 Chirici ⁷, ⁸, Martin Schwartz ², Ibrahim Fayad ², Sarah Brood ¹, Agnes Pellissier-tanon ², Ke Yu ², Haotian

7 Chen ⁹, ³, Songchao Chen ¹⁰, Alexandre d'Aspremont ¹, Philippe Ciais ²

8

9 Affiliations

10 ¹ CNRS & Département d'Informatique, École Normale Supérieure – PSL, 45 Rue d'Ulm, 75005 Paris,

11 France

12 ² Laboratoire des Sciences du Climat et de l'Environnement, CEA CNRS UVSQ Orme des Merisiers, 91190

13 Gif-sur-Yvette, France

14 ³ Université Paris-Saclay, AgroParisTech, INRAE, UMR ECOSYS, 91120 Palaiseau, France

15 ⁴ IGN, ENSG, Laboratoire d'inventaire forestier (LIF), 54000 Nancy, France

16 ⁵ Department of Science and Technology of Agriculture and Environment (DISTAL), University of Bologna,

17 40126 Bologna, Italy

18 ⁶ Department of Agriculture, Food, Environment and Forestry, Università degli Studi di Firenze, Via San

19 Bonaventura, 13, 50145 Firenze, Italy

20 ⁷ Fondazione per il Futuro delle Città, Florence, Italy

21 ⁸ CNR Istituto dei Sistemi Complessi, Sesto Fiorentino (Florence), Italy

22 ⁹ Institute of Crop Sciences, Chinese Academy of Agricultural Sciences / Key Laboratory of Crop

23 Physiology and Ecology, Ministry of Agriculture and Rural Affairs of China, Beijing 100081, China

24 ¹⁰ College of Environmental and Resource Sciences, Zhejiang University, 310058 Hangzhou, China

25

26 Corresponding Author

27 Yang Su yang.su@ens.fr +33 1 89 10 07 67 École Normale Supérieure – PSL

28

29 Abstract

30 Forests are vital for the carbon sequestration, biodiversity conservation, and climate regulation, making the

31 precise and continuous monitoring of forest structure attributes such as canopy height essential. Here we

32 present a two decades long (2004–2024), 30m resolution annual canopy height dataset for Italy, developed

33 using a time-series deep learning framework that integrates Landsat optical imagery with LiDAR

34 observations. Two UNET models were independently trained using canopy height reference data from

35 airborne laser scanning (ALS) and NASA's Global Ecosystem Dynamics Investigation (GEDI) spaceborne

36 LiDAR mission. Annual canopy height predictions from each model were fused using Bayesian Model

37 Averaging (BMA) to enhance spatial consistency and temporal continuity. Validation against ground-based

38 measurements from the Italian National Forest Inventory (NFI) demonstrated high predictive accuracy

39 (mean absolute error = 3.98 m). To further evaluate the utility of our dataset, we derived a canopy height

40 change-based disturbance product and validated it against observed events (mean precision = 0.64 for 2005–
<figure><img src="image_3.png" /></figure>

2016). In addition, we assessed post-disturbance recovery by monitoring canopy height regrowth in areas affected during 2004–2005, tracking changes annually through 2024 across various Italian biomes. Our results highlight the importance of integrating multi-source remote sensing data with deep learning and Bayesian data fusion for monitoring forest structural dynamics. The final dataset is publicly available via Zenodo and provides a reproducible and scalable resource to support forest research, ecological monitoring, and climate-related policy-making.

47

## 1. Introduction

Forests play a central role in regulating the Earth system through carbon sequestration, climate moderation, water cycling, and biodiversity conservation (Shukla et al., n.d.). Quantifying forest structure—particularly canopy height—is critical for estimating forest biomass, understanding ecosystem dynamics, and informing forest management and climate mitigation strategies (Chave et al., 2014). Despite its importance, accurate and long-term canopy height and height change data remain scarce, particularly at high spatial resolutions. In Italy, where forests range from Mediterranean woodlands to alpine coniferous ecosystems (Selvi et al., 2023), spatially detailed and temporally continuous monitoring is crucial for assessing ecological change and supporting sustainable forest management. Yet, most existing canopy height datasets over Italy represent only single-year snapshots, lacking the temporal depth necessary to capture forest dynamics over time (Lang et al., 2023; Pauls et al., 2024; Potapov et al., 2021).

Recent advances in remote sensing and deep learning have substantially enhanced our capacity to monitor forest structural attributes such as canopy height (Lang et al., 2023; Pauls et al., 2024). Increasingly, studies have combined discrete LiDAR measurements with continuous satellite imagery to generate spatially explicit canopy height maps over large areas (Fayad et al., 2024; Hansen et al., 2013a; Lang et al., 2023; Pauls et al., 2024; Potapov et al., 2021; Rajab Pourrahmati et al., 2024; Schwartz et al., 2023; Su et al., 2025). Deep learning techniques, such as convolutional neural networks, have further improved the accuracy and scalability of these approaches by learning complex, multi-scale spatial patterns from large training datasets (Fayad et al., 2024; Lang et al., 2023; Pauls et al., 2024; Rajab Pourrahmati et al., 2024; Schwartz et al., 2023; Su et al., 2025).

Among the available satellite platforms, the Landsat missions are the only ones providing four decades of 30m-resolution optical multispectral imagery with global coverage and frequent temporal observations (Crawford et al., 2023; Kovalskyy and Roy, 2013; Roy et al., 2014; Wulder et al., 2022). For LiDAR sources, airborne laser scanning (ALS) is widely regarded as one of the most accurate methods for estimating canopy height, providing high spatial resolution and accurate structural information (Wulder et al., 2012) (with data available from 2004 to 2017 in Italy (Montaghi et al., 2013)). However, the operational costs of ALS campaigns limit their spatial and temporal availability, especially at national or continental scales (Coops et al., 2021). In contrast, the Global Ecosystem Dynamics Investigation (GEDI), a spaceborne LiDAR mission, offers a more recent (with operational data available from 2019 in Italy) and globally distributed dataset, capturing canopy height at a large scale (Dubayah et al., 2020; Potapov et al., 2021). While GEDI offers greater spatial coverage than ALS, its temporal coverage is shorter and the footprint-based measurements remain sparse.

Critically, the capability of large-scale monitoring the vertical structure of forests opens new possibilities for both forest disturbance detection and post-disturbance recovery monitoring (Francini et al., 2022), which is more critical than ever undercurrent climate change scenarios (Palahí et al., 2021). While forest disturbances have traditionally been monitored using optical satellite data (Hansen et al., 2013b; Hermosilla et al., 2015; Kennedy et al., 2010), and several classification products have recently been developed for Italy at varying temporal resolutions (Francini et al., 2020, 2021, 2023), there remains a significant gap in the availability of maps that track changes in canopy height over time. Large-scale, long-term predictions of forest canopy height have the potential to fill this gap by enabling the creation of height change maps, which reveal structural transformations that are often invisible to optical sensors. These maps provide
<figure><img src="image_4.png" /></figure>

89

90

91

valuable complementary information to traditional disturbance products, offering enhanced insights into disturbance events, post-disturbance recovery trajectories, and forest resilience—critical components for understanding forest dynamics under changing environmental conditions.

92

93

94

95

96

97

98

99

100

101

102

103

104

105

In this study, we present a long-term (2004–2024), 30m resolution annual map of maximum canopy height for Italy. Our workflow coupled the ALS and GEDI data with Landsat multispectral data in a deep-learning framework: one UNET model (Ronneberger et al., 2015) was trained with ALS–Landsat pairs from 2004–2017, and a second with GEDI–Landsat pairs from 2019–2023. Annual canopy height predictions from 2004 to 2024 were generated independently using each model. To enhance consistency and accuracy across space and time, we applied Bayesian Model Averaging (BMA) (Hoeting et al., 1999) to fuse these ALS–Landsat and GEDI-Landsat based canopy height maps, with the measurements from the Italian National Forest Inventory (NFI) guiding the spatially adaptive weights. The resulting dataset offers a harmonized 30m-resolution time series record of canopy height across Italy, serving as a valuable tool for forest monitoring, carbon stock assessment, and land-use change analysis. We demonstrate its applicability by deriving a disturbance product based on canopy height change and tracking post-disturbance regrowth in areas affected during 2004–2005 across different Italian biomes (Figure S1). Our results highlight the utility of this dataset as a robust baseline for assessing forest structural dynamics in response to climate change and anthropogenic influences.

106

107

## 2. Data

108

### 2.1 Remote Sensing Data

109

110

111

112

113

114

115

116

117

Landsat surface reflectance products from the Thematic Mapper (TM, Landsat 5) (Crawford et al., 2023; Kovalskyy and Roy, 2013; USGS Landsat 5 Level 2, Collection 2, Tier 1, 2025; Wulder et al., 2022), Enhanced Thematic Mapper Plus (ETM+, Landsat 7) (Kovalskyy and Roy, 2013; USGS Landsat 7 Level 2, Collection 2, Tier 1, 2025), and Operational Land Imager (OLI, Landsat 8) (Roy et al., 2014; USGS Landsat 8 Level 2, Collection 2, Tier 1, 2025) were used as the primary optical inputs for canopy height modeling over the 2004–2024 period. All imagery was sourced from the USGS Collection 2, Tier 1 Level-2 Surface Reflectance datasets, which include atmospheric correction and quality assurance bands for cloud and shadow masking (Landsat Algorithms, 2025; USGS Landsat 5 Level 2, Collection 2, Tier 1, 2025; USGS Landsat 7 Level 2, Collection 2, Tier 1, 2025; USGS Landsat 8 Level 2, Collection 2, Tier 1, 2025).

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

Given the evolution of the Landsat program and sensor transitions (Vermote et al., 2016; Wulder et al., 2016), we split the data processing pipeline into two distinct periods: 2004–2012 and 2013–2024. For the first period (2004–2012), annual median composites were generated from Landsat 5 and 7. Since both sensors have spectral response characteristics that differ from the more recent Landsat 8, a harmonization step was necessary to ensure consistency across the full time series (Roy et al., 2016). We applied a linear transformation to normalize TM and ETM+ surface reflectance values to the OLI spectral space using the available slope and intercept coefficients (Landsat ETM+ to OLI Harmonization, 2025; Roy et al., 2016). The same transformation was applied to both Landsat 5 and 7, a common approach in long-term forest monitoring studies due to their comparable spectral bands and the need for alignment with the OLI reference system (Landsat ETM+ to OLI Harmonization, 2025; KC et al., 2021; Roy et al., 2014, 2016; Savage et al., 2018; Vogeler et al., 2018). Quality assessment and masking were conducted using the QA bands provided in Collection 2 Level-2 products (USGS Landsat 5 Level 2, Collection 2, Tier 1, 2025; USGS Landsat 7 Level 2, Collection 2, Tier 1, 2025; USGS Landsat 8 Level 2, Collection 2, Tier 1, 2025). From the QAPIXEL band, we excluded pixels flagged as Fill (bit 0), Cirrus (bit 2; populated for OLI only), Dilated Cloud (bit 1), Cloud (bit 3) and Cloud Shadow (bit 4). In addition, pixels flagged as radiometrically saturated in any band (QA\_RADSAT ≠ 0) were removed. After QA filtering, we then computed vegetation and water-related indices, including Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Land Surface Water Index (LSWI) and Normalized Difference Water Index (NDWI).
<figure><img src="image_5.png" /></figure>

137 For the second period (2013–2024), data availability was extended with the launch of Landsat 8, which provides improved radiometric performance and refined spectral band definitions (Roy et al., 2014). OLI imagery was processed directly with the above-mentioned QA-based masking, and no reflectance normalization was required. Data from Landsat 7 continued to be harmonized to match OLI reflectance values, ensuring spectral consistency across the entire time series (Landsat ETM+ to OLI Harmonization, 2025; Roy et al., 2016). For each year, an annual median composite was generated from all valid cloud-free images from both Landsat 7 and 8. Spectral indices were calculated using the same formulas as in the earlier period to ensure comparability across time. The script for downloading the Landsat data is accessible via Figshare (Code used in the study - "A fused canopy height map of Italy (2004–2024) from spaceborne and airborne LiDAR, and Landsat via deep learning and Bayesian averaging ", 2025). Although Landsat 9, launched in late 2021, provides data with nearly identical spectral and radiometric properties to Landsat 8 (Kaewmanee et al., 2023; Xu et al., 2023a), it was not incorporated in this study as the strong sensor similarity between OLI and OLI-2, the inclusion of Landsat 9 would not materially change the annual median composites or downstream canopy height predictions.

151 To characterize topographic variation across Italy, we used the Advanced Land Observing Satellite (ALOS) World 3D 30 m Digital Surface Model (AW3D30), provided by the Japan Aerospace Exploration Agency (JAXA) (Tadono et al., 2014, 2016; Takaku et al., 2014, 2016). The AW3D30 dataset offers global elevation data at a 30 m spatial resolution, derived from PRISM stereo imagery. We extracted the “DSM” band, which represents surface elevation including vegetation and built structures. And to support terrain-related analysis of canopy height patterns, we derived slope information using the ee.Terrain.slope() function (ee.Terrain.slope, 2025) applied to the AW3D30 mosaic. The code for downloading digital terrain information is available via Figshare (Code used in the study - "A fused canopy height map of Italy (2004–2024) from spaceborne and airborne LiDAR, and Landsat via deep learning and Bayesian averaging ", 2025).

161 **2.2 LiDAR Reference Data**

162 **2.2.1 Airborne Laser Scanning Data (ALS)**

163 ALS data used in this study were sourced from various national and regional airborne survey campaigns conducted between 2004 and 2017 across Italy (D'Amico et al., 2021; Montaghi et al., 2013). The datasets were acquired using different airborne lidar sensors, most commonly Optech ALTM systems (Gemini, 3032, 3033, and 3100; Teledyne Optech, Canada), Riegl LMS sensors (Q560, Q680i; RIEGL Laser Measurement Systems, Austria), Leica ALS50-II (Leica Geosystems, Switzerland), and TopoSys Falcon II (Germany) (Table S1). These sensors are full-waveform or discrete-return lidar systems operating in the near-infrared (NIR) spectrum (Table S1), which is standard for vegetation and topographic mapping (Budei et al., 2018; Wagner et al., 2008).

171 Flight characteristics varied widely across surveys. Flight altitudes ranged from 180 m to 4700 m above ground level, while point density ranged from 0.4 to over 5 pulses per square meter, and spatial resolution from 1 m to 5 m (Table S1). However, additional details on flight speed and overlaps were not consistently available for all surveys. When available, the lidar system settings used a forward overlap of 50–60% and a side overlap of 20–30%, with aircraft speeds generally between 150 and 250 km/h depending on altitude and swath width (Pirotti, 2010). These acquisition parameters are designed to ensure adequate point density and minimize occlusion effects in forested areas. Vertical accuracy for these ALS datasets typically ranges from 10 to 30 cm RMSE, and horizontal accuracy from 20 to 50 cm, depending on sensor type, flight altitude, and GPS/IMU quality (Optech, n.d.).

180 The temporal distribution of the datasets reflects various institutional objectives: some surveys were created for hydraulic risk assessment (e.g., MATTM acquisitions along rivers and coastlines), others for regional environmental monitoring (e.g., Piemonte, Trentino, Valle d'Aosta), and some for local forest or ecological studies (e.g., Bosco Fontana). The seasonal timing of the acquisitions is not consistently reported across providers. However, available metadata and project documentation show that most campaigns took place
<figure><img src="image_6.png" /></figure>

185 during late spring to summer months, when cloud-free and leaf-on conditions support high-quality data
186 collection (Montaghi et al., 2013; Pirotti, 2010). Leaf-on conditions maximize canopy detection but may
187 reduce ground return availability, while occasional leaf-off campaigns improve ground model accuracy at
188 the cost of slightly lower canopy height estimates. In several regions, multitemporal acquisitions are
189 available due to overlapping campaigns, such as MATTM and regional surveys, especially in Liguria, Valle
d'Aosta, Molise, Piemonte, Trentino-Alto Adige, and Basilicata. Depending on the provider, the delivered
190 products ranged from raw point clouds to pre-processed Digital Terrain Models (DTMs) and Digital Surface
191 Models (DSMs) (Table S1).
192

193 Despite their high precision, the use of ALS data is constrained by the high cost of acquisition and the lack
194 of systematic temporal and spatial coverage. As a result, ALS availability is fragmented across both time
195 and geography. ALS data were available for a total area of approximately 16.86 million hectares, covering
196 about 56% of Italy's national territory (D'Amico et al., 2021) (Figure S2b). To ensure compatibility with
197 the Landsat-based modeling framework, the ALS canopy height measurements were aggregated to a 30 m
198 spatial resolution. Within each 30 m pixel, the maximum canopy height was extracted from the underlying
199 1 m ALS measurements, yielding a consistent reference layer for model training and validation. Given the
200 high vertical accuracy of ALS (~10–30 cm RMSE (Optech, n.d.)), resampling to 30 m is expected to
201 contribute only minor additional error, although local variability in acquisition settings (altitude, density,
202 season) may explain part of the residual uncertainty in model validation.

### 2.2.2 Global Ecosystem Dynamics Investigation Data (GEDI)

204 GEDI canopy height data were retrieved via GEE from the collection labeled
205 “LARSE/GEDI/GEDI02_A_002_MONTHLY”, which corresponds to the GEDI Version 2 Level 2A
206 Geolocated Elevation and Height Metrics Product. Version 2 includes improved waveform geolocation
207 algorithms and calibration relative to Version 1, leading to enhanced vertical accuracy and reduced
208 geolocation bias (Dubayah et al., 2021). Reported accuracy assessments indicate an average horizontal
209 geolocation error of 10.2 m and an elevation error of approximately 17.8 cm over gentle terrain (Xu et al.,
210 2023b).

211 We extracted data spanning 2019 to 2023, using the Relative Height at 98% energy return (RH98) metric
212 to represent maximum canopy height, which match better with the NFI maximum heights (Besic et al.,
213 2025; Su et al., 2025). To ensure high data quality, we applied a series of rigorous filters (Schwartz et al.,
214 2023; Su et al., 2025). First, to minimize errors associated with solar contamination, we included only
215 nighttime acquisitions, defined as observations where solar elevation was below 0 degrees. Second, to
216 reduce the influence of terrain-induced uncertainty—known to affect LiDAR signal accuracy on steep
217 slopes (Kutchartt et al., 2022), we limited the inclusion of GEDI data to the regions with a slope of no more
218 than 10 degrees. To further avoid geolocation artifacts near forest–non-forest boundaries, we excluded GEDI
219 footprints within 25 m of forest edges, based on a morphological edge filter (focal_max) applied to the ESA
220 WorldCover 2020 tree cover classification (ESA WorldCover 10 m 2021 v200, 2024). This additional buffer
221 helped mitigate misclassification and spatial mismatch errors (Schleich et al., 2023). After all filtering steps,
222 a total of 1,485,449 GEDI footprints were retained for analysis (Figure S2a). Although GEDI data are
223 natively available at a 25 m spatial resolution, we resampled the RH98 values to 30 m to align with the
224 Landsat data used for model input. The code for downloading GEDI canopy height data is available via
225 Figshare (Code used in the study - "A fused canopy height map of Italy (2004–2024) from spaceborne and
226 airborne LiDAR, and Landsat via deep learning and Bayesian averaging ", 2025).

## 2.3 Ground Validation Data

### 2.3.1 Canopy Height Data from National Forest Inventory (NFI)

229 The NFI data used in this study were obtained from the Italian National Forest Inventory 2015 campaign
230 (INFC2015) (Gasparini and Di Cosmo, 2015; Gasparini and Floris, 2022; National Inventory of Forests
231 and forest Carbon pools - INFC, 2025), which was designed to provide consistent, spatially detailed
<figure><img src="image_7.png" /></figure>

232 information on forest structure and biomass across the country (Gasparini and Di Cosmo, 2015). INFC2015
233 employed a three-phase stratified sampling approach (Gasparini and Di Cosmo, 2015; Gasparini and Floris,
234 2022). In the first phase, a systematic grid of 1 km × 1 km cells is established over the entire national
235 territory. Within each cell, a single point is randomly selected and classified using high-resolution imagery
236 in land cover types. The second phase, carried out on a subsample of forested points identified in the first
237 phase, involves a qualitative survey that refines forest type classification. In the third phase, a random
238 subsample of 6,894 second phase points is selected for field survey. Each point was surveyed in the field
239 to measure forest attributes within nested circular design plots, centered on the selected points. Trees with
240 diameter at breast height (DBH) ≥ 4.5 cm were measured within a 4m radius subplot (approximately 50.27
241 m2), while trees with DBH ≥ 9.5 cm were recorded within a larger 13 m radius subplot (approximately 530
242 m2). The 13m subplot was also used to collect measurements of tree height, increment cores, stumps, and
243 lying deadwood. For each tree, species identity, canopy height, and health status were recorded, enabling
244 estimation of volume and biomass through species-specific allometric models. Beyond tree-level attributes,
245 the campaign also assessed forest regeneration, dead organic matter, and land use change, with a focus on
246 maintaining methodological consistency with the earlier INFC2005 inventory and aligning with
247 international forest monitoring standards (Gasparini and Floris, 2022).

248 In our analysis, we initially collected 217,421 records of canopy height for individual living trees from
249 forested sites surveyed in the INFC2015 dataset. Since each site included multiple tree measurements, we
250 selected the highest canopy height recorded at each location (based on geographic coordinates) to represent
251 the site-level maximum After processing, a total of 6,894 field observations were retained for model
252 validation (Figure S2c).

253 **2.3.2 Reference Data for Forest Disturbance**

254 The reference data (Monitoring clearcutting and subsequent rapid recovery in Mediterranean coppice
255 forests with Landsat time series. Zenodo repository. [dataset]., 2025; Chirici et al., 2020; Francini et al.,
256 2021, 2022) were derived through manual interpretation of Landsat imagery since 1999, following an
257 approach comparable to the TimeSync protocol (Cohen et al., 2010). Interpreters delineated the spatial
258 boundaries of clearcuts and recorded the corresponding year of harvest, using a minimum mapping unit of
259 0.1 ha. In addition to Landsat imagery, high-resolution aerial photographs accessed via a Web Map Service
260 (WMS) were used to support interpretation. Field visits were conducted in coordination with local forest
261 authorities to confirm harvest timing and extent. After quality control procedures (Chirici et al., 2020), the
262 resulting geodatabase was deemed an accurate record of clearcut location and timing from 1999 to 2016
263 (Monitoring clearcutting and subsequent rapid recovery in Mediterranean coppice forests with Landsat time
264 series. Zenodo repository. [dataset]., 2025; Chirici et al., 2020).

265

266 **3. Method**

267 Our workflow began by developing two UNET-Landsat based canopy height models using ALS and GEDI
268 data, respectively (Figure S3). Each model was used to generate annual 30 m resolution canopy height
269 maps across Italy from 2004 to 2024 (Figure 1b). To integrate their complementary strengths, we applied
270 Bayesian Model Averaging (BMA) to fuse the predictions (Figure 1a), producing a harmonized, spatially
271 consistent, and temporally continuous canopy height dataset. To demonstrate its utility, we used the annual
272 maps to identify canopy disturbances based on interannual height losses and analyzed post-disturbance
273 recovery patterns across Italian biomes over the two-decade period. A detailed description of the methods
274 is provided in the subsequent sections.
<figure><img src="image_8.png" /><figcaption>(a) BMA algorithm</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Canopy height map fusion with BMA</figcaption></figure>

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

Figure 1 | Algorithm of Bayesian Model Averaging. Plot a is the Flow chart of to search for the best local weights that can be used in BMA. Landsat imagery from 2015, combined with elevation and slope data, was input into both the GEDI-Landsat and ALS-Landsat UNET models to generate canopy height predictions across Italy. Predicted heights at NFI plot locations were extracted from both models and compared with observed NFI canopy heights. The BMA procedure (right panel) then optimized local weights for each province by minimizing the residuals between model predictions and NFI reference values. These weights were subsequently used to fuse the two UNET model outputs, ensuring accurate and region-specific canopy height estimates. Plot b is the workflow for generating the final canopy height product through BMA fusion from 2004 to 2024. Annual Landsat imagery, along with elevation and slope data, was used as input to two trained UNET models: one based on ALS reference data and the other on GEDI. Each model produced annual canopy height predictions across Italy. These predictions were then fused using spatially optimized BMA weights—previously derived from comparison with NFI reference data—to produce a temporally consistent, high-resolution canopy height dataset for the full 2004–2024 period.

275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
3.1 Canopy height mapping model

3.1.1 Model workflow

A UNET deep learning architecture, a fully convolutional neural network originally developed for biomedical image segmentation (Ronneberger et al., 2015), was employed to estimate annual canopy height across Italy from 2004 to 2024 at 30 m resolution. Characterized by its distinctive U-shaped structure, the UNET architecture is able to capture both global context and fine-grained spatial detail, making it particularly well-suited for remote sensing tasks. It facilitates the analysis of images at different scales, effectively capturing the spatial hierarchy and the relationship between adjacent pixels, even with limited
<figure><img src="image_11.png" /></figure>

297 data (Solórzano et al., 2021). This enables UNET to effectively identify and segment complex land cover
298 features such as forests, croplands, grasslands, human-built environments and lakes (Zhu et al., 2017). This
299 feature is particularly valuable for detecting canopy heights in various types of landscapes (Schwartz et al.,
300 2023).

301 To effectively model both spatial and temporal patterns of forest structure, we integrated ALS- or GEDI-
302 derived canopy height measurements from all available years with temporally corresponding Landsat
303 imagery during the training phase. Each canopy height reference point was matched with satellite input data
304 from the same acquisition year, ensuring that the model learned from consistent and synchronous
305 observations. This temporal alignment allowed the UNET models not only to recognize spatial patterns—
306 such as texture, structure, and contextual cues related to vegetation—but also to learn interannual changes
307 in canopy height. By incorporating this temporal dimension into the training process, the models gained the
308 ability to detect and generalize patterns of canopy growth, degradation, and regeneration over time, making
309 the approach well-suited for generating long-term, consistent canopy height estimates for forest monitoring.

310 The model architecture consisted of five down sampling and four up sampling blocks with skip connections, 311 batch normalization, and ReLU activation functions (Ronneberger et al., 2015). Each convolutional layer 312 used a 3×3 kernel size, and max-pooling was applied during down sampling to reduce spatial dimensions 313 (Ronneberger et al., 2015). The final output layer consisted of a single channel with a linear activation 314 function to predict continuous canopy height values.

315 **3.1.2 Training, validation and testing for ALS-Landsat and GEDI-Landsat based UNET models**

316 In our implementation, the UNET takes as input a multi-channel image stack composed of 6 Landsat surface
317 reflectance bands, 2 digital terrain bands, and 4 vegetation indices including NDVI, EVI, LSWI and NDWI.
318 The inclusion of these indices aims to improve the interannual stability of the model outputs and mitigate
319 the effects of noise caused by year-to-year quality variability in satellite observations, which can be
320 influenced by factors such as cloud cover, atmospheric conditions, and shadow effects in time-series
321 imagery (Zhu and Woodcock, 2012). Each input tile corresponds to a 1000×1000 pixels (30km×30km)
322 patch with those 12 bands. To increase data diversity and improve model robustness, each tile was randomly
323 cropped into smaller 256×256-pixel patches prior to being fed into the model (Zheng et al., 2020). To
324 enhance generalization and prevent overfitting, batch normalization and dropout regularization were
325 applied throughout the network. Dropout randomly disables a subset of neurons during training, effectively
326 performing model averaging and reducing overfitting risk (Srivastava et al., 2014). And the mean absolute
327 error (MAE) loss function was used to optimize the model by minimizing the difference between predicted
328 and observed canopy height values (Schwartz et al., 2023). The detailed procedure for training, validating,
329 and testing our canopy height models was outlined in Figure S3. The detailed model structure and model
330 parameters can be found in Table S1 and Table S2.

331 In this work, the UNET models were trained separately on ALS-derived or GEDI-derived canopy height
332 labels depending on the available time period. While the structure of both models remains identical, the
333 difference lies in the reference data employed. GEDI data provide extensive spatial coverage each year, yet
334 it does not provide spatially continuous height measurements due to its configuration and spanning from
335 2019 to 2023. Conversely, ALS data offer smaller yearly spatial coverage, but with a significantly higher
336 data density, and spanning from 2004 to 2017. Each dataset was randomly split into 3 parts at the input
337 patch level, 75% were assigned to the training set, 5% were utilized to the validation set for monitoring the
338 training progress of the UNET models (with training stopping once the validation loss converges), and the
339 remaining 20% were reserved as the out-of-box testing set to evaluate the final models' performances
340 (Figure S4). Model performance was quantified using three metrics: MAE, Root Mean Squared Error
341 (RMSE), and Coefficient of Determination (R²).

342 **3.2 Model fusion with Bayesian Model Averaging (BMA)**

343 **3.2.1 BMA description**
<figure><img src="image_12.png" /></figure>

344 The two UNET models developed in this study were trained on distinct reference datasets—ALS (2004–2017) and GEDI (2019–2023)—each characterized by different spatial and temporal properties. To harness the complementary strengths of both models and generate a spatially complete and temporally consistent canopy height dataset for Italy from 2004 to 2024, we adopted a Bayesian Model Averaging (BMA) fusion framework (Besic et al., 2025; Duan et al., 2007; Hoeting et al., 1999) (Figure 1a).

349 BMA is a probabilistic ensemble technique that combines predictions from multiple models by assigning weights based on their relative predictive performance (Hoeting et al., 1999). Under the BMA framework, it is assumed that each model's predictive error follows a Gaussian distribution centered around the true value. The posterior probability of each model—used as its weight—is computed based on how well the model's predictions align with a reference dataset (in our case, the Italian NFI data in 2015). This allows BMA to effectively quantify uncertainty and balance the strengths of different models (Besic et al., 2025; Chen et al., 2020; Duan et al., 2007; Hoeting et al., 1999).

356 In our implementation, BMA was used to compute a weighted average of ALS-Landsat and GEDI-Landsat based canopy height predictions at the provincial level. Italy comprises 107 provinces, and the optimal weights for each province were derived by comparing the model outputs at NFI plot locations within that specific province for the year 2015 against the corresponding observed canopy height measurements (Figure S5). This approach enables spatially adaptive weighting: in regions where one model underperforms, such as areas without ALS coverage or mountainous regions where GEDI accuracy may degrade (Kutchartt et al., 2022), the other model receives greater weight in the final fusion.

363 Mathematically, the fused canopy height $H_{BMA}$ can be written as:

$$H_{BMA} = w_{ALS, province i} \times H_{ALS} + w_{GEDI, province i} \times H_{GEDI}$$

365 Where $H_{ALS}$ and $H_{GEDI}$ are the predicted canopy heights from ALS-Landsat and GEDI-Landsat based model, respectively. The weight $w_{ALS, province i}$ and $w_{GEDI, province i}$ are the corresponding optimal weights in province $i$. The weights are derived from province-level agreement with NFI plots, and constrained by $w_{ALS, province i} + w_{GEDI, province i} = 1$. This fusion strategy leverages the strengths of both models while compensating for their individual limitations, producing a more reliable and consistent canopy height dataset across Italy. By dynamically combining the two models in a location-specific manner, the BMA fusion approach improves the robustness of the final canopy height estimates.

### 3.2.2 Canopy height map fusion with BMA

373 Once the optimal weights for each province were computed using 2015 NFI data, they were used to fuse annual predictions from both UNET models over the 2004–2024 period (Figure 1b). Specifically, first, canopy height maps for each year were produced independently using both the ALS-Landsat based and GEDI-Landsat based models. These were then fused using the BMA weights at the provincial level to create a harmonized annual time series.

### 3.3 Disturbance detection

379 Following the generation of the BMA-fused annual canopy height maps, we derived a canopy height reduction-based disturbance product by identifying substantial interannual reductions in canopy height of at least 5 meters, signaling potential forest disturbances. This threshold was chosen based on the model’s RMSE (Figure 2) to balance detection sensitivity and reliability. While applying a larger threshold enhances the confidence in detected disturbances, it may also result in fewer events being identified.

384 A disturbance was defined as a 30 × 30 m pixel—or cluster of adjacent pixels—exhibiting a canopy height loss greater than 5 meters between a reference year (t) and the following year (t+1). To reduce false detections caused by Landsat image artifacts such as cloud cover, seasonal snow, or shadow effects, the height reduction was required to persist into the second subsequent year (t+2). This persistence criterion assumes that a true maximum canopy height reduction of more than 5 meters is unlikely to regenerate within one year, thereby enhancing the robustness of disturbance detection.
<figure><img src="image_13.png" /></figure>

390 It is important to note that disturbances smaller than a Landsat pixel (30 m) or partial canopy removals that do not reduce the maximum canopy height by more than 5 m cannot be consistently detected by our approach. As a result, selective logging, thinning, or low-intensity disturbances may be underrepresented in our disturbance product. This limitation is inherent to the spatial resolution of Landsat imagery and our maximum-height-based detection strategy. Nevertheless, the method captures stand-replacing disturbances and large-scale partial events that are most relevant for national forest monitoring and carbon accounting. For finer-scale management applications, integration with higher-resolution data (e.g., Sentinel-2) would be required.

398 **3.4 External validation**

399 **3.4.1 External Validation of canopy height product**

400 To assess the generalizability of the ALS-Landsat based and GEDI-Landsat based UNET models and the BMA fusion model, we conducted external validation using canopy height observations from the last available NFI data of Italy (Gasparini and Di Cosmo, 2015; Gasparini and Floris, 2022; National Inventory of Forests and forest Carbon pools - INFC, 2025). Each model was used to generate a canopy height map of Italy for the year 2015 (Figure S6). Model predictions were extracted at the geographic coordinates of NFI plots and compared to the corresponding maximum canopy height values reported in the NFI dataset (Figure 2). This independent assessment enabled us to evaluate how well the models generalize to field-based measurements, thereby supporting their applicability for operational forest monitoring and ecosystem modeling.

409 To evaluate the accuracy of BMA fusion model across other years, we randomly sampled maximum 10,000 non-zero pixels from the ALS and GEDI testing datasets for each year. The corresponding BMA predictions were extracted and compared with sampled ALS or GEDI values.

412 We further evaluated the BMA canopy height product by comparing the BMA estimates against the global canopy height map for 2020 by Lang et al. (Lang et al., 2023), and the European-scale products for 2019 by Turubanova et al. (Turubanova et al., 2023) and by Liu et al. (Liu et al., 2023). To evaluate agreement and consistency, we randomly sampled 50,000 non-zero pixels across Italy and compared canopy height values across corresponding locations in each product.

<figure><img src="image_14.png" /><figcaption>(b) GEDI-Landsat model performance</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) BMA fusion model performance</figcaption></figure>

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

418 **Figure 2 | External validation of ALS-Landsat, GEDI-Landsat based UNET model, and BMA fusion model using NFI data in 2015.** Plot a is scatterplot of predicted canopy heights of our ALS-Landsat based UNET canopy height model versus the canopy height of NFI observations in 2015. Plot b is scatterplot of predicted canopy heights of our GEDI-Landsat based UNET canopy height model versus the canopy height of NFI observations in 2015. Plot c is scatterplot of predicted canopy heights of our BMA fusion canopy height model versus the canopy height of NFI observations in 2015.
<figure><img src="image_17.png" /></figure>

424 **3.4.2 External Validation of height derived disturbance product**

425 To evaluate the accuracy of the canopy height change-based disturbance detection, we compared our
426 disturbance product against reference disturbance datasets (Monitoring clearcutting and subsequent rapid
427 recovery in Mediterranean coppice forests with Landsat time series. Zenodo repository. [dataset]., 2025;
428 Chirici et al., 2020; Francini et al., 2021, 2022) within the survey area. We assessed performance using
429 precision, defined as the proportion of detected disturbances that were confirmed by the reference data. A
430 disturbance identified by our product was considered a true positive (TP) if it also appeared in the reference
dataset; otherwise, it was classified as a false positive (FP). Precision was then calculated as: Precision =
431 TP / (TP + FP).
432

433 We also compared our disturbance product with those developed by Saverio et al. (Francini and Chirici,
434 2022), as well as the product presented in the studies by Viana-Soto et al. (Viana-Soto and Senf, 2024). To
435 assess the precision, a random sample of 1,000 discrete disturbance patches (each >0.5 ha) from our map
436 was compared with those identified by reference products (Francini and Chirici, 2022; Viana-Soto and Senf,
437 2024) over a comparable time window. Similarly, we performed a reciprocal comparison by random
438 sampling 1,000 disturbances from reference products and evaluating their presence in our results. We
439 calculated the proportion of events jointly detected, uniquely detected by our product, and uniquely detected
440 by the reference datasets. This multi-source validation helps demonstrate the reliability of our disturbance
441 mapping approach, while avoiding dependence on a single reference dataset.
442

442

443 **4 Results**

444 **4.1 Italy canopy height map at 30m resolution covering the last two decades**

445 We produced 30 m resolution canopy height maps of Italy over 2004-2024 using BMA to fuse predictions
446 from two UNET models trained on ALS-Landsat and GEDI-Landsat data. The resulting national-scale map
447 of 2024 is shown in **Figure 3a**, while plots b-g illustrate the map's ability to capture detailed forest
448 landscape features across diverse land use types, as verified against Landsat imagery. These examples
449 highlight the model's capacity to detect canopy height variation in a range of environments.

450 Using this canopy height map, we also derived a tree line map (Figure S7), defined as the highest elevation
451 at which trees are present. Our results indicate that tree lines in Italy typically range between 1800 m and
452 2500 m above sea level, with higher elevations observed in the northern Alps compared to the southern
453 Apennines and volcanic slopes of southern Italy, such as Mount Etna in Sicily. These spatial variations
454 likely reflect a combination of climatic gradients, species composition, and local environmental conditions,
455 consistent with established knowledge on tree line ecology (Körner, 2012).
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Figure 3 | Canopy height map of Italy at 30m resolution for the year 2024. Plot a is the nationwide spatial distribution of canopy height. Plots b-g provide detailed comparisons between Landsat satellite imagery and predicted canopy height across different environments.</figcaption></figure>

To evaluate the performance of the canopy height estimation, we performed external validation using independent canopy height measurements from NFI. Predictions in 2015 from the ALS-Landsat based, GEDI-Landsat based, and BMA-fused models (Figure S6) were compared against maximum tree height recorded at NFI plots. The ALS-Landsat model achieved an MAE of 4.17m and R² of 0.4 (Figure 2a). The GEDI-Landsat model showed slightly lower performance, with an MAE of 4.33m and a R² of 0.37 (Figure 2b). The BMA-fused model outperformed both individual models, with an MAE of 3.98 m and R² of 0.46 (Figure 2c). We further assessed the performance of the BMA product by comparing it directly against the random sampled non-zero observations in the ALS and GEDI testing datasets. For each year, up to 10,000 observations were sampled, depending on data availability. Against ALS data in 2004-2017, the BMA model achieved an MAE of 4.09m and R² of 0.55 (Figure S8a). When compared to GEDI data in 2019-2023, the model produced a MAE of 5.11m and R² of 0.40 (Figure S8b). These evaluations confirm that the BMA fusion strategy provides balanced accuracy across multiple lidar sources.

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471
<figure><img src="image_20.png" /></figure>

472 To investigate how terrain complexity influences prediction performance, we analyzed the BMA model's errors relative to slope and elevation classes, using the random sampled ALS samples (Figure S9).
473 Prediction errors were generally low on gentle slopes, but increased progressively in steeper terrain, e.g., slope beyond 35°, reflecting reduced model accuracy in rugged mountainous areas (Figure S9a). A similar pattern was observed for elevation: errors remained low and stable up to 2000 m, but increased at higher altitudes (Figure S9b). These findings underscore the difficulty of accurately estimating canopy height in complex high-mountain environments, especially above 2000 m or on slopes exceeding 35°.

479 We further assessed our BMA canopy height estimates by comparing them with three independently developed large-scale datasets: the European-scale products from Turubanova et al. (Turubanova et al., 2023) and Liu et al. (Liu et al., 2023) in 2019, and the global canopy height map from Lang et al. (Lang et al., 2023) in 2020. Our BMA map in 2019 showed an overall good alignment with Liu's dataset (Figure 4a, b, e, h), with a small negative bias of -0.92 m. In contrast, when compared to Turubanova's product (Figure 4a, b, d, g), our BMA estimates were consistently higher, with an average positive bias of 5.16 m—indicating that Turubanova's model tends to underestimate canopy height relative to ours. Similarly, comparison with the 2020 Lang's product revealed that their model generally overestimates canopy height for trees taller than 10 m, resulting in an average negative bias of 4.19 m relative to our BMA estimates (Figure 4h).

489 To further examine these discrepancies, we compared each external product against independent GEDI canopy height observations (Figure S10). We randomly sampled 10,000 non-zero GEDI footprints from 2019 and 2020 and extracted corresponding canopy height estimates from Liu and Turubanova (2019) and Lang (2020). Liu's product showed a good agreement with GEDI, with a minor bias of 0.67 m. In contrast, Turubanova's estimates exhibited a substantial negative bias of -5.37 m, indicating consistent underestimation, while Lang's product showed a positive bias of +4.50 m, suggesting overestimation. Additionally, since Turubanova's dataset also includes 2015 canopy height estimates, we compared them with field measurements from the Italian National Forest Inventory (NFI). The results mirrored those from the GEDI comparison, showing an overall bias of -5.75 m (Figure S11), further indicating underestimation. Together, these results indicate that Liu's regional dataset aligns with GEDI and our BMA product, whereas Turubanova and Lang's models appear less accurate in capturing forest structure in Italy. Overall, the comparison supports the robustness of our BMA canopy height estimates at the national scale.
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 4** | Comparison between BMA products with external products from Turubanova et al. (Turubanova et al., 2023), Liu et al. (Liu et al., 2023) and Lang et al. (Lang et al., 2023). Plot a is the Google Maps (© Google Maps 2025) reference image in forested region in southern Italy; Plots b and c are the predicted canopy height from our BMA at this area in 2019 and 2020, respectively. Plot d is the predicted canopy height of Turubanova et al. at this area. Plot e is the predicted canopy height of Liu et al. at this area. Plot f is the predicted canopy height of Lang et al. at this area. Plot g is the scatter plot of BMA predicted canopy heights versus the canopy height from Turubanova in 2019. Plot h is the scatter plot of BMA predicted canopy heights versus the canopy height from Liu in 2019. Plot i is the scatter plot of BMA predicted canopy heights versus the canopy height from Lang in 2019.</figcaption></figure>

501

502

503

504

505

506

507

508

509

510

511

## 4.2 Disturbance map

512

513

To further validate our canopy height maps and demonstrate their practical application, we generated a disturbance map based on interannual canopy height changes. The accuracy of this product was assessed
<figure><img src="image_23.png" /></figure>

514 by comparing it with observed reference dataset (Monitoring clearcutting and subsequent rapid recovery in
515 Mediterranean coppice forests with Landsat time series. Zenodo repository. [dataset]., 2025; Chirici et al.,
516 2020; Francini et al., 2021, 2022), as shown in **Figure 5**. We showed that the range of precision is from
517 0.33 to 0.87 from 2005 to 2016, with an average precision of 0.64. The precision of disturbance detection
518 ranged from 0.33 to 0.87 between 2005 and 2016, with an average precision of 0.64. Variability in precision
519 is likely influenced by fluctuations in Landsat image quality, particularly due to the presence of clouds,
520 shadows, and snow (Zhu and Woodcock, 2014).

<figure><img src="image_24.png" /><figcaption>521
**Figure 5 | Validation of canopy height change derived forest disturbance detection using independent**
522 observations (2005–2016). Plot a is spatial distribution of the reference observation areas across Italy used
523 for evaluating canopy disturbance detection. Plot b and c are the comparison of BMA-predicted
524 disturbances (in red) with reference disturbances (in yellow) for the years 2006 and 2016 within selected
525 observation regions. Plot d is the annual precision of BMA-based disturbance detection from 2005 to 2016,
526 calculated as the proportion of correctly identified disturbances relative to the reference data.
527</figcaption></figure>

528 We further compared our canopy height change-based disturbance product with two independently
529 developed datasets (Francini and Chirici, 2022; Viana-Soto and Senf, 2024). **Figure S12** shows our
530 disturbance map in 2022-2023 alongside the Landsat-based product from Viana-Soto (Viana-Soto and Senf,
531 2024), while **Figure S13** presents a comparison for 2018–2019 with a Sentinel-2-based map from Saverio
532 (Francini and Chirici, 2022). Although all three maps rely on satellite imagery, they differ in methodology:
<figure><img src="image_25.png" /></figure>

533 our approach detects structural changes via canopy height reduction (>5 m), whereas the methods used by Viana-Soto and Saverio rely on spectral-temporal segmentation to classify areas as disturbed or undisturbed, with Viana-Soto further attributing disturbances to specific events such as fire, windthrow, insect outbreaks, or harvesting. Visual assessments across multiple sites (Figure Ss 12c-i, 13c-i) show broad agreement.

537 In quantitative evaluations, our disturbance map showed 59.9% overlap with Viana-Soto's product across 2,000 randomly sampled disturbance patches (≥0.5 ha), while 14.0% were uniquely detected by us and 26.1% only by Viana-Soto (Figure S14b). Similarly, comparison with Saverio's Sentinel-2-based product showed 54.0% agreement, with 24.4% detected only by us and 21.6% only by Saverio (Figure S14a). These discrepancies are probably linked to a combination of commission (false positives) and omission (false negatives) errors in all datasets. Additionally, it is important to note that not all disturbance types—such as low-intensity fires, selective logging, insect outbreaks, or windthrow—necessarily result in detectable canopy height reductions exceeding 5 m. As a result, some disturbances captured by spectral changes may not produce sufficient structural changes to be identified by our height-based approach. Conversely, our method may capture height changes from subtle or mixed disturbances that lack strong spectral signals. Integrating both structural (e.g., canopy height variation) and spectral information could enable the development of an ensemble or multi-sensor disturbance detection framework, potentially improving detection accuracy and robustness across diverse disturbance types.

550 Across each forest biome in Italy (Figure S1), we tracked post-disturbance recovery in 100 random disturbance patches from 2004–2005. Results show that canopy height recovery is most rapid in the initial years following disturbance, with the rate of regrowth gradually slowing over time. Among the biomes, Mediterranean forests exhibited the fastest recovery, followed by temperate broadleaf and mixed forests, while temperate conifer forests showed the slowest regrowth (Figure 6). The ability to detect disturbance and quantify subsequent regrowth underscores the value of our annual canopy-height series for assessing forest resilience and informing long-term monitoring and management.

<figure><img src="image_26.png" /><figcaption>Figure 6 | Post-disturbance canopy height recovery trajectories (2004–2024) in different forest biomes in Italy. The figure illustrates the recovery of canopy height over time in three forest biomes following disturbance events that occurred between 2004 and 2005. For each biome—(plot a) Mediterranean Forests, Woodlands & Scrub, (plot b) Temperate Broadleaf & Mixed Forests, and (plot c) Temperate Conifer Forests—100 disturbance patches were randomly sampled, and their 5-year average canopy height was tracked through 2024.</figcaption></figure>

564

565 **5 Limitations and further improvements**

566 While this study presents a novel, long-term, high-resolution canopy height and disturbance dataset for Italy,
567 several limitations should be acknowledged:
<figure><img src="image_27.png" /></figure>

568 **5.1 Uncertainty in extreme terrain conditions:**

569 Despite incorporating elevation and slope data into the model, prediction errors increased in high-slope and

570 high-altitude regions. This reflects persistent challenges in accurately modeling canopy height in rugged

571 terrain, where terrain-induced distortions and LiDAR signal occlusions affect both training and validation.

572 **5.2 Dependence on Landsat image quality:**

573 Although cloud and artifact filtering were applied using QA bands and additional masking techniques,

574 residual cloud contamination and seasonal snow may still influence Landsat reflectance values. These

575 factors can lead to erroneous canopy height predictions or false disturbance detections in isolated cases.

576 **5.3 BMA fusion relies on 2015 NFI reference only:**

577 The BMA fusion weights were calibrated using NFI plot data from 2015. While this provides a reliable

578 spatial benchmark, the assumption that relative model performance is stable across all years may introduce

579 bias in time periods far from the reference year.

580 **5.4 Data limitations**

581 Post-disturbance canopy height recovery trajectories (e.g., Figure 6) could not be validated against

582 remeasured NFI plots, as repeated inventory measurements capturing regrowth over time were unavailable.

583 Consequently, the recovery curves in this study should be regarded as model-based trajectories rather than

584 empirically validated growth signals. This is an important limitation that future work could address by

585 incorporating multi-temporal ground inventory data or harmonized European networks to provide direct

586 benchmarks for post-disturbance recovery.

587 **5.5 Resolution limitations**

588 The disturbance detection framework is constrained by the 30 m spatial resolution of Landsat imagery and

589 the requirement of a ≥5 m canopy height reduction. As a result, sub-pixel or low-intensity disturbances—

590 such as selective logging or thinning—cannot be reliably detected (see Section 3.3).

591

592 **6 Data availability**

593 The maps from 2004 to 2024 (https://doi.org/10.5281/zenodo.15627897), along with the covering 2005 to

594 2023 (https://doi.org/10.5281/zenodo.15627927), are publicly available on Zenodo (Canopy height map in

595 Italy - 2004-2024 [Dataset]. Zenodo., 2025; Canopy height change derived disturbance map in Italy - 2005-

596 2023 [Dataset]. Zenodo., 2025).

597

598 **7 Conclusions**

599 We present a long-term (2004–2024), 30 m resolution canopy height dataset for Italy, developed by fusing

600 ALS-Landsat and GEDI-Landsat based UNET models through a Bayesian model averaging framework.

601 The fusion approach balances the strengths of each data source, producing spatially continuous and

602 temporally consistent canopy height estimates. Using interannual canopy height change, we also generated

603 annual disturbance maps from 2005 to 2023. Validation against NFI data shows robust accuracy (R² = 0.46,

604 MAE = 3.98 m in 2015), while comparison with observed events indicates reliable detection of stand-

605 replacing disturbances (≥5 m canopy height loss, average precision = 0.64). The dataset also reveals biome-

606 specific recovery trajectories, with faster regrowth in Mediterranean than in temperate coniferous forests.

607 This dataset provides the first spatially consistent, two-decade record of canopy height dynamics in Italy,

608 supporting national forest monitoring, carbon accounting, climate reporting, and ecological studies on

609 disturbance and resilience. Two main limitations remain: (i) sub-pixel or low-intensity disturbances (e.g.,

610 selective logging, thinning) are not detectable at 30 m resolution, and (ii) post-disturbance recovery
<figure><img src="image_28.png" /></figure>

611 trajectories are model-based and lack validation with repeated NFI plots. Addressing these gaps will require

612 higher-resolution remote sensing (e.g., Sentinel-2, UAV lidar) and multi-temporal ground inventories.

613 The framework is scalable to continental or global domains by leveraging upcoming lidar missions (e.g.,

614 GEDI follow-on, ICESat-2 synergy) and dense optical archives. All datasets are openly available on Zenodo,

615 offering a valuable resource for forest monitoring, carbon accounting, and land-use change studies.

616

617 Code Availability

618 The code that used in this study can be found in Figshare (https://doi.org/10.6084/m9.f figshare.29416658)

619 (Code used in the study - "A fused canopy height map of Italy (2004–2024) from spaceborne and airborne

620 LiDAR, and Landsat via deep learning and Bayesian averaging ", 2025).

621

622 Disclaimer

623 Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the

624 text, published maps, institutional affiliations, or any other geographical representation in this paper. While

625 Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies

626 with the authors.

627

628 Acknowledgement

629 This work was supported by the ANR under the AI4Forests program with the reference ANR-22-FAI1-

630 0002. YS, PC, MS, IF and AD are supported by the French German project AI4FOREST (ANR-22-FAI1-

631 0002-01) funded by ANR and DLR. This research includes data from the Global Ecosystem Dynamics

632 Investigation (GEDI), distributed by NASA’s Land Processes Distributed Active Archive Center (LP

633 DAAC).

634 This study was also partially supported via the following projects: 1. MULTIFOR “Multi-scale observations

635 to predict Forest response to pollution and climate change” PRIN 2020 Research Project of National

636 Relevance funded by the Italian Ministry of University and Research (prot. 2020E52THS). 2. SUPERB

637 "Systemic solutions for upscaling of urgent ecosystem restoration for forest related biodiversity and

638 ecosystem services" H2020 project funded by the European Commission, number 101036849 call LC-GD-

639 7-1-2020. 3. FORWARDS "The ForestWard Observatory to Secure Resilience of European Forests",

640 HORIZON project funded by the European Commission, number 101084481 call HORIZON-CL6-2022-

641 CLIMATE-01. 4. MONIFUN. H2020 project funded by the European Commission, number 101134991

642 call HORIZON-CL6-2023-CircBio-01-14. 5. NextGenCarbon. H2020 project funded by the European

643 Commission, number 101184989 call HORIZON-CL5-2024-D1-01-07. NBFC “National Biodiversity

644 Future Center” project funded under the National Recovery and Resilience Plan (NRRP),

645 NextGenerationEU; Project code CN_00000033, adopted by the Italian Ministry of University and

646 Research, CUP B83C22002930006. 6. Space It Up. Project funded by the Italian Space Agency and the

647 Ministry of University and Research - Contract No. 2024-5-E.0 - CUP No. 153D24000060005. 7.

648 Interdisciplinary program ARTEMIS of Lorraine Université d'Excellence (ANR-15-IDEX-04-LUE).

649 This research uses Landsat 5 Thematic Mapper, Landsat 7 Enhanced Thematic Mapper Plus, and Landsat

650 8 Operational Land Imager data distributed by the U.S. Geological Survey (USGS), and Global Ecosystem

651 Dynamics Investigation (GEDI) data distributed by NASA’s Land Processes Distributed Active Archive

652 Center (LP DAAC). All datasets were accessed and processed via Google Earth Engine. The ALS and NFI

653 data were acquired under the help of Dr. Saverio Francini and Dr. Giovanni Forzieri. Other datasets used in

654 this study are open-access and freely available, including those from GEDI, Landsat, Liu, Lang, Turubanova,

655 Saverio, and Viana-Soto.
<figure><img src="image_29.png" /></figure>

656 We acknowledge the use of Google Maps satellite imagery (© Google Maps 2025), accessed via QGIS,

657 which was used solely for academic, non-commercial purposes to support geographic visualization and

658 spatial referencing in this study.

659

660 **Author Contributions**

661 Conceptualization: YS, MS, IF, PC

662 Methodology: YS, MS, IF, PC, NB

663 Investigation: YS, AP, SB

664 Visualization: YS, XZ, YX, KY, HC, SC

665 Data acquisition: YS, SF, GF

666 Funding acquisition: PC, AA

667 Supervision: YS, PC, AA

668 Writing—original draft: YS, PC

669 Writing—review & editing: YS, MS, IF, PC, NB, XZ, SC, HC, AA, SB, AP, KY, SF, GF

670

671 **Competing Interests.**

672 None declared.
<figure><img src="image_30.png" /></figure>

Reference

Besic, N., Picard, N., Vega, C., Bontemps, J.-D., Hertzog, L., Renaud, J.-P., Fogel, F., Schwartz, M., Pellissier-Tanon, A., Destouet, G., Mortier, F., Planells-Rodriguez, M., and Ciais, P.: Remote-sensing-based forest canopy height mapping: some models are useful, but might they provide us with even more insights when combined?, Geosci Model Dev, 18, 337–359, https://doi.org/10.5194/gmd-18-337-2025, 2025.

Budei, B. C., St-Onge, B., Hopkinson, C., and Audet, F.-A.: Identifying the genus or species of individual trees using a three-wavelength airborne lidar system, Remote Sens Environ, 204, 632–647, https://doi.org/10.1016/j.rse.2017.09.037, 2018.

Canopy height change derived disturbance map in Italy - 2005-2023 [Dataset]. Zenodo.: https://doi.org/10.5281/zenodo.15627927, last access: 30 June 2025.

Canopy height map in Italy - 2004-2024 [Dataset]. Zenodo.: https://doi.org/10.5281/zenodo.15627897, last access: 30 June 2025.

Chave, J., Réjou-Méchain, M., Búrquez, A., Chidumayo, E., Colgan, M. S., Delitti, W. B. C., Duque, A., Eid, T., Fearnside, P. M., Goodman, R. C., Henry, M., Martínez-Yrizar, A., Mugasha, W. A., Muller-Landau, H. C., Mencuccini, M., Nelson, B. W., Ngomanda, A., Nogueira, E. M., Ortiz-Malavassi, E., Pélissier, R., Ploton, P., Ryan, C. M., Saldarriaga, J. G., and Vieilledent, G.: Improved allometric models to estimate the aboveground biomass of tropical trees, Glob Chang Biol, 20, 3177–3190, https://doi.org/10.1111/gcb.12629, 2014.

Chen, S., Mulder, V. L., Heuvelink, G. B. M., Poggio, L., Caubet, M., Román Dobarco, M., Walter, C., and Arrouays, D.: Model averaging for mapping topsoil organic carbon in France, Geoderma, 366, 114237, https://doi.org/10.1016/j.geoderma.2020.114237, 2020.

Chirici, G., Giannetti, F., Mazza, E., Francini, S., Travaglini, D., Pegna, R., and White, J. C.: Monitoring clearcutting and subsequent rapid recovery in Mediterranean coppice forests with Landsat time series, Ann For Sci, 77, 40, https://doi.org/10.1007/s13595-020-00936-2, 2020.

Code used in the study - "A fused canopy height map of Italy (2004–2024) from spaceborne and airborne LiDAR, and Landsat via deep learning and Bayesian averaging": https://doi.org/10.6084/m9 figshare.29416658, last access: 3 July 2025.

Cohen, W. B., Yang, Z., and Kennedy, R.: Detecting trends in forest disturbance and recovery using yearly Landsat time series: 2. TimeSync — Tools for calibration and validation, Remote Sens Environ, 114, 2911–2924, https://doi.org/10.1016/j.rse.2010.07.010, 2010.

Coops, N. C., Tompalski, P., Goodbody, T. R. H., Queinnec, M., Luther, J. E., Bolton, D. K., White, J. C., Wulder, M. A., van Lier, O. R., and Hermosilla, T.: Modelling lidar-derived estimates of forest attributes over space and time: A review of approaches and future trends, Remote Sens Environ, 260, 112477, https://doi.org/10.1016/j.rse.2021.112477, 2021.

Crawford, C. J., Roy, D. P., Arab, S., Barnes, C., Vermote, E., Hulley, G., Gerace, A., Choate, M., Engebretson, C., Micijevic, E., Schmidt, G., Anderson, C., Anderson, M., Bouchard, M., Cook, B., Dittmeier, R., Howard, D., Jenkerson, C., Kim, M., Kleyians, T., Maiersperger, T., Mueller, C., Neigh, C., Owen, L., Page, B., Pahlevan, N., Rengarajan, R., Roger, J.-C., Sayler, K., Scaramuzza, P., Skakun, S., Yan, L., Zhang, H. K., Zhu, Z., and Zahn, S.: The 50-year Landsat collection 2 archive, Science of Remote Sensing, 8, 100103, https://doi.org/10.1016/j.srs.2023.100103, 2023.

D'Amico, G., Vangi, E., Francini, S., Giannetti, F., Nicolaci, A., Travaglini, D., Massai, L., Giambastiani, Y., Terranova, C., and Chirici, G.: Are we ready for a National Forest Information System? State of the art of forest maps and airborne laser scanning data availability in Italy, IForest, 14, 144–154, https://doi.org/10.3832/ifor3648-014, 2021.
<figure><img src="image_31.png" /></figure>

Duan, Q., Ajami, N. K., Gao, X., and Sorooshian, S.: Multi-model ensemble hydrologic prediction using Bayesian model averaging, Adv Water Resour, 30, 1371–1386, https://doi.org/10.1016/j.advwatres.2006.11.014, 2007.

Dubayah, R., Blair, J. B., Goetz, S., Fatoyinbo, L., Hansen, M., Healey, S., Hofton, M., Hurtt, G., Kellner, J., Luthcke, S., Armston, J., Tang, H., Duncanson, L., Hancock, S., Jantz, P., Marselis, S., Patterson, P. L., Qi, W., and Silva, C.: The Global Ecosystem Dynamics Investigation: High-resolution laser ranging of the Earth's forests and topography, Science of Remote Sensing, 1, 100002, https://doi.org/10.1016/j.srs.2020.100002, 2020.

Dubayah, R., Hofton, M., Blair, J., Armston, J., Tang, H., and Luthcke, S.: GEDI L2A Elevation and Height Metrics Data Global Footprint Level V002 [Data set], NASA Land Processes Distributed Active Archive Center, 2021.

Fayad, I., Ciais, P., Schwartz, M., Wigneron, J.-P., Baghdadi, N., de Truchis, A., d'Aspremont, A., Frappart, F., Saatchi, S., Sean, E., Pellissier-Tanon, A., and Bazzi, H.: Hy-TeC: a hybrid vision transformer model for high-resolution and large-scale mapping of canopy height, Remote Sens Environ, 302, 113945, https://doi.org/10.1016/j.rse.2023.113945, 2024.

Francini, S. and Chirici, G.: A Sentinel-2 derived dataset of forest disturbances occurred in Italy between 2017 and 2020, Data Brief, 42, 108297, https://doi.org/10.1016/j.dib.2022.108297, 2022.

Francini, S., McRoberts, R. E., Giannetti, F., Mencucci, M., Marchetti, M., Scarascia Mugnozza, G., and Chirici, G.: Near-real time forest change detection using PlanetScope imagery, Eur J Remote Sens, 53, 233–244, https://doi.org/10.1080/22797254.2020.1806734, 2020.

Francini, S., McRoberts, R. E., Giannetti, F., Marchetti, M., Scarascia Mugnozza, G., and Chirici, G.: The Three Indices Three Dimensions (313D) algorithm: a new method for forest disturbance mapping and area estimation based on optical remotely sensed imagery, Int J Remote Sens, 42, 4693–4711, https://doi.org/10.1080/01431161.2021.1899334, 2021.

Francini, S., D'Amico, G., Vangi, E., Borghi, C., and Chirici, G.: Integrating GEDI and Landsat: Spaceborne Lidar and Four Decades of Optical Imagery for the Analysis of Forest Disturbances and Biomass Changes in Italy, Sensors, 22, 2015, https://doi.org/10.3390/s22052015, 2022.

Francini, S., Hermosilla, T., Coops, N. C., Wulder, M. A., White, J. C., and Chirici, G.: An assessment approach for pixel-based image composites, ISPRS Journal of Photogrammetry and Remote Sensing, 202, 1–12, https://doi.org/10.1016/j.isprsjprs.2023.06.002, 2023.

Gasparini, P. and Di Cosmo, L.: Forest carbon in Italian forests: Stocks, inherent variability and predictability using NFI data, For Ecol Manage, 337, 186–195, https://doi.org/10.1016/j.foreco.2014.11.012, 2015.

Gasparini, P. and Floris, A.: Definitions and Sampling Design, 17–48, https://doi.org/10.1007/978-3-030-98678-0_2, 2022.

ee.Terrain.slope: https://developers.google.com/earth-engine/apidocs/ee-terrain-slope, last access: 8 June 2025.

Landsat Algorithms: https://developers.google.com/earth-engine/guides/landsat#landsat-processing-methods, last access: 8 June 2025.

Landsat ETM+ to OLI Harmonization: https://developers.google.com/earth-engine/tutorials/community/landsat-etm-to-oli-harmonization, last access: 8 June 2025.

Hansen, M. C., Potapov, P. V., Moore, R., Hancher, M., Turubanova, S. A., Tyukavina, A., Thau, D., Stehman, S. V., Goetz, S. J., Loveland, T. R., Kommareddy, A., Egorov, A., Chini, L., Justice, C. O., and
<figure><img src="image_32.png" /></figure>

Townshend, J. R. G.: High-Resolution Global Maps of 21st-Century Forest Cover Change, Science (1979), 342, 850–853, https://doi.org/10.1126/science.1244693, 2013a.

Hansen, M. C., Potapov, P. V., Moore, R., Hancher, M., Turubanova, S. A., Tyukavina, A., Thau, D., Stehman, S. V., Goetz, S. J., Loveland, T. R., Kommareddy, A., Egorov, A., Chini, L., Justice, C. O., and Townshend, J. R. G.: High-Resolution Global Maps of 21st-Century Forest Cover Change, Science (1979), 342, 850–853, https://doi.org/10.1126/science.1244693, 2013b.

Hermosilla, T., Wulder, M. A., White, J. C., Coops, N. C., and Hobart, G. W.: An integrated Landsat time series protocol for change detection and generation of annual gap-free surface reflectance composites, Remote Sens Environ, 158, 220–234, https://doi.org/10.1016/j.rse.2014.11.005, 2015.

Hoeting, J. A., Madigan, D., Raftery, A. E., and Volinsky, C. T.: Bayesian Model Averaging: A Tutorial, Statistical Science, 14, 382–401, 1999.

Kaewmanee, M., Leigh, L., Shah, R., and Gross, G.: Inter-Comparison of Landsat-8 and Landsat-9 during On-Orbit Initialization and Verification (OIV) Using Extended Pseudo Invariant Calibration Sites (EPICS): Advanced Methods, Remote Sens (Basel), 15, 2330, https://doi.org/10.3390/rs15092330, 2023.

KC, K., Zhao, K., Romanko, M., and Khanal, S.: Assessment of the Spatial and Temporal Patterns of Cover Crops Using Remote Sensing, Remote Sens (Basel), 13, 2689, https://doi.org/10.3390/rs13142689, 2021.

Kennedy, R. E., Yang, Z., and Cohen, W. B.: Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrend — Temporal segmentation algorithms, Remote Sens Environ, 114, 2897–2910, https://doi.org/10.1016/j.rse.2010.07.008, 2010.

Körner, C.: Alpine Treelines, Springer Basel, Basel, https://doi.org/10.1007/978-3-0348-0396-0, 2012.

Kovalskyy, V. and Roy, D. P.: The global availability of Landsat 5 TM and Landsat 7 ETM+ land surface observations and implications for global 30m Landsat data product generation, Remote Sens Environ, 130, 280–293, https://doi.org/10.1016/j.rse.2012.12.003, 2013.

Kutchartt, E., Pedron, M., and Pirotti, F.: ASSESSMENT OF CANOPY AND GROUND HEIGHT ACCURACY FROM GEDI LIDAR OVER STEEP MOUNTAIN AREAS, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, V-3-2022, 431–438, https://doi.org/10.5194/isprs-annals-V-3-2022-431-2022, 2022.

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

Liu, S., Brandt, M., Nord-Larsen, T., Chave, J., Reiner, F., Lang, N., Tong, X., Ciais, P., Igel, C., Pascual, A., Guerra-Hernandez, J., Li, S., Mugabowindekwe, M., Saatchi, S., Yue, Y., Chen, Z., and Fensholt, R.: The overlooked contribution of trees outside forests to tree cover and woody biomass across Europe, Sci Adv, 9, https://doi.org/10.1126/sciadv.adh4097, 2023.

Monitoring clearcutting and subsequent rapid recovery in Mediterranean coppice forests with Landsat time series. Zenodo repository. [dataset]. https://doi.org/10.5281/zenodo.3689194, 2020.

Montaghi, A., Corona, P., Dalponte, M., Gianelle, D., Chirici, G., and Olsson, H.: Airborne laser scanning of forest resources: An overview of research in Italy as a commentary case study, International Journal of Applied Earth Observation and Geoinformation, 23, 288–300, https://doi.org/10.1016/j.jag.2012.10.002, 2013.

National Inventory of Forests and forest Carbon pools - INFC: https://www.inventarioforestale.org/en/, last access: 8 June 2025.
<figure><img src="image_33.png" /></figure>

Palahí, M., Valbuena, R., Senf, C., Acil, N., Pugh, T. A. M., Sadler, J., Seidl, R., Potapov, P., Gardiner, B., Hetemäki, L., Chirici, G., Francini, S., Hlásny, T., Lerink, B. J. W., Olsson, H., González Olabarria, J. R., Ascoli, D., Asikainen, A., Bauhus, J., Berndes, G., Donis, J., Fridman, J., Hanewinkel, M., Jactel, H., Lindner, M., Marchetti, M., Marušák, R., Sheil, D., Tomé, M., Trasobares, A., Verkerk, P. J., Korhonen, M., and Nabuurs, G.-J.: Concerns about reported harvests in European forests, Nature, 592, E15–E17, https://doi.org/10.1038/s41586-021-03292-x, 2021.

Pauls, J., Zimmer, M., Kelly, U. M., Schwartz, M., Saatchi, S., Ciais, P., Pokutta, S., Brandt, M., and Gieseke, F.: Estimating Canopy Height at Scale, in: Proceedings of the 41st International Conference on Machine Learning, 39972–39988, 2024.

Pirotti, F.: Assessing a Template Matching Approach for Tree Height and Position Extraction from Lidar-Derived Canopy Height Models of Pinus Pinaster Stands, Forests, 1, 194–208, https://doi.org/10.3390/f1040194, 2010.

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

Rajab Pourrahmati, M., Baghdadi, N., Scolforo, H. F., Alvares, C. A., Stape, J. L., Fayad, I., and le Maire, G.: Integration of very high-resolution stereo satellite images and airborne or satellite Lidar for Eucalyptus canopy height estimation, Science of Remote Sensing, 10, 100170, https://doi.org/10.1016/j.srs.2024.100170, 2024.

Ronneberger, O., Fischer, P., and Brox, T.: U-Net: Convolutional Networks for Biomedical Image Segmentation, 234–241, https://doi.org/10.1007/978-3-319-24574-4_28, 2015.

Roy, D. P., Wulder, M. A., Loveland, T. R., C.E., W., Allen, R. G., Anderson, M. C., Helder, D., Irons, J. R., Johnson, D. M., Kennedy, R., Scambos, T. A., Schaaf, C. B., Schott, J. R., Sheng, Y., Vermote, E. F., Belward, A. S., Bindschadler, R., Cohen, W. B., Gao, F., Hipple, J. D., Hostert, P., Huntington, J., Justice, C. O., Kilic, A., Kovalskyy, V., Lee, Z. P., Lymburner, L., Masek, J. G., McCorkel, J., Shuai, Y., Trezza, R., Vogelmann, J., Wynne, R. H., and Zhu, Z.: Landsat-8: Science and product vision for terrestrial global change research, Remote Sens Environ, 145, 154–172, https://doi.org/10.1016/j.rse.2014.02.001, 2014.

Roy, D. P., Kovalskyy, V., Zhang, H. K., Vermote, E. F., Yan, L., Kumar, S. S., and Egorov, A.: Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity, Remote Sens Environ, 185, 57–70, https://doi.org/10.1016/j.rse.2015.12.024, 2016.

Savage, S., Lawrence, R., Squires, J., Holbrook, J., Olson, L., Braaten, J., and Cohen, W.: Shifts in Forest Structure in Northwest Montana from 1972 to 2015 Using the Landsat Archive from Multispectral Scanner to Operational Land Imager, Forests, 9, 157, https://doi.org/10.3390/f9040157, 2018.

Schleich, A., Durrieu, S., Soma, M., and Vega, C.: Improving GEDI Footprint Geolocation Using a High-Resolution Digital Elevation Model, IEEE J Sel Top Appl Earth Obs Remote Sens, 16, 7718–7732, https://doi.org/10.1109/JSTARS.2023.3298991, 2023.

Schwartz, M., Ciais, P., De Truchis, A., Chave, J., Ottlé, C., Vega, C., Wigneron, J.-P., Nicolas, M., Jouaber, S., Liu, S., Brandt, M., and Fayad, I.: FORMS: Forest Multiple Source height, wood volume, and biomass maps in France at 10 to 30 m resolution based on Sentinel-1, Sentinel-2, and Global Ecosystem Dynamics Investigation (GEDI) data with a deep learning approach, Earth Syst Sci Data, 15, 4927–4945, https://doi.org/10.5194/essd-15-4927-2023, 2023.

Selvi, F., Campetella, G., Canullo, R., Chelli, S., Domina, G., Farris, E., Gasperini, C., Rosati, L., Wellstein, C., and Carrari, E.: The Italian endemic forest plants: an annotated inventory and synthesis of knowledge, Plant Ecol Evol, 156, 29–45, https://doi.org/10.5091/plecevo.95929, 2023.
<figure><img src="image_34.png" /></figure>

Shukla, P. R., Skea, J., Buendia, E. C., Masson-Delmotte, V., Pörtner, H.-O., Roberts, D. C., Zhai, P., Slade, R., Connors, S., Diemen, R. van, Ferrat, M., Haughey, E., Luz, S., Neogi, S., Pathak, M., Petzold, J., Pereira, J. P., Vyas, P., Huntley, E., Kissick, K., Belkacemi, M., and Malley, J.: IPCC, 2019: Summary for Policymakers. In: Climate Change and Land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems, n.d.

Solórzano, J. V., Mas, J. F., Gao, Y., and Gallardo-Cruz, J. A.: Land Use Land Cover Classification with U-Net: Advantages of Combining Sentinel-1 and Sentinel-2 Imagery, Remote Sens (Basel), 13, 3600, https://doi.org/10.3390/rs13183600, 2021.

Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R.: Dropout: A Simple Way to Prevent Neural Networks from Overfitting, Journal of Machine Learning Research, 15, 1929–1958, 2014.

Su, Y., Schwartz, M., Fayad, I., García, M., Zavala, M. A., Tijerín-Triviño, J., Astigarraga, J., Cruz-Alonso, V., Liu, S., Zhang, X., Chen, S., Ritter, F., Besic, N., d'Aspremont, A., and Ciais, P.: Canopy height and biomass distribution across the forests of Iberian Peninsula, Sci Data, 12, 678, https://doi.org/10.1038/s41597-025-05021-9, 2025.

Tadono, T., Ishida, H., Oda, F., Naito, S., Minakawa, K., and Iwamoto, H.: Precise Global DEM Generation by ALOS PRISM, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, II–4, 71–76, https://doi.org/10.5194/isprsannals-II-4-71-2014, 2014.

Tadono, T., Nagai, H., Ishida, H., Oda, F., Naito, S., Minakawa, K., and Iwamoto, H.: GENERATION OF THE 30 M-MESH GLOBAL DIGITAL SURFACE MODEL BY ALOS PRISM, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XLI-B4, 157–162, https://doi.org/10.5194/isprs-archives-XLI-B4-157-2016, 2016.

Takaku, J., Tadono, T., and Tsutsui, K.: Generation of High Resolution Global DSM from ALOS PRISM, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XL–4, 243–248, https://doi.org/10.5194/isprsarchives-XL-4-243-2014, 2014.

Takaku, J., Tadono, T., Tsutsui, K., and Ichikawa, M.: VALIDATION OF “AW3D” GLOBAL DSM GENERATED FROM ALOS PRISM, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, III–4, 25–31, https://doi.org/10.5194/isprs-annals-III-4-25-2016, 2016.

Turubanova, S., Potapov, P., Hansen, M. C., Li, X., Tyukavina, A., Pickens, A. H., Hernandez-Serna, A., Arranz, A. P., Guerra-Hernandez, J., Senf, C., Häme, T., Valbuena, R., Eklundh, L., Brovkina, O., Navrátilová, B., Novotný, J., Harris, N., and Stolle, F.: Tree canopy extent and height change in Europe, 2001–2021, quantified using Landsat data archive, Remote Sens Environ, 298, 113797, https://doi.org/10.1016/j.rse.2023.113797, 2023.

USGS Landsat 5 Level 2, Collection 2, Tier 1: https://developers.google.com/earth-engine/datasets/catalog/LANDSAT LT05_C02_T1_L2#description, last access: 8 June 2025.

USGS Landsat 7 Level 2, Collection 2, Tier 1: https://developers.google.com/earth-engine/datasets/catalog/LANDSAT LE07_C02_T1_L2#description, last access: 8 June 2025.

USGS Landsat 8 Level 2, Collection 2, Tier 1: https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2#description, last access: 8 June 2025.

Vermote, E., Justice, C., Claverie, M., and Franch, B.: Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product, Remote Sens Environ, 185, 46–56, https://doi.org/10.1016/j.rse.2016.04.008, 2016.

Viana-Soto, A. and Senf, C.: The European Forest Disturbance Atlas: a forest disturbance monitoring system using the Landsat archive, https://doi.org/10.5194/essd-2024-361, 21 November 2024.
<figure><img src="image_35.png" /></figure>

Vogeler, J. C., Braaten, J. D., Slesak, R. A., and Falkowski, M. J.: Extracting the full value of the Landsat archive: Inter-sensor harmonization for the mapping of Minnesota forest canopy cover (1973–2015), Remote Sens Environ, 209, 363–374, https://doi.org/10.1016/j.rse.2018.02.046, 2018.

Wulder, M. A., White, J. C., Nelson, R. F., Næsset, E., Ørka, H. O., Coops, N. C., Hilker, T., Bater, C. W., and Gobakken, T.: Lidar sampling for large-area forest characterization: A review, Remote Sens Environ, 121, 196–209, https://doi.org/10.1016/j.rse.2012.02.001, 2012.

Wulder, M. A., White, J. C., Loveland, T. R., Woodcock, C. E., Belward, A. S., Cohen, W. B., Fosnight, E. A., Shaw, J., Masek, J. G., and Roy, D. P.: The global Landsat archive: Status, consolidation, and direction, Remote Sens Environ, 185, 271–283, https://doi.org/10.1016/j.rse.2015.11.032, 2016.

Xu, H., Ren, M., and Yang, L.: Evaluating the consistency of surface brightness, greenness, and wetness observations between Landsat-8 OLI and Landsat-9 OLI2 through underfly images, International Journal of Applied Earth Observation and Geoinformation, 124, 103546, https://doi.org/10.1016/j.jag.2023.103546, 2023a.

Xu, Y., Ding, S., Chen, P., Tang, H., Ren, H., and Huang, H.: Horizontal Geolocation Error Evaluation and Correction on Full-Waveform LiDAR Footprints via Waveform Matching, Remote Sens (Basel), 15, 776, https://doi.org/10.3390/rs15030776, 2023b.

ESA WorldCover 10 m 2021 v200: https://doi.org/10.5281/zenodo.7254221, last access: 11 June 2024.

Zheng, Q., Yang, M., Tian, X., Jiang, N., and Wang, D.: A Full Stage Data Augmentation Method in Deep Convolutional Neural Network for Natural Image Classification, Discrete Dyn Nat Soc, 2020, 1–11, https://doi.org/10.1155/2020/4706576, 2020.

Zhu, X. X., Tuia, D., Mou, L., Xia, G.-S., Zhang, L., Xu, F., and Fraundorfer, F.: Deep Learning in Remote Sensing: A Comprehensive Review and List of Resources, IEEE Geosci Remote Sens Mag, 5, 8–36, https://doi.org/10.1109/MGRS.2017.2762307, 2017.

Zhu, Z. and Woodcock, C. E.: Object-based cloud and cloud shadow detection in Landsat imagery, Remote Sens Environ, 118, 83–94, https://doi.org/10.1016/j.rse.2011.10.028, 2012.

Zhu, Z. and Woodcock, C. E.: Automated cloud, cloud shadow, and snow detection in multitemporal Landsat data: An algorithm designed specifically for monitoring land cover change, Remote Sens Environ, 152, 217–234, https://doi.org/10.1016/j.rse.2014.06.012, 2014.