Responding to Large-Scale Forest Damage in an Alpine Environment with Remote Sensing, Machine Learning, and Web-GIS

: This paper reports a semi-automated workﬂow for detection and quantiﬁcation of forest damage from windthrow in an Alpine region, in particular from the Vaia storm in October 2018. A web-GIS platform allows to select the damaged area by drawing polygons; several vegetation indices (VIs) are automatically calculated using remote sensing data (Sentinel-2A) and tested to identify the more suitable ones for quantifying forest damage using cross-validation with ground-truth data. Results show that the mean value of NDVI and NDMI decreased in the damaged areas, and have a strong negative correlation with severity. RGI has an opposite behavior in contrast with NDVI and NDMI, as it highlights the red component of the land surface. In all cases, variance of the VI increases after the event between 0.03 and 0.15. Understorey not damaged from the windthrow, if consisting of 40% or more of the total cover in the area, undermines signiﬁcantly the sensibility of the VIs to detecting and predicting severity. Using aggregational statistics (average and standard deviation) of VIs over polygons as input to a machine learning algorithm, i.e., Random Forest, results in severity prediction with regression reaching a root mean square error (RMSE) of 9.96, on a severity scale of 0–100, using an ensemble of area averages and standard deviations of NDVI, NDMI, and RGI indices. The results show that combining more than one VI can signiﬁcantly improve the estimation of severity, and web-GIS tools can support decisions with selected VIs. The reported results prove that Sentinel-2 imagery can be deployed and analysed via web-tools to estimate forest damage severity and that VIs can be used via machine learning for predicting severity of damage, with careful evaluation of the effect of understorey in each situation. topography [20] and are applied to analyze seasonal and annual trends [21]. It must be noted that VIs depend strongly also on the phenology, i.e., the growing stage, not only on disturbance and stress factors induced by climate change or pathogen attacks. NDVI has a non-linear response, in sparsely vegetated areas where the background reﬂectance of the red spectral band increases [22–24]. Over densely vegetated areas, NDVI enhances the low ratio values while compressing higher ratio values. The enhanced vegetation index (EVI) was created to correct the background effect using a soil adjustment factor, L, and correcting the red band with blue band coefﬁcients, C1 and C2 [23]. In the windthrow-damaged area, the forestry patch is sparse and many trunks lie on the ground. Hence, EVI has been tested to reduce the noise effect of the canopy background. A damaged forest, with weakened, stressed, and fallen trees as in Figure 1, can be the optimal habitat for forest pests like the European spruce bark beetle ( Ips typhographus ) [25].


