# Forest attribute maps: a support for small area estimation of forest disturbances

Ankit Sagar, Cédric Vega, Olivier Bouriaud, Christian Piedallu, Jean-Pierre Renaud

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

To cite this version:

Ankit Sagar, Cédric Vega, Olivier Bouriaud, Christian Piedallu, Jean-Pierre Renaud. Forest attribute maps: a support for small area estimation of forest disturbances. Annals of Forest Science, 2025, 82 (1), pp.23. 10.1186/s13595-025-01293-8. hal-05124727

HAL Id: hal-05124727

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

Submitted on 23 Jun 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_2.png" /></figure>

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

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

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

# Forest attribute maps: a support for small area estimation of forest disturbances

Ankit Sagar¹,²,³,⁴* , Cédric Vega² ID, Olivier Bouriaud²,⁵ ID, Christian Piedallu³ and Jean-Pierre Renaud²,⁶ ID

**Abstract**

**Key Message** Forest attribute maps generated using national forest inventory and remote sensing data can be used to quantify forest attributes (i.e., growing stock volume, basal area) affected by disturbance events (like bark-beetle induced damage) through small area estimation techniques. Our results showed that mature stands were more sensitive to bark-beetle infestation than younger ones.

**Context** Recurrent drought episodes in Europe weakened trees and promoted biotic disturbances such as bark-beetle outbreaks, heavily affecting spruce (Picea abies (L.) H. Karst.) and fir stands (Abies alba Mill.). While remote sensing data make it possible to assess and monitor the spatial extent of outbreaks, the assessment of the actual impact on forest attributes is still lacking.

**Aims** We propose a small area estimation approach to assess the amount of forest attributes impacted in outbreak areas during the period 2021–2022 using forest attribute maps.

**Methods** We produced forest attribute maps using National Forest Inventory (NFI) plots, and auxiliary data from Light Detection and Ranging (LiDAR) and Sentinel-2, acquired in 2021. The maps were generated using a general k-nearest neighbour model for all forest species (general model) and a model specific to spruce and fir only (specific model). Small estimation domains, defined as bark-beetle outbreak areas, were obtained from an independent source using time series of Sentinel-2 data (2018–2022). Small area estimations were generated using bootstrapped model-predictions with reliability assessment.

**Results** In total, 1137 ha of forest area was affected by bark-beetle outbreaks during 2021–2022, corresponding to 494,000 m³ (general model) and 505,000 m³ (specific model) of growing stock volume, respectively. The associated errors were 2.2% for the general model and 2.6% for the specific model. In the specific model, the estimated mean growing stock volume in outbreak areas for 2021 and 2022 was 441.5 m³ ha⁻¹, which was 64.5 m³ ha⁻¹ higher than the estimated mean volume of spruce and fir for the whole area (377 m³ ha⁻¹). This suggests that mature stands were more sensitive to infestation.

**Conclusion** Regular updates to forest attribute maps, possibly every third year using photogrammetric methods, could improve monitoring of forest attributes impacted by disturbances such as drought and bark-beetle outbreaks and contribute to forest decision making and planning.

Handling editor: Andreas Bolte.

This article is part of the topical collection on "Impacts of natural disturbances on forest ecosystems under climate change: large-scale analysis and modelling".

*Correspondence:
Ankit Sagar
sagar1@nus.edu.sg

Full list of author information is available at the end of the article

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

BMC

© The Author(s) 2025. **Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
**Keywords** Multisource forest inventories, Model-based inference, Forest attribute maps, Uncertainty analysis, Forest disturbances, LiDAR

# 1 Introduction

With the ongoing climate change, the frequency of forest disturbance events, such as windstorms, fires, drought episodes, and insect outbreaks, is increasing worldwide (Canelles et al. 2021). These events are important drivers of local and global change in forest ecosystems, as they disrupt their composition, structure and functioning, and alter forest resource availability (Seidl et al. 2017). In order to inform decision-making processes in the aftermath or during such events, forest managers and policy makers require precise estimates (e.g. population mean and total) of impacted forest attributes (e.g. GSV—growing stock volume and BA – basal area), preferably in the form of maps that pinpoint the location and quantity of impacted attributes (Bright et al. 2013; Renaud et al. 2017).

While methods are available to identify and map the extent of forest disturbance areas (Abdullah et al. 2019; Lausch et al. 2013; Myroniuk et al. 2020), techniques for estimating the impacted GSV and BA are rarely explored. This can be attributed to the discrete nature of these disturbances, characterized by varying spatial resolution (extent), to which National Forest Inventories (NFIs) are not adapted (Coulston et al. 2020). In such cases, Multi-source National Forest Inventories (MSNFIs) become an attractive option. MSNFIs rely on a statistical combination of NFI data with spatially continuous auxiliary data (e.g. thematic maps, remote sensing data) (Tomppo et al. 2008), and produce landscape-level estimates as well as wall-to-wall maps of forest attributes, irrespective of spatial resolution constraints (Esteban et al. 2019; Saarela et al. 2020; Sagar et al. 2022). The resulting MSNFI forest attribute maps could be used for estimating the amount of impacted forest attributes by intersecting them with disturbance maps. The estimation framework in such small disturbed areas can be further supported by small area estimation (SAE) techniques (Rao and Molina 2015).

SAE techniques have gained popularity in forest inventory applications due to their ability to provide precise estimates of forest attributes in smaller geographic domains (Gaines and Affleck 2021; Stanke et al., 2022). The technique has been mainly used to improve estimation precision at administrative (Breidenbach et al. 2010; Vega et al. 2021) or management unit scales (Hill et al. 2018), but it has been marginally used in cases of forest fires (Gaines and Affleck 2021) or forest harvesting (Coulston et al. 2021). The precision improvement is achieved by "borrowing strength" of

NFI and auxiliary data from other similar areas through statistical models (Dettmann et al. 2022). Because disturbed areas can vary greatly in shape, size, and spatial distribution (Seidl et al. 2017) and there may not be sufficient NFI data to produce reliable estimates for such areas, dedicated estimation frameworks should be used instead.

The choice of the statistical model or estimation framework depends on the characteristics of the disturbance areas and the availability of field data within those areas. In the case of large disturbance patches, such as those resulting from large fires or extensive windstorms, the number of NFI plots available before the disturbance event may be sufficient for both design-based direct estimators (Gaines and Affleck 2021) and design-based model-assisted estimators (Bell et al. 2022; McConville et al. 2020) to produce unbiased estimates (Gregoire 1998). Conversely, disturbances occurring at smaller scales or exhibiting heterogeneous spatial patterns (Bright et al. 2012; Li et al. 2021) do not bear enough NFI plots, making traditional design-based approaches unfeasible. Addressing this limitation would require either additional field sampling (which would be costly as it require extra field work and adjustments to both the sampling and estimation framework), or adoption of alternative estimation methods better suited for small-area inference. Model-based estimators can provide precise estimates for such areas when well-specified models are used, even in the absence of local NFI plots (Frescino et al. 2022). For example, when estimating forest harvesting in the Southeastern United States, Coulston et al. (2021) found that model-based estimators achieved greater precision compared to design-based ones. However, if the model is not well specified (auxiliary data may not fully capture the variability present in NFI plot data) there is a risk of bias (Sagar et al. 2022). Additionally, estimating the uncertainty of model-based approaches remains challenging, especially at the pixel level (Kangas et al. 2023). A possible way to improve model performance and reduce prediction biases involves developing site-specific models (Deo et al. 2017) or constraining model development and application to small scales (Westfall 2015).

Another challenge in estimating attributes impacted by disturbances, relates to the availability of appropriate auxiliary data to minimize both the impact of the disturbance on the data and the time lag between data collection and the occurrence of the disturbance. As
mentioned earlier, the former depends on the nature
of the disturbance and its spatial-temporal properties,
while the latter depends mainly on the availability of
auxiliary data. The 3D auxiliary data from Light Detec-
tion and Ranging (LiDAR) or photogrammetry are the
most efficient remote sensing products for character-
izing forest attributes related to structure, and their
availability at large scales is well established (Guerra-
Hernández et al. 2024; Nilsson et al. 2017). However,
despite the growing popularity of LiDAR surveys, large-
scale renewals remain limited to a handful of countries
and typically occur only at decennial intervals (Kangas
et al. 2018). In countries that conduct regular aerial
surveys, photogrammetric data could be available with
a renewal frequency of a few years (Irulappa-Pillai-
Vijayakumar et al. 2019), allowing for forest attribute
estimation with minimal time lags to the disturbance
event. This aligns with the expectations of forest prac-
titioners for a more frequent renewal of forest attribute
maps (Fassnacht et al. 2024). High-resolution (~ 10–30
m) and high-frequency (~ 1 year) height maps at land-
scape level are now emerging as a result of development
of 3D datasets from multiple spaceborne LiDAR instru-
ments and their integration with optical and/or radar
imagery through deep learning techniques (Lang et al.
2022; Potapov et al. 2021). These products will gradu-
ally become key to support the development of moni-
toring systems of forest structure change at large scale.

With the development of remote-sensing based forest attribute maps, the assessment of prediction uncertainty at pixel level becomes crucial to identify areas with low prediction performance (Kangas et al. 2023; Saarela et al. 2020; Sagar et al. 2021). The calibration data from NFI may not capture the full variability present in the auxiliary data (due to low sampling intensity), causing an increased amount of extrapolation in pixel predictions (Sagar et al. 2021), thereby raising concerns about the reliability of these products. The spatial arrangement of uncertainty can be used to identify forest areas with low reliability or to highlight small area estimates that have higher uncertainty above the mean. This could be used as supplementary diagnostic information to evaluate the reliability of estimates in disturbed areas compared to those for the whole area.

The objective of this research was to develop a model-
based small area estimation (SAE) approach for estimat-
ing the amount of forest attributes impacted by forest
disturbance events, not only in the classical form of an
area-based estimation (i.e. area statistics) but also as
maps. The proposed approach relies on the combina-
tion of MSNFI forest attribute maps with an indepen-
dently generated disturbance map. The approach aimed
to model two major forest attributes, namely GSV and

