## 3 Methods
### Characterization of Forest Ecosystems
#### 3.1.3 Forest Density Mapping  
Sample drawing and interpretation  
The first step consists in establishing a reference value of Forest Density for the sample units, based on the interpretation of VHSR imagery (Ikonos, Pleiades). The prerequisite is the co-registration of VHSR data to HR data, i.e., Ikonos 2007 and Landsat 5 nearest epoch 2005, Pleiades 2014 and Landsat 8 nearest epoch 2015, Pleiades 2018 and Sentinel-2 nearest epoch 2019. The sampling design is based on a stratified random sampling of 100 squared Primary Samples Units (PSU) of 100 \(\times\) 100 m that intersect the forest stratum. Each PSU is then subdivided into 100 Secondary Sample Unit (SSU) of 10 \(\times\) 10 m spacing. Each SSU is then visually interpreted as forest or non-forest. Zonal statistics are then computed for each PSU of 100 m\({}^{2}\) as the sum of SSUs interpreted as forest from which is then derived the percentage of forest density (see Figure 5 below). Figure 5: Visual interpretation of SSU points inside a PSU for computing the percentage of Forest Density reference value. 2. Computing vegetation indices The second step consists in computing several vegetation indexes within each PSU over the HR times-series, in order to establish a correlation between those indexes and the Forest Density value obtained through the interpretation of samples. Common indexes for the assessment of Forest and the characterization of vegetation and leaf properties are the Normalized Difference Vegetation Index (NDVI), Brightness Index (BI), Normalized Burnt Ratio (NBR), Green Normalized Difference Vegetation Index (GNDVI), Modified Normalized Difference Vegetation Index (GNDVI), Modified Normalized Difference Water Index (MNDWI), as well as the mean value of the spectral bands available, such as Blue, Green, Red, Near Infrared (NIR) and Short Wave Infrared (SWIR) [15]. Normalized Difference Vegetation Index (NDVI) \[\text{NDVI}=\frac{\text{NIR}-\text{RED}}{\text{NIR}+\text{RED}}\] (1) Brightness Index (BI) \[\text{BI}=\frac{(\text{RED}\times\text{RED})/(\text{GREEN}\times\text{GREEN})}{ 2}\] (2) Normalized Burn Ratio (NBR) \[\text{NBR}=\frac{\text{NIR}-\text{SWIR}}{\text{NIR}+\text{SWIR}}\] (3) Green Normalized Difference Vegetation Index (GNDVI) \[\text{GNDVI}=\frac{\text{NIR}-\text{GREEN}}{\text{NIR}+\text{GREEN}}\] (4) Modified Normalized Difference Water Index (MNDWI) \[\text{MNDWI}=\frac{\text{RED}-\text{SWIR}}{\text{RED}+\text{SWIR}}\] (5) (iii) Multiple linear regression analysis The third step lies in defining the model that will predict the Forest Density value, through an analysis of the best multiple linear regression model, according to the following equation: \[\text{Y}=\text{a}_{1}\text{x}_{1}+\text{a}_{2}\text{x}_{2}+\text{a}_{3}\text {x}_{3}+\ldots\text{+a}_{n}\text{x}_{n}+\text{b}\] (6) \[\text{Y}=\text{estimated Forest Density value}\] \[\text{a}_{n}=\text{coefficients}\] \[\text{x}_{n}=\text{ explanatory variables (vegetation indexes)}\] \[\text{NDVI}=\frac{\text{NIR}-\text{RED}}{\text{NIR}+\text{RED}}\]0 Candidate models are assessed with the use of the Bravais-Pearson correlation coefficient, and the model with the higher coefficient was selected. \[\text{NDVI}=\frac{\text{NIR}-\text{RED}}{\text{NIR}+\text{RED}}\]1 (7) \[\text{NDVI}=\frac{\text{NIR}-\text{RED}}{\text{NIR}+\text{RED}}\]2 \[\text{NDVI}=\frac{\text{NIR}-\text{RED}}{\text{NIR}+\text{RED}}\]3 Finally, the most accurate model was applied over the complete area of interest in order to derive the Forest Density values. An example of Forest Density for the most recent 2019 epoch is shown in \[\text{NDVI}=\frac{\text{NIR}-\text{RED}}{\text{NIR}+\text{RED}}\]4 (iv) Quality Check of the consistency of the Forest Density between all epochsThe Computed Forest Density value over a stable area can present change over 1991-2019 due to several factors: for instance, EO imagery from various sensors or acquisition time of EO image reflecting seasonality.