Introduction 

According to Food and Agriculture Organization (FAO) definition, the forest is any “territory with arboreal coverage greater than 10% compared to an extension greater than 0.5 ha, where the trees reach a minimum height of 5 m when ripe”, including also: “[…] windbreak barriers and wooded strips of a width exceeding 20m, rubber tree plantations, and cork trees” ([10]). Thirty-one percent of the World’s total surface is occupied by forests (4.06 billion hectares) and Europe occupies a prominent position in this context since European forests account for 1,017,461 ha ([11]). Forests play a strategic role in terms of landscape, history, culture, and economy ([40]). They contribute significantly to human well-being, providing a multitude of ecosystem services among which wood production, carbon sink, food production, soil protection, climate control, biodiversity conservation, environmental decontamination, and hydrological cycle regulation ([47]).

The Italian forest heritage is one of the most important in Europe - 11.8 million ha or 40% of the Italian land ([34]) - but it is exposed to many threats, in the majority of cases directly or indirectly related to anthropic activities. Indeed, great pressure is currently given by climate change, which undoubtedly led to an increase in the occurrence of extreme events such as storms, fires, and droughts ([13]). Therefore, it results of primary importance to efficiently monitor forest ecosystems and their status, as well as to detect change processes and provide new tools for their conservation.

Forest change processes are territorial transformations of primary importance, and have a central role in the definition of policies aimed to preserve forest resources and its ecosystem services provision. International organizations and other institutions at different levels have undertaken forest changes issues, giving definitions to frame this phenomenon.

The United Nations Food and Agriculture Organization Forestry Department established a series of definitions related to changes between forest and other land use classes, including the loss of forest cover. According to these definitions, it is possible to distinguish two different processes of tree cover growth: (i) afforestation, i.e., the conversion from other land uses into the forest, or the increase of the canopy cover to above 10%; and (ii) reforestation, i.e., the re-establishment of forest formations after less than 10 years with less than 10% canopy cover due to human-induced or natural perturbations. In this sense, reforestation is presented as a short-term process that, if longer than 10 years, falls within the definition of afforestation ([10]).

Forest changes occur heterogeneously around the globe, mainly because of different socio-economic backgrounds. While in Africa and South America and Eastern Europe deforestation trends up, in the European Mediterranean areas afforestation, generally resulting from agricultural to forest land-use changes, represents broadly the dominant process ([33], [21], [37]). Regarding the Italian territory, historical cartographic and inventory data confirms this assumption, showing that forest surface has been constantly increasing through the years, with an average annual growth rate of 0.3% ([2]).

For the above reason, afforestation represents a process of increasing importance, which certainly needs to be characterized to make the most of its potential and to correctly drive its development in line with forest sustainable management principles. As on one hand, afforestation certainly involves all the benefits connected to forest ecosystem services, on the other hand, if not well managed, it can cause land fragmentation and other related issues ([50], [8]). In other words, most of the issues and challenges associated with afforestation depend on its location ([33]).

Since the afforestation process has a significant effect on the land cover dynamics, it is necessary to monitor its evolution, and to this purpose many different independent data in national and European environments are available, both inventory and cartographic. One of the most important is the Land Use Inventory of Italy (IUTI), an inventory data based on a sampling system and the classification in 6 macro-categories ([32]). Although IUTI is a valuable source of information, it is seldom updated (1990, 2000, 2008, 2013, and 2016) and provides information on the predominant land cover in the area of analysis (in Italy 1,217,032 random points distributed in as many quadrants of 25 ha, with a minimum mapping unit of 5000 m² and a minimum width of 20 m).

In this context, Remote Sensing (RS) certainly represents a high potential instrument ([29]. [12]), which allows to analyse and classify large areas, and enables to obtain frequently updatable results at different spatial scales ([1]). Remote sensing data analysis can be improved thanks to cloud platforms as Google Earth Engine^® (GEE - [19]), which allow managing large amounts of data and information even with small computational capabilities. Google Earth Engine^® provides numerous satellite missions as Landsat and the more recent European Sentinel Copernicus program which are by far the most used, being open-source data and covering many years of observations at different spatial resolution ([27]).

