A Statistical Approach to Detect Land Cover Changes in Mediterranean Ecosystems Using Multi-Temporal Landsat Data: The Case Study of Pianosa Island, Italy

The normalized difference vegetation index (NDVI) is commonly used to detect spatiotemporal changes of vegetation cover. This study modeled the spatiotemporal changes of land cover on Pianosa Island, Italy, in the period 1999–2015, using the multi-temporal Landsat images. Since the end of the 1990s, the natural vegetation has been re-colonizing an area of abandoned agricultural land and the island is undergoing a process of re-naturalization in harsh (drought and hot) environmental conditions. Hence, it is an ideal test site to monitor the effects of anthropogenic and climatic stressors on vegetation dynamics under Mediterranean climate. In this work, we proposed a new statistical approach based on a pixel-by-pixel analysis of multi-temporal Landsat images. Mean (µ) and standard deviation (σ) values of the NDVI images taken in 2015 were used for the determination of the pixel thresholds (µ ± 3σ). The evaluation of land cover change was carried out by comparing the µ value of a single NDVI pixel for 2015 with the same pixel of different years of the study period. The results indicate that surface reflectance (SR) Landsat images are more suitable in detecting the vegetation dynamics on the island than the top of atmosphere (TOA) ones and highlight an increasing trend of vegetation cover on Pianosa Island, mainly during the early seven years following the land abandonment in all the main land cover classes: abandoned crops and pastures, Mediterranean macchia, and woodland. However, the abandoned agricultural and pasture areas showed a higher increase in the vegetation cover and a shift in the shape of the normalized frequency distribution of the SR NDVI data during the study period, suggesting that a colonization process from other vegetation classes is occurring (i.e., Mediterranean macchia and trees are colonizing the abandoned land, partly replacing herbaceous species). Our data highlight that the statistical approach applied in this study is suitable for detecting vegetation cover changes associated with anthropogenic and climatic drivers in a typical Mediterranean environment and could be proposed as a new methodological approach in several other land monitoring studies.