Introduction
High severity natural disturbances can significantly alter the structure and morphology of forest landscape [1][2][3]. The damage assessment, generally evaluated using severity classes, has been usually carried out through field surveys.
Wind is the main disturbance agent affecting European forests, being responsible for more than 50% of all damage by volume [4]. The damage from fire and windstorms have been increasing over the past centuries, and they are likely to continue to increase in the next future [5,6]. According to Gregow et al. [7], windstorms that affected the Central and Northern European forests in the last decade were stronger than the ones that occurred before 1990. Moreover, in the period between 1980-2010, the impact on forest areas, considering the same amount of biomass, was 3-4 times higher than that registered in the period between 1950-1980. Part of this increase might be due to climate change [8], but most of it is high likely due to factors related to the current stand conditions of the European forests [9]. The stages of the pest attack are the green-attack, red-attack, and grey-attack, where stressed vegetation loses foliage, increases the visible red spectral component, and decreases the near-infrared response [26]. In the green-attack, the foliage looks unchanged. Thus, it is difficult to identify, and the moisture decreases in the sapwood [27]. The redstage happens the next year and the foliage turns yellow and red. Finally, the tree loses needles, and turns to grey. The grey stage has been studied with NDVI [28], but several studies have focused on the red stage applying the RGI index [29][30][31].
The relative greenness index (RGI) is an index of anthocyanin content [32], and it takes into consideration the response of the red band over the green band. Study results on this index bring evidence of a high correlation coefficient, from 0.7 to 0.9, between pixels covered with red-stage pest attack and the tree crown using a very high spatial resolution imagery e.g., ~1 m resolution or less. The epidemic red stage was studied in Canada The stages of the pest attack are the green-attack, red-attack, and grey-attack, where stressed vegetation loses foliage, increases the visible red spectral component, and decreases the near-infrared response [26]. In the green-attack, the foliage looks unchanged. Thus, it is difficult to identify, and the moisture decreases in the sapwood [27]. The red-stage happens the next year and the foliage turns yellow and red. Finally, the tree loses needles, and turns to grey. The grey stage has been studied with NDVI [28], but several studies have focused on the red stage applying the RGI index [29][30][31].
The relative greenness index (RGI) is an index of anthocyanin content [32], and it takes into consideration the response of the red band over the green band. Study results on this index bring evidence of a high correlation coefficient, from 0.7 to 0.9, between pixels covered with red-stage pest attack and the tree crown using a very high spatial resolution imagery e.g.,~1 m resolution or less. The epidemic red stage was studied in Canada using Landsat-TM imagery and applying the enhanced wetness difference index (EWDI) in a mixed forest [33]. The EWDI is calculated by subtracting the tasseled cap wetness index [34] of two images in different years. Again, the red-stage attack was studied by [35] using the normalized difference moisture index (NDMI). Wilson and Sader [36] found higher accuracies of NDMI over NDVI in change detection of forest harvest. Biotic and abiotic disturbances were studied considering the red-edge chlorophyll indices (CI) and NDVI for estimating the temporal change in leaf chlorosis and defoliators in [37].
Previous scientific literature supports further investigation to understand the accuracy of modelling damage from windthrow using vegetation indices, to identify those that provide optimal results over mountain regions characterized by an Alpine environment. The specific goal of this paper is to report on results of testing NDVI, EVI, RGI, EWDI, NDMI, and CI indices through machine learning and to identify which indices are more suitable for quantifying damage in areas in the Agordino forest hit by the Vaia storm, considering also the disturbance and recovery stage. The results have been integrated into a web-GIS application to support the public administration and stakeholders involved in forestry management.

Study Area
The study area is located in the upper and middle Cordevole Valley, named Agordino. Here, the forests were severely damaged by the Vaia storm, which affected almost 4000 ha (6% of the area) with more than 600,000 m 3 of windthrown timber. Consequently, 22 testing areas ( Figure 2), representative of a gradient of severity, aspect, and slope, were inspected between June and September 2019. The size of the testing areas ranges between 1 ha and 61.58 ha. using the normalized difference moisture index (NDMI). Wilson and Sader [36] found higher accuracies of NDMI over NDVI in change detection of forest harvest. Biotic and abiotic disturbances were studied considering the red-edge chlorophyll indices (CI) and NDVI for estimating the temporal change in leaf chlorosis and defoliators in [37].
Previous scientific literature supports further investigation to understand the accuracy of modelling damage from windthrow using vegetation indices, to identify those that provide optimal results over mountain regions characterized by an Alpine environment. The specific goal of this paper is to report on results of testing NDVI, EVI, RGI, EWDI, NDMI, and CI indices through machine learning and to identify which indices are more suitable for quantifying damage in areas in the Agordino forest hit by the Vaia storm, considering also the disturbance and recovery stage. The results have been integrated into a web-GIS application to support the public administration and stakeholders involved in forestry management.

Study Area
The study area is located in the upper and middle Cordevole Valley, named Agordino. Here, the forests were severely damaged by the Vaia storm, which affected almost 4000 ha (6% of the area) with more than 600,000 m 3 of windthrown timber. Consequently, 22 testing areas (Figure 2), representative of a gradient of severity, aspect, and slope, were inspected between June and September 2019. The size of the testing areas ranges between 1 ha and 61.58 ha.

Field Survey
The severity was assessed by determining the canopy cover loss by visual interpre-

