Yes

**Competing Interests:** No competing interests were disclosed.

**Reviewer Expertise:** My areas of expertise are in lidar processing, individual-based modelling, as well as the creation of simulated forest stands (cf. my 2020 paper on this topic, mentioned in the review), which is very close to what the authors have been working on.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

Author Response 24 Nov 2023

Raphaël Aussenac

We would like to thank the reviewers for their useful comments which have helped significantly improve the manuscript. A special effort has been made to provide an in-depth evaluation of the generated landscapes, to clarify the functioning of our downscaling algorithm and to consolidate the rational of our work.

This revision has led us to optimise our code (better memory management, more effective parallelisation, etc.). It is now much faster, taking less than 1 hour to generate the 42 million trees of our three landscapes, while it used to take "less than 5 hours" to generate the 35 million trees of the Bauges landscape. We updated the introduction accordingly.

**Please find below our replies to your comments.**

The article by Aussenac *et al.* describes a statistical procedure to generate a large data set of individual trees from airborne laser scanning (ALS) and inventory data. The variables include trunk diameter, tree height and species identity, and are provided across three European landscapes. The result is an impressive number of simulated/potential trees, which is a useful data set in forest ecology. As applications, the authors mention studies of scale and (more vaguely) forest management/ecosystem prediction, but one could easily think of a number of other concrete applications, such as input/validation of individual-based models of forest dynamics, or comparisons with automatically mapped tree crowns from airborne imagery, e.g. as in Weinstein *et al.* 2021¹, Ball *et al.* 2022², or spaceborne imagery, as in Tucker, Brandt, Hiernaux, *et al.* 2023³.

I also found the paper generally well-written and with a well-thought through methodology for the mapping. The authors carefully tune their models to obtain optimal performance at every step and clearly have spent considerable amounts of time and effort to improve the prediction of stand attributes. In particular, I found the idea of matching predicted basal area to real stands and then filling in/removing trees until the basal area matches intriguing. This bears similarities with model-based estimations of forest attributes/tree attributes from lidar (Hurtt *et al.* 2004⁴, Taubert *et al.* 2015⁵, Rödig *et al.* 2017⁶, Fischer *et al.*
Yes

**Competing Interests:** No competing interests were disclosed.

**Reviewer Expertise:** My areas of expertise are in lidar processing, individual-based modelling, as well as the creation of simulated forest stands (cf. my 2020 paper on this topic, mentioned in the review), which is very close to what the authors have been working on.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

Author Response 24 Nov 2023

Raphaël Aussenac

We would like to thank the reviewers for their useful comments which have helped significantly improve the manuscript. A special effort has been made to provide an in-depth evaluation of the generated landscapes, to clarify the functioning of our downscaling algorithm and to consolidate the rational of our work.

This revision has led us to optimise our code (better memory management, more effective parallelisation, etc.). It is now much faster, taking less than 1 hour to generate the 42 million trees of our three landscapes, while it used to take "less than 5 hours" to generate the 35 million trees of the Bauges landscape. We updated the introduction accordingly.

**Please find below our replies to your comments.**

The article by Aussenac *et al.* describes a statistical procedure to generate a large data set of individual trees from airborne laser scanning (ALS) and inventory data. The variables include trunk diameter, tree height and species identity, and are provided across three European landscapes. The result is an impressive number of simulated/potential trees, which is a useful data set in forest ecology. As applications, the authors mention studies of scale and (more vaguely) forest management/ecosystem prediction, but one could easily think of a number of other concrete applications, such as input/validation of individual-based models of forest dynamics, or comparisons with automatically mapped tree crowns from airborne imagery, e.g. as in Weinstein *et al.* 2021¹, Ball *et al.* 2022², or spaceborne imagery, as in Tucker, Brandt, Hiernaux, *et al.* 2023³.

I also found the paper generally well-written and with a well-thought through methodology for the mapping. The authors carefully tune their models to obtain optimal performance at every step and clearly have spent considerable amounts of time and effort to improve the prediction of stand attributes. In particular, I found the idea of matching predicted basal area to real stands and then filling in/removing trees until the basal area matches intriguing. This bears similarities with model-based estimations of forest attributes/tree attributes from lidar (Hurtt *et al.* 2004⁴, Taubert *et al.* 2015⁵, Rödig *et al.* 2017⁶, Fischer *et al.*
2020<sup>7)</sup> and shares some of these models' advantages (e.g. more fine-scale distribution of biomass, no shrinking to the mean).

However, like these models, the authors' method also involves a lot of complex modelling steps, and it is in the validation step of the procedure that I see deficiencies that need to be addressed. I see two main issues:

a) the robustness of the models to extrapolation issues and spatial autocorrelation is not evaluated, so it is hard to assess how good the models are outside their calibration range and how much we can trust the predictions across the landscape. **The leave-one-out cross validation (LOOCV) performed on our entire workflow presented in this revised version of the article gives an indication on the robustness of our approach (see responses to the comments on model validation).**