Information about land cover characteristics can be evaluated using the different bands of multispectral data, which can be summarised in single values, called Spectral Indices, or, when they are related to vegetated land cover, Spectral Vegetation Indices (SVI - [48]). The most used SVI in literature are the Normalized Difference Vegetation Index (NDVI), the Enhanced Vegetation Index (EVI) and the Normalized Burn Ratio (NBR). Besides the SVI, other indices known as Tasseled Cap (TC) represent a valuable tool for this kind of analysis. These indices resume multispectral data into four dimensions of the land cover considered: Brightness (B), Wetness (W), Greenness (G) and Angle (A). SVI and TC acquire even more value when referred to different periods, allowing to build temporal trends which represent important additional information about the phenology of a given territory ([38]). Hence, to elaborate a land cover classification and to produce forest cartography, an efficient method consists in analysing multitemporal series of SVI or TC in order to attain synthetic though very complete information. In doing so, it is possible to establish classification rules by imposing thresholds to discriminate different land cover classes and recognising the ones of interest ([46], [47], [31], [9]).

Landsat and Sentinel images have different applications related to their intrinsic characteristics. Landsat data are available since the 1970s, thus providing very long time data series at medium spatial resolution (30 m) which are useful to analyse land cover trends and evolution. Instead, Sentinel is available since 2015, providing images with high spatial resolution (10 meters) and short revisiting time (in Italy 2/3 days), which make it a tool of great utility for applications where high precision is needed. Sentinel-1 data and Sentinel-2 images, thanks to their spatial and temporal resolutions, have been recently used in many studies on the analysis of forest disturbances and land cover change detection. In Laurin et al. ([30]), Sentinel-1 data were used to detect the damages caused by the Vaia wind storm in 2018 in Northern Italy, while Sentinel-2 data were used to map clearcuts in Tuscany ([17]), and to estimate areas and relative confidence intervals of forest fires, wind damages and clearcuts in Italy for 2018 (Francini et al., in review).

Other research activities have taken advantage of the high temporal resolution of some satellite images to monitor land cover and land cover changes at the annual or sub-annual level. De Fioravante et al. ([9]) integrated the Sentinel-2 and Sentinel-1 data to assess land cover (2018) and land cover changes (2017-2018) in Italy, while Francini et al. ([15]) used PlanetScope imagery to develop a methodology to detect near-real-time forest disturbances in Tuscany (Italy). While the land cover and its changes monitored by remote sensing is a widely analyzed topic, afforestation assessment remains a theme mainly unrepresented, being treated as a central subject only in a few studies ([41], [44]) in which relatively low accuracy was reached (less than 75%). For example, Huang et al. ([26]) and Qiu et al. ([39]) used MODIS images (500 m spatial resolution) to map forest growth areas in China. However, products with a spatial resolution of 500 meters are not useful for monitoring afforestation in fragmented and morphologically complex areas typical of the Italian territory, where afforestation processes are characterised by patches of limited size.

The present research aims to illustrate a new methodology that exploits the Landsat time series to demonstrate that RS can provide valuable support not only for monitoring stable land cover classes ([9]) and forest disturbances ([18]), but also for predicting afforestation areas.

In the first section of this paper, we describe the training data and Landsat images (1984-2020). Then we explain the use of the input data to assess the temporal trend difference of photosynthetic activity of the classification periods (non-forest, afforestation and forest). Afterwards, we focus on the use of the Random Forest™ (RF) algorithm to automatically distinguish the classes and identify the combination of predictors which allows reaching the largest accuracy.

In the last part of the paper, we discuss the results and how they support the use of remote sensing data in the distinction between non-forest, afforestation and forest.

  Materials 

Study area and reference data