Field Survey
The severity was assessed by determining the canopy cover loss by visual interpretation and validation by experts. First, interpretation was done using high-resolution orthoimagery from aerial surveys (acquired in 2019), determining the boundaries of the area considered damaged and providing a value of severity (canopy loss) on a scale from 0 to 100. The second step consisted of a visual inspection directly on site, validating or changing the severity class assigned in the first step.
The field survey consisted of visual inspection of the 22 plots reporting (i) the severity in the percentage of fallen trees; a description of the damage with percentages of (ii) residual canopy cover, (iii) understorey, (iv) overall health index, and (v) crown status of standing trees. Visual inspection means that expert personnel physically went in the area and validated the initial severity assessment for each plot. The survey is necessary to assess the forest condition, and to correlate the ground truth with the values of VIs that were derived from Sentinel-2 imagery. The fieldwork schedule was coordinated with the Sentinel-2 acquisition dates, also considering weather conditions (cloud-free days). Plot names and values of severity, residual canopy cover, and understorey are reported in Table 1 in the results.

Satellite Imagery Processing: Vegetation Indices
Level-2A Sentinel-2 images over the summer season in 2018 and 2019 have not been clipped to obtain six cloud-free images. Sentinel-2 Level-2A provides the bottom of surface reflectance, i.e., correction of images for atmospheric effects, through the sen2cor algorithm [38]. Two images were sensed in August and September 2018 (before the event) and four were sensed between June and October 2019 (after the event). Using these images, we calculated the following indices NDVI, EVI1, EVI2, RGI, NDMI, CI, and EWDI.
Below are the specific formulas, where "nir" is the near infrared band, specifically band 8 having 10 m spatial resolution for Sentinel-2; and "red", which is band 4 in the Sentinel-2 specifications. The NDVI formula is: The EVI index was calculated in two versions called EVI1 and EVI2. The EVI1 formula is: The EVI2 is: RGI is a simple ratio of the red and green bands, as stated in (4) RGI = red green (4) NDMI measures the response to the NIR spectral band of the leaf structure in contrast with the mid-infrared (MIR) absorbed by the water in the vegetation foliage [33,34], as stated in (5): The CI red-edge chlorophyll index is based on the ratio R750/R710, as defined by [35] as stated in (6), and the corresponding Sentinel-2 bands are band 6 and band 5 (respectively, having 740 nm and 705 nm). The stressed vegetation has a decrease in chlorophyll content, which absorbs the blue and red light at 430-660 nm. Consequently, the visible red spectral component increases, whereas the near-infrared response decreases: CI = red 750 red 705 (6) EWDI is the difference between the most recent tasselled cap image and older image. The tasselled cap formula is stated in (7): T. cap wetnes = 0.1509 * blue + 0.1973 * red + 0.3279 * red + 0.3406 * nir8 − 0.7112 * nir11 − 0.4572 * nir12 The indices were extracted over 19 damaged areas and 3 control zones that were not damaged. The severity ranges from 0 (LCL_NW1, LCL_NW2) to 100% (RP_01), and it has been evaluated with field survey. Average and standard deviation of the index values inside the polygon were aggregated for each date. Spearman's rank correlation between the average values and severity has been calculated.

Satellite Imagery Processing: Machine Learning
Vegetation indices can be used for predicting the severity of damage using machine learning. Assessment of accuracy was carried out for the following methods: Support Vector Machine (SVM), Random Forest (RF), and K-nearest neighbours (KNN).
The SVM split the dataset into two groups using a separating hyperplane. The class is assigned to all class membership using a Kernel function [39]. The RF creates trees and applies a predictor to all branches. The tree is built using a Bootstrap statistical technique for data aggregation [40]. The KNN uses a distance metric, calculating the k nearest neighbors of the sample. The rminer library for R CRAN was used for classification and accuracy metrics in this task [41], which in turn uses the kknn library. Default values were used, i.e., optimal kernel and the number of neighbours k = 7.
The dataset consists of the VIs from the 6 images over the 22 areas, for a total of 132 observations with VIs as predictors and severity as the target variable to predict. All severity values before the event were set to zero. The dataset was split in 60% as a training set and 40% as a testing set. The training consisted of using 10-k fold cross-validation, a Remote Sens. 2021, 13, 1541 7 of 23 robust methodology to prevent overfitting issues whereas the testing set was used for the final prediction. The metrics to assess the accuracy of the regression are the following: mean absolute error (MAE), root mean square error (RMSE), correlation (COR), and adjusted R 2 (adj. R2), whereas the metric used to evaluate the classification is Kappa metric (K). The importance of the VIs is calculated over sensitivity analysis, where the input sensitivity is calculated for each input VIs using a One-dimensional SA (1D-SA) with seven levels described in [42].