Regarding lidar modelling, cross-validation (such as used in this study) remains a common method for accuracy assessment of landscape predictions, although some limitations have been highlighted in the literature. Accuracy assessment remains an important research topic, with no consensus on the best practices, as field data are usually scarce and modeling approaches complex (see paragraph 2.5 of the review in https://doi.org/10.1016/j.rse.2021.112477).

We acknowledge that spatial autocorrelation is necessarily present both in the lidar acquisition and in the forest structure. Hence spatial autocorrelation can be expected in the lidar model predictions and errors. Unfortunately, calibration plots were sampled on a regularly spaced grid in each study area which does not make it possible to quantify spatial autocorrelation at distances smaller than the grid step.

b) two of the key attributes of the data set (tree diameter and species identity) are not validated at all, despite featuring prominently in the title and in the results section (Figure 5). This should be a priority in a revised version. **In this revised version, we provide an in-depth evaluation of the generated landscapes (see responses to the specific comments below).** We changed the abstract accordingly: "We carried out an in-depth evaluation of our workflow including, among other analyses, a leave-one-out cross validation. Overall, the landscapes we generated are in good agreement with the landscapes they aim to reproduce. In the most favourable conditions, the root mean square error (RMSE) of stand basal area (BA) and mean quadratic diameter (Dg) predictions were respectively 5.4 m<sup>2</sup> and 3.9 cm, and the generated main species corresponded to the observed main species in 76.2% of cases."

We also modified the paragraph presenting the evaluation in the General Approach section: "We evaluated the overall reliability of our workflow, i.e. its ability to produce virtual landscapes as close as possible to the real ones (see Dataset validation). In particular we carried out a leave-one-out cross validation (LOOCV) on our entire workflow. This analysis consisted in: • comparing the observed and predicted values of BA, Dg, Bab and the quantiles of tree height and diameter; • comparing the observed and predicted values of species abundance at the landscape level; • calculating the frequency at which the most abundant species was correctly predicted at the cell
2020<sup>7)</sup> and shares some of these models' advantages (e.g. more fine-scale distribution of biomass, no shrinking to the mean).

However, like these models, the authors' method also involves a lot of complex modelling steps, and it is in the validation step of the procedure that I see deficiencies that need to be addressed. I see two main issues:

a) the robustness of the models to extrapolation issues and spatial autocorrelation is not evaluated, so it is hard to assess how good the models are outside their calibration range and how much we can trust the predictions across the landscape. **The leave-one-out cross validation (LOOCV) performed on our entire workflow presented in this revised version of the article gives an indication on the robustness of our approach (see responses to the comments on model validation).**

Regarding lidar modelling, cross-validation (such as used in this study) remains a common method for accuracy assessment of landscape predictions, although some limitations have been highlighted in the literature. Accuracy assessment remains an important research topic, with no consensus on the best practices, as field data are usually scarce and modeling approaches complex (see paragraph 2.5 of the review in https://doi.org/10.1016/j.rse.2021.112477).

We acknowledge that spatial autocorrelation is necessarily present both in the lidar acquisition and in the forest structure. Hence spatial autocorrelation can be expected in the lidar model predictions and errors. Unfortunately, calibration plots were sampled on a regularly spaced grid in each study area which does not make it possible to quantify spatial autocorrelation at distances smaller than the grid step.

b) two of the key attributes of the data set (tree diameter and species identity) are not validated at all, despite featuring prominently in the title and in the results section (Figure 5). This should be a priority in a revised version. **In this revised version, we provide an in-depth evaluation of the generated landscapes (see responses to the specific comments below).** We changed the abstract accordingly: "We carried out an in-depth evaluation of our workflow including, among other analyses, a leave-one-out cross validation. Overall, the landscapes we generated are in good agreement with the landscapes they aim to reproduce. In the most favourable conditions, the root mean square error (RMSE) of stand basal area (BA) and mean quadratic diameter (Dg) predictions were respectively 5.4 m<sup>2</sup> and 3.9 cm, and the generated main species corresponded to the observed main species in 76.2% of cases."

We also modified the paragraph presenting the evaluation in the General Approach section: "We evaluated the overall reliability of our workflow, i.e. its ability to produce virtual landscapes as close as possible to the real ones (see Dataset validation). In particular we carried out a leave-one-out cross validation (LOOCV) on our entire workflow. This analysis consisted in: • comparing the observed and predicted values of BA, Dg, Bab and the quantiles of tree height and diameter; • comparing the observed and predicted values of species abundance at the landscape level; • calculating the frequency at which the most abundant species was correctly predicted at the cell
level; As a complement, we also compared the stands dominant heights measured by ALS (HdomALS) to those calculated from the trees we generated (HdomT). Finally, we compared the spatial distribution of species to current expert knowledge.”

In the following I will provide a few comments on the article following roughly the overall structure, and give suggestions on how to improve the model validation.

## Justification for the data set