We selected sixty-one areas distributed along the Italian peninsula, especially on the mountains range (Alps and Appennines) where afforestation occurred in the period 1988-2020. Forty-two sites are located in areas with elevation higher than 600 m a.s.l., and twenty-six of them are located in areas above 1000 m a.s.l; the slope of almost 3/4 of the points is higher than 10% (Fig. 1). These areas cover a surface of 269.197 ha, and each polygon is on average 4 ha large. The areas were defined using the IUTI dataset and other images, as explained below.

Fig. 1 - Study area and distribution of reference dataset. On the right, an example of afforestation occurred between 1990-2000 is displayed: in the upper image (1988) there was no forest cover, in the centre (2000) the forest partially cover the area, in the last image (2019) the forest cover is complete.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Land Use Inventory of Italy (IUTI)

Land Use Inventory of Italy IUTI was realized by the Italian Ministry of Environment and Protection of Land and Sea in 1990 (updated in 2000, 2008, 2013, and 2016) in the framework of the Extraordinary Plan of Environmental Remote Sensing. It is composed of 1,217,032 randomly selected points, which have been classified through photointerpretation, considering a minimum mapping unit of 5000 m² and a minimum width of 20 m. The classification system is based on the Intergovernal Panel on Climate Change guidelines, and it is composed of three hierarchical levels, which aim to identify the land cover categories which are important for the Kyoto Protocol and to integrate the National Inventory on Forest and Carbon Pools (INFC) results on woods and other woody areas categories, based on the FAO definition ([42]). Analysing the evolution of the dataset through time makes it possible to identify those areas where afforestation occurred.

Landsat data and Best Available Pixel Composite

In this study, the Landsat surface reflectance data available (for Landsat mission 5-7-8) in the Google Earth Engine^® (GEE) archive were used. We used Landsat missions 5, 7 and 8, because the product of the previous missions had a lower spatial and temporal resolution. Moreover, we analysed data starting from 1985 because they are comparable with the training dataset. Landsat images are composed of six bands: three visible (blue, green, and red), one near-infrared (nir), and two short-wave infrared (swir1 and swir2). All of them are already processed to orthorectify surface and brightness (the thermal infrared) reflectance, atmospherically corrected, and cloud masked. Landsat images were used to create Italian annual composites from 1984 to 2020, from June 1^st until August 31^st, using the Best Available Pixel (BAP) procedure ([51]).

The BAP aims to fill the final image mosaic with the composite best available pixel surface reflectance value. The selection of the best pixel of the images collection takes into account four different criteria: (i) sensor score, to penalize Landsat 7 images where the Scan Line Corrector malfunction (SLC-off) is present; (ii) target day score, to preferably select the images acquired close to a defined acquisition day (in this case 15^th of August); (iii) distance to cloud/cloud shadow score, to decrease the scores of those pixels which are in the proximity the cloud cover; (iv) opacity score (calculated using opacity band produced by LEDAPS - [43]), to prefer the pixel with low opacity. Scores (i) and (ii) are applied to the whole image, while scores (iii) and (iv) are applied to each pixel. BAP composites calculation was performed using the GEE application developed by Francini et al. ([16]) which allows calculating and downloading BAP composites ([24], [25], [51]). For more information on the criteria involved in the BAP pixels selection see the BAP-GEE documentation ([14]). For more info on the BAP algorithm see White et al. ([51]) and Griffiths et al. ([20]).

  Methods 

The methodology presented here employed Landsat images to define the feasibility of the assessment of afforestation areas in the last thirty-six years through remote sensing. The analysis of many years allows verifying the evolution of the classification periods, considering the afforestation, which can have a different duration depending on the species.

Using 61 training areas consisting of polygons defined using IUTI elements, orthophotos and very high resolution images to identify the afforestation, a preliminary analysis of photosynthetic activity in non-forest, afforestation and forest areas was developed. Then RF algorithm, calibrated and validated, was applied to automatically distinguish the non-forest, afforestation and forest land cover (Fig. 2). The analyses were carried out using QGIS, GEE^® and R softwares, which are useful in processing spatial and statistical data.