Web-GIS Interface
The structure of the web-GIS is composed of the back end, the analysis module, and the front end, as described in Figure 3. The satellite imagery and metadata are available through the Sentinel API Hub, so updated Level-2 images over the region of interest are automatically downloaded and stored as PostGIS elements in the PostgreSQL database. The data is kept in the original jp2 file format, and PostGIS stores meta-information (out_db option). The analysis module is based on R CRAN language, which connects the database with the front end. Therefore, the query can be launched by the user through the front-end interface that has been developed using Shiny and Leaflet packages. The user sets the parameters, which are the vegetation index, the areas drawn on the map, and the temporal window. In contrast, for the querying of the database, the threshold for the cloud coverage for each image is automatically set below 20%, and the cloud masking is automatically applied. The results are a colour-coded index map and a time-series plot (see Appendix A).
set and 40% as a testing set. The training consisted of using 10-k fold cross-validation robust methodology to prevent overfitting issues whereas the testing set was used for t final prediction. The metrics to assess the accuracy of the regression are the followin mean absolute error (MAE), root mean square error (RMSE), correlation (COR), and a justed R 2 (adj. R2), whereas the metric used to evaluate the classification is Kappa met (K). The importance of the VIs is calculated over sensitivity analysis, where the input se sitivity is calculated for each input VIs using a One-dimensional SA (1D-SA) with sev levels described in [42].

Web-GIS Interface
The structure of the web-GIS is composed of the back end, the analysis module, a the front end, as described in Figure 3. The satellite imagery and metadata are availa through the Sentinel API Hub, so updated Level-2 images over the region of interest a automatically downloaded and stored as PostGIS elements in the PostgreSQL databa The data is kept in the original jp2 file format, and PostGIS stores meta-informati (out_db option). The analysis module is based on R CRAN language, which connects t database with the front end. Therefore, the query can be launched by the user through t front-end interface that has been developed using Shiny and Leaflet packages. The u sets the parameters, which are the vegetation index, the areas drawn on the map, and t temporal window. In contrast, for the querying of the database, the threshold for the clo coverage for each image is automatically set below 20%, and the cloud masking is au matically applied. The results are a colour-coded index map and a time-series plot (s Appendix A).

Results and Discussion
This section is composed of three parts: (i) the temporal profile of each VI, (ii) t weight of each VI in a prediction model using machine learning, and (iii) the integrati in the web-GIS portal.

Temporal Analysis
The results of the field survey and the analysis are reported in Tables 1 and 2. Ta 1 shows the names of the surveyed areas, with values of estimated severity, percent post-event canopy cover and understorey with respective average and standard deviati of VI value in post-event imagery (September 2019). The field survey highlighted th some areas have a similar severity of damage, but different understorey, such as VOA_

Results and Discussion
This section is composed of three parts: (i) the temporal profile of each VI, (ii) the weight of each VI in a prediction model using machine learning, and (iii) the integration in the web-GIS portal.