I see the value of a fine-grain large-scale data set, and having such a data set is indeed rare, but it would be helpful to mention concrete applications. At the moment, the only justification given is the sentence: "Yet, this type of data could help address the scaling issues in ecology and could prove useful for testing forest management strategies and accurately predicting the dynamics of ecosystem services". This is the sentence from the abstract, but the same point is made at the end of the first paragraph. Could the authors rephrase and add literature references in the main text? The vast majority of data sets can be useful for the testing of forest management strategies or predicting dynamics of ecosystem services. What is unique to your data set? Why do we need detailed, tree-based data at large scales? **As suggested, we rephrased this part and added some references:** "Yet, fine-grain descriptions of large forest areas could help address the pervasive scaling issues in forest ecology, modelling and management. In practice, such data could help better understand at which spatial scale ecological processes emerge in forest ecosystems [Craven et al., 2020; With, 2019]. They could also be extremely valuable to compare forest dynamics models operating at different scales (organ, tree, stand, landscapes) and evaluate their validity across scales [Papaik et al., 2010]. They could ultimately help develop and test management strategies at different spatial scales [Seidl et al., 2013]." Added references: Craven, Dylan, et al. "A cross-scale assessment of productivity-diversity relationships." Global Ecology and Biogeography 29.11 (2020): 1940-1955. Papaik, Michael J., et al. "Forest processes from stands to landscapes: exploring model forecast uncertainties using cross-scale model comparison." Canadian journal of forest research 40.12 (2010): 2345-2359. Seidl, Rupert, et al. "deepl" European Journal of Forest Research 132 (2013): 653-666.

## Model for mapping of tree attributes

ALS metrics: which metrics precisely did you use? **Candidate metrics for selection are those computed with the aba_METRICS and std_tree_METRICS functions from the lidaRtRee R package, as stated in section ALS mapping. Selected metrics result from an automated procedure applied separately on each study area, forest parameter, stratum (in case a specific model is calibrated for each stratum). The following sentence was added to the Stratification" subsection of the ALS mapping section: "The metrics selected in the 32 models for BA and Dg (which include at most six independent variables) are presented in Table S1 of the Extended data."**

The seven stratum-specific logistic regression models obtained for BA<sub>b</sub> with stepAIC variable selection include 20 to 30 metrics and are not presented.
level; As a complement, we also compared the stands dominant heights measured by ALS (HdomALS) to those calculated from the trees we generated (HdomT). Finally, we compared the spatial distribution of species to current expert knowledge.”

In the following I will provide a few comments on the article following roughly the overall structure, and give suggestions on how to improve the model validation.

## Justification for the data set

I see the value of a fine-grain large-scale data set, and having such a data set is indeed rare, but it would be helpful to mention concrete applications. At the moment, the only justification given is the sentence: "Yet, this type of data could help address the scaling issues in ecology and could prove useful for testing forest management strategies and accurately predicting the dynamics of ecosystem services". This is the sentence from the abstract, but the same point is made at the end of the first paragraph. Could the authors rephrase and add literature references in the main text? The vast majority of data sets can be useful for the testing of forest management strategies or predicting dynamics of ecosystem services. What is unique to your data set? Why do we need detailed, tree-based data at large scales? **As suggested, we rephrased this part and added some references:** "Yet, fine-grain descriptions of large forest areas could help address the pervasive scaling issues in forest ecology, modelling and management. In practice, such data could help better understand at which spatial scale ecological processes emerge in forest ecosystems [Craven et al., 2020; With, 2019]. They could also be extremely valuable to compare forest dynamics models operating at different scales (organ, tree, stand, landscapes) and evaluate their validity across scales [Papaik et al., 2010]. They could ultimately help develop and test management strategies at different spatial scales [Seidl et al., 2013]." Added references: Craven, Dylan, et al. "A cross-scale assessment of productivity-diversity relationships." Global Ecology and Biogeography 29.11 (2020): 1940-1955. Papaik, Michael J., et al. "Forest processes from stands to landscapes: exploring model forecast uncertainties using cross-scale model comparison." Canadian journal of forest research 40.12 (2010): 2345-2359. Seidl, Rupert, et al. "deepl" European Journal of Forest Research 132 (2013): 653-666.

## Model for mapping of tree attributes

ALS metrics: which metrics precisely did you use? **Candidate metrics for selection are those computed with the aba_METRICS and std_tree_METRICS functions from the lidaRtRee R package, as stated in section ALS mapping. Selected metrics result from an automated procedure applied separately on each study area, forest parameter, stratum (in case a specific model is calibrated for each stratum). The following sentence was added to the Stratification" subsection of the ALS mapping section: "The metrics selected in the 32 models for BA and Dg (which include at most six independent variables) are presented in Table S1 of the Extended data."**

The seven stratum-specific logistic regression models obtained for BA<sub>b</sub> with stepAIC variable selection include 20 to 30 metrics and are not presented.
Point cloud properties: Could the authors add information on/discussion of the sensitivity of their point metrics to scanner acquisitions? **Visual inspection of metrics maps was performed to exclude metrics exhibiting spatial patterns which might be linked to acquisition. We hypothesise that the selected metrics and models based on them are robust to variations in acquisition settings, as a systematic sampling of field plots should ensure that the heterogeneities of lidar acquisition and forest structure on the whole landscape are represented in the data. In the case of Milicz, several metrics seemed quite sensitive to points acquired with a large scan angle. Scan angle values above 21 degrees were excluded in order to achieve a trade-off between metrics robustness, point density and comprehensive coverage of the study area. In the case of Sneznik, some intensity-related metrics where also found to be influenced by acquisition pattern and where thus removed from the analysis ("imean", "imax", "isd", "iskew", "ikurt"). The following sentences were added to the manuscript: "Each metric map was visually checked for spatial patterns potentially linked to acquisition patterns, which eventually led to: - discard some intensity-related metrics in Sneznik study area; - remove ALS points acquired with a scan angle larger than 21 degrees in Milicz study area, in order to achieve a trade-off between metrics robustness, point density and comprehensive coverage of the study area."