Fig. 2 - Workflow of the applied methodology. The steps to distinguish among non-forest, afforestation and forest areas are briefly described.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Reference data acquisition

The orthophotos used for the photointerpretation process are available as Web Map Services (WMS) of the Italian Ministry of Environment, Land and Sea, and they are referred to 1988-1989, 1994-1998 (black and white) and to 2000, 2006, 2012 (colour), with a spatial resolution of 50 cm. The more recent Very High Resolution images (30 cm spatial resolution, available since early 2000s) are freely available thanks to the web service of QGIS and Google Earth Pro^®.

The reference dataset (Fig. 1) was produced through photointepretation of orthophotos and Very High Resolution images, available from 1988 to 2020; in addition to known areas, the IUTI dataset has been used to establish the afforestation areas by considering the points which identify a forest change (from non forest to forest area). Sixty-one training areas have been collected and finally integrated with other information on observation date and polygon dimension.

Orthophotos and IUTI datasets were used to define the temporal range in which the afforestation process occurred. More specifically, two years were identified: (i) the year corresponding to the transition between non-forest and afforestation; and (ii) the year corresponding to the transition between afforestation and forest.

Temporal predictors calculation and analysis

The indices calculation for each composite BAP and the analysis of the resulting temporal series were carried out to calculate the temporal statistics. More specifically, we augmented the six Landsat bands of our BAP composites (see above) by calculating seven additional indices: (i) Normalized Difference Vegetation Index (NDVI); (ii) Normalized Burnt Ratio (NBR); (iii) Enhanced Vegetation Index (EVI); (iv) Tasseled Cap brightness (B); (v) wetness (W); (vi) greenness (G); and (vii) angle (A). As a result of this step, we obtained 13 temporal series of 36 years, hereinafter called Time Series (TS - Fig. 3).

Fig. 3 - Example of NDVI series in an afforestation pixel. For each classification period the trend line was obtained and used to calculate the Pearson’s correlation coefficient (r).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

For each polygon and for each of the non-forest, afforestation, and forest temporal ranges, we calculated from the TS four temporal statistics: (i) the mean and (ii) the standard deviation of the TS; (iii) the slope; and (iv) the Pearson’s correlation coefficient (r) of the regression line obtained by linear interpolation of the TS over time. As a result of this step, a set of 52 temporal predictors (13 TS × 4 temporal statistics) were obtained to distinguish photosynthetic activity between non-forest / afforestation / forest areas and to investigate the differences in photosynthetic trends using density plots.

Random Forest

Random Forest model and hyperparameters

To classify forest, non-forest, and afforestation areas we used a random forest model and the temporal predictors as variables. Random Forest™ (RF) is a decision tree algorithm often applied to forest classification ([23]), where a subset of predictors is used to build a set of regression trees, applying the Out-Of-Bag sample (OOB) procedure ([6]). RF is also useful as it allows to assess the variable importance in the model. RF algorithm was chosen because it achieves large accuracies when compared with other machine learning approaches ([35]) and it is not influenced by overfitting ([4]). Random Forest hyperparameters are model parameters that can be calibrated to obtain the minimum error rate ([30]). In this study, we tuned the following two: max features (the number of features considered) and max depth (the maximum depth of the tree, the longest path between the root node and the leaf node). The number of trees in the algorithm was set to 300 and the minimum number of samples for splitting a node was set to one.

Max depth and max features calibration

Max features and max depth hyperparameters were calibrated using the procedure named random search ([5], [30]). Specifically, we defined a grid of max features (1-52) and max depth (1-40) ranges, from which a subset of 100 combinations were randomly selected. To assess the performance of each tested combination we used the k-fold Cross-Validation (CV) with k = 5.