Introduction
The Mediterranean Basin is considered as one of the most vulnerable regions to future climate change due to the foreseen increase in temperature to approximately double compared to the average global warming and decrease in precipitation, leading to potentially harmful impacts on its hydrological and biological resources [1]. Although representing only 0.8% of the world's ocean surface, the Mediterranean Sea is one of only 34 global biodiversity hotspots, testimony of its incredible importance for the biological wealth of our planet [2]. On the other hand, this region is currently subject to considerable spatial and temporal changes in land use and vegetation cover, with consequent relevant alterations in ecosystem composition and productivity [3]. As an example, the increasing abandonment of agricultural practices in the Mediterranean rural landscape over the past five decades has profoundly affected the vegetation dynamics in this region [4]. In addition, land use changes between agriculture and forest can significantly affect the carbon budget and the sequestration capacity of ecosystems. Hence, the detection of land use/cover changes is one of the major goals of future research and monitoring for evaluating the effects of climatic and anthropogenic pressures on primary productivity and carbon budget, biodiversity, and resilience of terrestrial ecosystems in the Mediterranean region. Here, the ecological consequences of global changes are highly complex and insular ecosystems can represent relevant natural study systems to investigate acclimation/adaptation mechanisms to climatic and anthropogenic pressures of taxa or communities. Such islands can represent potential "sentinels" for studying the effects and interactions of climate change and vegetation dynamics in the Mediterranean region [5]. In particular, the use of such "sentinels" might allow an improved understanding of the rate of environmental and vegetation changes under global warming and hence elaborate management planning and mitigate biodiversity loss.
Nowadays, optical remote sensing techniques are widely used for vegetation change detection [6] and as an effective temporal monitoring tool for protected areas [7], where it is very difficult to obtain time series of in situ information that is needed for ecological preservation planning [8]. The regional and global land cover products are made possible by the synoptic and periodic surveys of Earth Observation (EO) satellite missions. Civilian global land-cover EO started with the launch of Landsat 1 on 23 July 1972 [9]. The increased availability of EO data, especially from the Landsat archive, coupled with improved computing and storage capacity with novel image compositing approaches, have resulted in the availability of large-area and gap-free annual remote sensing data [10]. Several vegetation indices are commonly used to monitor the Earth's photosynthetic vegetation performance in support of phenological observations and biophysical interpretations [11]. In this framework, vegetation change detection is one of the most important topics in environmental monitoring [6]. Remote sensing data provide an estimate of the vegetation status and its spatiotemporal changes using the traditional normalized difference vegetation index (NDVI) [12,13]. Site-based and remote sensing methods are effective tools for mapping and monitoring vegetation and ecosystem health conditions from local to regional scale [14], with possible up-scaling integration to global level [15]. Although widely used, the interpretation of satellite images and derived thematic maps is strongly limited by the presence of clouds, shadows, and hazes, which all reduce the number of usable observations for a detailed analysis, both within and across years [16]. For this reason, there is a strong need to develop new methodological approaches to correctly interpret and homologate the temporal variations of satellite data using the existing databases. Since 2008, the Landsat data archive has become freely available to the scientific community [17]. This offers a unique data source to monitor the global land surface at a high spatial resolution (30 m). Hence, with more than 40 years of data, the Landsat archive is a highly valuable dataset for analyzing long-term land surface dynamics [18].
In a review paper, Galidaki et al. [19] discussed optical and active remote sensing at various spatial scales for estimating forest biomass in Mediterranean climate regions. They highlighted that literature on the Mediterranean biome is limited with a significant research gap in knowledge and, therefore, there is a need for reducing uncertainties of remote-sensing-based biomass estimates over these areas. Further development of image processing methods is needed for operational regional and national forest biomass assessment over Mediterranean areas.
In the study of land use change and vegetation cover, several threshold identification methods including fuzzy, neural networks, expert systems, change vector analysis, and Bayes theory have been applied (e.g., [10,[20][21][22][23][24][25][26][27][28][29]). Nevertheless, due to its simplicity, the image thresholding methods based on standard deviation are the most commonly used tools for monitoring vegetation changes [3]. The land cover detection methods typically employ medium spatial resolution Landsat imagery [26]. The most common approach is a bi-temporal change analysis through a comparison of pairs of images or characterizations [30] or multi-temporal change detection using more imagery.
The present study aimed to detect the spatiotemporal changes in vegetation cover in Pianosa Island (Tuscan Archipelago National Park, Italy, Figure 1a) during the period 1999-2015 using data from multi-temporal Landsat Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+), and Operational Land Imager (OLI). Pianosa is an ideal Mediterranean test site due to its recent history of land use change. Indeed, this island is a typical Mediterranean ecosystem where, after the abandonment of agricultural practices in previously cultivated land, natural vegetation is currently re-colonizing and the land is undergoing a process of re-naturalization or secondary succession [5,30]. Hence, this insular site may be considered a typical "sentinel" of the impact of climate change and anthropogenic disturbances on plant communities in a Mediterranean environment, also providing useful information on their interactions and feedbacks in terrestrial ecosystems [5,31]. Finally, a novel land cover change detection approach, based on a pixel-by-pixel statistical analysis of multi-temporal Landsat images, is developed and evaluated.