Lidar scans often exhibit considerable variation in pulse density even within a single acquisition (e.g. scan line centre vs. overlapping scan lines). What is each scan's standard deviation of point/pulse density? **The mean +- sd point density (/m2) for each study site, computed at 25 m resolution for pixels with at least 1 point above 2 m height is :** Bauges (southern part) : 5.9 +- 3.1 Bauges (northern part) : 27.6 +- 13.3 Milicz : 16.5 +- 7.1 Sneznik : 18.4 +- 10.1 Those values are now indicated in the article. Mean values slightly differ from the previous ones which were computed for pixels of the whole landscapes.

Could you include that as a variable in stratification? Could this improve your models (e.g. stratify by pulse densities between 5 and 10, 10 and 15, 15 and 20, etc., or even smaller step sizes)? **In order to obtain robust models and valid inferences, [https://doi.org/10.5589/m12-052] suggests to use at least 55 observations in calibration, which limits the number of strata that the whole dataset can be partitioned in. We hypothesised that a stratification based on forest structure defined by ancillary data would help in reducing the prediction error. For the Bauges study area where ALS data originate from two different ALS campaigns, campaign was also tested as stratification criterion. Several stratifications where tested, and the one retained for the final map production is the one with the best RMSE improvement in cross validation compared to the single model. Parcimony was also considered when selecting the model. We have not tested a stratification based on pulse densities, as most computed metrics are a priori robust to density and as we have checked that metrics maps did not display any pattern due to acquisition conditions, such as density variations.**

**Descriptions:** I appreciate that the paper is already quite dense, but quite a few steps in the methods section remain unclear to me, particularly in step 3. E.g., in the matching of BA and BAb, why do you need a correction value alpha? Can you explain the weighting better and
Point cloud properties: Could the authors add information on/discussion of the sensitivity of their point metrics to scanner acquisitions? **Visual inspection of metrics maps was performed to exclude metrics exhibiting spatial patterns which might be linked to acquisition. We hypothesise that the selected metrics and models based on them are robust to variations in acquisition settings, as a systematic sampling of field plots should ensure that the heterogeneities of lidar acquisition and forest structure on the whole landscape are represented in the data. In the case of Milicz, several metrics seemed quite sensitive to points acquired with a large scan angle. Scan angle values above 21 degrees were excluded in order to achieve a trade-off between metrics robustness, point density and comprehensive coverage of the study area. In the case of Sneznik, some intensity-related metrics where also found to be influenced by acquisition pattern and where thus removed from the analysis ("imean", "imax", "isd", "iskew", "ikurt"). The following sentences were added to the manuscript: "Each metric map was visually checked for spatial patterns potentially linked to acquisition patterns, which eventually led to: - discard some intensity-related metrics in Sneznik study area; - remove ALS points acquired with a scan angle larger than 21 degrees in Milicz study area, in order to achieve a trade-off between metrics robustness, point density and comprehensive coverage of the study area."

Lidar scans often exhibit considerable variation in pulse density even within a single acquisition (e.g. scan line centre vs. overlapping scan lines). What is each scan's standard deviation of point/pulse density? **The mean +- sd point density (/m2) for each study site, computed at 25 m resolution for pixels with at least 1 point above 2 m height is :** Bauges (southern part) : 5.9 +- 3.1 Bauges (northern part) : 27.6 +- 13.3 Milicz : 16.5 +- 7.1 Sneznik : 18.4 +- 10.1 Those values are now indicated in the article. Mean values slightly differ from the previous ones which were computed for pixels of the whole landscapes.

Could you include that as a variable in stratification? Could this improve your models (e.g. stratify by pulse densities between 5 and 10, 10 and 15, 15 and 20, etc., or even smaller step sizes)? **In order to obtain robust models and valid inferences, [https://doi.org/10.5589/m12-052] suggests to use at least 55 observations in calibration, which limits the number of strata that the whole dataset can be partitioned in. We hypothesised that a stratification based on forest structure defined by ancillary data would help in reducing the prediction error. For the Bauges study area where ALS data originate from two different ALS campaigns, campaign was also tested as stratification criterion. Several stratifications where tested, and the one retained for the final map production is the one with the best RMSE improvement in cross validation compared to the single model. Parcimony was also considered when selecting the model. We have not tested a stratification based on pulse densities, as most computed metrics are a priori robust to density and as we have checked that metrics maps did not display any pattern due to acquisition conditions, such as density variations.**

**Descriptions:** I appreciate that the paper is already quite dense, but quite a few steps in the methods section remain unclear to me, particularly in step 3. E.g., in the matching of BA and BAb, why do you need a correction value alpha? Can you explain the weighting better and
why it is divided by 16? Maybe this is more exhaustively explained in the Extended Data, but
this needs to be clear from the main text already. **We have gone into more detail on
points 3a and 3b. They now read as follows:**

o "(a) For this, we first calculated α, a multiplier correction coefficient to be applied to all tree diameters of a field plot. The idea is to increase or decrease tree diameters so that their Dg reaches the Dg value of the cell to which they are associated. α is given by:"

o "(b) Thereafter, we calculated the weight (ω in n.ha⁻¹) of these trees with corrected diameters, so that the generated stand matches the BA and BA b values of the cell to which it is associated. ω is given by:"

As for the weighting, it is now explained (in a foot note and in point 3c), but see
answer to Reviewer 2's comment on scale factors.

Model validation

As pointed out above, this is the point of the paper that needs to be more comprehensive. At the moment, the authors validate their approach by comparing dominant height, as obtained from lidar (mean height of six highest local maxima), to dominant height of the simulated stands, obtained via local allometries (mean height of six highest trees). It is definitely useful to do this comparison and good to see that the results are broadly consistent, so I would keep it in the paper. However, there are issues with circularity, as the authors first use a number of lidar metrics that involve height / basal area-to-height relationships to create the maps and then compare the inferred results (+ independently derived height allometries) again to lidar-derived height metrics. Furthermore, height of the dominant trees may be related to basal area, but it cannot be used to evaluate basal area/tree diameter predictions as such, nor does it validate predicted species composition - both are key features of the data set. **We agree with the reviewer, the comparison of Hdom only partially validates our virtual landscapes. The revised version of the article includes a much more comprehensive evaluation of the virtual landscapes.**

We added this point on circularity in the section presenting the comparation of Hdom
as follows: "There is some circularity in comparing Hdom ALS and Hdom T as models
predicting BA, Dg and BA b from ALS point clouds may include ALS derived height
metrics or more generally metrics which are correlated with the dominant height
estimated from ALS point clouds. The results of this comparison must therefore be
interpreted with caution."

Given that the author's simulation approach seems fast (only ca. 5 hours on a modern
laptop, amazing!), another approach suggests itself, namely within-site cross-validation,
ideally in the form proposed by Ploton et al. 2020⁸. Since a spatially explicit leave-one-out
cross-validation, as suggested in Ploton et al. 2020⁸, may be too computationally intensive, I
would recommend the simpler approach proposed in the same paper: for each of the
European landscapes, I would recommend the authors to split their field data sets into, e.g.,
5 spatially aggregated folds (i.e., spatial clusters), and run their model 5 times, each times
using 4 folds to train the model and 1 separate geographic fold of plots to validate the
model. In this 1 fold, the authors could directly compare predictions of tree values to actual
data according to some simple standard metrics (total basal area, mean quadratic diameter,
why it is divided by 16? Maybe this is more exhaustively explained in the Extended Data, but
this needs to be clear from the main text already. **We have gone into more detail on
points 3a and 3b. They now read as follows:**

o "(a) For this, we first calculated α, a multiplier correction coefficient to be applied to all tree diameters of a field plot. The idea is to increase or decrease tree diameters so that their Dg reaches the Dg value of the cell to which they are associated. α is given by:"

o "(b) Thereafter, we calculated the weight (ω in n.ha⁻¹) of these trees with corrected diameters, so that the generated stand matches the BA and BA b values of the cell to which it is associated. ω is given by:"

As for the weighting, it is now explained (in a foot note and in point 3c), but see
answer to Reviewer 2's comment on scale factors.

Model validation

As pointed out above, this is the point of the paper that needs to be more comprehensive. At the moment, the authors validate their approach by comparing dominant height, as obtained from lidar (mean height of six highest local maxima), to dominant height of the simulated stands, obtained via local allometries (mean height of six highest trees). It is definitely useful to do this comparison and good to see that the results are broadly consistent, so I would keep it in the paper. However, there are issues with circularity, as the authors first use a number of lidar metrics that involve height / basal area-to-height relationships to create the maps and then compare the inferred results (+ independently derived height allometries) again to lidar-derived height metrics. Furthermore, height of the dominant trees may be related to basal area, but it cannot be used to evaluate basal area/tree diameter predictions as such, nor does it validate predicted species composition - both are key features of the data set. **We agree with the reviewer, the comparison of Hdom only partially validates our virtual landscapes. The revised version of the article includes a much more comprehensive evaluation of the virtual landscapes.**

We added this point on circularity in the section presenting the comparation of Hdom
as follows: "There is some circularity in comparing Hdom ALS and Hdom T as models
predicting BA, Dg and BA b from ALS point clouds may include ALS derived height
metrics or more generally metrics which are correlated with the dominant height
estimated from ALS point clouds. The results of this comparison must therefore be
interpreted with caution."

Given that the author's simulation approach seems fast (only ca. 5 hours on a modern
laptop, amazing!), another approach suggests itself, namely within-site cross-validation,
ideally in the form proposed by Ploton et al. 2020⁸. Since a spatially explicit leave-one-out
cross-validation, as suggested in Ploton et al. 2020⁸, may be too computationally intensive, I
would recommend the simpler approach proposed in the same paper: for each of the
European landscapes, I would recommend the authors to split their field data sets into, e.g.,
5 spatially aggregated folds (i.e., spatial clusters), and run their model 5 times, each times
using 4 folds to train the model and 1 separate geographic fold of plots to validate the
model. In this 1 fold, the authors could directly compare predictions of tree values to actual
data according to some simple standard metrics (total basal area, mean quadratic diameter,
95th percentile of diameter, percentage of species xyz, 95th percentile of height, mean height, dominant height). For comparison and to broadly assess whether spatial autocorrelation makes a difference, the authors could do the same validation procedure also with 5 folds containing plots randomly distributed in space (so no spatial clusters). This would only take 25 hours for each validation and give a good impression of how easy it is to accurately map individual trees and species at landscape scale and how realistic the produced inventories are. It would likely also increase interest in the data set, as it would give potential users higher confidence in the results. **We carried out a leave-one-out cross-validation (LOOCV). This validation is actually not too computationally intensive because it has to run only on a few hundred field plots and not on the whole landscape. We also preferred this LOOCV to a k-fold cross-validation because the latter would require to check that sample size in each stratum would remain relevant in each calibration fold, while performing repetitions of the k-fold samples to assess the reliability of the statistics. Besides we are not certain that spatially correlated folds will help in understanding the robustness of the map. One can expect a model calibrated on a spatially-selected subset to perform poorly on another spatial subset as the predicted area may be different from the calibration area both for the lidar and forest. Spatial k-folds do not allow to assess the local spatial variations of performance of a model calibrated with systematic samples (which is what our models are) in a better way than what a LOOCV does. In both cases, spatial correlation at distances smaller than the grid step, which is interesting e.g. for small area estimations, is impossible to evaluate.**

We added the following paragraph to the method subsection of the validation section: "We carried out a leave-one-out cross validation (LOOCV) to evaluate the realism of the virtual landscapes we generated. This consisted in excluding a field plot from our entire workflow and comparing the predicted values obtained to the observed values. This operation was repeated within each landscape for all field plots. We calculated the root mean square error (RMSE) of the predictions of BA, Dg, Bab and the quantiles of tree height and diameter. As part of the LOOCV, we also compared the observed and predicted values of species abundance at the landscape level (in BA) and calculated the frequency at which the most abundant species was correctly predicted at the cell level."

We presented the results of this new evaluation procedure in the results subsection of the validation section: "Overall, the virtual landscapes are in good agreement with the landscapes they aim to reproduce. The generated stand structures and compositions are consistent with the observations and make it possible to distinguish stands at different stages of development and with different compositions.

At Milicz, predictions are the most accurate. The RMSE of all evaluated variables are the lowest in comparison with the other landscapes (Table T1, Figure 4). Species abundance at the landscape level is also better reproduced (Figure F1). Finally, in 76.2% of cases, the generated main species corresponds to the observed main species. This higher quality of predictions can be explained by the fact that Milicz has the highest density of inventory plots and the least complex landscape, with a predominance of even-aged monospecific stands and the lowest species diversity
95th percentile of diameter, percentage of species xyz, 95th percentile of height, mean height, dominant height). For comparison and to broadly assess whether spatial autocorrelation makes a difference, the authors could do the same validation procedure also with 5 folds containing plots randomly distributed in space (so no spatial clusters). This would only take 25 hours for each validation and give a good impression of how easy it is to accurately map individual trees and species at landscape scale and how realistic the produced inventories are. It would likely also increase interest in the data set, as it would give potential users higher confidence in the results. **We carried out a leave-one-out cross-validation (LOOCV). This validation is actually not too computationally intensive because it has to run only on a few hundred field plots and not on the whole landscape. We also preferred this LOOCV to a k-fold cross-validation because the latter would require to check that sample size in each stratum would remain relevant in each calibration fold, while performing repetitions of the k-fold samples to assess the reliability of the statistics. Besides we are not certain that spatially correlated folds will help in understanding the robustness of the map. One can expect a model calibrated on a spatially-selected subset to perform poorly on another spatial subset as the predicted area may be different from the calibration area both for the lidar and forest. Spatial k-folds do not allow to assess the local spatial variations of performance of a model calibrated with systematic samples (which is what our models are) in a better way than what a LOOCV does. In both cases, spatial correlation at distances smaller than the grid step, which is interesting e.g. for small area estimations, is impossible to evaluate.**