Cross Validation serves to train and validate the model on different data and it is useful to reduce the overfitting ([52]) and to test the robustness of the model ([30]). Using k-fold CV, the dataset is divided into k groups, which are composed of the same quantity of data: one of these folds is defined as the validation group, while the remaining k-1 groups are considered as training data. The procedure is then repeated k times, each time using a different fold as a validation group. Finally, the validation data is merged and the k-fold CV Overall Accuracy (OA_cv) is calculated. As a result of this step, we obtained the max features - max depth combination that allowed to obtain the larger OA_cv.

Performance assessment

Although the k-fold CV procedure (see above) allowed to rigorously estimate the accuracy (OA_cv), we further assessed the accuracy of the model exploiting the best max features - max depth combination using Out Of Bag (OOB) data, which allows getting an unbiased estimate of the error ([6]). During random forest training, two independent sets are created: (i) the bootstrap sample, i.e., the data chosen to be “in-the-bag” by sampling with replacement; and (ii) the out-of-bag (OOB) sample, i.e., all data not chosen in the sampling process. During model training, this process is repeated and many OOB samples are created. Since they are never used to train the model, they are never-seen-before data and can thus be used to assess the performance of the model. Finally, all OOB samples are aggregated and OOB predictions are compared with reference data by constructing a confusion matrix from which we calculated the OOB. Overall Accuracy (OA_oob), omission errors (omissions_oob), and commission errors (commissions_oob). The procedure was repeated four times, to compare the classification models as follows: (i) non-forest/afforestation, (ii) afforestation/forest, (iii) non-forest/forest, and (iv) non-forest/afforestation/forest. For each of them, classification omission and commission error were calculated, to respectively identify the pixel non-classification and the misclassification.

Variables importance calculation

To identify temporal predictors that mostly contributed to increasing the accuracy of the model we analysed the variable importance ranking outputted by random forest. The variable importance was expressed in terms of the Mean Decrease of the Gini index (MDG - [36]). MDG is the measure of the variable contribution to the homogeneity of the random forest trees component. Large MDG values correspond to variables that relevantly contributed to increasing the performance of the model and thus to the accuracy achieved. Conversely, small MDG values correspond to variables that did not contribute to increasing the accuracy of the model.

  Results 

Classification accuracies

Random forest validation using Out-Of-Bag data (OA_oob) was performed on all the temporal predictors of the best iteration of 5-fold CV (Tab. 1). The largest OA_oob was reported by the model classifying between afforestation and forest (80%), followed by the classification between non-forest and forest (77%). Non-forest and afforestation and non-forest/afforestation/forest classification model reached an overall accuracy smaller than 70%.

For each class of the four classification models, omission (omissions_oob) and commission (commissions_oob) have been calculated (Tab. 2). Both omission and commission errors are larger than 10%, with the maximum value (40%) registered by non forest class in the classification between non-forest, afforestation and forest; the smallest error values are registered in the afforestation and forest classification, where omission error of the forest class is 11% and the afforestation one is 27%; as regard commission error, the afforestation registered the smallest value (14%), while the forest one is 24%.

The large random forest OA_oob values were reached thanks to the use of the hyperparameters with the largest OA_cv, calculated using random search and 5-fold CV on the classification models, to identify the couple of hyperparameters which allow to reach the best results. Fig. 4 underlines that the OA_cv is larger than 65% in every hyperparameters combination considered. Regarding the best iteration (value “max” in Fig. 5) the smallest OA_cv is registered in the non-forest, afforestation and forest (72%) and non-forest and afforestation (75%) classification models. The model performance is larger in the classification of non-forest/forest and afforestation/forest, where the best iteration OA_cv reaches 83% in the non-forest/forest classification and 87% in afforestation/forest one.

Fig. 4 - Five-fold cross validation overall accuracy (Oa_cv). The value “max” is referred to the best iteration. Each OA_cv of every couple of hyperparameters for each classification models is displayed, grouped by accuracy.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 5 - Importance ranking of temporal predictors. The importance value (Mean Decrease Gini) of each temporal predictors in every classification model is reported.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Random forest temporal predictors importance and classification periods differences analysis