Temporal Analysis
The results of the field survey and the analysis are reported in Tables 1 and 2. Table 1 shows the names of the surveyed areas, with values of estimated severity, percent of post-event canopy cover and understorey with respective average and standard deviation of VI value in post-event imagery (September 2019). The field survey highlighted that Remote Sens. 2021, 13, 1541 8 of 23 some areas have a similar severity of damage, but different understorey, such as VOA_02, LCL_04, TA_02 and TA_01. The severity in VOA_02 is 80% and understorey cover is 70, whereas the severity in LCL_04 is 80% and understorey cover is 35%. TA_02 and TA_01 have 60% of severity and understorey cover, respectively, of 40% and 10%. Table 1 also shows the values of the VIs in September 2019. The photosynthesis rate in the Alpine regions is higher in September, i.e., tree growth. Understorey affects NDVI value in the areas with similar severity. VOA_02 and TA_02 are examples of high NDVI (~0.7) and high understorey cover in contrast with LCL_04 and TA_01. The last three rows represent control areas nearby, which were healthy stands not affected by wind damage.
The NDVI temporal series for each plot area are reported in Figure 4. The control areas with no damage are named LCL_NW1, LCL_NW2, RP_NW (red outline), and clearly have stable NDVI values, as expected. The NDVI values in the damaged areas decrease clearly and the standard deviation increase e.g., the areas RP_03, RP_04. It is worth noting that in some damaged areas, such as the VOA_02, the NDVI values decrease to a lesser extent than other damaged areas, and the post-event values seem stable over 2019. Again, TA_01 and TA_02 both have 60% severity, but TA_01 has 10% understorey and NDVI of 0.44, TA_02 has 70% understorey and NDVI of 0.69. LCL_04 and CLS_03 are two forestry harvest area where the trunks have been removed. The soil is covered by woodchips mixed with herbaceous recovery; NDVI values in these two cases range from 0.33 to 0.59.
Detailed information about the effect of residual shrub cover and vegetation regeneration over NDVI at pixel scale can be observed in Figure 5. The figure shows red/nir spectral space with values of each pixel in the imagery. The points show the significant effect of post-event residual (shrub) and regenerated vegetation with the scatterplot and the NDVI baselines in September 2019. On one hand, the points of RA_02 and VOA_02 are compressed along the Y (red) axis. Both have the same shrub cover of 70%, but RA_02 site, with 30% severity, results in 0.81 value of NDVI, whereas VOA_02 has 80% severity, resulting in lower NDVI value of 0.74. In both cases, no information can be extracted about the recovery, because the NDVI value reflects vegetation that can be either from existing surviving shrub (understorey) or from recovery.
In contrast, in TA_02, the understorey is 70%, and the pixels' NDVI values range between the baselines of 0.6 and 0.8. Each pixel in area TA_02 positions itself in the red/nir scatterplot, depending on the vegetation/damage fraction mixture in the pixel. Scenarios range from having a significant fraction of vegetation (0.8 NDVI value), likely due to heavy shrub or even trees resisting the storm, which cannot be differentiated in the scatterplot from regeneration. The harvesting sites show similar behavior-CLS_03 has 90% severity, but 45% of the area is recovering with herbaceous vegetation and shrubs. The points in the scatterplot in Figure 5 range between the NDVI baselines, representing 0.4 and 0.8 values. LCL_04 has 80% severity and 35% understorey, but the points in the scatterplot cover an area below the NDVI baseline 0.4. Consequently, the background reflectance of the soil can contribute to defining severity significantly when the sparse vegetation coverage (remaining shrubs or recovery) is less than 40%. Above 40% residual vegetation, the NDVI does not allow to easily distinguish a damaged area. It is worth therefore stating that it is important to record the initial post-event NDVI value, thus defining the initial per-pixel NDVI "residual" from shrubs or other resisting vegetation. Finally, a strong negative Spearman's rank correlation between NDVI and severity was found for 2019, whereas for 2018, the correlation was not significant, as can be expected.   RGI is a simple ratio of the red and green band. In the pre-event and into the co zones, the values range between 0.65 and 0.7 (see Figure 4), which means the green r tance is greater than the red reflectance. In the harvesting sites, the ratio becomes g than 1, so the reflectance of red band exceeds the reflectance of the green band. I Agordino windthrown forests, the trunks are lying on the ground since the major trees were uprooted. In this position, the canopy does not completely cover the trun the branches; therefore, the bark and soil are directly exposed and visible from th agery. Consequently, the red spectral component increases. A strong correlation RGI is a simple ratio of the red and green band. In the pre-event and into the control zones, the values range between 0.65 and 0.7 (see Figure 4), which means the green reflectance is greater than the red reflectance. In the harvesting sites, the ratio becomes greater than 1, so the reflectance of red band exceeds the reflectance of the green band. In the Agordino windthrown forests, the trunks are lying on the ground since the majority of trees were uprooted. In this position, the canopy does not completely cover the trunk and the branches; therefore, the bark and soil are directly exposed and visible from the imagery. Consequently, the red spectral component increases. A strong correlation of R2 greater than 0.7 between RGI and the severity has been found (see Table 2).
NDMI has a similar trend to NDVI response, and it has a strong negative correlation with the severity. NDMI values range from 0.07 to 0.7 (see Figure 4). Healthy vegetation ranges between 0.6 and 0.8, whereas damaged areas reach a minimum of 0.2 such as RP_01. In VOA_02 that had 80% of the area damaged and the remaining understorey of 70%., average NDMI value is 0.48. Likewise, for the site RA_02, which has exactly the same understorey of 70%, the NDMI value is a bit greater, 0.55, because the severity is less than half (30%). In areas TA_01 and TA_02, both with severity of 60%, the NDMI values range between 0.07 and 0.37 due to the different remaining shrub coverage that is 10% and 40%, respectively, so the index helps to identify the fraction of bare soil and shrubs.
EVI1 and EVI2 have been tested to study the effect of the canopy fraction. Both indices have a similar trend, and they range from 0.19 to 0.52. In the pre-event, the values are quite constant around 0.5, and in the after-event, the value of the indices reach a peak at 0.6. Despite the Spearman's correlation being negative moderate, the difference pre-post storm is not very clear, especially in the plot. These indices will not be integrated into the web application.
EWDI is calculated as the difference between the tasseled cap wetness in two sequential temporal dates. When the difference was calculated in the same year e.g., in September and August 2018, obviously no correlation with forestry damage was found. In contrast, on subtracting September 2018 to September 2019 and August 2018 to August 2019, a strong negative correlation between the index and the severity was found.
The CI index studies the response of the chlorophyll content that decreases in damaged areas. The content of chlorophyll in healthy vegetation has a value around 3 [37]. In areas with sparse green cover such as TOA_01 and TOA_02, CI ranges between 1.45 and 2.10, and the absolute minimum value from all our areas is 1.34 in LCL_04. Therefore, CI shows a strong negative correlation with the percentage of damage. The temporal variation of the CI and NDVI is reported in the scatterplots in Figure 6. In the control areas, there is no variation in chlorophyll content and NDVI response, so the scatterplot contains sparse points that do not reflect any change. In contrast, when an area is damaged, some pixels are aligned with the healthy baseline and other pixels moves from the black baseline (healthy) to the red baseline (stress), reflecting coherent results with what is reported in [37].
In conclusion, all indices except EVI1 and EVI2 can identify quite well the damaged areas. This is in line with results that report that EVIs give significant information on green biomass in non-degraded forests [43]. The visual plot of NDVI, RGI, and NDMI are easier to read and understand, and will be used for the web application. The VIs have been used for other disturbances such as bark-beetle infestation [25,33], which is a different case, but present a similar drop in photosynthetic material.