We added the following paragraph to the method subsection of the validation section: "We carried out a leave-one-out cross validation (LOOCV) to evaluate the realism of the virtual landscapes we generated. This consisted in excluding a field plot from our entire workflow and comparing the predicted values obtained to the observed values. This operation was repeated within each landscape for all field plots. We calculated the root mean square error (RMSE) of the predictions of BA, Dg, Bab and the quantiles of tree height and diameter. As part of the LOOCV, we also compared the observed and predicted values of species abundance at the landscape level (in BA) and calculated the frequency at which the most abundant species was correctly predicted at the cell level."

We presented the results of this new evaluation procedure in the results subsection of the validation section: "Overall, the virtual landscapes are in good agreement with the landscapes they aim to reproduce. The generated stand structures and compositions are consistent with the observations and make it possible to distinguish stands at different stages of development and with different compositions.

At Milicz, predictions are the most accurate. The RMSE of all evaluated variables are the lowest in comparison with the other landscapes (Table T1, Figure 4). Species abundance at the landscape level is also better reproduced (Figure F1). Finally, in 76.2% of cases, the generated main species corresponds to the observed main species. This higher quality of predictions can be explained by the fact that Milicz has the highest density of inventory plots and the least complex landscape, with a predominance of even-aged monospecific stands and the lowest species diversity
among our three landscapes.