The importance ranking of temporal predictors is illustrated in Fig. 5. The most important predictors for the model afforestation/forest were the median value of the green (MDG = 2.62) and the blue (MDG = 2.62) bands, the Tasseled Cap Greenness R (MDG = 3.04), and the slope value of the Tasseled Cap Brightness (MDG = 1.76) and the green band (MDG = 1.72). The other predictors had smaller importance scores. In the non-forest and forest models, only two predictors were indicated as more important, that is the median value of the green (MDG = 10.35) and the blue (MDG = 19.36) bands; the same ranking was registered in the non-forest, afforestation, and forest models, even if the other predictors had a larger score than the ones of forest, non-forest classification. Finally, the first predictor in the importance ranking of non-forest and afforestation models is the R value of the EVI (MDG = 9.90) and, to a lesser extent, the slope value of the blue band (MDG = 2.69).

Fig. 6 displays the temporal predictors of the non-forest, afforestation and forest, related to the variables which were considered as most important in the random forest. Tasseled Cap Greeness R value (TCG - Fig. 6d), registered the largest difference between data, in particular the non-forest reaches the largest values, as in the EVI (Fig. 6c). As regard the median values of blue (Fig. 6a) and green (Fig. 6b) bands and the slope values of green (Fig. 6e) band and TCB (Fig. 6f), the forest reaches the largest values, and the afforestation and the non-forest are quite similar.

Fig. 6 - Density distribution plots. The most important temporal predictors of non-forest, afforestation and forest are displayed; plots allow to estimate the differences of the classification periods.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Discussion 

Afforestation monitoring is important because it influences ecosystem services but it can also cause some issues, as land fragmentation and habitat loss. In this research we aimed to demonstrate the efficacy of remote sensing data in the assessment of the difference in photosynthetic activity and in the automatic classification of non-forest, afforestation and forest areas. To classify them, we analysed 36 years of Landsat images using sixty-one training areas and a random forest algorithm. The large overall accuracy reached with the classification algorithm supports the approach chosen, confirming the efficacy of the remote sensing in the automatic classification between non-forest, afforestation and forest.

The methodology by using Landsat images allows analysing long time series of images, which is essential to assess forest evolution, which is generally slow. Moreover, both the first analysis on temporal predictors, and the random forest algorithm are easily replicable, also using different input images. To limit the overfitting due to the use of neighbouring pixels as training and validation data the data were elaborated at polygon level, considering the training data boundaries. Furthermore, random search and k-fold cross validation have been useful to identify the combination of random forest hyperparameters with the highest accuracy, to ensure the maximum performance of the algorithm. Considering the results of the random forest validation, the difference in OA_oob values can be due to the similarity in photosynthetic activity in non-forest and afforestation land covers, which registered smaller accuracy values, while forest and non-forest or forest and afforestation classification models gives better results because there is a stronger difference in the spectral response of the two classes considered. It is notable that, even if the k-fold cross validation and the Out-Of-Bag accuracy results are not comparable, the OA_cv and OA_oob ranking is the same, with the largest accuracy value registered in distinguishing between afforestation and forest and the smallest ones in distinguishing among non-forest, afforestation and forest. This is important to support the results and the methodology, because the same result is obtained from two independent validation models.

The analysis of the temporal predictors importance allowed us to identify those variables which were more significant to separate the three classification periods. The median value of the green and the blue bands was part of the main predictors in almost every distinction considered; this probably concerned the fact that the blue band is associated with soil moisture, while the green band is associated with the reflectance value of the vegetation: the different water content and the different reflectance value are therefore two important factors to examine to classify non-forest, afforestation and forest. Considering the other predictors in the importance ranking, the Vegetation Spectral Indices and the Tasseled Cap in the first positions were those ones which are related to the blue band (as TCB, EVI), to the green band (TCG) and, to a lesser extent, those related to nir and red bands, as the Normalized Difference Vegetation Index.

