Capturing Temporal Dynamics in Large-Scale Canopy Tree Height Estimation

Thank you for your thorough review and for acknowledging our work. Let us address your remarks and questions in detail below.

However, the analysis section needs further improvement. [...]

We agree that this is an important aspect. Our primary goal was to develop an openly available model that others can use for their own analyses. To demonstrate the model's temporal capabilities, we analyzed potential deforestation events by tracking pixels for which the corresponding height decreased from above 8m to below 5m between years. The affected area increased from 7,729.1 km² in 2020-2021 to 15,942.5 km² in 2021-2022, also see our answer to reviewer 1141. We have updated our manuscript accordingly

However, a key weakness is that the study primarily applies existing models and algorithms [...]

[1] If the improved results are merely due to [...].

We would like to emphasize that our work follows the application-driven track. While pinpointing the exact source of improvement is challenging due to significant architectural differences, the results from Table 2 demonstrate three key improvements:

1. Utilizing 12 monthly images rather than a single composite image (-Composite- vs -Stack-) reduces MAE by 3.6%
2. Processing temporal data with 3D convolutions instead of 2D convolutions (2D-Stack- vs 3D-Stack-) further reduces MAE by 1.4%
3. Training on multiple years of data (2019-2022) versus a single year (2020) (-Year vs -MultiYear) reduces MAE by 6%

[2] The model structure mentions, "12 monthly Sentinel-2 images concatenated with an aggregated Sentinel-1 composite." [...]

The model receives 12 monthly images, each containing 16 channels - the Sentinel-2 bands for that month concatenated with the yearly median of Sentinel-1 data. So technically, the model overall receives 12*16=192 channels as input, although we do not concatenate these, but rather process them in the 3D-U-Net as (12, 16, 256, 256) image input. Regarding your second question, we do not explicitly filter clouds from the input images, but the model demonstrates robust performance even with partially cloudy inputs, maintaining high prediction accuracy. We utilize the L2A product, which includes atmospheric correction and other preprocessing steps.

[3] Why is the 'Sentinel-1 data aggregated into a single median-based image over the year' ? Is it because the data lacks clear periodicity?

While Sentinel-1 has similar temporal resolution to Sentinel-2, we opted to use composite radar data since it reduces noise and speckle artifacts through temporal aggregation, while keeping computational costs manageable. In contrast, we found using the full temporal sequence beneficial for optical Sentinel-2 data.

[4] The model's training configuration, duration, and computational resources should be provided in detail.

In Section 3.3 of the paper, we provide detailed information about the experimental setup, including the considered optimizer, learning rate, batch size, and more. Regarding duration and computational resources, we trained the model for a duration of roughly four days on two A100 GPUs. We will add detailed information about the computational resourced used to the revision. Thank you very much for this remark.

[5] If the reported MAE in this paper is 4.64m, what is the accuracy of GEDI? What is GEDI's temporal accuracy? Is the accuracy upper bound for time-series forest height inversion determined by GEDI's accuracy?

GEDI has a measurement resolution of 15cm but varying accuracy by surface type - overall MAE of 0.98m, with higher error in tree-covered areas (1.67m) versus grassland (0.79m) and cropland (0.57m) [PEL]. GEDI cannot reliably measure heights 0-3m and has reduced accuracy below 5m. Regarding temporal accuracy - GEDI rarely measures the same location twice, so it lacks a temporal dimension. For accuracy bounds - as GEDI provides our training labels, this likely sets an upper limit on our model's achievable accuracy. Is this what you meant?

[PEL] Pronk, M., Eleveld, M., & Ledoux, H. (2024). Assessing vertical accuracy and spatial coverage of ICESat-2 and GEDI spaceborne lidar for creating global terrain models. Remote Sensing, 16(13), 2259.

[6] Can you provide some failed cases and explain why the current method cannot handle them? Additionally, for low vegetation (below 5 meters), does this study face challenges in height estimation and time-series measurement?

We provide two failure cases: 1) GEDI measurements on slopes, where terrain changes are mistaken for vegetation height differences, and 2) checkerboard artifacts in harbor areas. We will expand this analysis in the appendix. Additionally, as noted earlier, both GEDI and our model struggle with accurate measurements of low vegetation.

We hope to have addressed all your concerns, please let us know if further clarification is needed. Thank you!

Thank you for your detailed review. In our manuscript, we proposed seasonal variation and geolocation shifts as possible factors influencing the model's performance. We agree that some of our statements have been to explicit from a environmental perspective (e.g., to make use of the geolocation shifts given in Sentinel-2 time series data to improve model performance or the use of time series data to capture seasonal patterns). We regret these unprecise statements. We have revised the manuscript to present these as hypotheses rather than definitive claims. Below, we include further experiments and analyses that explore these factors in more detail. We hope that our comments and additional experiments address your concerns.

Seasonal Pattern Analysis: How does the model specifically utilize seasonal information? Could you provide analysis showing:

a) Model attention/activation patterns across different seasons

The figure in https://ibb.co/1f8zh2gR shows activation patterns across months for different patches using Guided Attention [1]. We observe varying activation strengths across months and patches, suggesting the model processes temporal information differently by location. However, further research would be needed to confirm these hypotheses.

[1] Striving for simplicity: The all convolutional net.

b) Performance comparison between different types of forests that show different seasonal patterns

We evaluated our model separately on broadleaf and coniferous forests using the Copernicus Land Monitoring Service Forest Type Map (2018). These forest types show different seasonal patterns - broadleaf forests have distinct leaf-on/off periods while coniferous forests maintain constant canopy. The metrics below show our findings:

c) Ablation study isolating seasonal pattern benefits from general data volume benefits

Note that the number of training labels remains identical for all variants. To investigate the benefits of seasonal information independently from data volume, we conducted an ablation study comparing three models trained on different 4-month subsets:

• Winter (Nov-Feb)
• Summer (Jun-Sep)
• Mixed (Jan-Feb, Aug-Sep)

The results below show that using a mix of winter and summer months yields better validation performance.

Thank you for this comment. We will provide these additional findings in the updated version of our manuscript.

Geolocation Shift Benefits: Can you provide direct evidence that the model leverages Sentinel-2 geolocation shifts? Specifically:

a) Comparative analysis with models that cannot use shift information c) Quantification of improvement specifically attributable to shift utilization

Direct analysis of geolocation shifts is challenging as they are inherently embedded in raw satellite data. However, our ablation study (Table 3) shows superior performance using raw vs composite data. Note that although the amount of input data changes, we keep the number of training labels constant.

b) Examples showing improved edge detection or spatial detail

We provide a comparative figure (https://ibb.co/WpHqStqM) demonstrating enhanced edge detection across three model variants: 2D-Composite, 2D-Stack, and 3D-Stack. Additional examples of improved edge detection can be found in Figures 1, 4, and 8, as well as in our interactive Google Earth Engine application.

Growth Detection Validation: How reliable is the model at detecting forest growth? Please provide:

a) Validation using known growth areas

While our model's ability to detect forest growth is constrained by the 4-year observation period, we observe clear growth signals in forest plantation regions like the Le Landes forest plantation in France (search of "Garein" in our GEE app).

b) Comparison with ground measurements over time

We collected LiDAR data from the Vosges forest in France for 2020 and 2022 and analysed the measured vs predicted growth (https://ibb.co/YTbkXfpf). Our model is able to predict the growth of the trees, however the model uncertainty is high, which leads to "heavier tails" than what was measured. We have added these results to our manuscript.

c) Analysis of minimum detectable growth rate

Given our model's MAE of 4.76m and GEDI's uncertainty of 0.98m, we estimate that reliable growth detection requires changes of at least 5m over the observation period, varying by region and forest type.

We have also revised the paper about computational requirements, model performance and robustness across different forest types. Please let us know if you need any additional clarification.

Thank you for your thorough review and for acknowledging the contributions of our work. Let us address your concerns and questions one by one.

More references to remote sensing timeseries understanding literature such as [1] can be discussed.

Thank you for the suggestion. We have added additional references to our manuscript.

Although this paper mainly focuses on methodology and dataset [...]

We agree this is important. While our main focus is providing an openly available model for others to analyze, we did conduct some initial analysis: By tracking pixels where height decreased from >8m to <5m between years (indicating potential deforestation), we found affected areas of 9,747.9 km² (2019-2020), 7,729.1 km² (2020-2021), and 15,942.5 km² (2021-2022). We have added these findings to our manuscript.

Line 152: What’s the reason for including coastal aerosol (B01) and water vapour (B09) bands?

We decided to include all L2A (bottom-of-the-atmospohere) information, which includes all bands except for B10. We agree that B01 and B09 might not necessarily be relevant for the task at hand. To analyze their impact, we ran additional experiments to assess their importance. The following table reports the results on the validation part of our dataset (L1>15m refering to the L1 loss for all labels that exceed 15m):

As it can be seen, B01 and B09 only have little impact on the results and are, hence, candidates to be removed from the set of input channels. We will discuss these findings in the updated version of our manuscript.

Line 185: Could the authors clarify the sparsity of the labels? What percentage of the pixels have a canopy height label derived from GEDI?

GEDI measures roughly 4% of Earth's surface. For our 256x256 pixel training samples, only about 100 pixels (0.15%) have usable GEDI labels from the same year, due to noise and the need to match measurements temporally with satellite imagery.

Line 197: How are 2.56 km × 2.56 km patches created? Are they created by buffering GEDI point measurements?

We first make sure to only take image patches from our training areas to not have an overlap with the validation patches. We then randomly select a 2.56km x 2.56km area and load all GEDI measurements within that area. Since we only have the coordinates of the GEDI measurements, we assign them to the closest Sentinel-pixel. We hope that this clarifies your questions. Please let us know if that is not the case.