At the Bauges and Sneznik, the RMSE of the evaluated variables are comparable (Table T1., Figure 4). In contrast, predictions of species abundance at the landscape level are more accurate at Sneznik (Figure F1). The same applies to the compositions predicted at the plot level: the predicted main species corresponds to the observed main species in 63.1% of cases at Sneznik and in 37.2% of cases at the Bauges. However, two datasets were used in the Bauges. In the local forest inventory (LFI) not all trees were identified at the species level and trees with a dbh between 7.5 and 17.5 cm were not measured but counted by diameter classes and grouped in two categories (coniferous and broadleaf). This led us to use a local subset of the NFI from which composition is derived in our downscaling algorithm. The poorer composition predictions in the Bauges might therefore partly be an artefact arising from the evaluation itself, as the LFI may not be suitable to serve as a field reference.

The fact that the RMSE values obtained from the LOOCV carried out on our entire workflow are almost similar to the RMSE values obtained from the LOOCV of ALS models shows that the downscaling algorithm hardly adds any error (Table 1, Table T1). The main way of increasing the realism of our virtual landscapes would therefore be to improve the ALS models."

Finally, we added the RMSE values associated to the three scatter plots of our Figure 4 for greater consistency in the document and modified the paragraph where this figure is presented: "With $R^2$ values ranging from 0.61 to 0.83 (Figure 4) and RMSE values below 5 m, $H_{dom\_ALS}$ and $H_{dom\_T}$ are consistent with one another. This provides a general validation of our workflow. As discussed above, the better predictions obtained at Milicz might steam from the higher density of inventory plots and the lower complexity of the landscape. At Sneznik, $H_{dom\_T}$ tends to be overestimated as $H_{dom\_ALS}$ decreases. This divergence could be due to the ice storm that occurred between the field inventory and the ALS acquisition and that might have biased the ALS models."