These are important information to set the best combination of parameters to obtain the maximum accuracy in the distinction of non-forest, afforestation and forest, using a small number of input data. The results of the analysis of the temporal predictors confirm that the remote sensing data allow to assess the differences in the trend of photosynthetic activity in relation to non-forest, afforestation and forest. The median value of blue and green bands, r value of EVI, TCG, TCB and green slope values, which were the most important variables in the random forest ranking, showed the most significant dissimilarity in the photosynthetic activity among classes. This can be due to several reasons, not only to the spectral response of the different vegetation type, but also to their phenology; in fact, a non-forest land cover, like pastures or other grasslands, could be more influenced by the microclimate or the soil moisture than the afforestation or forest land cover. Similarly, shrublands could be more susceptible to local weather condition than the forest, which is the most stable land cover.

There are few studies on the use of remote sensing to assess afforestation areas, and most of them evaluate the forest gain by comparing different land cover maps. For example, Shen et al. ([44]) estimated afforestation areas at a regional scale using annual forest mask and combining them to extract the afforestation areas. Brovelli et al. ([7]) and Townshend et al. ([49]) assessed forest cover change using machine learning and remote sensing data through the comparison of two forest cover maps. Other researches examine the forest gain at global level using remote sensing, comparing forest and non-forest land cover, but they do not consider the afforestation land cover ([22], [28], [45]).

One of the limits of the method proposed is the medium spatial resolution of images: in fact, Landsat images resolution is 30 meters, and this could cause the omission of the smaller elements. This could be overcome by integrating other data, like radar, lidar, or other multispectral images. Another issue could occur by increasing the area of study outside the training areas; in fact, it will be necessary to consider potential commission or omission errors which could be due for example to the orography, which could reduce the accuracy because of slope and shadows.

The results of the present research offer the possibility to deepen the existing analysis approaches, directly highlighting the afforestation, a land cover class which is not stable through time, but represent a transition from non-forest to forest cover. Furthermore, it will be possible to elaborate an afforestation areas map, at national or European level. The new cartographic data could satisfy the requirements proposed by European legislation for monitoring activities, being aligned to the EIONET Action Group on Land monitoring in Europe (EAGLE) Concept. EAGLE is based on a classification system which considers the distinction between land cover and land use and it is composed of three descriptors (land cover components, land use attributes and further characteristics), combined to define a classification system suitable to be integrated to existing classes, maintaining the three components independent ([3]). The classification structure of the EAGLE matrix specifically includes the afforestation land cover class. Considering this new classification method, the map will be aligned not only to the existing cartographic data, but also to the future steps of the land cover monitoring and classification planned by Europe and Italy, being easily updated in all the descriptors.

  Conclusions 

Non-forest, afforestation and forest can be automatically predicted using Landsat data and classification models aimed to evaluate the different responses of the photosynthetic activity of non-forest, afforestation and forest areas. Since this approach is suitable for application at different scales of analysis and thanks to the results obtained in this study, the methodology will be applied to the entire Italian national area, to locate the potential afforestation areas outside the training areas. This step will allow to create a map of afforestation areas across Italy. This product follows the indications given in article 15 of the Italian national legislation on forestry and forestry supply chain (legislative decree 34/2018), which encourages the creation of a geo-referenced forest mapping and contributes advantageously to forest planning and management activities. In fact, the expansion of forest areas is often due to the agricultural and pasture surfaces loss ([34], [40]), and this could bring to an uncontrolled vegetation growth, which needs to be located and monitored, to avoid negative consequences, such as fires, that could be more destructive if the vegetation is not managed.

The afforestation area assessment can be important to analyze the forest growth trend, allowing to make forecasts on the forest evolution, both in quantitative and qualitative terms. Moreover, the integration between the map and other data will allow to study the evolution of the forest areas to assess the forest ecosystem services.

  Acknowledgements 

Conceptualization: SF, AC, MMa, GSM, GC, MMu; methodology: SF; statistical analysis and results elaborations: AC and SF; writing-original draft preparation: SF, AC, VF, GC, GLS, CC, LC, MMa, GC, GSM, MMu. All authors have read and agreed to the published version of the manuscript.

  References