Machine Learning Analysis
Regression with machine learning can be used to predict severity. The model parameters were tuned automatically using the Rminer package and a number of search equal to 10. Results were assessed using k-fold cross-validation, with k = 10. The regression error characteristic curve (REC), which describes the accuracy of the models, shows very similar results comparing RF, KNN, and SVM (Figure 7) with the former two having slightly better results. The differences between the first two algorithms are the importance of input variables. The most important input in both cases is NDMI, but RF takes advantage of EVI2 information to a greater degree. In contrast, KNN considers RGI as an important variable (Figure 8). Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 23 Figure 6. The plot compares NDVI/NDVI+1 and CI/CI+1. In the area with medium and high severity, the points move from the healthy black baseline to the right. The areas with low or null severity have a compact set of points. Red and black baselines are from [37]. ter results. The differences between the first two algorithms are the impo variables. The most important input in both cases is NDMI, but RF take EVI2 information to a greater degree. In contrast, KNN considers RGI a variable (Figure 8).   eters were tuned automatically using the Rminer package and a number of search e to 10. Results were assessed using k-fold cross-validation, with k = 10. The regression e characteristic curve (REC), which describes the accuracy of the models, shows very sim results comparing RF, KNN, and SVM (Figure 7) with the former two having slightly ter results. The differences between the first two algorithms are the importance of i variables. The most important input in both cases is NDMI, but RF takes advantag EVI2 information to a greater degree. In contrast, KNN considers RGI as an impo variable ( Figure 8).   The severity was then predicted using KNN, RF, SVM algorithms, and a combination of indices over the testing (Figure 9). Results for the RF algorithm, calculated using the Rminer package, are reported in Table 3. NDMI is among the better predictors in our study for severity. It must be noted that NDMI is used in the literature for detecting water stress that can be caused by abiotic or biotic factors, e.g., water availability and pathogen attacks, respectively, and damaged trees reflect a similar spectral signature. It is therefore important to use NDMI when the causality can be identified specifically.
emote Sens. 2021, 13, x FOR PEER REVIEW 17 of 2 The severity was then predicted using KNN, RF, SVM algorithms, and a combination of indices over the testing (Figure 9). Results for the RF algorithm, calculated using the Rminer package, are reported in Table 3. NDMI is among the better predictors in ou study for severity. It must be noted that NDMI is used in the literature for detecting wate stress that can be caused by abiotic or biotic factors, e.g., water availability and pathogen attacks, respectively, and damaged trees reflect a similar spectral signature. It is therefore important to use NDMI when the causality can be identified specifically. The index performing less well is EVI2 with a minimum of 0.38 of adjusted R 2 . The NDVI has a average adjusted R 2 value of 0.66, but the value increases using NDVI togethe with the standard deviation of NDVI; in this case the adjusted R 2 reach 0.84.

