Analysis of Normalized Difference Vegetation Index (NDVI) multi-temporal series for the production of forest cartography

Abstract

This study aims to produce a forest raster cartography, valid for the Italian territory and characterized by a spatial resolution of 10 m, using remote sensing techniques. In detail, Sentinel-2 multispectral images have been used in order to obtain Normalized Difference Vegetation Index (NDVI) values, useful to distinguish vegetation among all the different other kinds of land cover.

An average monthly value of NDVI has been calculated using the cloud service Google Earth Engine (GEE), by processing all the Sentinel-2 images acquired in a time range of 18 months, from March 2018 to August 2019. Analysing these multi-temporal NDVI series, forest cover has been classified through decision rules aimed to distinguish the different phenology of broad-leaved trees from coniferous and from other types of vegetation, such as meadows, shrubs and vineyards, and understand how external factors such as elevation and latitude have an impact on these patterns. In order to assess the performance of the methodology, it has been tested in 3 study areas of Italy, characterized by different vegetation cover and latitude.

This study has conducted to satisfying results with the chance of giving both a qualitative and quantitative characterization to the different temporal NDVI patterns useful to define an index, called NDVI summer max, for the classification of the forest land cover. It has been possible so, just with the use of this index, to obtain cartographic products defined by values of global accuracy around 60%. As well, the study has highlighted a series of criticalities whose resolution can be a good starting point for further researches such as the necessity of a better comprehension of disturbance element (shrubs, meadows and vineyards) to reduce commission errors and the use of ancillary data to link the NDVI responses to geo-specific factors and meteorological conditions.

Keywords: Remote sensing, Forest monitoring, NDVI.

1. Introduction

Forests are an essential resource, providing several services (known as “Ecosystem Services”) such as goods supply, soil-water conservation, climate regulation, generation of natural habitats and decontamination of the environment (Tao et al., 2012; Nink et al., 2019). In order to protect this important natural resource, the monitoring and mapping of forest land is needed.

In this study, we use the Food and Agriculture Organization (FAO) definition of forest that is “any land with a tree cover major than 10% and with an extension of minimum 0.5 ha in which trees have to reach at least 5 m of highness to maturity” (SIAN, 2015). Considering this definition, forest land in Italy extends for 11,8 millions of hectares, corresponding to 39,4% of the national surface (Munafò et al., 2018). Italian forest patrimony is incredibly various and it comprises very different arboreal species depending on altitude, latitude and climatic conditions (Aleffi et al., 2007).

There are several tools, inventory and cartographic, designed to map forest in Italy. The most accurate, in terms of spatial resolution, is the High Resolution Layer Forests (HRL Forests), included in Copernicus project, which is a raster layer characterized by a pixel size of 20 m. It is available for two different reference years, 2012 and 2015, and it distinguishes the land cover classes of Coniferous forest and Broadleaved forest from all other non-forest areas (European Environment Agency Eea, 2017).

The project in which this study fits aims to elaborate cartographic products defined by a spatial resolution of 10 m and updatable every year.

In this context, Remote Sensing techniques find an interesting application, if considered as “the scientific discipline that permits to obtain qualitative and quantitative information about distant objects and on the surrounding environment, on the bases of measures of electromagnetic energy that, emitted, reflected or diffused, interacts with the surfaces of interest” (Brivio et al., 2006).

For this study, the images used have been acquired by satellites Sentinel-2, developed by European Space Agency (ESA) in the frame of land monitoring service “Copernicus”, whose sensors, produce multi-spectral images characterized by thirteen different spectral bands (Drusch et al., 2012). By combining the different values of reflectance related to the different bands, parameters can be obtained that summarize the entire multispectral information. These parameters are known as spectral indices (Oberti, 2009). Among these, there are the Spectral Vegetation Indices (SVI) that allow to grasp the presence of vegetation on Earth’s surface and its variation (Spinsi et al., 2012).

Several studies have aimed at monitoring forest and forest disturbance, with different approaches such as: the use of long time series of Landsat satellite data for the detection of forest changes for more decades (Alencar et al., 2020), the combined use of Sentinel-2 images and Spectral Angle Mapper (SAM) to map forest species in mountain regions (Mohajane et al., 2017), the elaboration of Sentinel-2 images to get NDVI time series to map Mangrove species (Li et al., 2019).