Since the paper puts its focus on the value of individual trees, there should, in my opinion, also be one result/validation graph that shows individual trees in some way. It could be, for example, a zoomed-in image of lidar-derived canopy height models + a predicted distribution of trees. If the 5-fold cross-validation is carried out, as above, the authors could simply show sample lidar canopy height models on top of plots, and the diameter distributions for the simulated and the inferred plots. **Here, we do not exactly agree with the reviewer. We do not generate individual trees, per se, but rather a list of trees for each cell, and those generated trees are not spatialised within the cells. For such a purpose different methods and denser ALS data would be required. Also, it is not clear to us how the comparison of canopy height models and diameter distributions would constitute an evaluation of our algorithm. The comparison of diameter and height quantiles carried out in response to the previous comment seems to us a better way to evaluate the quality of the generated tree lists.**

Overall, it would also be interesting to readers to understand in how far the predicted
among our three landscapes.

At the Bauges and Sneznik, the RMSE of the evaluated variables are comparable (Table T1., Figure 4). In contrast, predictions of species abundance at the landscape level are more accurate at Sneznik (Figure F1). The same applies to the compositions predicted at the plot level: the predicted main species corresponds to the observed main species in 63.1% of cases at Sneznik and in 37.2% of cases at the Bauges. However, two datasets were used in the Bauges. In the local forest inventory (LFI) not all trees were identified at the species level and trees with a dbh between 7.5 and 17.5 cm were not measured but counted by diameter classes and grouped in two categories (coniferous and broadleaf). This led us to use a local subset of the NFI from which composition is derived in our downscaling algorithm. The poorer composition predictions in the Bauges might therefore partly be an artefact arising from the evaluation itself, as the LFI may not be suitable to serve as a field reference.

The fact that the RMSE values obtained from the LOOCV carried out on our entire workflow are almost similar to the RMSE values obtained from the LOOCV of ALS models shows that the downscaling algorithm hardly adds any error (Table 1, Table T1). The main way of increasing the realism of our virtual landscapes would therefore be to improve the ALS models."

Finally, we added the RMSE values associated to the three scatter plots of our Figure 4 for greater consistency in the document and modified the paragraph where this figure is presented: "With R² values ranging from 0.61 to 0.83 (Figure 4) and RMSE values below 5 m, Hdom ALS and Hdom T are consistent with one another. This provides a general validation of our workflow. As discussed above, the better predictions obtained at Milicz might steam from the higher density of inventory plots and the lower complexity of the landscape. At Sneznik, Hdom T tends to be overestimated as Hdom ALS decreases. This divergence could be due to the ice storm that occurred between the field inventory and the ALS acquisition and that might have biased the ALS models."

Since the paper puts its focus on the value of individual trees, there should, in my opinion, also be one result/validation graph that shows individual trees in some way. It could be, for example, a zoomed-in image of lidar-derived canopy height models + a predicted distribution of trees. If the 5-fold cross-validation is carried out, as above, the authors could simply show sample lidar canopy height models on top of plots, and the diameter distributions for the simulated and the inferred plots. **Here, we do not exactly agree with the reviewer. We do not generate individual trees, per se, but rather a list of trees for each cell, and those generated trees are not spatialised within the cells. For such a purpose different methods and denser ALS data would be required. Also, it is not clear to us how the comparison of canopy height models and diameter distributions would constitute an evaluation of our algorithm. The comparison of diameter and height quantiles carried out in response to the previous comment seems to us a better way to evaluate the quality of the generated tree lists.**

Overall, it would also be interesting to readers to understand in how far the predicted
species distributions reflect current expert knowledge, but this is not a necessity. **In the method subsection of the validation section, we added: "Finally, we examined the spatial distribution of species at each site and compared it to current expert knowledge."**