Conclusions
The goal of this investigation is to test the potential of VIs that have proven useful in previous applications, i.e., NDVI, EVI, RGI, EWDI, NDMI, CI vs. NDVI for their potentia to predict severity of forest areas damaged by windthrow in an alpine environment. The practical application of these results is to find the most suitable VIs for integrating thei use in a web-GIS application, which automates the typical basic tasks in remote sensing Remote sensing analysis requires a piece of specific knowledge and the application aims to simplify the remote sensing workflow as a support to public administration and stake holder, which are not experts in remote sensing. The architecture is based on the Post greSQL database and cloud computing, and allows easy and fast on-demand analysis. The analysis focuses on windthrow damage and forestry recovery. Furthermore, the indice  The index performing less well is EVI2 with a minimum of 0.38 of adjusted R 2 . The NDVI has a average adjusted R 2 value of 0.66, but the value increases using NDVI together with the standard deviation of NDVI; in this case the adjusted R 2 reach 0.84.

Conclusions
The goal of this investigation is to test the potential of VIs that have proven useful in previous applications, i.e., NDVI, EVI, RGI, EWDI, NDMI, CI vs. NDVI for their potential to predict severity of forest areas damaged by windthrow in an alpine environment. The practical application of these results is to find the most suitable VIs for integrating their use in a web-GIS application, which automates the typical basic tasks in remote sensing. Remote sensing analysis requires a piece of specific knowledge and the application aims to simplify the remote sensing workflow as a support to public administration and stakeholder, which are not experts in remote sensing. The architecture is based on the PostgreSQL database and cloud computing, and allows easy and fast on-demand analysis. The analysis focuses on windthrow damage and forestry recovery. Furthermore, the indices take into consideration disturbances that can occur in the future, such as parasite attack. Studying the damaged areas using Sentinel-2 imagery, the results show a strong correlation of NDVI, RGI, and NDMI. NDVI and NDMI decrease in the damaged areas, and they have a strong negative correlation. RGI has a trend in contrast with the NDVI, and it highlights the red component of the image, such as bark, branches, and soil. It is a simple ratio, easy to calculate using a red green blue (RGB) camera, which is the base equipment of commercial drones. If detailed inspection is necessary, UAVs imagery can integrate satellite analysis [35]. NDMI has a behaviour similar to NDVI, but it is specific for water stress. Finally, the plots in Figure 4 summarize the temporal behavior of the VIs at the sites. The control areas have a stable trend in 2018 and 2019 and the damaged areas show a variation in the index with an increase in standard deviation. Hence, damaged areas can be recognized, but green coverage greater than 40% influences the indices' response. Above this percentage, the apparent forest recovery can be due to shrubs, herbs, or regenerating trees, so field survey is required. KNN and RF assign different importance to the VI. KNN uses NDMI, RGI, and NVI whereas RF focus on NDVI, EVI2, and RGI. RF gives less importance to NDVI than KNN. However, severity prediction using RF produced good results using NDMI or NDVI and the NDVI standard deviation. In this study, the surveyed areas describe the various case histories in terms of exposure, slope, and severity in the Agordino area. However, further analysis is necessary for integrating Sentinel-1 and Sentinel-2, or applying the model to other Alpine or Pre-Alpine areas affected by the Vaia storm. After two years from the storm, more investigation can focus on disturbances due to parasites near damaged areas, such as to study forest health by applying machine learning.
Future work is expected, both related to the analytical part and to the web-GIS. Regarding the analytical part, tests will be done adding Sentinel-1 information, which might provide added and non-collinear information with respect to the optical information from Sentinel-2. In addition, the regeneration will be monitored with the VIs that have been assessed, to see the sensibility of such information to recovery of the vegetation in the area. The web-GIS will continue development, to provide automated visualization and analysis of Sentinel-2 data online. In particular, it is planned to add specific per-pixel cloud probability that can be calculated with more accurate methods (e.g., s2cloudless [44]) in order to have a finer definition of reliable pixel values. Users will also be able to define their own threshold for cloud probability.
management. In the current stage, the software integrates three indices, determined using previous investigation on correlation and variable importance using machine learning regression. The indices used are NDVI, RGI, and NDMI. Figure A1 shows an example of the interface. On the left, a toolbar allows to select the layer displayed, to activate a measuring tool and to draw a polygon where the VI analysis will be done. Then, VI and the temporal domain are chosen in the setting panel. All available Sentinel-2 images, which spatially overlap the testing area, are clipped and masked. Masking consists of removing pixels which belong to "cloud (high probability)", "cloud (medium probability)", and "cloud-shadow" in the "Scene Classification" (SCL) product of the Level 2C Sentinel-2 product. Masking is an important step, which removes gross errors in the time-series data that are created by the user in real-time and might be used for further processing.
A downloadable zonal statistic plot with average and standard deviation of the selected VI is produced. The color-scaled index map is automatically displayed, and the time bar allows to change the VI images that are created on-the-fly and projected in the map. The following items provide an overview of the service provided by the web platform to end-users: • draws, over the view area, a styled raster of the chosen VI, automatically clipping and masking; • creates a downloadable zonal statistic plot reporting average and standard deviation values of pixels inside the selected area for all images in the temporal scale; • provides a downloadable report with raw values with id, timestamp, average, and standard deviation calculated from all images over the user-selected area (the data used for creating the plot described in the previous point).
Available Sentinel-2 Level-2C images are automatically downloaded through a script that checks for available new images each day. This is done only over one specific tile, (TQS32) as the project had a specific study area to be investigated.