BA, using auxiliary data from forest maps, 3D LiDAR
point cloud data from airborne laser scanning (ALS), and
metrics from Sentinel-2 remote sensing images. In addi-
tion, the prediction maps of impacted forest attributes
are complemented with reliability maps, following the
method suggested by Sagar et al. (2022).

The specific focus of this study was on the ongoing
bark beetle (Ips typographus L.) outbreak in the Vosges
mountain forest (North East France), which has been
triggered by prolonged drought episodes (Piedallu et al.
2023). Bark beetles are responsible for around 8% of all
tree mortality induced by natural disturbances over
Europe since 1850 (Schelhaas et al. 2003). The beetles
mainly attack spruce and fir species, and rapidly spread
to the neighbouring forest stands, resulting in the prema-
ture death of healthy trees and widespread damage. This
not only reduces the current availability of timber but
also affects its future supply, leading to significant social
and economic impacts (Hlásny et al. 2021). Since the
bark beetle primarily attacks spruce and fir in European
temperate forests (Morris et al. 2017), a secondary objec-
tive was to compare the performance of a general model,
based on all NFI plots, to a specific model, based only
on field plots dominated by the above two species. This
comparison was conducted based on the hypothesis that,
while the general model benefits from a larger dataset it
could also be more prone to prediction errors because
the entire field information is used to predict a subset of
it (Thapa et al. 2015). This raises a critical question: does
the advantage of having more data in a general model
outweigh the potential benefit of reduced noise in a more
tailored, specific model when estimating bark beetle-
impacted attributes? Clarifying this balance is important
for determining which modelling approach provides the
most reliable estimates.

**2 Materials and methods**

2.1 Study site

The study was carried out in the northern sector of
the Vosges mountain forest, located in North-Eastern
France (Grand Est region) (Fig. 1). The area covers 3,612
km² and extends from 6.57° to 7.50° east and 48.01° to
48.89° north. The area is predominantly mountainous,
with elevations ranging from 167 to 1299 m and slopes
sometimes exceeding 30° (Giacona et al. 2018). The area
is characterized by temperate semi-continental climate
(Minářová et al. 2018) with mean annual temperature
between 6° C and 11° C, and mean annual precipitation
between 500 and 2000 mm. About 64% (2,305 km²) of
the study area is covered with forests. Forested areas are
dominated by coniferous stands (51%) mainly of spruce
and fir (i.e., *Picea abies* (L.) H. Karst. and *Abies alba*
Mill.) found at higher altitudes. Pure broadleaved stands
<figure><img src="image_7.png" /><figcaption>**Fig. 1** Location of the study site and spatial distribution of the 430 NFI plots inside the forest types. Subfigures (a and b) are showing patches of bark beetle-impacted forest area</figcaption></figure>

on low altitude account for 22% of the forest area and are dominated by beech and oak (*Fagus sylvatica* L. and *Quercus petraea* (Matt.) Liebl.). The remaining forested area is covered by mixed stands (27%) of coniferous and broadleaved species.

## 2.2 NFI data

The study area includes 430 NFI field plots from the aggregation of five annual samples acquired between 2016 and 2020. This 5-year time interval aligns with the standard practice of the French NFI for generating official statistics (Robert et al. 2009).

Field information was collected within four concentric circular plots of radii 6 m, 9 m, 15 m, and 25 m. Tree measurements were conducted within the first three concentric plots, while stand cover descriptions were assessed over the larger one. Trees were categorized based on their diameter at breast height (DBH at 1.3 m).