Also, different studies have been carried on considering different SVI, such as Enhanced Vegetation Index (EVI) (Bajocco et al., 2018), Plant Phenology Index (PPI) (Karkauskaite et al., 2017) and Normalized Difference Vegetation Index (NDVI) (Jonsson et al., 2018).

NDVI index estimates the photosynthetic activity by a combination of Red (R) and Near Infrared (NIR) bands that are generally conditioned by the presence of chlorophyll. NDVI index varies in a gap included between −1 and +1 and increases at the increasing of vegetation land cover. Threshold values such as 0.2 can be considered to distinguish vegetation from other kinds of cover (Spinsi et al., 2012).

The information that can be deduced from this index are much more relevant if considered with their temporal trend, in order to have a clearer differentiation between the different land covers, being aware of the leave-on/leave-off periods of the trees and of the other kinds of vegetation (Zeng et al., 2019; Yang et al., 2018). In some other researches, this approach has been considered but with the utilization of MODIS images, characterized by a temporal accuracy of 2 days and so potentially very useful to get a clean and complete data series (Beck et al., 2005). The issue connected with the use of MODIS images stands for its very low spatial resolution (from 250 m to 1000 m, depending on the spectral band considered). That’s the reason why our choice fell on the use of Sentinel-2 images, characterized by a revisiting time of 5 days, that still allows to get a satisfactory amount of clean images, though by a spatial resolution of 10 m.

This method has been supported, in the last years, by several studies which reached good results defined by a high resolution (Yang et al., 2018). This approach though finds it difficult to distinguish trees from other kind of vegetation forms (Bayle et al., 2019). This issue is particularly important for the Italian territory which is rich of shrubland, especially in proximity of the coast, natural meadows, vineyards and other no arboreal vegetation forms (Blasi, Biondi, 2017).

Considering this, this study was structured in order to achieve the following objectives. Firstly, the one to extract and analyse temporal series of NDVI to get some rules useful to improve the classification and, in particular, the distinction between arboreal and no arboreal land cover. Secondly, to apply these rules and obtain forest cartographic products and evaluating its accuracy.

2. Materials and methods

The workflow shown in Fig. 1, represents the different steps in which the materials and methods section is subdivided.

2.1. The study areas

The study has been carried on for three different zones of the Italian territory, each one representative of the different landscapes of the North, centre and South of the Country, in order to analyse as well the influence of latitude on the NDVI patterns. The zones have taken in correspondence of the Sentinel-2 images tiles format, being so squared areas of 100 × 100 km and identified by an alpha-numeric code of two numbers and three letters.

The first study area (identified as 33TUL) is situated in the North-East of Italy and it partially covers the two regions of Friuli-Venezia Giulia and Veneto (some territory of Slovenia and Croatia are included in this tile but have not been considered in this study). The second study area (identified as 33TTG) is situated on the west coast of the Italian peninsula and it extends over a relevant part of the region of Lazio, including the city of Rome. The third study area (identified as 33SWC) covers the southernmost part of Calabria region and a portion of the North-East side of Sicily region.

2.2. Definition of the ROIs

With a view to a semi-automatic classification method, used for the production of forest cartography, it has been necessary to define some Regions Of Interest (ROI) (or training areas) to be used as samples to analyse, in order to define the right classification rules to apply then automatically to the entire study areas (Jiang et al., 2012). With the use of the program QGIS, it has been possible to draw approximately 150 polygons representatives of the ROIs for each study area. To each polygon a land cover class has been assigned, through photo-interpretation methods, considering the ISPRA (Italian Superior Institute for the Protection and the Research of the Environment) classification code and the following categories: broadleaves (code 211), conifers (code 212), meadows (code 221) (Munafo et al., 2019). Particular attention has been put in choosing homogeneous training areas, namely formed by the same typology of vegetation. The three study areas and the training areas taken inside those are represented in the following Fig. 2.

2.3. Extraction of multi-temporal series of NDVI for the ROIs

After several comparisons among the most used SVI, the Normalized Difference Vegetation Index (NDVI) appeared to be the most appropriate for the purpose of this study.