FOR PEER REVIEW
19 of 23 previous investigation on correlation and variable importance using machine learning regression. The indices used are NDVI, RGI, and NDMI. Figure A1 shows an example of the interface. On the left, a toolbar allows to select the layer displayed, to activate a measuring tool and to draw a polygon where the VI analysis will be done. Then, VI and the temporal domain are chosen in the setting panel. All available Sentinel-2 images, which spatially overlap the testing area, are clipped and masked. Masking consists of removing pixels which belong to "cloud (high probability)", "cloud (medium probability)", and "cloudshadow" in the "Scene Classification" (SCL) product of the Level 2C Sentinel-2 product.
Masking is an important step, which removes gross errors in the time-series data that are created by the user in real-time and might be used for further processing. A downloadable zonal statistic plot with average and standard deviation of the selected VI is produced. The color-scaled index map is automatically displayed, and the time bar allows to change the VI images that are created on-the-fly and projected in the map. The following items provide an overview of the service provided by the web platform to end-users: • draws, over the view area, a styled raster of the chosen VI, automatically clipping and masking; • creates a downloadable zonal statistic plot reporting average and standard deviation values of pixels inside the selected area for all images in the temporal scale; • provides a downloadable report with raw values with id, timestamp, average, and standard deviation calculated from all images over the user-selected area (the data used for creating the plot described in the previous point).
Available Sentinel-2 Level-2C images are automatically downloaded through a script that checks for available new images each day. This is done only over one specific tile, (TQS32) as the project had a specific study area to be investigated.