Small trees with a DBH in the range [7.5 cm – 22.5 cm [were measured on the 6 m radius plot, medium-sized trees with a DBH in the range [22.5 cm – 47.5 cm [were measured on the 9 m radius plot, and larger trees on the 15 m radius plot. Field measurements included species, circumference, and status (i.e. dead or alive) among others. Tree volume was computed from tree data using species-specific equations (Deleuze et al. 2014). Plot level information was aggregated from tree level data and corresponding inclusion probabilities (Hervé, 2016; Hervé et al. 2014). For this study, two forest attributes were considered: growing stock volume (GSV, m³ ha⁻¹) and basal area (BA, m² ha⁻¹).

## 2.3 Auxiliary data

Along with the NFI data, four freely available auxiliary datasets were used. These data include: (i) a forest map, (ii) ALS data, (iii) Sentinel-2 optical data, and (iv) bark
beetle outbreak maps. These datasets were collected
before and after the bark beetle outbreak. The timeline of
these datasets is presented in Fig. 2.

2.3.1 **Forest map**

The forest map (BD-Forêt® version 2.0) is a freely avail-
able vector map provided by the Institut National de
l'Information Géographique et Forestière (IGN, https://
geoservices.ign.fr/, last consulted August 24, 2024). The
map provides geographical location, quantitative (area),
and qualitative (composition, 32 classes) descriptions of
forest stands of at least 0.5 hectares.

2.3.2 Airborne laser scanning data

The airborne laser scanning (ALS) data was acquired by
IGN in January 2021 using a Leica ALS70 laser scanner.
The scanner operated at an altitude of 1,797 m above
ground level, with a pulse repetition frequency of 309
kHz and a scanning frequency of 35 Hz. The maximum
half scan angle was 25°, and the mean point density at
ground level was 8.3 pts m⁻². The preliminary processing
of the ALS point cloud was done at IGN using Terrasol-
id's Terrascan software (www.terrasolid.com) and deliv-
ered as 1 × 1 km tiles.

Ground classified points were used to generate a 1 m
resolution digital terrain model through inverse distance
weighted interpolation. The elevation of the points above
the ground was converted into height by subtracting the
corresponding values from the digital terrain model. Sub-
sequently, the normalized point cloud data was used to
generate a canopy height model of 1 m pixel resolution
by selecting the maximum point height per pixel. In cases
where pixels did not have values, interpolation was per-
formed by averaging the values of the eight neighbouring
pixels (Vega et al., 2016).

2.3.3 Sentinel-2 data

During the growing seasons of 2020 and 2021, no cloud-
free images were available from the Copernicus Senti-
nel-2 satellite. Consequently, two images acquired on

June 29, 2019, were used. The images were obtained from
the Theia portal (www.theia-land.fr) with the processing
level 2A products. This processing level included ortho-
rectification, atmospheric correction, and cloud detec-
tion with a shadow mask (Hagolle et al. 2015). The data
included spectral bands with spatial resolutions of 10 m
(blue, green, red, and near-infrared), and 20 m (red edge,
three near-infrared and two short-wave infrared bands).

2.3.4 Bark-beetle outbreak maps

Bark-beetle outbreak maps from 2018 to 2022 were
obtained from Dutrieux (2023) and generated using the
ForDead method (Dutrieux et al. 2020) on Coperni-
cus Sentinel-2 images. These maps assess tree health by
detecting dieback in stands dominated by pure fir, pure
spruce, and two mixed conifer stand types, with or with-
out pine trees in the mixture according to BD Forêt®. The
method uses time series of CR<sub>SWIR</sub> (Continuum Removal
SWIR; Fig. 3) vegetation indices to detect vegetation anom-
alies. CRswir relies on a convex hull fit and accounts for
the differences between SWIR1 and SWIR 2 bands with
respect to the NIRa band, providing information about
water stress.

Within the spruce and fir stands mask, bare ground
and clouds were identified to filter out non-vegetation
pixels. Then, a harmonic model was fitted to 2 years of
Sentinel-2 data collected prior to the bark-beetle event,
and anomalies were detected by comparing the actual
vegetation index with the predicted index from the
harmonic model (Dutrieux et al. 2020). These anoma-
lies were accumulated over the year and classified into
four classes. Class 1 and 2 consist of pixels with one or
two successive anomalies detected, respectively. Class
3 includes pixels affected by at least three consecutive
anomalies and was considered as confirmed outbreak
pixels. Lastly, class 4 consists of pixels where anoma-
lies have been detected before being identified as bare
ground, indicating sanitary logging. It should be noted
that anomalies can be reversed (i.e., from outbreak to a
healthy state). For example, a pixel previously labelled as

<figure><img src="image_8.png" /><figcaption>**Fig. 2** A timeline showing the various datasets used in the study. The distinct colours on the timeline correspond to their respective acquisition year or date</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 3** Illustration of the calculation of the CR<sub>SWIR</sub> index on a spectral profile corresponding to a "healthy" stand and an "impacted" stand profile (adapted from Dutrieux et al. 2020)</figcaption></figure>

confirmed outbreak can return to a healthy state if three consecutive non-detections of anomalies were observed. More detailed information about the method may be found in Dutrieux et al. (2020). The source code of the method is available at the following address: https://fordead.gitlab.io/fordead_package/ (last consulted, February 12, 2025). The method's accuracy was evaluated through field verification. A total of 397 detected anomalies over France were visited and less than 20% error was reported for Classes 3 and 4 (Bernard and Doridant 2024).

## 2.4 NFI data filtering

NFI data were collected up to 4 years before the ALS data acquisition (i.e. 2016—2020). To ensure that the models constructed to estimate GSV and BA remained unaffected by disturbances, only NFI data outside of areas affected by bark beetles were used to calibrate the models. The NFI plot positions were intersected with ForDead maps from 2018 to 2020, and plots located within outbreak areas were discarded from the analysis to

maintain consistency between the field data and auxiliary data used for generating attribute maps. For the specific model, NFI plots including fir and spruce, and located inside the four BD-Forêt® classes used by ForDead, were chosen. The plot statistics are presented in Table 1.

## 2.5 Feature-level auxiliary data fusion

Auxiliary data processing was performed at a spatial resolution of 30 m to match the diameter of the NFI plots. A computation grid with a 30 m resolution was set according to the Lambert conic projection using Lambert-93 geodetic system (EPSG: 2154). The origin of the grid was defined as the origin of the coordinate system, and the extent of the grid fixed according to the bounding box of the area of interest. Auxiliary variables/features derived from the auxiliary data were fused within this grid for modelling purposes. The following paragraphs provide a detailed description of auxiliary data processing.

BD Forêt® was used to classify each grid cell as forest or non-forest. Cells whose centre was located within

<table><caption>**Table 1** Number of NFI plots (N), mean and standard deviation (SD) of growing stock volume (GSV, m³ ha⁻¹) and basal area (BA, m² ha⁻¹) for all, training and testing plots in general (all NFI plots) and specific (only spruce and fir plots) model. Training and testing plots are described in Sect. 2.6</caption><thead><tr><th rowspan="2">Forest attribute</th><th rowspan="2">Model</th><th colspan="3">All plots</th><th colspan="3">Training plots</th><th colspan="3">Testing plots</th></tr><tr><th>N</th><th>Mean</th><th>SD</th><th>N</th><th>Mean</th><th>SD</th><th>N</th><th>Mean</th><th>SD</th></tr></thead><tbody><tr><td rowspan="2">GSV</td><td>General</td><td>430</td><td>308.5</td><td>208.3</td><td>344</td><td>317.1</td><td>211.2</td><td>86</td><td>274.1</td><td>193.4</td></tr><tr><td>Specific</td><td>176</td><td>369.2</td><td>215.7</td><td>141</td><td>360.6</td><td>209.0</td><td>35</td><td>403.8</td><td>240.6</td></tr><tr><td rowspan="2">BA</td><td>General</td><td>430</td><td>30.7</td><td>15.2</td><td>344</td><td>31.4</td><td>15.3</td><td>86</td><td>27.9</td><td>14.5</td></tr><tr><td>Specific</td><td>176</td><td>34.9</td><td>14.7</td><td>141</td><td>35.0</td><td>14.3</td><td>35</td><td>34.6</td><td>16.5</td></tr></tbody></table>
the boundaries of BD Forêt® polygons were classified as forest, otherwise non-forest. Non-forest cells and cells partially inside the BD Forêt® were discarded from the computation domain. Furthermore, the forest type (conifer, hardwood, or mixed) of each cell was determined using the attributes of the BD Forêt® polygons, with the cell centroid serving as the reference point (Irulappa-Pilai-Vijayakumar et al. 2019).

Two types of features were extracted from the ALS data: point-based and raster-based. Point-based features were computed for each grid cell using the first return points above 5 m above ground. These features included: minimum, maximum, and mean height, along with the standard deviation, skewness, kurtosis, and entropy of the height distribution. Additionally, the percentage of returns above the mean height, the percentile of the height distribution, and the cumulative percentage of return were also computed (Roussel et al. 2020).

On the other hand, given that the study area represents a mountainous region, a set of raster-based topographic features was derived using the digital terrain model. These features included the mean and standard deviation of altitude and slope of the terrain. Another set of raster-based features was computed using the canopy height model, which included: Rumple index, gap area, gap ratio, total inner volume, total outer volume, canopy inner volume, canopy outer volume, gap inner volume, and gap outer volume (Véga et al. 2016). Computational details are provided in Table 2.

All spectral bands with a spatial resolution above 10 m were resampled to a resolution of 10 m. Along with 10 spectral bands, 17 vegetation indices (Table 3) were

<table><caption>**Table 2** Canopy height model (CHM) metrics. *P* represents a pixel height of the CHM, and *s* its spatial resolution. *pMax* represents the maximum height of the portion of the CHM considered, and *i* and *j* represents line and column indices of *p*. *CHM<sub>g</sub>* and *CHM<sub>pa</sub>* represented respectively the area and the projected area of the CHM (Véga et al. 2016)</caption><thead><tr><th>Metric name</th><th>Formulae</th></tr></thead><tbody><tr><td>Total inner volume</td><td>TIV = &Sigma;<sub>p&gt;0</sub> p &times; s<sup>2</sup></td></tr><tr><td>Total outer volume</td><td>TOV = &Sigma;<sub>p&gt;0</sub> (p<sub>Max</sub> &times; s<sup>2</sup>) &minus; TIV</td></tr><tr><td>Canopy inner volume</td><td>CIV = &Sigma;<sub>p&gt;Th</sub> p &times; s<sup>2</sup></td></tr><tr><td>Canopy outer volume</td><td>COV = &Sigma;<sub>p&gt;Th</sub> (p<sub>Max</sub> &times; s<sup>2</sup>) &minus; ICV</td></tr><tr><td>Gap Area</td><td>GA = &Sigma;<sub>p&lt;Th</sub> s<sup>2</sup></td></tr><tr><td>Gap inner volume</td><td>GIV = &Sigma;<sub>p&lt;Th</sub> p &times; s<sup>2</sup></td></tr><tr><td>Gap outer volume</td><td>GOV = (p<sub>Max</sub> &times; GA) &minus; IGV</td></tr><tr><td>Rumple area</td><td>RA = CHM<sub>g</sub> / CHM<sub>pa</sub><br>CHM<sub>pa</sub> = &Sigma;<sub>p&gt;Th</sub> s<sup>2</sup><br>CHM<sub>g</sub> = &Sigma;<sub>p&gt;Th</sub> || <sup>(pMax+pMin)/2</sup>&frasl;<sub>2</sub> + <sup>(pMax+pMin)/2</sup>&frasl;<sub>2</sub> || <sup>(pMax+pMin)/2</sup>&frasl;<sub>2</sub> || <sup>s<sup>2</sup></sup>&frasl;<sub>s</sub></td></tr></tbody></table>

<table><caption>**Table 3** Vegetation indices calculated from Sentinel-2 data (Sonobe et al. 2018; www.custom-scripts.sentinel-hub.com). "b" refers to band number</caption><thead><tr><th>Vegetation indices</th><th>Formulae</th></tr></thead><tbody><tr><td>Atmospherically resistant vegetation index</td><td>b8-b4-1*(b4-b2)<br>b8-b4-1*(b4-b2)</td></tr><tr><td>Advanced vegetation index</td><td>√b8 * (1 - b4) * (b8 - b4)</td></tr><tr><td>Canopy chlorophyll content index</td><td>(b8-b5)/(b8+b5)<br>(b8-b4)/(b8+b4)</td></tr><tr><td>Chlorophyll index green</td><td>b8/b3 - 1</td></tr><tr><td>Disease-water stress index</td><td>b3/b4</td></tr><tr><td>Enhanced vegetation index</td><td>2.5*(b8-b4)<br>(b8+6.0*b4-7.5*b2)+1.0</td></tr><tr><td>Green normalized difference vegetation index</td><td>b8-b3<br>b8+b3</td></tr><tr><td>Modified chlorophyll absorption reflectance index</td><td>(b5 - b4) - 0.2 * (b5 - b3) * <sup>b5</sup>&frasl;<sub>b4</sub></td></tr><tr><td>Modified soil adjusted vegetation index</td><td>2.0*b8+1.0 - √<sup>(2.0*b8+1.0)<sup>2</sup> - 8.0*(b8-b4)</sup>&frasl;<sub>2</sub></td></tr><tr><td>Moisture stress index</td><td>b11<br>b8</td></tr><tr><td>Normalized difference moisture index</td><td>b8-b11<br>b8-b11</td></tr><tr><td>Normalized difference red-edge index</td><td>b8-b5<br>b8+b5</td></tr><tr><td>Normalized difference vegetation index</td><td>b8-b4<br>b8-b4</td></tr><tr><td>Normalized difference water index</td><td>b3-b8<br>b3+b8</td></tr><tr><td>Pigment specific simple ratio</td><td>b7/b4</td></tr><tr><td>Simple ratio</td><td>b8/b4</td></tr><tr><td>Visible atmospherically resistant index green</td><td>b3-b4<br>b3+b4-b2</td></tr></tbody></table>

computed as auxiliary variables (Sonobe et al. 2018). The resampled images and vegetation indices were then mosaicked and aggregated to a 30 m resolution.

For the bark-beetle analysis, only infestations that occurred after the period following ALS acquisition were considered for estimating the bark-beetle impacted forest attributes (i.e., 2021 and 2022). This procedure was used to mitigate potential bias in predictions due to alteration of auxiliary variables by bark-beetle activity. The estimation domain was defined as the area formed by ForDead's cells belonging to classes 3 and 4, which correspond to confirmed outbreak detection. A computation grid cell was classified as disturbed if it contained at least one outbreak pixel. Furthermore, for each cell, the fraction of the surface area impacted was computed using an area-weighted approach (Eq. 1) (Hill et al. 2021):

$$prop_{outbreak} = \frac{area_{outbreak}}{area_{grid}} \qquad (1)$$

where, *prop*<sub>outbreak</sub> is the proportion of outbreak area in the grid cell, *area*<sub>outbreak</sub> is the total outbreak area in one grid cell and *area*<sub>grid</sub> is the area of one grid cell.
The same set of features was extracted for the cells centred on each NFI plot for the model development process.

## 2.6 Model selection

The workflow of model selection and assessment is presented in Fig. 4.

The relationship between NFI forests attributes of interest and auxiliary data was investigated using a non-parametric k-nearest neighbour (kNN) model. Such models are widely used in MSNFI because they do not rely on assumptions regarding the distribution of the variables, can predict multiple forest attributes at once, and preserve the covariance structure of the variables (Kim and Tomppo 2006). Also kNN model relies on a limited set of parameters (e.g., number of nearest neighbours (k), distance metric, weighting function, and set of auxiliary variables) making them attractive for modelling purposes (McRoberts 2012).

Following earlier studies, the models were trained using k= 5 (Chirici et al. 2016; Irulappa-Pillai-Vijayakumar et al. 2019; Sagar et al. 2022), and two distance metrics were considered: Most Similar Neighbours (msn; Eq. 2) and Euclidean (euc; Eq. 3) (Crookston and Finley 2008), along with an inverse distance weighting function (Eq. 4) (Altman 1992):

$$d_{ij}(\text{msn}) = \sqrt{(x_i - x_j)^2 \Gamma \Lambda^2 \Gamma^T (x_i - x_j)^T} \quad (2)$$

$$d_{ij}(\text{euc}) = \sqrt{(x_i - x_j)^T (x_i - x_j)} \quad (3)$$

$$w_{ij} = \frac{1}{1 + d_{ij}} \quad (4)$$

where, $d_{ij}(\text{msn})$ and $d_{ij}(\text{euc})$ are most similar neighbour and Euclidean distances, $x_i$ and $x_j$ are the vector of target and reference auxiliary variables, $\Lambda$ is a matrix of canonical correlation, $\Gamma$ is a matrix of canonical vectors of auxiliary data, $T$ is the transpose of the matrix, and $w_{ij}$ is the weighted distance.

Auxiliary variables were selected using an iterative feature selection algorithm applied in both addition and deletion modes (Crookston and Finley 2008). For each distance metric and variable selection approach, the best models made of 1 up to 7 variables were selected. Such a limit was set because kNN models are sensitive to the dimensionality of auxiliary variables which may lead to model overfitting problems (Breidenbach et al. 2012). Prior to variable selection, the set of auxiliary variables was filtered to avoid multicollinearity. The filtering process relied on the correlation matrix between the response variables and the auxiliary variables. The auxiliary data were ranked in descending order with respect to their correlation with the response variables. A vector of candidate auxiliary variables was initiated by selecting the auxiliary variable having the highest correlation with the response variables. The remaining auxiliary variables were iteratively tested against the vector in descending order and included in the vector if their correlation with the auxiliary variables contributing to it was lower than a threshold fixed here to 0.7 (Sagar et al. 2022).

Among the 28 models generated (2 distance metrics × 2 selection approaches × 7 sets of auxiliary data) for both the general and specific models, the best one was selected and evaluated using a train-test approach. The NFI data were split into train and test sets using a geographical partitioning method (Meyer and Pebesma 2021). In short, according to this method, the bounding box of the coordinates of the NFI plots was first generated, and one of its corners was

<figure><img src="image_10.png" /><figcaption>**Fig. 4** Workflow of the plot-level model selection. Field plot data were intersected with auxiliary data, and plots falling inside the outbreak areas were discarded. The resulting plots were split into train and test sets, and the train set was used to select the best model using a leave-one-out cross-validation approach by minimizing an error criterion. The performance of the selected model was evaluated using the independent test set</figcaption></figure>
randomly selected. Then, 20% of the closest plots from the selected corner were kept as an independent testing set and the remaining 80% of the plots as a training set. The same corner was used for generating the training-testing sets for both the general and specific models. The training set was used for model development and selection, and the testing set was used to validate the model against independent external data (Karasiak et al. 2022).

Model selection and assessment were performed on the training set using a leave one out cross-validation (LOOCV) method. Performance metrics included the calibration and the validation relative root mean square error (RMSE%, Eq. 5), the relative mean error (ME%, Eq. 6), the relative total error (TE%, Eq. 7), and a measure of overfitting (OFI, Eq. 8, Valbuena et al. 2017):

$$ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}; \quad \text{RMSE}\% = \frac{\text{RMSE}}{\bar{y}} \times 100 \quad (5) $$

$$ \text{ME} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i); \quad \text{ME}\% = \frac{\text{ME}}{\bar{y}} \times 100 \quad (6) $$

$$ \text{TE} = \sqrt{\text{RMSE}^2 + \text{ME}^2}; \quad \text{TE}\% = \frac{\text{TE}}{\bar{y}} \times 100 \quad (7) $$

$$ \text{OFI} = \frac{1}{n} \sum_{i=1}^{n} (y_i^{cv} - \hat{y}_i) / \frac{1}{n} \sum_{i=1}^{n} (y_i^{train} - \hat{y}_i) \quad (8) $$

where $y_i$ is the observed value, $\hat{y}_i$ is predicted value, $\bar{y}$ is the mean of observed values, n is the number of NFI plots, OFI is the overfitting index, $y_i^{cv}$ are the cross-validated predictions and $y_i^{train}$ are the training predictions.

The best general and specific models were selected in a two-stage process. First, the models with an OFI above 1.1 (Valbuena et al. 2017) were considered as overfitted and were filtered out. Then, the models with TE% within 1% of the minimum TE% were selected, of which the more parsimonious model (Pestov 2013), (i.e. using the least number of auxiliary variables), was considered as the best model.

Finally, the resulting models were calibrated again using the entire training set, and their performance was assessed on the test set, serving as an external assessment.

## 2.7 Model propagation

Forest attribute and reliability maps were generated for the 2 selected models (general and specific) using bootstrap (Fig. 5). For each bootstrap sample ($n_{boot} = 100$), a subset of NFI plots was randomly selected with replacement from the entire NFI data. The NFI sample was then used to train the two kNN models and to predict the grid cells. The final prediction for each forest attribute over

the entire study site was computed by taking the overall mean of bootstrapped predictions and produce predictive raster maps at 30 m spatial resolution.

Reliability maps were generated (Sagar et al. 2022) using classes of precision and accuracy. Precision classes (low, intermediate, high) were derived from the kNN model's calibration-validation RMSE% and coefficient of variation (CV) of the bootstrapped prediction at pixel level. For each pixel, the precision was considered high when CV was below calibration RMSE%, intermediate when CV was in between calibration and validation RMSE%, and low otherwise. Accuracy was defined through a measure of extrapolation (Renaud et al. 2022). It consisted first in generating the convex hull of the auxiliary data used in the models at plot level, and then, to test if the auxiliary data describing each cell of the computation grid were included or not in the convex hull. Cells outside the convex hull were considered extrapolated and classified as having low accuracy, as their validity domain was unknown. On the opposite, cells inside the convex hull were considered as having a high accuracy. Finally, the reliability classes were computed by intersecting precision and accuracy, and assigning the quality class to each cell.

## 2.8 Statistical inference

Population parameters (mean, variance and relative error) were derived from the kNN predictions using a bootstrapped model-based estimator proposed by McRoberts et al. (2022) (Fig. 5):

$$ \hat{\mu}_{\text{boot}} = \frac{1}{n_{\text{boot}}} \sum_{b=1}^{n_{\text{boot}}} \hat{\mu}_{\text{boot}}^b \quad (9) $$

$$ \hat{V}_{\text{boot}}(\hat{\mu}) = \frac{1}{n_{\text{boot}} - 1} \sum_{b=1}^{n_{\text{boot}}} \left( \hat{\mu}_{\text{boot}}^{b} - \hat{\mu}_{\text{boot}} \right)^2 \quad (10) $$

$$ \text{RE}(\hat{\mu}_{\text{boot}}) = \frac{\sqrt{\hat{V}_{\text{boot}}(\hat{\mu})}}{\hat{\mu}_{\text{boot}}} \times 100 \quad (11) $$

where $n_{boot}$ is the number of bootstrap samples, $\hat{\mu}_{boot}^b$ is the mean of the boot b, $\hat{\mu}_{boot}$ and $\hat{V}_{boot}(\hat{\mu})$ are the mean and variance of the mean of the bootstrap samples, respectively. RE($\hat{\mu}_{boot}$) is the relative error of the bootstrap samples, which is the error associated with the estimated mean by expressing the standard error as a percentage of the mean estimate.

The estimation of the GSV and BA impacted by the bark-beetle attacks was achieved by sub-setting the grid cells using the outbreak map. Predicted GSV and BA were converted from per hectare values to per cell
<figure><img src="image_11.png" /><figcaption>**Fig. 5** Workflow of the grid-level model propagation and estimation framework. Forest attribute and reliability maps were generated using bootstrapped predictions and a convex hull approach to identify pixels outside the model calibration domain. The resulting maps were combined with outbreak maps to estimate the amount of forest attributes impacted</figcaption></figure>

values and weighted using the fraction of affected grid cells. From these weighted values, estimates of total, variance of total and relative errors were estimated as follows:

$$ \hat{\bar{t}}_{\text{boot}} = \frac{1}{n_{\text{boot}}} \sum_{b=1}^{n_{\text{boot}}} \sum_{i=1}^{N} \hat{y}_{i} \quad (12) $$

$$ \hat{V}_{\text{boot}}(\hat{t}_{\text{boot}}) = \frac{1}{n_{\text{boot}} - 1} \sum_{b=1}^{n_{\text{boot}}} (\hat{t}_{\text{boot}}^{b} - \hat{t}_{\text{boot}})^2 \quad (13) $$

$$ \text{RE}(\hat{\bar{t}}_{\text{boot}}) = \frac{\sqrt{\hat{V}_{\text{boot}}(\hat{t}_{\text{boot}})}}{\hat{\bar{t}}_{\text{boot}}} \times 100 \quad (14) $$

where, $\hat{y}_i$ is the predicted value, N is the number of outbreak samples, $n_{\text{boot}}$ is the number of bootstrap samples, $\hat{t}_{\text{boot}}$ and $\hat{V}_{\text{boot}}(\hat{t}_{\text{boot}})$ are the total and variance of the total of the bootstrap samples, respectively. RE($\hat{\bar{t}}_{\text{boot}}$) is the relative error of the bootstrap samples, which is the error associated with the estimated total by expressing the standard error as a percentage of the total estimate.

## 3 Results

Since the beginning of the bark-beetle crises in North-Eastern France in 2018, the outbreak covered a cumulative area of 2837 ha (ForDead classes 3 and 4). This represents 3.0% of the area covered by spruce and fir dominated forests within the study site, according to the BD-Forêt map (Table 4). For the period 2021–2022, the total disturbed area reached 1137 ha, corresponding to 1.2% of the surveyed coniferous areas.

<table><caption>**Table 4** Forest surfaces detected by the ForDead method in classes 3, 4 (confirmed outbreak), and the total of the two classes. Values are presented in hectares. Numbers in parenthesis correspond to percent to the total area covered by spruce and fir dominated forests within the study site</caption><thead><tr><th rowspan="2">Period</th><th colspan="3">ForDead classes (ha)</th></tr><tr><th>3</th><th>4</th><th>Total</th></tr></thead><tbody><tr><td>2018–2022</td><td>925 (1.0)</td><td>1912 (2.0)</td><td>2837 (3.0)</td></tr><tr><td>2021</td><td>52 (0.1)</td><td>224 (0.2)</td><td>277 (0.3)</td></tr><tr><td>2022</td><td>642 (0.7)</td><td>219 (0.2)</td><td>861 (0.9)</td></tr><tr><td>2021–2022</td><td>694 (0.7)</td><td>443 (0.5)</td><td>1137 (1.2)</td></tr></tbody></table>
The two selected kNN models (general and specific) for estimating the amount of forest attributes impacted were built with a Euclidean distance and variables selected using a forward (addition) variable selection approach. The general model of forest attributes included 4 auxiliary variables, of which 3 were derived from the ALS data (i.e. total inner canopy volume, gap area, 2nd quantile cumulative percentage of returns) and one from Sentinel-2 data (i.e. mean of near infrared values). The specific model included only two auxiliary variables extracted from the ALS data: the total inner canopy volume and the gap area. These two variables were also included in the general model. Both the NFI and auxiliary data with selected auxiliary variables are publicly available online (Vega et al. 2025; https://doi.org/10.5281/zenodo.15276127.).

For GSV, the general model cross-validated RMSE% ranged from 28.4% (calibration) to 41.7% (test) and from 25.1% (calibration) to 33.1% (test) for the specific model (Table 5). The mean error (bias) for both models ranged from -1.1% to 0.4% in calibration and validation, but exhibited larger values with the independent test set, reaching -5.9% for the general model. For BA, the general model cross-validated RMSE% ranged from 26.8% (calibration) to 42% (test) and from 22.9% (calibration) to 30.8% (test) for the specific model (Table 5). The mean error for both models was below 0.8% in calibration and validation but exhibited greater and negative values with the independent test set, reaching -1.8% for the general model.

General model residuals computed for the independent test set (Fig. 6) did not show significant departure from normality (p-value = 0.14) but exhibited a slight underestimation of GSV above 600 m³ ha⁻¹ while specific model residuals showed a random distribution.

Table 6 presents estimate of forest attributes for the whole study area, and the subset of spruce and fir forests. The overall mean GSV was estimated at 312 m³ ha⁻¹ with an error of 1.78%. The mean BA was estimated at 30.8 m² ha⁻¹ with a relative error of 1.66%.

For the sub-population of spruce and fir forests, the general model estimated a GSV of 366 m³ ha⁻¹ and a BA of 34.2 m² ha⁻¹, with relative errors of 1.92% and 1.77%, respectively. With the specific model, GSV was estimated at 377 m³ ha⁻¹ with a relative error of 2.36%, and the mean BA was 34.4 m² ha⁻¹ with an error of 2.09%.

The estimated forest attributes impacted by bark-beetles with the two models is presented in Table 7, and an illustration of the underlying prediction maps is presented in Fig. 7. In 2021, the GSV of outbreak areas (classes 3–4), covering 277 ha, is estimated at 117,646 m³ and 120,889 m³ with the general and specific models, respectively. These estimates are supported by respectively 75.9% (GSV, general model) and 79.6% (GSV, specific model) of reliable predictions. The corresponding BA totals were respectively 10,187 m² and 10,038 m². In 2022, the 861 ha of outbreaks account for 376,543 m³ and 384,487 m³ of GSV for the general and specific models, respectively. The percentages of reliable pixel were of the same magnitude, ranging from 78.2% (GSV, general model) to 82.6% (GSV, specific model). The relative error between the general and specific models was also of the same magnitude, with a maximum error obtained of 2.6% (2022—GSV, specific model). The reliability of bark-beetle affected pixels was on average 6.5% lower than the one achieved in areas of healthy forest.

Overall, higher frequencies of pixels affected by bark-beetle outbreaks were found for GSV classes below 450 m³ ha⁻¹ (Fig. 8a), although relative values indicated a slight trend towards an increasing proportion of outbreaks occurring in areas with greater GSV (Fig. 8b). Notably in the specific model, the mean of estimated mean growing stock volume in outbreak areas for 2021 and 2022 was 441.5 m³ ha⁻¹ (Table 7), which is 64.5 m³ ha⁻¹ higher than the estimated mean volume of spruce and fir for the whole area (Table 6: 377 m³ ha⁻¹). This suggests that mature stands were more susceptible to bark beetle outbreaks than younger ones in the study area.

<table><caption>**Table 5** Model prediction mean, training (calibration-validation) and testing relative errors of growing stock volume (GSV, m³ ha⁻¹) and basal area (BA, m² ha⁻¹) for the general and specific models</caption><thead><tr><th rowspan="2">Forest attribute</th><th rowspan="2">Model</th><th colspan="3">Train—calibration</th><th colspan="3">Train—validation</th><th colspan="3">Test—external validation</th></tr><tr><th>Prediction Mean</th><th>RMSE%</th><th>ME%</th><th>Prediction Mean</th><th>RMSE%</th><th>ME%</th><th>Prediction Mean</th><th>RMSE%</th><th>ME%</th></tr></thead><tbody><tr><td rowspan="2">GSV</td><td>General</td><td>318</td><td>28.4</td><td>0.4</td><td>317</td><td>37.7</td><td>0.1</td><td>258</td><td>41.7</td><td>-5.9</td></tr><tr><td>Specific</td><td>357</td><td>25.1</td><td>-1.1</td><td>356</td><td>32.8</td><td>-1.1</td><td>411</td><td>33.1</td><td>1.7</td></tr><tr><td rowspan="2">BA</td><td>General</td><td>31.6</td><td>26.8</td><td>0.7</td><td>31.7</td><td>36.0</td><td>0.8</td><td>27.5</td><td>42.0</td><td>-1.8</td></tr><tr><td>Specific</td><td>35.0</td><td>22.9</td><td>0.3</td><td>35.1</td><td>30.1</td><td>0.4</td><td>36.7</td><td>30.8</td><td>1.5</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Fig. 6** Observation and residuals versus prediction of growing stock volume (GSV) over the test set for both general and specific models. A similar pattern was observed for basal area (BA)</figcaption></figure>

<table><caption>**Table 6** Estimation means and relative errors for growing stock volume (GSV, m³ ha⁻¹) and basal area (BA, m² ha⁻¹) for the entire study area, obtained using both the general and specific models. *To test model transferability, the general model was applied to all NFI plots and to conifer plots only</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Forest type</th><th colspan="2">GSV</th><th colspan="2">BA</th></tr><tr><th>μ̂<sub>boot</sub></th><th>RE(μ̂<sub>boot</sub>) (%)</th><th>μ̂<sub>boot</sub></th><th>RE(μ̂<sub>boot</sub>) (%)</th></tr></thead><tbody><tr><td rowspan="2">General</td><td>*All</td><td>312</td><td>1.78</td><td>30.8</td><td>1.66</td></tr><tr><td>*Conifer</td><td>366</td><td>1.92</td><td>34.2</td><td>1.77</td></tr><tr><td>Specific</td><td>Conifer</td><td>377</td><td>2.36</td><td>34.4</td><td>2.09</td></tr></tbody></table>

## 4 Discussion

Our results showed that forest attribute maps originally generated for a general purpose of small area estimation could be combined with polygons of bark-beetle

outbreaks to estimate the amount of impacted forest attributes and to produce maps of the impact intensity. In model-based approaches, model specifications are key to limit biases, and using a site-specific model could be considered as a way to reduce local bias (Westfall 2015). Furthermore, at the pixel level, estimating the uncertainty associated with model predictions remains challenging. Various approaches have been proposed to address this issue in the field of statistics (Kangas et al. 2023) and deep learning (Baumann et al. 2023), or to assess the validity domain of predictions (Sagar et al. 2022), thus contributing to evaluate prediction bias, and avoiding to draw inappropriate conclusions from the estimates.

With a mean GSV of 366 m³ ha⁻¹, the general model estimated a lower amount of overall forest attributes than the specific model (377 m³ ha⁻¹) and exhibited slightly lower relative errors (Table 6). A similar trend was observed in the estimation total, mean, and relative
<table><caption>**Table 7** Estimation totals, means and percentage of high reliable pixels of growing stock volume (GSV) and basal area (BA) in confirmed outbreak areas (classes 3–4). The numbers in parentheses represent the relative error of the total and the percentage of highly reliable undisturbed pixels, respectively</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Year</th><th colspan="3">GSV</th><th colspan="3">BA</th></tr><tr><th>Total (m³)</th><th>Mean (m³ ha⁻¹)</th><th>High reliability (%)</th><th>Total (m²)</th><th>Mean (m²ha⁻¹)</th><th>High reliability (%)</th></tr></thead><tbody><tr><td rowspan="2">General</td><td>2021</td><td>117,646 (± 2.2%)</td><td>425</td><td>75.9 (82.8)</td><td>10,187 (± 2.1%)</td><td>36.8</td><td>77.5 (85.5)</td></tr><tr><td>2022</td><td>376,543 (± 2.2%)</td><td>437</td><td>78.2 (82.9)</td><td>32,245 (± 2.2%)</td><td>37.5</td><td>79.5 (85.7)</td></tr><tr><td rowspan="2">Specific</td><td>2021</td><td>120,889 (± 2.5%)</td><td>436</td><td>79.6 (86.3)</td><td>10,038 (± 2.3%)</td><td>36.2</td><td>78.3 (86.7)</td></tr><tr><td>2022</td><td>384,487 (± 2.6%)</td><td>447</td><td>82.6 (86.3)</td><td>31,564 (± 2.3%)</td><td>36.7</td><td>81.3 (86.8)</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Fig. 7** Map of growing stock volume (GSV) classes over the area of interest (a), with a focus on a subset of outbreak area (d). Corresponding GSV and reliability maps are shown for the general (b, c) and specific (e, f) models</figcaption></figure>

error of total for GSV in dieback areas during 2021 and 2022 (Table 7). However, general model showed a lower amount of highly reliable pixels compared to the specific model, though the difference did not exceed 4.4%

(Table 7). Several factors may explain these results. On one side, the lower amount of field plots in the specific model might generate a higher variability in bootstrap predictions than in the general model. Also, by including
<figure><img src="image_14.png" /><figcaption>**Fig. 8** Frequencies (a) and percentages (b) of pixels affected by bark-beetle outbreaks ventilated by growing stock volume classes of 50 m³ intervals. A similar pattern was observed for basal area</figcaption></figure>

the variability of all forest plots, the general model may benefit from "borrowing strength" from mixed forest stands, which can be similar to conifer stands in terms of composition and structure (Ghosh and Rao 1994). On the other side, the larger variability in field data increases the amount of extrapolation in the general model (11%) compared to the specific model (5%), thus decreasing the overall reliability at pixel level (Sagar et al. 2022). Notably, the differences in the proportions of reliable pixels (i.e., those associated within confidence of model predictions) were consistent between models when comparing pixels for disturbed and undisturbed areas. Overall, our results suggest that a general resource model can be efficiently used to derive small area estimates for specific parts of the domain. It is thus possible to compute predictions in a more effective way if attribute maps are readily available and produced regularly, to support information needs for managing sanitary crises.

While LiDAR is commonly regarded as the most effective source of 3D data for modelling forest attributes, there are still limitations in its capacity to produce timely renewed information. Despite technological advancements allowing broader coverage at lower costs, LiDAR data remains expensive for a frequent renewal (Barrett et al. 2016). However, when a LiDAR digital terrain model is available, photogrammetric measurements can be used to renew the forest information layer (Price et al. 2020), leverage time-series estimates to assess flux variables (Irulappa-Pillai-Vijayakumar et al. 2019), or quantify

forest changes (Renaud et al. 2017). In many countries, aerial photographs with photogrammetric capabilities (i.e., sufficient overlap between images along and across tracks) are routinely acquired for surveying purposes. For instance, in France, such coverage is renewed every 3 years at the departmental level (administrative subdivisions). This acquisition frequency could be adequate to estimate the amount of forest attributes in disturbed areas, with a maximum time difference of two years. Estimations may pose challenges for small areas distributed across different departments that may have been surveyed in different years. For larger scale disturbances, design-based direct estimators may be sufficient, as the number of NFI points could also be sufficient. However, if estimation precision is insufficient, or specific estimates are needed, data from spaceborne sensors such as GEDI or ICESat-2, in combination with optical and/or radar data like Sentinel, SPOT or Landsat data, could support estimation (Duncanson et al. 2022; Lang et al. 2022; Potapov et al. 2021; Schwartz et al. 2023), and enhance the capabilities of countries to monitor forest attributes at national level on a yearly basis. Due to the lower resolution (~ 10–30 m) and precision of the later sources of information compared to LiDAR, it is hypothesized that the resulting attribute maps are more appropriate for design-based model-assisted SAE. For smaller domain estimation under model-based inference, their capacity to provide unbiased estimates of the forest attributes still have to be investigated.
The outbreak maps used as auxiliary data in our study were considered errorless while an independent field assessment indicated a detection error rate of about 20% (Bernard and Doridant 2024). Taking this error into account was out of the scope of the study. Despite those possible errors, the detected areas affected by severe outbreaks corresponded to 1.2% of the total surface (Table 4), which is in agreement with the 1.4% reported by Piedallu et al., (2023) for the same study site using another approach. Considering the proportion of basal area of fir and spruce, and under the hypothesis that both species are affected in the same proportions, this represents 2.5% and 4.2% of the fir and spruce areas, respectively. The result for spruce is in accordance with the one reported by Gilles et al. (2023) who estimated that 5.6% of the spruce stands were lost during the 2017–2024 period in the Vosges mountains. Overall, our results suggested that the mature stands are more impacted than the younger ones (Fig. 8, right). This can be explained by the fact that the development of bark beetle populations is associated with the presence of dead trees which are mainly found in the mature stands (Bennett et al. 2015). The lower reliability obtained in dieback areas (~ 6.5%) compared to healthy areas may also be attributed to this phenomenon. Mature stands are more likely to be in the periphery of the distribution of auxiliary data and thus to be in extrapolation. Another factor contributing to the lower reliability could be early"green"infestation (Hollaus and Vreugdenhil 2019; Huo et al. 2021), which was not detected in optical data but impacted LiDAR data.

Considering the large variability of abiotic and biotic factors influencing forest health, and the equally large variability of forest responses to these stresses, monitoring forest health is challenging. The development of remote sensing products and modelling techniques is central toward the development of such monitoring systems (Massey et al. 2023; Seidl et al. 2017), and should, in the end, allow for the detection of a range of disturbances, the identification of their nature, and provide insights into their dynamics with the appropriate spatial and temporal resolution (Francini et al. 2020). In the context of the bark-beetle crisis, differentiating between clear-cut areas and vegetation anomalies associated with sanitary logging is expected. The accelerated development of time series data from remote sensing, both optical and LiDAR, offers excellent data sources for training dedicated deep learning semantic change detection approaches (Zheng et al. 2022). As indicated above, a promising approach involves leveraging high-temporal-frequency height models trained with 3D data from airborne or spaceborne LiDAR sensors (Ge et al., 2022; Boucher et al. 2020) to support simultaneous surface detection, classification, and estimation of corresponding

attributes. Incorporating such high-resolution, high-temporal-frequency height data into the proposed method could enhance both the renewal frequency and spatial extent of monitored areas, thereby regularly informing policies and supporting decision-making.

## 5 Conclusion

Forest attribute maps generated for assessing the overall forest attribute over an area can be combined with independent disturbance maps to generate SAE of attributes impacted by bark beetles. Pixel predictions in disturbed areas were consistent with those in undisturbed areas, with a loss of reliability limited to 6.5%. This loss is primarily attributed to extrapolation issues, as the mature stands affected by bark beetles are located at the edges of the auxiliary data space (i.e. at the borders of the model calibration domain) that are less common or poorly represented in the training data, making it harder for the model to predict them accurately.

With the development of remote sensing databases enabling to both predict forest attributes and detect growing conditions at both high spatial resolution and high temporal frequency, the estimation of forest attributes over large areas and its monitoring at local scale is gaining in maturity and will contribute to forest decision making and planning.

### Acknowledgements

Not applicable.

### Code availability

Not applicable.

### Authors' contributions

Conceptualization: Ankit Sagar, Cédric Vega, Jean-Pierre Renaud; Methodology: Ankit Sagar, Cédric Vega, Jean-Pierre Renaud; Formal analysis and investigation: Ankit Sagar, Cédric Vega, Jean-Pierre Renaud; Writing - original draft preparation: Ankit Sagar, Cédric Vega; Writing - review and editing: Ankit Sagar, Cédric Vega, Jean-Pierre Renaud, Olivier Bouriaud, Christian Piedallu; Funding acquisition: Cédric Vega, Olivier Bouriaud; Resources: Cédric Vega, Jean-Pierre Renaud, Christian Piedallu; Supervision: Cédric Vega, Jean-Pierre Renaud. The authors read and approved the final manuscript.

### Funding

This work was supported by the French PIA project "Lorraine Université d'Excellence", reference ANR-15-IDEX-04-LUE, through the project Impact DeepSurf and the Labex ARBRE (ANR-11-LABX-0002-01). The authors also received the support of CNES-TOSCA programme through the project TOSCA Sylvosansat (0010160963). Olivier Bouriaud was supported by Romania National Council for Higher Education Funding, CNFIS, project number CNFIS-FDI-2024-F-0155.

### Data availability

The data that support the findings of this study are available online (Vega et al. 2025; https://doi.org/10.5281/zenodo.15276127).

### Declarations

#### Ethics approval and consent to participate

Not applicable.
**Consent for publication**

All authors gave their informed consent to this publication and its content.

**Competing interests**

The authors declare that they have no conflict of interest.

**Author details**

¹Faculté des Sciences et Technologies - Campus Aiguillettes, Université de Lorraine, Vandœuvre Les Nancy 54506, France. ²Laboratoire d'Inventaire Forestier, Univ Gustave Eiffel, ENSG, IGN, 14 Rue Girardet, Nancy F-54042, France. ³Université de Lorraine, AgroParisTech, INRAE, SILVA, Nancy F-54000, France. ⁴Department of Geography, Centre for Nature-based Climate Solutions, National University of Singapore, Singapore 117546, Singapore. ⁵Ştefan cel Mare University of Suceava, 13 University Street, Suceava 720229, Romania. ⁶Office National des Forêts, Pôle Recherche Développement Innovation, Villers-Lès-Nancy 54600, France.

Received: 25 October 2024 Accepted: 23 May 2025

Published online: 16 June 2025

**References**

Abdullah H, Skidmore AK, Darvishzadeh R, Heurich M (2019) Sentinel-2 accurately maps green-attack stage of European spruce bark beetle (Ips typographus, L.) compared with Landsat-8. Remote Sens Ecol Conserv 5:87-106. https://doi.org/10.1002/rse2.93

Altman NS (1992) An introduction to kernel and nearest-neighbor nonparametric regression. Am Stat 46:175. https://doi.org/10.2307/2685209

Barrett F, McRoberts RE, Tomppo E, Cienciala E, Waser LT (2016) A questionnaire-based review of the operational use of remotely sensed data by national forest inventories. Remote Sens Environ 174:279-289. https://doi.org/10.1016/j.rse.2015.08.029

Baumann A, Roßberg T, Schmitt M (2023) Probabilistic MIMO U-Net: efficient and accurate uncertainty estimation for pixel-wise regression. In: 2023 IEEE/CVF International Conference on Computer Vision Workshops (ICCVW). IEEE, Paris, France. pp. 4500-4508. https://doi.org/10.1109/ICCVW60793.2023.00484

Bell DM, Wilson BT, Werstak CE, Oswalt CM, Perry CH (2022) Examining k-nearest neighbor small area estimation across scales using national forest inventory data. Front for Glob Change 5:763422. https://doi.org/10.3389/ffgc.2022.763422

Bennett AC, McDowell NG, Allen CD, Anderson-Teixeira KJ (2015) Larger trees suffer most during drought in forests worldwide. Nat Plants 1:15139. https://doi.org/10.1038/nplants.2015.139

Bernard D, Doridant JB (2024) Méthode FORDEAD - analyse de la validité des détections d'anomalies de végétation dans le cas des résineux par contrôle sur le terrain - Synthèse des principaux résultats. DRAAF Auvergne-Rhône-Alpes Dir. Régionale Aliment. Agric. For. URL https://draaf.auvergne-rhone-alpes.agriculture.gouv.fr/methode-fordead-analyse-de-la-validité-des-detections-d-anomalies-de-vegetation-a5597.html. Accessed 25 Apr 2025

Boucher P, Hancock S, Orwig D, Duncanson L, Armston J, Tang H et al (2020) Detecting change in forest structure with simulated GEDI LiDAR waveforms: a case study of the Hemlock Woolly Adelgid (HWA; *Adelges tsugae*) infestation. Remote Sens 12:1304. https://doi.org/10.3390/rs12081304

Breidenbach J, Nothdurft A, Kändler G (2010) Comparison of nearest neighbour approaches for small area estimation of tree species-specific forest inventory attributes in central Europe using airborne laser scanner data. Eur J Forest Res 129:833-846. https://doi.org/10.1007/s10342-010-0384-1

Breidenbach J, Næsset E, Gobakken T (2012) Improving k-nearest neighbor predictions in forest inventories by combining high and low density airborne laser scanning data. Remote Sens Environ 117:358-365. https://doi.org/10.1016/j.rse.2011.10.010

Bright BC, Hicke JA, Hudak AT (2012) Estimating aboveground carbon stocks of a forest affected by mountain pine beetle in Idaho using LiDAR and multispectral imagery. Remote Sens Environ 124:270-281. https://doi.org/10.1016/j.rse.2012.05.016

Bright BC, Hudak AT, McGaughey R, Andersen HE, Negrón J (2013) Predicting live and dead tree basal area of bark beetle affected forests from discrete-return LiDAR. Can J Remote Sens 39:599-S111. https://doi.org/10.5589/m13-027

Canelles Q, Aquilué N, James PMA, Lawler J, Brotons L (2021) Global review on interactions between insect pests and other forest disturbances. Landsc Ecol 36:945-972. https://doi.org/10.1007/s10980-021-01209-7

Chirici G, Mura M, McInerney D, Py N, Tomppo EO, Waser LT et al (2016) A meta-analysis and review of the literature on the k-nearest neighbors technique for forestry applications that use remotely sensed data. Remote Sens Environ 176:282-294. https://doi.org/10.1016/j.rse.2016.02.001

Coulston JW, Edgar CB, Westfall JA, Taylor ME (2020) Estimation of forest disturbance from retrospective observations in a broad-scale inventory. Forests 11:1298. https://doi.org/10.3390/f11121298

Coulston JW, Green PC, Radtke PJ, Prisley SP, Brooks EB, Thomas VA et al (2021) Enhancing the precision of broad-scale forestland removals estimates with small area estimation techniques. Forestry 94:427-441. https://doi.org/10.1093/forestry/cpa045

Crookston NL, Finley AO (2008) yalmpute: an R package for kNN imputation. J Stat Softw 23. https://doi.org/10.18637/jss.v023.i10

Deleuze C, Morneau F, Renaud JP, Vivien Y, Rivoire M, Santenoise P et al (2014) Estimer le volume total d'un arbre, quelles que soient l'essence, la taille, la sylviculture, la station. RDV Tech ONF 44:22-32

Deo R, Russell M, Domke G, Andersen H-E, Cohen W, Woodall C (2017) Evaluating site-specific and generic spatial models of aboveground forest biomass based on Landsat time-series and LiDAR strip samples in the eastern USA. Remote Sens 9:598. https://doi.org/10.3390/rs9060598

Dettmann GT, Radtke PJ, Coulston JW, Green PC, Wilson BT, Moisen GG (2022) Review and synthesis of estimation strategies to meet small area needs in forest inventory. Front for Glob Change 5:813569. https://doi.org/10.3389/ffgc.2022.813569

Duncanson L, Kellner JR, Armston J, Dubayah R, Minor DM, Hancock S (2022) Aboveground biomass density models for NASA's Global Ecosystem Dynamics Investigation (GEDI) LiDAR mission. Remote Sens Environ 270:112845. https://doi.org/10.1016/j.rse.2021.112845

Dutrieux R (2023) Cartographie des scolytes par la méthode fordead - Juin 2023. [Data set]. V1. https://doi.org/10.57745/8LJGQ9

Dutrieux R, Férét JB, Ose K (2020) Mise au point d'une méthode reproductible pour le suivi généralisé des dégâts de scolytes par télédétection satel- litaire. RDV Tech ONF 69-70:37-44

Esteban J, McRoberts R, Fernández-Landa A, Tomé J, Næsset E (2019) Estimating forest volume and biomass and their changes using random forests and remotely sensed data. Remote Sens 11:1944. https://doi.org/10.3390/rs11161944

Fassnacht FE, Mager C, Waser LT, Kanjir U, Schäfer J, Buhvald AP et al (2024) Forest practitioners' requirements for remote sensing-based canopy height, wood-volume, tree species, and disturbance products. Forestry 98:233-252. https://doi.org/10.1093/forestry/cpa021

Francini S, McRoberts RE, Giannetti F, Mencucci M, Marchetti M, Scarascia Mugnozza G et al (2020) Near-real time forest change detection using PlanetScope imagery. Eur J Remote Sens 53:233-244. https://doi.org/10.1080/22797254.2020.1806734

Frescino TS, McConville KS, White GW, Toney JC, Moisen GG (2022) Small area estimates for national applications: a database to dashboard strategy using FIESTA. Front for Glob Change 5:779446. https://doi.org/10.3389/ffgc.2022.779446

Gaines GC, Affleck DLR (2021) Small area estimation of postfire tree density using continuous forest inventory data. Front For Glob Change 4:761509. https://doi.org/10.3389/ffgc.2021.761509

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

Ghosh M, Rao JNK (1994) Small area estimation: an appraisal. Stat Sci 9(1):55-76. https://doi.org/10.1214/ss/1177010647

Giacona F, Eckert N, Mainieri R, Martin B, Corona C, Lopez-Saez J et al (2018) Avalanche activity and socio-environmental changes leave strong footprints in forested landscapes: a case study in the Vosges medium-high mountain range. Ann Glaciol 59:111-133. https://doi.org/10.1017/aog.2018.26

Gilles A, Lisein J, Cansell J, Latte N, Piedallu C, Claessens H (2023) Evolution of the bark beetle crisis in spruce (*Picea abies*) forests: a spatial and temporal
remote sensing analysis in Belgium and North- eastern France. https://doi.org/10.21203/rs.3.rs-3051830/v1

Gregoire TG (1998) Design-based and model-based inference in survey sampling: appreciating the difference. Can J for Res 28:1429–1447. https://doi.org/10.1139/x98-166

Guerra-Hernández J, Pascual A, Tupinambá-Simões F, Godinho S,otequim B, Jurado-Varela A et al (2024) Using bi-temporal ALS and NFI-based time-series data to account for large-scale aboveground carbon dynamics: the showcase of Mediterranean forests. Eur J Remote Sens 57:2315413. https://doi.org/10.1080/22797254.2024.2315413

Hagolle O, Huc M, Villa Pascual D, Dedieu G (2015) A multi-temporal and multi-spectral method to estimate aerosol optical thickness over land, for the atmospheric correction of FormoSat-2, LandSat, VENUS and Sentinel-2 images. Remote Sens 7:2668–2691. https://doi.org/10.3390/rs70302668

Hervé JC (2016) France, in: Vidal, C., Alberdi, I.A., Hernández Mateo, L., Redmond, J.J. (Eds.), National forest inventories - assessment of wood availability and use. Springer, Cham, pp. 385–404. https://doi.org/10.1007/978-3-319-44015-6_20

Hervé JC, Wurpillot S, Vidal C, Roman-Amat B (2014) L'inventaire des ressources forestières en France: un nouveau regard sur de nouvelles forêts. Revue Forestière Française 247–260. https://doi.org/10.4267/2042/56055

Hill A, Mandallaz D, Langshausen J (2018) A double-sampling extension of the German national forest inventory for design-based small area estimation on forest district levels. Remote Sens 10:1052. https://doi.org/10.3390/rs10071052

Hill A, Massey A, Mandallaz D (2021) The R package forestinventory: design-based global and small area estimations for multiphase forest inventories. J Stat Softw 97(4):1–40. https://doi.org/10.18637/jss.v097.i04

Hlásny T, König L, Krokene P, Lindner M, Montagné-Huck C, Müller J et al (2021) Bark beetle outbreaks in Europe: state of knowledge and ways forward for management. Curr Forestry Rep 7:138–165. https://doi.org/10.1007/s40725-021-00142-x

Hollaus M, Vreugdenhil M (2019) Radar satellite imagery for detecting bark beetle outbreaks in forests. Curr Forestry Rep 5:240–250. https://doi.org/10.1007/s40725-019-00098-z

Huo L, Persson HJ, Lindberg E (2021) Early detection of forest stress from European spruce bark beetle attack, and a new vegetation index: Normalized distance red & SWIR (NDRS). Remote Sens Environ 255:112240. https://doi.org/10.1016/j.rse.2020.112240

Irulappa-Pillai-Vijayakumar DB, Renaud JP, Morneau F, McRoberts RE, Vega C (2019) Increasing precision for French forest inventory estimates using the k-NN technique with optical and photogrammetric data and model-assisted estimators. Remote Sens 11:991. https://doi.org/10.3390/rs11080991

Kangas A, Astrup R, Breidenbach J, Fridman J, Gobakken T, Korhonen KT et al (2018) Remote sensing and forest inventories in Nordic countries – road-map for the future. Scand J for Res 33:397–412. https://doi.org/10.1080/02827581.2017.1416666

Kangas A, Myllymäki M, Mehtätalo L (2023) Understanding uncertainty in forest resources maps. Silva Fenn 57(2). https://doi.org/10.14214/sf.22026

Karasiak N, Dejoux JF, Monteil C, Sheeren D (2022) Spatial dependence between training and test sets: another pitfall of classification accuracy assessment in remote sensing. Mach Learn 111:2715–2740. https://doi.org/10.1007/s10994-021-05972-1

Kim H-J, Tomppo E (2006) Model-based prediction error uncertainty estimation for kNN method. Remote Sens Environ 104:257–263. https://doi.org/10.1016/j.rse.2006.04.009

Lang N, Kalischek N, Armston J, Schindler K, Dubayah R, Wegner JD (2022) Global canopy height regression and uncertainty estimation from GEDI LIDAR waveforms with deep ensembles. Remote Sens Environ 268:112760. https://doi.org/10.1016/j.rse.2021.112760

Lausch A, Heurich M, Gordalla D, Dobner HJ, Gwillym-Margianto S, Salbach C (2013) Forecasting potential bark beetle outbreaks based on spruce forest vitality using hyperspectral remote-sensing techniques at different scales. For Ecol Manag 308:76–89. https://doi.org/10.1016/j.foreco.2013.07.043

Li Y, Wu Z, Xu X, Fan H, Tong X, Liu J (2021) Forest disturbances and the attribution derived from yearly Landsat time series over 1990–2020 in the Hengduan mountains region of Southwest China. For Ecosyst 8:73. https://doi.org/10.1186/s40663-021-00352-6

Massey R, Berner LT, Foster AC, Goetz SJ, Vepakomma U (2023) Remote sensing tools for monitoring forests and tracking their dynamics, in: Girona, M.M., Morin, H., Gauthier, S., Bergeron, Y. (Eds.), Boreal forests in the face of climate change, Advances in Global Change Research. Springer International Publishing, Cham, pp. 637–655. https://doi.org/10.1007/978-3-031-15988-6_26

McConville KS, Moisen GG, Frescino TS (2020) A tutorial on model-assisted estimation with application to forest inventory. Forests 11:244. https://doi.org/10.3390/f11020244

McRoberts RE (2012) Estimating forest attribute parameters for small areas using nearest neighbors techniques. For Ecol Manag 272:3–12. https://doi.org/10.1016/j.foreco.2011.06.039

McRoberts RE, Næsset E, Saatchi S, Quegan S (2022) Statistically rigorous, model-based inferences from maps. Remote Sens Environ 279:113028. https://doi.org/10.1016/j.rse.2022.113028

Meyer H, Pebesma E (2021) Predicting into unknown space? Estimating the area of applicability of spatial prediction models. Methods Ecol Evol 12:1620–1633. https://doi.org/10.1111/2041-210X.13650

Minárová V, Müller M, Clappier A, Kašpar M (2018) Comparison of extreme precipitation characteristics between the Ore Mountains and the Vosges Mountains (Europe). Theor Appl Climatol 133:1249–1268. https://doi.org/10.1007/s00704-017-2247-x

Morris JL, Cottrell S, Fettig CJ, Hansen WD, Sherriff RL, Carter VA (2017) Managing bark beetle impacts on ecosystems and society: priority questions to motivate future research. J Appl Ecol 54:750–760. https://doi.org/10.1111/1365-2664.12782

Myroniuk V, Bilous A, Khan Y, Terentiev A, Kravets P, Kovalevskyi S et al (2020) Tracking rates of forest disturbance and associated carbon loss in areas of illegal amber mining in Ukraine using Landsat time series. Remote Sens 12:2235. https://doi.org/10.3390/rs12142235

Nilsson M, Nordkvist K, Jonzén J, Lindgren N, Axensten P, Wallerman J et al (2017) A nationwide forest attribute map of Sweden predicted using airborne laser scanning data and field data from the national forest inventory. Remote Sens Environ 194:447–454. https://doi.org/10.1016/j.rse.2016.10.022

Pestov V (2013) Is the k-NN classifier in high dimensions affected by the curse of dimensionality? Comput Math Appl 65:1427–1437. https://doi.org/10.1016/j.camwa.2012.09.011

Piedallu C, Dallery D, Bresson C, Legay M, Gégout JC, Pierrat R (2023) Spatial vulnerability assessment of silver fir and Norway spruce dieback driven by climate warming. Landsc Ecol 38:341–361. https://doi.org/10.1007/s10980-022-01570-1

Potapov P, Li X, Hernandez-Serna A, Tyukavina A, Hansen MC, Kommareddy A et al (2021) 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

Price B, Waser LT, Wang Z, Marty M, Ginzler C, Zellweger F (2020) Predicting biomass dynamics at the national extent from digital aerial photogrammetry. Int J Appl Earth Obs Geoinf 90:102116. https://doi.org/10.1016/j.jag.2020.102116

Rao JNK, Molina I (2015) Small area estimation, 2nd edn. John Wiley & Sons, Hoboken, New Jersey

Renaud JP, Vega C, Durrieu S, Lisein J, Magnussen S, Lejeune P et al (2017) Stand-level wind damage can be assessed using diachronic photogrammetric canopy height models. Ann for Sci 74:74. https://doi.org/10.1007/s13595-017-0669-3

Renaud JP, Sagar A, Barbillon P, Bouriaud O, Deleuze C, Vega C (2022) Characterizing the calibration domain of remote sensing models using convex hulls. Int J Appl Earth Obs Geoinf 112:102939. https://doi.org/10.1016/j.jag.2022.102939

Robert N, Vidal C, Colin A, Hervé J-C, Hamza N, Cluzeau C (2009) Developement of France's national forest inventory. In: Tomppo E, Gschwantner T, Lawrence M, McRoberts RE (eds) National forest inventories – pathways for common reporting. Springer International Publishing, Cham, pp 207–221

Roussel J-R, Auty D, Coops NC, Tompalski P, Goodbody TRH, Meador AS et al (2020) lidR: An R package for analysis of Airborne Laser Scanning (ALS) data. Remote Sens Environ 251:112061. https://doi.org/10.1016/j.rse.2020.112061

Saarela S, Wästlund A, Holmström E, Mensah AA, Holm S, Nilsson M et al (2020) Mapping aboveground biomass and its prediction uncertainty
using LiDAR and field data, accounting for tree-level allometric
and LiDAR model errors. For Ecosyst 7:43. https://doi.org/10.1186/
s40663-020-00245-0

Sagar A, Véga C, Piédallu C, Bouriaud O, Renaud JP (2021) High resolution mapping of forest resources and prediction uncertainty using multi-source inventory approach. In: Proceedings of the SilviLaser Conference 2021, pp. 219–221. https://doi.org/10.34726/wim.1986

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

Schelhaas MJ, Nabuurs GJ, Schuck A (2003) Natural disturbances in the European forests in the 19th and 20th centuries. Glob Chang Biol 9:1620–1633. https://doi.org/10.1046/j.1365-2486.2003.00684.x

Schwartz SEO, Benoit L, Clayton S, Parnes MF, Swenson L, Lowe SR (2023) Climate change anxiety and mental health: environmental activism as buffer. Curr Psychol 42:16708–16721. https://doi.org/10.1007/s12144-022-02735-6

Seidl R, Thom D, Kautz M, Martin-Benito D, Peltoniemi M, Vacchiano G et al (2017) Forest disturbances under climate change. Nat Clim Chang 7:395–402. https://doi.org/10.1038/nclimate3303

Sonobe R, Yamaya Y, Tani H, Wang X, Kobayashi N, Mochizuki K (2018) Crop classification from Sentinel-2-derived vegetation indices using ensemble learning. J Appl Remote Sens 12:1. https://doi.org/10.1117/1JRS.12.026019

Stanke H, Finley AO, Domke GM (2022) Simplifying small area estimation with rFIA: a demonstration of tools and techniques. Front for Glob Change 5:745874. https://doi.org/10.3389/ffgc.2022.745874

Thapa RB, Watanabe M, Motohka T, Shiraishi T, Shimada M (2015) Calibration of aboveground forest carbon stock models for major tropical forests in central Sumatra using airborne LiDAR and field measurement data. IEEE J. Sel Top Appl Earth Obs Remote Sens 8:661–673. https://doi.org/10.1109/JSTARS.2014.2328656

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

Valbuena R, Hernando A, Manzanera JA, Górgens EB, Almeida DRA, Mauro F et al (2017) Enhancing of accuracy assessment for forest above-ground biomass estimates obtained from remote sensing via hypothesis testing and overfitting evaluation. Ecol Model 366:15–26. https://doi.org/10.1016/j.ecolmodel.2017.10.009

Véga C, Renaud JP, Durieu S, Bouvier M (2016) On the interest of penetration depth, canopy area and volume metrics to improve LiDAR-based models of forest parameters. Remote Sens Environ 175:32–42. https://doi.org/10.1016/j.rse.2015.12.039

Véga C, Renaud JP, Sagar A, Bouriaud O (2021) A new small area estimation algorithm to balance between statistical precision and scale. Int J Appl Earth Obs Geoinf 97:102303. https://doi.org/10.1016/j.jag.2021.102303

Véga C, Sagar A, Bouriaud O, Renaud JP, Piedallu C (2025) Raw data for the manuscript entitled "Forest attribute maps: a support for small area estimation of forest disturbances". [Data set]. V1. Zenodo. https://doi.org/10.5281/zenodo.15276127

Westfall J (2015) Spatial-scale considerations for a large-area forest inventory regression model. Forestry 88:267–274. https://doi.org/10.1093/forestry/cpv001

Zheng Z, Zhong Y, Tian S, Ma A, Zhang L (2022) ChangeMask: Deep multi-task encoder-transformer-decoder architecture for semantic change detection. ISPRS J Photogramm Remote Sens 183:228–239. https://doi.org/10.1016/j.isprsjprs.2021.10.015

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.