Study Area
Pianosa Island (Figure 1) has an area of 10.2 km 2 and a coastal perimeter of approximately 20 km. The highest elevation is 29 m above sea level (a.s.l.), with an average of about 18 m a.s.l. The flat morphology of Pianosa prevents condensation of moist air and thus the mean annual rainfall is very low. Moreover, due to the high permeability of the soils, the rain is usually quickly drained. Based on a historical meteorological dataset (1951-2009), the mean air temperature is 15.8°C and the mean annual rainfall is 496.6 mm, ranging between a minimum of 176.2 mm in 1999 and a maximum of 716.2 mm in 1984. Precipitation shows the typical seasonal pattern of the Mediterranean environment, with a maximum between October and December and a minimum between July and August. The meteorological data for the years under study were obtained from the regional meteorological service of Tuscany, Italy (Environmental Modeling and Monitoring Laboratory for Sustainable Development -LaMMA-Consortium). Values of daily precipitation, cumulative precipitation, and temperatures are reported in Figure S1 (Supplementary Materials). The island was a penal colony from the mid 19 th century to the late 1990s. The agricultural activities of the penal colony strongly affected the land use and the vegetation composition of Pianosa [31]. All agricultural and pastoral activities were then interrupted in the late 90s, following the transfer of the jail to another island. Hence, for more than 20 years, the anthropogenic pressure on land management has ceased or at least has been significantly reduced. The land previously used for crops and pastures is no longer managed and is now colonized by a plant community typical of the degraded soil of abandoned agricultural land in the Mediterranean region [30]. Dominant species are herbaceous therophytes. The Mediterranean shrubland vegetation covering other portions of the island is typical of calcareous soil and is dominated by Rosmarinus officinalis l., Cistus spp., and Juniperus phoenicea l. [32]. Additionally, there are also small woodlands with Pinus halepensis mill., Quercus ilex l., Arbutus unedo l., and a few Eucalyptus globulus labill trees. Therefore, three main types of ecosystems have been identified at Pianosa: (i) abandoned crops and pastures that are currently undergoing a process of renaturalization; (ii) Mediterranean macchia; and (iii) woodlands. The in situ land cover information is as of 2004 [33]. In total, the three types of ecosystems account for more than 90% of the entire island area (Figure 1b). The abandoned crops and pastures are the main ecosystems in terms of area (52%), followed by macchia (33%), and woodland (15%), while the Mediterranean macchia is the main ecosystem in terms of biomass (61%), followed by abandoned crops and pastures (26%) and woodland (13%). Among these ecosystems, Mediterranean macchia showed the highest aboveground net primary productivity per surface unit.

Acquisition of Long-Term Multi-Temporal Image Data
A flowchart of the methodological approach used in this study is shown in Figure S2 (Supplementary Materials). The first step was the identification of the time period for dataset collection in the study area. With this aim, a preliminary study on vegetation cover changes in the three main ecosystem types of Pianosa Island was carried out at a coarse spatial resolution of 250 m using the Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI data for the period 2000-2016 (MODIS product type Mod44/Myd44 [34]; Figures 1c,d). Figure 1d shows the seasonal variation of the interannual average values of the MODIS NDVI data during the 2000-2016 period. Each data point is therefore the average among the 17 years under study. Subsequently, a higher spatial resolution study (30 m) was performed by means of the Landsat top of atmosphere (TOA) and surface reflectance (SR) images. All the images were collected during the late spring-summer from day of year (DoY) 153 to DoY 281 (i.e., the seasonal period during which the MODIS NDVI data showed the lowest values and evidenced a clear distinction among the three studied ecosystems, Figure 1d). This period was also chosen because it is the less affected by cloud cover and simultaneously the most sensitive to vegetative activities in all ecosystems, as determined by our in situ studies.
Landsat TOA and SR orthorectified images (path 192-193, row 030-031) were acquired from the United States Geological Survey (USGS), Land Satellite Data Systems (LSDS) Science Research and Development (LSRD) [35]. To detect the spatiotemporal changes of land cover in the period 1999-2015, we acquired twenty-four data images: two from Landsat-5 TM, ten from Landsat-7 ETM+, and twelve from Landsat-8 OLI (see Table 1 for more details on the single data acquisition). All the data used for this study were standard Level 1 terrain-corrected (L1T) products [36,37].
Seasonal changes can induce strong differences in vegetated scenes due to variations in phenological stages (e.g., plant senescence), physiological activity (e.g., photosynthetic performance and stomatal conductance), and canopy architecture development. Diurnal and seasonal differences also affect sun azimuth and elevation [38]. For these reasons, all the data were acquired at the sun azimuth from 131.18° to 150.27° and sun elevation from 46.83° to 65.12° ( Table 1). The dominant weather conditions can affect atmospheric transmission and scattering [16]. Consistent differences in gross atmospheric conditions are often associated with seasonal changes. Hence, to reduce the influence of the atmospheric condition and the presence of clouds, the pixel quality assessment was carried out based on the CFMask (Function of Mask Algorithm) file and only the "clear" pixels (free of water, snow, clouds, cloud shadows, and cirrus) were used for the calculation of the seasonal variability of NDVI index. For this purpose, a decision tree based on the CFMask file was applied to each Landsat dataset used in this study (Tables 1 and S1) and only pixels equal to 0 (i.e., "clear" pixels) value were used for the subsequent analyses [39]. Moreover, due to a scan line corrector off (SLC-off) problem for the ETM+ data starting from the year 2003 [40], the decision tree was able to identify and remove pixels affected by the original data gaps in the primary SLC-off scene (i.e., pixels coded with number 255).