And in the results subsection of the validation section, we added: "Overall, species spatial distribution in the virtual landscapes is consistent with field observations. In the Bauges, pure and mixed stands of fir and spruce are more abundant at higher elevation while mixed stands of broadleaf species are found at lower elevation. At Milicz, pure stands of Scots pine are found at lower elevation while broadleaf species and mixed stands appear at higher elevation. Finally, at Sneznik, pure beech stands are found at higher elevation while fir is found at lower elevation in pure or mixed stands (a specific feature of the site)."

## Data set

I had a quick look at the data set. One variable I did not understand was the variable "n" or "number of trees". Could you explain it a bit better? Does this mean that the specific diameter exists n times in the specific data set? If this is true (and only in this case), I seem to get some cells (very few) of 25m by 25m (e.g. cellID25 = 2439821 in the "Bauges" data set) that contain more than 500 trees with dbh >= 9-10cm per 625m2 and a total basal area >= 6m2 (which would yield roughly 100m² per hectare, at densities of 8000 trees). These are outliers, and every model is allowed to have outliers (and nature is full of them too), but it would be interesting to get your take on that in terms of realism/stand type. It could also be part of the validation to assess the edges of the basal area distribution or to give readers a hint what to make of the most extreme values. **The variable "n" indeed means that a specific diameter exists n times. We clarified that in the presentation of our dataset on the zenodo website and in the Extended data section as follows: "- n: number of trees. n is an integer >= 1, meaning that a specific set of species "sp", diameter "dbh" and height "h" can be present multiple times in a cell."**

The extreme values you identified are indeed outliers. They are a direct consequence of the uncertainties associated with the models. The realism of the stands associated with these extreme values is, of course, open to question. However, it seems difficult to define thresholds separating realistic from unrealistic stands as extreme values can be locally observed. We believe it is up to the users of the dataset to decide whether or not to use all the stands depending on their objectives. We drew the reader's attention to the outliers at the end of the validation section as follows: "Our procedure is not free of flaws and some outliers are present in the generated data (i.e. stands with extreme values of BA, Dg, tree height or density). These outliers are a direct consequence of the uncertainties associated with the models we used. The realism of the stands associated with these extreme values is open to question. However, separating realistic from unrealistic stands seems difficult as extreme values can be locally observed. It is therefore up to the users of the dataset to decide whether or not to consider these stands depending on their objectives."

Finally, it is always tricky to assess the relevance of extreme values as the field data
species distributions reflect current expert knowledge, but this is not a necessity. **In the method subsection of the validation section, we added: "Finally, we examined the spatial distribution of species at each site and compared it to current expert knowledge."**

And in the results subsection of the validation section, we added: "Overall, species spatial distribution in the virtual landscapes is consistent with field observations. In the Bauges, pure and mixed stands of fir and spruce are more abundant at higher elevation while mixed stands of broadleaf species are found at lower elevation. At Milicz, pure stands of Scots pine are found at lower elevation while broadleaf species and mixed stands appear at higher elevation. Finally, at Sneznik, pure beech stands are found at higher elevation while fir is found at lower elevation in pure or mixed stands (a specific feature of the site)."

## Data set

I had a quick look at the data set. One variable I did not understand was the variable "n" or "number of trees". Could you explain it a bit better? Does this mean that the specific diameter exists n times in the specific data set? If this is true (and only in this case), I seem to get some cells (very few) of 25m by 25m (e.g. cellID25 = 2439821 in the "Bauges" data set) that contain more than 500 trees with dbh >= 9-10cm per 625m2 and a total basal area >= 6m2 (which would yield roughly 100m² per hectare, at densities of 8000 trees). These are outliers, and every model is allowed to have outliers (and nature is full of them too), but it would be interesting to get your take on that in terms of realism/stand type. It could also be part of the validation to assess the edges of the basal area distribution or to give readers a hint what to make of the most extreme values. **The variable "n" indeed means that a specific diameter exists n times. We clarified that in the presentation of our dataset on the zenodo website and in the Extended data section as follows: "- n: number of trees. n is an integer >= 1, meaning that a specific set of species "sp", diameter "dbh" and height "h" can be present multiple times in a cell."**

The extreme values you identified are indeed outliers. They are a direct consequence of the uncertainties associated with the models. The realism of the stands associated with these extreme values is, of course, open to question. However, it seems difficult to define thresholds separating realistic from unrealistic stands as extreme values can be locally observed. We believe it is up to the users of the dataset to decide whether or not to use all the stands depending on their objectives. We drew the reader's attention to the outliers at the end of the validation section as follows: "Our procedure is not free of flaws and some outliers are present in the generated data (i.e. stands with extreme values of BA, Dg, tree height or density). These outliers are a direct consequence of the uncertainties associated with the models we used. The realism of the stands associated with these extreme values is open to question. However, separating realistic from unrealistic stands seems difficult as extreme values can be locally observed. It is therefore up to the users of the dataset to decide whether or not to consider these stands depending on their objectives."

Finally, it is always tricky to assess the relevance of extreme values as the field data
itself does not allow to estimate them accurately. We believe such an analysis could hardly help in evaluating our approach and we would rather not include it as part of the validation.

**Competing Interests:** No competing interests were disclosed.
itself does not allow to estimate them accurately. We believe such an analysis could hardly help in evaluating our approach and we would rather not include it as part of the validation.

**Competing Interests:** No competing interests were disclosed.