In particular, at first, apart from NDVI, also the Enhanced Vegetation Index (EVI) and the Normalized Burn Ratio (NBR) had been taken into account. The entire methodology workflow has been repeated as well for these other two indices (NBR has been used in a complementary way, since everything it doesn't recognize as burn is therefore vegetation) to see which one could bring to the most accurate results. This test revealed that NDVI leads to the best result, since EVI appeared to be very conditioned by atmospheric effects and since NBR couldn't describe the seasonal variation of the trees phenology at a satisfactory level.

Multi-temporal series have been taken over a period of 18 months, from March 2018 to August 2019, resuming an average value of NDVI for each month.

To get a multi-temporal series of NDVI for each ROI, the cloud-based platform Google Earth Engine has been used, whose big computational power has allowed to obtain 18 values of NDVI (one per month) for each training area (Gorelick et al., 2017). Each single value is the result of an average done over the values of every pixel contained in the ROIs and previously obtained as median of all the NDVI values related to the single months. The statistic tool of median has been chosen so to cut off the possible presence of outliers (due to probability of having cloud coverage, snow, or other disturbances) (Bajocco et al., 2018). Furthermore, for computational costs reasons, original values of NDVI have been multiplied for 10, so to avoid decimal numbers and reduce the calculation time.

According to this approach, it is expected to have, as result, different NDVI patterns depending on the type of vegetation land cover. For instance, broadleaves are expected to be characterized by a seasonal trend, with maximum points reached during summer period (leave-on period) and minimum during winter (leave-off period), while coniferous are expected to have a more constant pattern, considering they don't have leave-on/leave-off periods (Eklundh et al., 2012).

2.4. Multi-temporal series analysis - extraction of phenological parameters

Several analysis on multi-temporal series have been carried on for this study to understand the phenological behaviour of the different land cover classes. In particular, some of these were focused on the characterization of different types of trees phenology, such as broadleaves and coniferous, while others have been done to understand how trees phenological behaviour changes under different conditions such as latitude, altitude and distance from the sea.

Here the two main ones are reported: the first on the comparison between NDVI patterns of broadleaves, coniferous and natural meadows within the same tile and the second on the comparison between broadleaves phenological patterns related to different tiles. They both have the aim of giving an analytical description to the patterns and to deduce a methodology of classification of arboreal land cover to be used in the next phase.

The main idea is the one of distinguishing trees from other kind of vegetation by considering the NDVI values relatively to the summer period, in which broadleaves and conifers reach their highest values while other types of vegetation reach their lowest (Bajocco et al., 2018).

In order to identify phenological patterns, different methods were tested.

A first one has been to identify certain NDVI thresholds and look at the periods of time during which the different patterns had values higher than these. To make this process more easily reproducible, this description has been done through a binary code 0-1 (where 0 stands for the period in which the pattern is lower that the threshold and 1 for the period in which is above it). Particular interest has been given to the crossing of the thresholds of values 6 and 7 of NDVI, useful to identify the periods in which broadleaves had their maximum values.

A second one has been the Dynamic Threshold Method (DTM), appositely projected to give an analytical description of seasonal pattern of vegetation index. After a previous data regularization, the method allows to extract "phenological parameters" known as SOS (Start Of Season), EOS (End Of Season) and L (length of season), graphically represented in Fig. 3, useful to represent the summer period of maximum NDVI values (Zeng et al., 2019). In particular, data regularization has been done with the Stavitzky-Golay filtering method (whose algorithm is below shown) by which original data are replaced by a linear combination of near values (Eklundh, Jonsson, 2017; Granero-Belinchon et al., 2020):

∑_j=-m^m c_j y_i+j

In which the different weights are given by:

c_j = 1 / (2m + 1)

In this case, m has been taken equal to 5, in order to get a moving average in a window of five consecutive values.

2.5. Arboreal land cover classification

The information previously calculated with the analysis on the multi-temporal series have been used to execute a forest land cover classification, with a particular focus on distinguishing trees from other vegetation forms. The main issue considered during this phase in fact was to recognize no arboreal vegetation form to exclude it from the classification (De Lillis, Fontanella, 1992). Elements like shrubland, vineyards and meadows have a phenological behaviour similar to the one of trees and that's the reason why analysis have been done comparing these two classes in order to get a classification rule. To accomplish this, as announced before, the different phenological behaviours have been considered during summer period since the analysis showed that's the period in which the differences between NDVI patterns are maximum (Tang et al., 2015; Bajocco et al., 2018). The workflow represented in Fig. 4 shows the main stages in which the process of classification has been subdivided.

Firstly, a raster layer with maximum values of NDVI during the summer period (renamed "NDVI summer max") has been elaborated. To define the summer period, the results of all the different analysis previously done have been considered, with the different methods of DTM and NDVI thresholds, separately for each of the three reference tiles. Having the NDVI summer max raster layers, it has been possible to define a threshold value considering the NDVI summer max values relatively to the ROI defined as tree. The threshold, in particular, has been chosen as the tenth percentile of the distribution of the ROIs values of NDVI summer max, so to exclude eventual outliers. This threshold value has been later applied to the NDVI summer max values of all the pixels of the three tiles so to separate the two land cover classes. This can be graphically represented by the histogram shown in Fig. 5, where the blue block of pixel values represents the pixel classified as tree cover and the green block represents the pixel associated to non-tree cover.

2.6. Validation of the results

Once obtained the forest cartography with the previous stage of classification, a phase of validation has been carried on to get the accuracy of the cartographic products and understand, for each tile, which analysis method (and corresponding period) has led to the best result.

To accomplish this, 200 points photo-interpreted for each tile have been used to represent the ground truth to calculate the confusion matrices and, consequently, the accuracy indices. Confusion matrices present the scheme shown in Table 1 and they compare the ground truth with the classification data. The index k represent the number of land cover classes considered in the classification and the index n is the total number of point used, hence 200 in this case. The items on the major diagonal (aii) represent the number of samples correctly classified, while the other items represent classification errors (Congalton and Green, 2009).

Table 1. Confusion matrix.

For this study three accuracy indices have been considered:

Global accuracy (GA), defined as the ratio between the number of samples correctly classified and the total number of samples (n):

GA = (∑_i=1^k a_ii) / n

User’s accuracy (UA), for each class, is expressed as the ratio between the number of samples correctly classified and the column total:

UA_i = a_ii / a_+i

Producer’s accuracy (PA), for each class, is defined as the ratio between the number of samples correctly classified and the raw total:

PA_i = a_ii / a_i+

The first one gives an overall statement over the goodness of the product. The other two are inversely proportional respectively to the presence of commission errors (false positives) and omission errors (false negatives) (Congalton and Green, 2009).

3. Results and discussion

3.1. Definition of the ROIs and extraction of the multi-temporal series

The ROIs have been taken homogeneously over the territory of the three different tiles. From each one of these a multi-temporal series has been extracted. Three different examples are shown in Fig. 6, relatively to the three main different types of vegetation land covers. As expected, broadleaves have a phenological pattern that reaches high values of NDVI during summer periods and low values during winter, coniferous have a constant pattern on high values of NDVI. Meadows present more irregularity but, in the majority of the cases, they have low values during summer periods. Besides, it has to be considered that every example has its own peculiarities and outliers due to site-specific conditions.

Table 2. Summer periods extracted with the three different methods and used for the classification of forest land cover.

3.3. Arboreal land cover classification and cartographic results

The classification has been carried on with the methods described in the previous chapter and with the summer period calculated during the analysis of the patterns and reported here in the following Table 2.

In Fig. 9 are reported the classification results for each study area, in which green represents the forest territory, yellow represents the non-tree cover and black represents the non-vegetation cover, compared with a satellite image representing the same zone.

Also, some zoom are given in Fig. 10 to show, in different scales, the high resolution that characterize the land cover cartography obtained. In particular, comparing Fig. 10a and b, the high level of accuracy in recognising border lines between the different classes can be appreciated. Observing Fig. 10c instead, it can be seen, at a pixel level scale, which is the level of detail reached by the classification and its better operating when classes are more compact.

3.4. Validation

Accuracy indices, reported in Table 3, have been considered relatively to each period, calculated by a different method, used in the classification stage for each tile. The aim was that to determine the reliability of the cartographic products generated and, at the same time, obtain a comparison between the different methods.

Global Accuracy values are averaged around 60% percent, a result that is acceptable thinking that it has been reached with the use of a single classification index. In the frame of these study, other indices are provided to work together with NDVI summer max index and determine so remarkable accuracy values. Observing the other two indices, it can be noticed that Producer Accuracy stands at high values (close or higher that 90% in the majority of the cases) while User Accuracy in maintained on low values. This fact shows that the classification operated works well in recognising trees (high PA stands for low omission errors), without confusing them with other vegetation forms. More difficulties have been encountered in distinguishing meadows, bushes, vineyards and other non-arboreal vegetation forms. This still remains a challenge to face (De Lillis, Fontanella, 1992). A purpose for a future research could be to focus specifically on the phenological characterization of these disturbance elements, possibly repeating for them the same analysis conducted for trees in this study, in order to exclude them more easily from the arboreal classification.

Consistent differences between the methods used haven't been found for the tiles 33TUL and 33TTG, while for the tile 33SWC the method NDVI threshold 6 shows a consistently worst result. This is probably due to the short temporal period, of the only month of May, that the method has extracted and that is less strong than the other two ones in excluding false positive, as the 56% of UA shows.

Considering Global Accuracy as the most descriptive index, it can be seen that NDVI threshold 6 brings to the best results in all the different tiles. In the tile 33TUL the best result is also obtained with the use of the method DTM.

Table 3. Accuracy indices relating to the different forest land cover cartography generated with the different methods.

4. Conclusions

The analysis executed on the NDVI multi-temporal series turned out to be useful for the characterization of the different vegetation land cover phenology and to obtain cartographic products very useful to the forest monitoring and forest change detection.

The methods used have allowed to obtain valuable results in terms of phenological parameters, summer periods and accuracy of classification. For this last point, in particular, it has been possible to get an averagely 60% accuracy of classification with the only use of the index NDVI summer max.

At the same time, this study allowed to identify some directions in which to move for a future research on this topic.

Firstly, the same approach could be used with time series longer than 18 months and defined for shorter intervals than the monthly ones. This way, more precisely defined summer periods would be obtained. Secondly, a more detailed study could be done on the disturbance elements such as meadows, bushes and vineyards in order to reduce commission errors and improve Global Accuracy of the cartographic products. Finally, ancillary data on arboreal species, meteorological conditions and, in general, on territorial variables could surely give a more clear comprehension of the NDVI responses. With these additional information in fact it would be possible to recognize how external factors such as the ground composition, the average steepness and the rainfall pattern affect the phenological trends and so achieve a more complete and accurate arboreal classification.

CRediT authorship contribution statement

Gian Luca Spadoni: Conceptualization, Methodology, Validation, Visualization, Writing - original draft, Writing - review & editing. Alice Cavalli: Investigation, Formal analysis, Writing - review & editing. Luca Congedo: Software, Data curation, Writing - review & editing. Michele Munafò: Resources, Project administration, Supervision, Writing - review & editing.

Declaration of competing interest

None.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.rsase.2020.100419.

References

1. Aleffi, M., et al., 2007. I Boschi Italiani Dalle Alpi Al Mediterraneo, Techno ed.
2. Alencar, A., et al., 2020. Mapping three decades of changes in the brazilian savanna native vegetation using Landsat data processed in the Google Earth Engine platform. Rem. Sens. 12, 924.
3. Bajocco, S., et al., 2018. Remotely-sensed phenology of Italian forests: going beyond the species. Int J Appl Earth Obs Geoinformation 74, 314.
4. Bayle, A., et al., 2019. Improved mapping of mountain shrublands using the Sentinel-2 red-edge band. Rem. Sens. 11, 2807.
5. Beck, P.S.A., et al., 2005. Improved monitoring of vegetation dynamics at very high latitudes: a new method using MODIS NDVI. Rem. Sens. Environ. 100, 321.
6. Blasi, C., Biondi, E., 2017. La flora in Italia. Ministero dell'Ambiente e della Tutela del Territorio e dell'Atare.
7. Brivio, P., Lechi, G., Zilioli, E., 2006. Principi e metodi di Telerilevamento, CittàStudi ed.
8. Congalton, R.G., Green, K., 2009. Assessing the accuracy of remotely sensed data. Int. J. Appl. Earth Obs. Geoinf. 11, 448.
9. De Lillis, M., Fontanella, A., 1992. Comparative phenology and growth in different species of the Mediterranean maquis of central Italy. Vegetatio 99, 83.
10. Drusch, M., et al., 2012. Sentinel-2: ESA's optical high-resolution mission for GMES operational services. Rem. Sens. Environ. 120, 25.
11. Eklundh, L., Jonsson, P., 2017. Timesat with Seasonal Trend Decomposition and Parallel Processing. In: Software Manual: Lund University.
12. Eklundh, L., et al., 2012. High Resolution Mapping of Vegetation Dynamics from Sentinel-2 to Sentinel-2 Preparatory Symposium Rome. European space agency, p. 7. April 23-27.
13. European Environment Agency Eea, 2017. Copernicus land monitoring service - high resolution layer forest. Product specifications document 1-38.
14. Gorelick, N., et al., 2017. Google Earth engine: Planetary-scale geospatial analysis for everyone. Rem. Sens. Environ. 202, 13.
15. Granero-Belinchon, C., et al., 2020. Phenological dynamics characterization of alignment trees with Sentinel-2 imagery: a vegetation indices time series reconstruction methodology adapted to urban areas. Rem. Sens. 12, 639.
16. Huang, X., et al., 2019. The optimal threshold and vegetation index time series for retrieving phenology based on a modified dynamic threshold method. Rem. Sens. 11, 2725.
17. Jiang, D., et al., 2012. A simple semi-automatic approach for land cover classification from multispectral remote sensing imagery. PloS One 7, 9.
18. Jonsson, P., et al., 2018. A method for robust estimation of vegetation seasonality from Landsat and Sentinel-2 time series data. Rem. Sens. 10, 635.
19. Karkauskaite, P., et al., 2017. Evaluation of the Plant Phenology Index (PPI), NDVI and EVI for Start-of-Season trend analysis of the northern hemisphere boreal zone. Rem. Sens. 9, 485.
20. Li, H., et al., 2019. Incorporating the plant phenological trajectory into mangrove species mapping with dense time series Sentinel-2 imagery and the Google Earth Engine platform. Rem. Sens. 11, 2479.
21. Mohajane, M., et al., 2017. Mapping forest species in the central middle Atlas of Morocco (Azrou forest) through remote sensing techniques. International Journal of Geoinformation 6, 275.
22. Munafò, M., et al., 2018. Territorio-processi e trasformazioni in Italia. Istituto Superiore per la Protezione e la ricerca Ambientale (ISPRA), Sistema Nazionale per la Protezione dell'Ambiente (SNPA).
23. Munafò, M., et al., 2019. Consumo di suolo, dinamiche territoriali e servizi ecosistemici. Istituto Superiore per la Protezione e la ricerca Ambientale (ISPRA), Sistema Nazionale per la Protezione dell'Ambiente (SNPA).
24. Nink, S., et al., 2019. Using Landsat and Sentinel-2 data for the generation of continuously updated forest type information layers in a cross-border region. Rem. Sens. 11, 2337.
25. Oberti, R., 2009. L'analisi dell'immagine per l'identificazione di stress nelle colture. Atti Workshop Elementi Tecnici di Visione Artificiale e Applicazioni negli Agroecosistemi D. Monterotondo.
26. Sistema Informativo Agricolo Nazionale (SIAN). http://www.sian.it/inventarioforestale/.
27. Spinsi, A., et al., 2012. Indici vegetazionali da satellite per il monitoraggio in continuo del territorio. Italian Journal of Agrometeorology 3, 49.
28. Tang, S., et al., 2015. Trends of climatic sensitivities of vegetation phenology in semiarid and arid ecosystem in the US Great Basin during 1982-2011. Biogeosciences 12, 6985.
29. Tao, Z., Yan, H., Zhan, J., 2012. Economic valuation of forest ecosystem services in Heshui watershed using Contingent Valuation Method. Procedia Environmental Sciences 13, 2445.
30. Yang, Y., et al., 2018. The NDVI-CV method for mapping evergreen trees in complex urban areas using reconstructed Landsat 8 time-series data. Forests 10, 139.
31. Zeng, L., et al., 2019. A review of vegetation phenological metrics extraction using time-series, multispectral satellite data. Rem. Sens. 237, 2020.

Received 10 June 2020; Received in revised form 27 August 2020; Accepted 23 September 2020; Available online 28 September 2020.