Calculation of the Normalized Difference Vegetation Index (NDVI)
The NDVI value for each TM, ETM+, and OLI image used in the present study was calculated as: where red and nir are the reflectance values acquired in the red and near-infrared regions, respectively, for both TOA and SR.

Image Pre-Processing
A pre-processing step of co-registration of the multi-temporal images was necessary. For this purpose, one Landsat 8 image (OLI 2015, DoY 180) was selected and used as a reference (i.e., master image) to register all the other images (i.e., slave images). A simple error-minimization routine was used to provide the relative geo-referencing, matching each image to the reference one by means of vertical and horizontal shifts. The subsequent processing steps were conducted on a region of interest (ROI) of 201  201 pixels.
All the data processing was carried out using the free software system for numerical computations Octave, under the terms of the GNU General Public License as published by the Free Software Foundation (https://www.gnu.org/software/octave/).

Landsat Data Continuity across Multi-Temporal NDVI Series
The Landsat data continuity across multi-temporal NDVI series was assured through specific quantitative spectral reflectance transformation functions [37] to compensate the relative spectral response (RSR) differences among the three Landsat sensors (TM, ETM+ and OLI; Figures 2 and S3).  The regression equation and determination coefficient (R 2 ) are shown. The black lines show a 1:1 relationship and the red and blue lines are the linear regression lines of best fit for TOA and SR data, respectively.

Seasonal Vegetation Change Detection
Seasonal changes of vegetation cover on the Island were detected by adopting a new statistical approach based on a pixel-by-pixel analysis of the multi-temporal time series. Mean (μ) and standard deviation () values of the OLI TOA and SR NDVI images (from 06/06/2015 to 20/09/2015; DoY: 157,  164, 173, 180, 189, 205, 212, 221, 228, 237, 244, 253) were used for the determination of the pixel thresholds (μ ± 3; Figure 3). For each pixel of the ROI, we calculated a mean (μij, with intervals i = 1, 201 and j = 1,201) and standard deviation (ij) on a minimum sample size of four DoY to a maximum of twelve DoY, where i and j identify the pixel inside the ROI. Only the pixels (ith,jth) following the normal distribution were considered for studying the vegetation change evolution during the time series. The Shapiro-Wilk test values were performed [41], resulting in 91.9% of the ROI following a normal distribution. It is well known that for random variables following a normal distribution, 99.7% of probability is within ±3 around μ [42].
The temporal evolution of the vegetation cover was studied through a pixel-by-pixel analysis, comparing the μ value of a single NDVI pixel for 2015 with the same pixel of a previous reference year of the selected time series (NDVI_year; where the reference year corresponds to 1999, 2003, 2007 and 2011, respectively). Each NDVI_year is a mean seasonal value of ROI for the reference year. For example, the NDVI_1999 was calculated as a mean value of three ETM+ images collected at DoY 192, 217, and 233. As previously explained, each pixel used for the calculation of NDVI_year (see Table 1) was clear and not affected by the SLC-off problem. The difference between the 2015 OLI NDVI and the NDVI of the reference years (NDVI_year) was evaluated for each pixel using a threshold value, defined as μ ± 3. In detail, because 2015 was the reference image to which we compared years prior, if the pixel showed a value of NDVI_year lower than μ − 3, it was classified as an area of increasing vegetation cover (gain in vegetation). Conversely, if the NDVI_year value was higher than μ + 3, the pixel was classified as an area of decreasing vegetation cover (loss in vegetation). Finally, if the NDVI_year was included within the μ ± 3 range, the pixel was classified as an area where the vegetation did not show any significant change during the time interval under study (stable in vegetation). .

Results and Discussion
The interannual average values of MODIS NDVI for the entire time series 2000-2016 (n = 17 years) in the three main types of ecosystems (woodland, Mediterranean macchia, and abandoned crops and pastures) showed a typical seasonal trend, with a decrease during the late spring-summer period (DoY 153-281, Figure 1d) due to the harsh (dry and hot) Mediterranean environmental conditions of the site [5]. During the summer period, all of the ecosystems under study showed a low phenological and environmental variability and reached the lowest NDVI values, most likely because of the limited water availability and high temperature leading to stress conditions (i.e., stomatal closure and possible molecular effects) and reduced photosynthetic activity [5]. The decrease in NDVI during this period was especially evident in the abandoned agricultural land because of the presence of herbaceous species (therophytes) that are more sensitive to environmental stresses compared to Mediterranean sclerophyll species of the macchia and woodland [5]. Therefore, during summer in the abandoned cropland, the photosynthetic activity of herbaceous species is decreasing more than that of Mediterranean shrubs, whose contribution to the spectral NDVI signal is, hence, increasing.
Hence, this index allows for distinguishing between the three types of ecosystems during the selected period, with the Mediterranean macchia showing the highest interannual average (±standard deviation) value of NDVI (0.61 ± 0.03), followed by woodland (0.57 ± 0.03) and abandoned crops and pastures (0.42 ± 0.04), in this order, in agreement with the differences observed in the above-ground net primary productivity per surface unit. Figure 2 shows the μ and  values of the SR NDVI OLI images for 2015. Figure 4 shows the multi-temporal SR NDVI images and the relative changes in vegetation cover during 1999-2015. The spatiotemporal changes in vegetation cover was estimated by a pixel-by-pixel comparison of the OLI NDVI 2015 with the NDVI images of the reference years (NDVI_years). In Figure 4, the pixels showing an increase (gain), a decrease (loss), or absence of changes (stable) in vegetation cover are marked in green, red, and blue, respectively.  (Figure 3) are shown. Similar but not equal results were also obtained for the TOA NDVI images (data not shown). Moreover, to validate our statistical approach, the ETM+ NDVI images for 2015 (Figure 4i) were compared with the OLI NDVI images relative to the same year (Figure 4l). The close agreement between the 2015 NDVI data from the two platforms (OLI and ETM+) supports the reliability of the implemented pre-processing procedures, assuring that the proposed approach is not affected by the different Landsat sensors and by possible noise-based (not vegetative) variations.
The quantitative results of this methodological approach are summarized in Tables 2 and 3   To investigate the spatiotemporal dynamic of vegetation cover on the island, a pixel-by-pixel analysis of the sequential change of NDVI was conducted during the study period, highlighting with different colors the last year during which a change was observed ( Figure 5). By means of these maps, it is possible to identify the areas subjected to natural vegetation expansion for use in the groundtruth validation. A comparison between the TOA and SR NDVI maps revealed a significant difference in the western area of Pianosa Island dominated by the Mediterranean macchia, where the TOA map highlighted a more consistent increase in vegetation cover than the SR NDVI map.
To further validate the appropriateness and reliability of our statistical threshold approach in detecting possible vegetation cover changes in the study area and to compare TOA and SR reflectance data, two sub-areas subject to vegetation expansion were identified and compared to the groundtruth observations obtained by the analysis of digital aerial orthophotos (2000-2013 by Regione Toscana-SITA: Cartoteca [43], Figures 6 and 7). These sub-areas were chosen because they are highly representative of the ecosystem types of the island and of the land-use changes occurring during the study period. In particular, in the first orthophoto image (Figure 6), a circular sub-area (marked in green) indicates a clear gain in vegetation during the 2000-2013 period, in agreement with the TOA and SR NDVI maps (Figures 5a,b). In the second orthophoto image (Figure 7), two circular sub-areas were selected: a sub-area (marked in blue) where the vegetation showed only limited changes during the 2000-2013 time interval, and another one (marked in green) showing an increase in vegetation cover. The comparison of these orthophotos with the TOA and SR NDVI images (Figures 5a,b), indicated a strong agreement with the SR NDVI, whereas the TOA NDVI was less related to the ground-truth observations. Hence, it can be concluded that in this study, the SR NDVI was more suitable to distinguish vegetation cover changes occurring in these two sub-areas, as suggested by Roy et al. [17]. In particular, these authors highlighted that SR reflectance is needed to derive consistent geophysical and biophysical products because the impact of atmospheric gases and aerosols on optical wavelength radiation is variable in space and time. For these reasons, SR reflectance can be considered a more suitable index in detecting vegetation dynamics than TOA. Table 4 shows the partitioning of the chronological gain/loss of vegetation for the different land cover classes (Mediterranean macchia, woodland, and abandoned crops and pastures) using the SR NDVI data. An increase in vegetation cover occurred following the year of crop abandonment in all three main ecosystem types on Pianosa Island. However, the abandoned crop and pasture areas, as expected, showed the highest vegetation gain. These variations in vegetation cover are mainly due to land use changes, although variability in precipitation and temperature regimes can also play a role in the observed inter-annual variations [44]. Accordingly, Maselli et al. [44] reported increasing fluxes of both carbon and water in the macchia ecosystem during the decade of 2001-2010, mainly depending on a similarly increasing spring rainfall. It is also worth noting that in 2015, a slight loss in vegetation cover compared to 2010 was recorded (Table 4), which was related to differences in thermo-pluviometric regimes between the two years ( Figure S1). In fact, 2010 showed much more frequent rainy events and distributed throughout the growing season compared to 2015. This likely prolonged the vegetative season, especially in the abandoned crop sites, leading to the increased NDVI values (Table 4). Indeed, these areas are mainly colonized by therophyte herbaceous species, with relatively low water-use efficiency, and therefore strongly affected by the precipitation regime during the vegetative season [5]. Finally, Figure 8 reports the temporal change in the frequency distribution of SR NDVI for the analyzed period in the three land cover classes. The results indicate a clear change in pixel distribution within each class from the years 1999-2007 to 2011, with a general upward shift of the mean SR NDVI values supporting the increasing trend of vegetation cover during the study period. However, both Mediterranean macchia and woodland classes did not show a substantial change in the shape of frequency distribution curves from 2011 to 2015, while in the abandoned crops and pastures, a progressive change in the frequency distribution shape occurred between 2007 and 2011. These data suggest that the increase in the mean SR NDVI in woodland and Mediterranean macchia is probably mainly due to a natural expansion of the preexisting vegetation during the study period, while the change in the shape of frequency distribution in the abandoned crops and pastures may be explained by the re-colonization process from other vegetation classes that is actually occurring. It is noteworthy that the shape and frequency distribution in the latter vegetation class in recent years have become more similar to that of woodland and/or Mediterranean macchia, which would indicate a recolonization of shrubs and bushes replacing, in part, herbaceous species.
The original vegetation of the island, dominated by a mixture of evergreen and deciduous trees, bushes, and grasslands [45], was greatly affected by the agricultural activities of the penal colony until its abandonment. Subsequently, a progressive secondary succession is currently undergoing after abandonment and the land previously used for agriculture is now covered by a plant community typical of degraded Mediterranean agricultural soil [5]. In particular, Scartazza et al. [5] showed an increase in the frequency of Cistus spp., a typical opportunistic semi-deciduous shrub species of Mediterranean environments that can play a relevant role in colonizing the degraded abandoned land of Pianosa. Furthermore, this secondary succession would favor the diffusion of Mediterranean species characterized by a higher water-use efficiency. Hence, it is possible that some species of the original Mediterranean macchia have an advantage with respect to others, leading to a novel plant association. Accordingly, the present data confirm that a natural expansion of vegetation is occurring on degraded and drought-prone soils of Pianosa. Abandonment of agriculture and pastoral practices in the Mediterranean rural landscape over the past five decades will have strong repercussions on the vegetation dynamics and biodiversity of these areas [4]. Due to its ecological interest and extreme drought environment, Pianosa Island could be a "sentinel" of climate change and anthropogenic pressure in the Mediterranean region, providing relevant information on biodiversity loss risks and vegetation dynamics as well as possible management strategies to preserve these extremely fragile ecosystems. Our statistical methodological approach can represent a useful tool for studying the vegetation dynamics in this and other unmanaged ecosystems. This methodology can take advantage of in-situ observations to study the ecological process of renaturalization and the species involved. In situ data are available in Pianosa for the years studied [5,33]. In future studies, the remote sensing technique using pixel-by-pixel analyses of multi-temporal Landsat data represents a potential tool for monitoring these processes in complex and variegated ecosystems such as those of Pianosa, in combination with accurate image pre-processing steps and ground-truth observations.

Conclusions
This study aimed to detect the spatiotemporal changes of vegetation cover of Pianosa Island (Italy) in the period 1999-2015 by a comparative analysis of TOA and SR NDVI Landsat maps. The results show that the SR NDVI data are more suitable to evaluate the spatiotemporal land cover changes, indicating an intense increasing trend in vegetation cover in the period 2011-2015. The total change in vegetation cover over the selected period was about 2.46 km 2 , more than one fifth of the overall island surface. Most of this change occurred during the early seven years following the agricultural abandonment and the sharp decrease in human presence, in the three main vegetation classes: Mediterranean macchia, woodland, and abandoned crops and pastures. However, in the latter vegetation class, the NDVI trend suggests that an increase in vegetation cover is currently in progress. If such an increasing trend continues, it is very likely that Mediterranean vegetation, as expected, will colonize the abandoned agricultural land and will undergo a process of renaturalization in the near future. This will have a significant impact on carbon budget and carbon sequestration, both in the above ground and in the soil organic C. Furthermore, a new methodological approach for vegetation change detection was applied by using a descriptive statistical analysis of the selected multi-temporal NDVI data series. This statistical approach employed a pixel-by-pixel image comparison analysis throughout the time series and was proven to be effective in monitoring the vegetation dynamics on the island. This insular ecosystem in the Tuscan Archipelago National Park has turned out to be a valuable test site, highly representative of several Mediterranean environments where the vegetation dynamic is strongly affected by anthropogenic land use and climatic drivers; hence, due to its recent history, environmental conditions, and geographical position, it could represent a "sentinel" for monitoring vegetation cover changes and biodiversity in the Mediterranean region. Accordingly, the method used in this study may prove useful in studying land cover changes in several other environmental conditions and world ecosystems as well as other aspects of land monitoring.
Supplementary Materials: The following are available online at www.mdpi.com/1999-4907/11/3/334/s1, Figure  S1: Meteorological data; Figure S2: Flowchart of the methodological approach; Figure S3: Linear relationships between spectral reflectance bands. Table S1: Number and percentage of pixel masked out in each datasets of Pianosa Island.