Line 199: Do “month images” refer to month composites or one image selected within a month? If the latter, could the authors clarify the reasons for not using monthly composites?

When creating the 12-months image stack, we select one of the images within each month, namely the one with the least amount of cloud cover. We decided not to use monthly composites, following Wolters et al., who showed that not using composites allows the model to learn finer details, possibly due to small geolocation shifts in the satellite images.

Line 207: Is smoothing applied before calculating validation metrics or just for the final mapping visualization?

Smoothing is only applied for the visualizations.

My only concern is the possible spatial autocorrelation in the validation dataset [...]

Thank you for raising this point. Indeed, there are spatial correlations between the different spatial areas. Note, however, that the learning scenario can be seen as a transductive learning setting, where one already has access to the (test) input images (but not the ground-truth labels). From that perspective, it is even valid to include the geolocation of an image stack as input. Note that the overall goal is to fill the remaining "gaps" that are not covered by GEDI groundtruth labels. We hope that this answers your question. Please let us know if there are still remaining concerns from your side.

Have the authors considered other architectures for remote sensing timeseries such as [1]?

We have considered other segmentation architectures and backbones, but we found that the 3D extension of the U-Net architecture to work surprisingly well. U-Nets are efficient and in our setting, yield good results. We will consider exploring more complex architectures in future work.

We hope to have clarified all concerns and questions, please let us know if further clarification is needed. Thanks!

Thank you for your positive and thorough review. Let us address your concerns and questions in detail.

How do you know that you are looking at trees (as opposed to buildings, for instance)?

That is a very good and true observation, indeed, given that GEDI measures the height of all objects, we cannot tell whether we are estimating a tree or a building: buildings are also measured at some height, therefore in cities one can see height prediction even in the absence of trees. This however is common practice in the remote sensing community. Canopy height maps are mainly used for forest monitoring and carbon stock estimation and both applications apply a forest mask before further use, which could be itself generated by a segmentation model; which is however a different problem we do not address here. While this is common practice, we will make this more clear in the paper to avoid confusion.

why is Sentinel-1 not mentioned in the abstract?

That was not intentional, thank you for pointing it out. We have updated our manuscript accordingly.

Please explain why sqrt{\text{MSE}} \neq \text{RMSE}?

We individually calculate the MSE and RMSE for each of the 1,500 validation patches and average them afterwards (weighted average based on the number of labels in each patch). This is why there can be a difference. Thank you for this remark. We will provide additional explainations in our manuscript.

Even though the authors do not provide their weights or datasets (they state that they will do so upon acceptance of their publication), the code provided appears reasonable.

Please note that the predictions can be viewed and downloaded from the Google Earth Engine (GEE) website linked in the paper. We are happy to share all the code to reproduce our results upon acceptance. Sharing terabytes of data in an anonymous way right now is challenging.

I would recommend training, validating, and testing in different years, in non-overlapping.

We tested this approach under the "-2020" setup, where we train only on 2020 data and evaluate on data from 2019, 2021 and 2022 (cf. e.g. Table 3 in the paper). We have made this more explicit in our manuscript, thank you.

A search yields the following articles related to the problem that are not mentioned: Satellite Image and Tree Canopy Height Analysis Using Machine Learning on Google Earth Engine with Carbon Stock Estimation

Thank you very much. We have added the reference to the background section of our manuscript.

what is the normalizing divisor in Table 1? Where does it come from?

We found it beneficial to rescale the data to a range that is more suitable for the model (i.e., we divide the input data by this divisor). To that end, we manually inspected the data to find out value ranges that contain valuable data, in particular because standard normalization did not work due to problem with atmospheric distortions and cloud cover.

provide R2 in Table 2. (no need to include both MSE and RMSE)

Thank you for the suggestion, we have added R2 to Table 2.

seperate, Appendix A

Fixed, thank you!

What are the sources of error when estimating tall trees?

Although we do not know the exact reason, we suspect that it has to do with the following two reasons:

1. Due to the natural distribution of trees in Europe, tall trees are less common, creating a skew in the label distribution.
2. Tall trees naturally have a higher canopy density (e.g. more leaves, branches, etc.), which leads to a higher fraction of LIDAR photons/measurements not penetrating the canopy. In that case, the photons do not reach the ground and we do not have a usable measure of the tree height (this filtering is already applied by GEDI, not by us).

Why not use stratified evaluation based on different forest types?

Thank you for this suggestion. While our paper focused on overall performance metrics, we have now included a detailed analysis comparing broadleaf and coniferous forests (see our response to reviewer qZ94). We welcome your input on additional forest categories that would be valuable to evaluate.

What post-processing steps were used to ensure temporal consistency?

We use a quadratic-spline approach to smooth the predictions over time, but only for visualization purposes.

We hope to have addressed all your remarks. Thank you again. If you have any further questions or concerns, please let us know.