Urban and Peri-Urban Vegetation Monitoring Using Satellite MODIS NDVI Time Series, Singular Spectrum Analysis, and Fisher–Shannon Statistical Method

Abstract

The purpose of this work was to evaluate the potential of Singular Spectrum Analysis (SSA) and the Fisher–Shannon method to analyse NDVI MODIS time series and to capture and estimate inner vegetation anomalies in forest covers. In particular, the Fisher–Shannon method allows to calculate two quantities, the Fisher Information Measure (FIM) and the Shannon entropy power (SEP), which are used to characterise the complexity of a time series in terms of organisation/disorder. Pilot sites located both in urban (Milano, Torino, and Roma) and peri-urban areas (Appia Park, Castel Porziano, and Castel Volturno) were selected. Among the six sites, Roma, Castel Porziano, and Castel Volturno are affected by the parasite Toumeyella parvicornis. The time series was analysed using the products available in Google Earth Engine. To explore and characterise long-term vegetation dynamics, the time series was analysed using a multistep processing chain based on the (i) normalisation of the satellite time series, (ii) removal of seasonality and any other periodical cycles using SSA, (iii) analysis of the de-trended data using the Fisher–Shannon statistical method, and (iv) validation through comparison with independent data and ancillary information. Our findings point out to a clear discrimination between healthy and unhealthy sites, being the first (Milano, Torino, Appia) characterised by a larger FIM (lower SEP) and the second (Roma, Castel Porziano, Castel Volturno) by a lower FIM (larger SEP). The results of the investigations showed that the use of the SSA and Fisher–Shannon statistical methods coupled with the NDVI time series of the MODIS satellite made it possible to effectively identify and characterise subtle but physically significant signals veiled by seasonality and annual cycles.

1. Introduction

Preserving the natural capital (CN) requires constant analysis and systematic monitoring, especially for the most critical environments, as those relating to forests in urban and peri-urban areas, useful in countering the effects of climate change and in improving the quality of the environment and life [1,2,3,4]. Any form of vegetation cover (e.g., cropland, grassland, and forest) provides numerous services to the ecosystems and humanity, but urban forests are of primary importance because they do improve life quality (strengthen social connections as well as physical and mental health) reduce air and water pollution, as well as heating and cooling costs, increase real estate values, and mitigate climate change impacts [5,6,7,8,9,10].

Therefore, considering the extremely important role played by “urban forests”, their management, monitoring, and preservation is mandatory [11,12,13,14,15].

Recently, ESA (https://eo4ea-2022.esa.int/agenda/, accessed on 6 June 2023) national and international institutions and agencies, such as, for example, the European Environmental Agency (EEA), emphasised the importance of earth observation (EO) and remote sensing (RS) data for monitoring and accounting for natural capital to support development and management policies. Undoubtedly, since the 1980s, EO and RS technologies have been gaining special attention in the monitoring of vegetation changes and dynamics and in the detection of plant diseases and pests [16,17,18,19,20]. Invasive pests and alien plant bacteria are considered as major threats worldwide because they can induce serious plant diseases with devastating impacts on both natural ecosystems and agriculture production with huge environmental (loss of biodiversity) and economic damage [21]. Remote sensing methods can capture the degradation signs induced by many pathogen infections [22]; nevertheless, reliable methodological approaches have to be devised and assessed, not only to identify and map damage but also, and above all, to devise effective systems for early warning.

Changing trends can be quantified using NDVI time series, but this estimation can significantly differ on the base of the (1) analysed satellite dataset, (2) the spatiotemporal resolution of the used satellite products, and (3) the applied statistical method [23,24,25,26,27].

The analysis is based on NDVI because it is regarded as one of the most significant parameters useful for the assessment of the vegetation conditions at multiple temporal and spatial scales [28,29,30,31,32,33,34]. NDVI has been widely adopted as a proxy indicator of vegetation status and degradation induced by fires, drought, desertification, soil erosion, and soil salinisation, for assessing soil organic carbon, and more recently for the identification of pest attacks, whose detection and monitoring has, therefore, been representing a much more recent issue compared to other forms of land degradation [33,34,35]. Thus, assessing the feasibility and effectiveness of MODIS satellite NDVI time series to provide early warning for mitigating the impact of parasites and contrasting their spread is truly challenging.

In fact, pathogens and parasites are alien species that, due to a lack of predators, are capable of spreading very fast-provoking and devastating effects [36,37]. This phenomenon has led the European community to adopt laws and regulations to prevent, or at least mitigate, the risk of spreading non-native pests [38]. The recording of early warnings in vegetation changes can be crucial for timely intervention and it is, therefore, necessary to consider a long time series, free of seasonal trends, in order to observe even the smallest changes in trends. Several methods exist in the literature for conducting multi-temporal satellite analyses, such as: (i) temporal segmentation algorithms, such as Continuous Change Detection and Classification (CCDC), Vegetation Regeneration and Disturbance Estimates through Time (VERDeT), and LandTrend [39,40,41,42,43,44]; and (ii) trend analyses that analyse changes in pixel values over time [45,46,47,48,49].

In this paper, the approach developed by Telesca et al. [22], based on the use of the Fisher–Shannon statistical method to analyse and characterise the temporal behaviour of healthy and unhealthy vegetation (trend analyses), was adopted. This approach was also used by Telesca et al. [50] to identify and map the Xylella fastidiosa in the Puglia region and the Toumayella in the Campania and Lazio test sites using the evapotranspiration (ET) products available in Google Earth Engine (GEE).

The aim of this paper was to assess the potential of SSA and the Fisher–Shannon statistical method applied to time series of NDVI data extracted from MODIS satellite to observe anomalies in peri-urban and urban forests, which generally refer to all trees within a densely populated area, including trees in parks, on street-ways, and on private property” (see, for example, https://www.fs.usda.gov/ccrc/topics/urban-forests, accessed on 6 June 2023). Pilot sites located both in urban (Milano, Torino, and Roma) and peri-urban areas (Appia Park, Castel Porziano, and Castel Volturno) were selected in North, Middle and South Italy.

These case studies were chosen because they provide an excellent model for assessing the possibility of understanding changes that have taken place in vegetation or are taking place in vegetation, over time, in case studies with different characteristics that may be subject to (i) sudden changes in the amount of greenery, as in the case of urban greenery, and (ii) gradual changes over time related to plague and pest contamination. In both cases, these are phenomena of great interest to the global scientific community. The choice was made for (i) Turin (Torino) and Milan (Milano), because they are two cities with urban green areas, whose change is mainly due to urban works activities as well as reforestation; (ii) Rome (Roma), because it is a city with a large urban natural heritage composed of pine trees that are more than one hundred years old and are now affected by pests, as frequently reported in newspaper articles and popular websites [51]; (iii) Appia Antica Park, Castel Porziano, and Castel Volturno because they represent three peri-urban parks, which in the cases of Castel Porziano and Castel Volturno present the same parasitic problems as Roma (see Section 2).

For the purposes of the investigation, the NDVI MODIS time series was used by exploiting the products available in Google Earth Engine (GEE) that, in recent years, has been widely used in many disciplines of RS and EO techniques, including the environmental monitoring and risk assessment studies [52].

GEE is widely recognised for its robust computing capabilities and efficient handling of large-scale datasets, making it an invaluable tool in the realms of remote sensing and big data analysis [53,54,55,56]. Over the years, GEE has gained substantial popularity across various disciplines, evident in the significant growth in the number of scientific papers focused on utilising GEE. Researchers have leveraged GEE’s capabilities in diverse fields, including forest and vegetation analysis [9,57,58], land use and land cover studies [59,60], hydrology [61], and ecosystems and climate research [62,63], as well as cultural heritage studies [64,65,66]. GEE’s accessibility has facilitated the development and sharing of numerous free tools, which can be readily accessed on the GEE website.

Some of the most renowned satellite imageries are available on the GEE including the archives of MODIS, Sentinels, and Landsat series. which are complemented by other ancillary data, such as DEMs (Digital Elevation Models), meteorological data, shape-files, and Land Cover maps [30,53,54,67].

2. Materials and Methods

2.1. Dataset and Study Areas

Investigations were conducted on six study areas selected as a combination of urban areas with natural parks in Italy (Figure 1). The study areas were (i) Milano, (ii) Torino, (iii) Roma, (iv) Appia Antica Park, (v) Castel Porziano, and (vi) Castel Volturno. They are situated in the north (i–ii), in the centre (iii–v), and in southern Italy (vi).

Table 1 and Figure 2 show the characteristics of the study areas.

• Torino is the Italian city richest in public green areas, with 320 km of tree-lined avenues and about 50 parks. Public parks are in the historic centre, as well as in each district where there is at least one “lung” with lawns, areas equipped for children, and spaces for running and playing. Torino’s green heritage is not only extensive, but also varied, complex, sometimes precious, and delicate (http://www.comune.torino.it/verdepubblico/il-verde-a-torino/, accessed on 6 June 2023).

• Milano has only 17 square metres of greenery per inhabitant well below the national average of 31 square metres as for other cities in the Lombardy Region where the availability of public urban green is 28.6 square metres per inhabitant (with significant variations from city-to-city). For this reason, recently a plan of urban reforestation, named the Forestami project, has been launched and promoted by the Metropolitan City of Milano, the Municipality of Milano, the Lombardy Region, Parco Nord Milano, Parco Agricolo Sud Milano, ERSAF, and the Milano Community Foundation. Born from research by the Milano Polytechnic, thanks to the support of the Falck Foundation and FS Sistemi Urbani. The aim of Forestami is the planting of 3 million new trees in the Metropolitan City by 2030. Moreover, in the framework of the PNRR, additional efforts for the urban reforestation have been planned (https://www.cittametropolitana.mi.it/export/sites/default/ambiente/news/2022/ALLEGATO_1_Piano_forestazione.pdf, accessed on 6 June 2023).

• Roma has been selected because, herein, the green areas are often part of the archaeological site (including Villa Ada, Villa Pamphilj, Villa Chigi, Villa Torlonia, Villa Borghese, Villa Glori, Nemorense Park, Castel Sant’Angelo, the Pincio, Colle Oppio, Villa Sciarra, the Gardens of Piazza Vittorio). As a whole, the public green areas within the urban fabric are made up of urban parks, historic villas, public gardens, flowerbeds, and green areas for a total of 3932 hectares (https://www.ansa.it/sito/notizie/magazine/numeri/2017/03/02/roma-la-citta-piu-verde-deuropa_80a6cd96-5ae7-4b3f-9742-c039cb38dfc6.html, accessed on 6 June 2023). Moreover, Roma has been selected because, in the last years, the urban green areas have been strongly affected by pests and parasite attacks, mainly Toumeylla Parvicornis (https://www.comune.roma.it/web/it/informazione-di-servizio.page?contentId=IDS908082, accessed on 6 June 2023), which has strongly devasted the urban areas generally made up of Pinus trees, which are an integral part of the historical landscape of Roma as for the whole peninsula (https://www.italianostra.org/news-nazionali/salviamo-il-pinus-pinea-salviamo-il-paesaggio-italiano/, accessed on 6 June 2023). The Toumeyella parvicornis is an alien parasite (therefore, without a phytosanitary protocol to draw on for its treatment) that arrived in the capital in 2018, probably from the Campania. The extreme prolificity of Toumeyella determined the rapid spread of the insect causing the pine trees to die within a few years [38].

• The Appia Antica Regional Park, with its 4580 hectares, is the largest urban protected area in Europe. A green wedge runs from the city centre towards the Castelli Romani. Here, history, archeology, and nature blend in a landscape and environment of exceptional interest that offers even the most intransigent visitors unexpected surprises (https://www.parcoappiaantica.it/wp-content/uploads/2021/04/13.LaFlora.pdf, accessed on 16 May 2023). At the centre is the Regina Viarum, an open-air museum that winds for over 16 kilometres, and is always accessible. On its sides extend the historic agricultural estates, interrupted by the spectacular remains of imperial villas and aqueducts.

• Castel Porziano is a Presidential Estate. It is 25 km approx. from Roma and covers an area of 69 km² approx., including some historic hunting estates such as “Trafusa, Trafusina, Riserve Nuove and Capocotta”. The land cover is: forest (2300 hectares), Mediterranean scrub (500 hectares), ilex grove (261 hectares), cork oak forest (460 hectares), and stone pine forests (750 hectares). The latter have been artificially grafted and are intended to protect the coast from erosion while simultaneously sheltering inland areas from sea winds. Although artificially introduced, stone pine is now a fundamental element of the Italian landscape, introduced by the ancient Romans.

• The Castel Volturno nature reserve covers about 268 hectares on the coast of the municipality of the same name. It occupies an area between the mouth of the Regi Lagni and the mouth of Lago Patria (Figure 1). Within it are the protected areas ZSC IT8010021 “Pineta di Patria” and the Regional Nature Reserve “Foce Volturno-Costa di Licola”, consisting of pine trees. Castel Volturno is a typical example of Mediterranean maquis, where the tree layer is mainly represented by umbrella pine (Pinus pinea), maritime pine (Pinus pinaster), holm oak (Quercus ilex); the shrubby one from Juniper (Juniperus sp.), Phillyrea (Phillyrea angustifoglia), Myrtle (Myrtus comunis), and lianas such as Asparagus (Asparago acutifolius). In the last years Castel Volturno has been severely affected by Toumeyella parvicornis [38].

2.2. Methodological Approach

The designed procedure has been based on two steps, which is to account for both the variability of the seasonal variations in the signal and to identify subtle multi-year trends or changes: (i) the first step is based on the removal of the annual/seasonal/sub-seasonal (or other potential) periodical cycle; (ii) the second step is based on the processing of de-trended signals, in order to identify the presence of subtle (but significant) changes. MODIS NDVI products were used because data from more than 20 years are available (Table 2).

For the purpose of this study, the MODIS product, referred to as MOD13Q1.006, Terra vegetation indices were used. The MODIS NDVI vegetation index is derived from the normalised difference of the atmospherically corrected reflectance of near-infrared (858 nm) and red bands (645 nm) [28], according to the Formula (1):

(1)$NDVI=\frac{({\rho }_{NIR}-{\rho }_{red})}{({\rho }_{NIR}+{\rho }_{red})}$

where ρ_NIR, ρ_red are the surface reflectance over the near-infrared (NIR) and red bands of MODIS. The NDVI is well correlated with vegetation cover, vegetation canopy, vegetation dynamics, biomass, and leaf area index, and often considered as the vegetation proxy [35,67,71,72,73,74].

The choice to use MODIS data was virtually a mandatory one. In fact, the MOD13Q1.006 data, as shown in Table 2, has:

(i) Freely available data in GEE;

(ii) A spatial resolution of 250 m/pixel attested to be useful in this type of analysis in other studies [50,73];

(iii) A temporal resolution or revisit time of 16 days, the result of which is an interpolation of multi-day acquisitions for which the best pixels are then selected as reported in [75,76,77];

(iv) A temporal coverage of more than twenty years.

In particular, points (ii)–(iv) were relevant in the choice of sensor, and in preferring MODIS to higher resolution sensors such as Landsat. Since the analyses conducted need a large and continuous (without null values) time series of observations in order to work at their best.

The spatially averaged NDVI series were collected for the study regions (Figure 1) using the GEE for the period of 2001–2020.

As a whole, the data processing can be summarised as follows (Figure 3):

• Computation of spatial average of NDVI time series for each investigated site;

• Normalisation of the satellite time series. The data are processed to have zero mean and unitary standard deviation;

• Application of SSA after normalisation of NDVI. First, the phase shift value M is selected. Taking into account the sampling time of the series (16 days), to identify at least the annual cycle, M was set equal to 24; this value is also in agreement with the criterion of Khan and Poskitt [78], considering that the size of each series is 503. The decomposition of a time series into independent components can be performed by using several methods; among these the Singular Spectrum Analysis (SSA) [79], based on phase-lagged copies of the series, is well known for its efficiency, even in case of noisy and short time series. The independent components that are derived by the application of the SSA behave generally as slowly changing trends, cyclic components and purely random noise [80]. This step will allow to identify trend and anomalies;

• Characterisation of the informational properties of NDVI time series through the Fisher Information Measure (FIM) and the Shannon entropy (SE), which in theory of are applied to measure the local and global smoothness of the distribution of a series, respectively. FIM and SE can be used to describe complex, nonstationary time series which are described in order and organisation (FIM) [81], and uncertainty or disorder (SE) [82];

• Validation step conducted by comparing the results with independent data sets and ancillary information. All of the test site results from the Fisher–Shannon statistical method were compared with: (i) imagery available in Google Earth; (ii) forest watch (https://www.globalforestwatch.org/, accessed on 17 May 2023), for which products are freely available on a global scale and well-documented. The forest watch products are obtained from Landast at 30 m and provide analyses and quantification on forest cover loss and forest gain, herein the TM 2000–2020 time series were processed using the Hansen et al. [83,84] Global Forest Change dataset; (iii) ancillary information including results and outputs from other projects; (iv) field survey.

SSA and FS Method

The time series decomposition into independent components can be performed by using several methods; among these the Singular Spectrum Analysis (SSA) [79], based on phase-lagged copies of the series, is well known for its efficiency, even in the case of noisy and short time series. The independent components that are derived by the application of the SSA behave generally as slowly changing trends, cyclic components and purely random noise [80].

Given the time series y_i (i = 1, …, N) its Toeplitz M-lagged correlation matrix is given by:

(2)$C_j=\frac{1}{N-j} \sum _{i=1}^{N-j}{y}_i{y}_{i+j}, 0 \le j \le M$

From the eigenvalues λ_k, sorted in decreasing order, and eigenvectors E_k,j, where j and k vary from 1 to M, of the Toeplitz M-lagged correlation matrix, the k-th principal component a_ik are calculated:

(3)$a_{ik}=\sum _{j=1}^M{y}_{i+j}{E}_{jk}, 0 \le i \le N-M$

The k-th is obtained as follow:

(4)$R_k=\frac{1}{M}\sum _{j=1}^M{a}_{i-j,k}{E}_{jk}, M \le i \le N-M+1$

The decreasing order of the eigenvalues λ_k corresponds to a decreasing order of the fraction of the total variance of the original series explicated by the reconstructed component R_k [85]. The only requirement of the SSA is the proper selection of the lag M. The maximum M = (log N)^c, 1.5 ≤ c ≤ 2.5 was calculated in [78].

The MDL (minimum description length) criterion is given by:

(5)$MDL\left(k\right)=-\mathrm{l}\mathrm{o}\mathrm{g}\left(\frac{\prod _{i=k+1}^p{\lambda }_i^{\frac{1}{p-k}}}{\frac{1}{p-k}\sum _{i=k+1}^p{\lambda }_i}\right){}_{ }{}^{\left(p-k\right)N}+\frac{1}{2}k\left(2p-k\right)logN$

where λ_k are the eigenvalues, p is the number of eigenvalues, identical to M, and N is the length of the original series. This criterion is used to separate the series into a trend and a detrended series. Calculating the value of k ∈ {0, 1, 2, …, p – 1} for which the MDL is minimised (k_min), the trend is given by summing up all of the reconstructed components up to k_min, while the detrended series is obtained by subtracting the trends from the original values.

Using the Fisher–Shannon method, it is possible to study the information content of a time series. In particular this method is based on the joint calculation of the Fisher Information Measure (FIM) and the Shannon entropy (SE), used to measure the local and global smoothness of the distribution of a series, respectively. FIM and SE can be used to describe complex and non-stationary time series as (i) order and organisation and (ii) uncertainty and disorder, respectively [73,81]. The FIM and SE are defined as follows:

(6)$FIM=\underset{-\infty }{\overset{+\infty }{\int }}{\left(\frac{\partial }{\partial x}f(x)\right)}^2\frac{dx}{f\left(x\right)}$

(7)$SE=\underset{-\infty }{\overset{+\infty }{\int }}{f}_x\left(x\right)log{f}_x\left(x\right)dx$

where f(x) is the distribution of the series x. Alternatively to SE, generally used is the Shannon entropy power N_X (8):

(8)$N_x=\frac{1}{2\pi e}{e}^{2H_x}$

to avoid dealing with negative quantities. FIM and N_X depend on each other by the isoperimetric inequality FIM·N_X ≥ D [50], where D is the dimension of the space (D = 1 for time series).

FIM and N_X depend on f(x), which has to be accurately estimated for obtaining reliable values. For the calculation of FIM and N_X, the kernel-based approach was employed [86] as shown in the following formula:

(9)$\widehat{f_M}(\mathrm{x})=\frac{1}{Mb}\sum _{i=1}^{M}K\left(\frac{x-x_i}{b}\right)$

where M and b denote the length of the series and the bandwidth respectively, while K(u) is the kernel that is a continuous, symmetric, and non-negative function, satisfying the two following constrains:

(10)$k\left(u\right)\ge 0 and\underset{-\infty }{\overset{+\infty }{\int }}k\left(u\right)du=1$

f(x) is estimated by means of an optimised integrated procedure [87] with a Gaussian kernel:

(11)$\widehat{f_M}(\mathrm{x})=\frac{1}{M\sqrt{{2\pi b}^2}}\sum _{i=1}^M{e}^{-\frac{{\left(x-x_i\right)}^2}{{2b}^2}}$

Due to the isoperimetric inequality, the time series can be represented in the so-called Fisher–Shannon (FS) information plane [88], whose x-axis and y-axis are N_X and FIM, respectively. The curve IN_X = 1 separates the FS information plane into two half-spaces, and a time series can be represented by a point located just in the half-space IN_X > 1.

3. Results and Discussion

NDVI values were extracted as the average value for the areas of interest (Figure 1) using the GEE for the period of 2001–2020.

Figure 4 shows, as an example, the application of the SSA to the normalised NDVI time series of Torino. Figure 4a shows the results of the eigenvalues of the SSA, where each equals a reconstructed component and represents the fraction of the total variance of the original series explained by that component. Figure 4b represents all reconstructed components, which have a signal ranging from oscillatory with amplitude modulation to apparently noisy.

Each time series was separated into a trend series and a detrending series by applying the MDL criterion (Figure 5). The minimum of the MDL curve is at k_min = 10. This value was used to extract the trend and obtain the detrended series.

The trend is the sum of the first 10 reconstructed components, while the detrended series is obtained by subtracting the trend from the normalised series (Figure 6).

The trend has an oscillatory behaviour that can be explained by the seasonal cyclicality of the series, probably due to climate variations. The detrended series would represent the internal temporal dynamics of the series unaffected by exogenous phenomena that, instead, are explained by the trend. Table 3 shows, for all investigated sites, the value of the minimum MDL criterion.

The Fisher–Shannon method was applied to each site to explore the information properties of the de-trended series, as the aim was to characterise the complexity of the temporal dynamics of vegetation within the study sites. Using the Fisher–Shannon information plane (FSIP), which has the Shannon entropy power (SEP) as the x-coordinate and the FIM as the y-coordinate, an overall representation of the six studied sites was obtained (Figure 7).

Figure 8 shows the quantification of forest loss and gain as available in Forest watch for all of the investigated areas (Torino, Milano, Roma, Appia park, Castel Porziano, and Castel Volturno) and re-computed by the author using the same algorithm.

The results obtained from the Fisher–Shannon statistical method applied to the MODIS times series were validated by comparison with the results from the TM time series 2000–2020 processed using the Hansen et al. [83,84] Global Forest Change dataset.

These comparisons confirmed the reliability of the Fisher–Shannon statistical method and its capability to assess and evaluate the net loss and gain for each study area. This is well-visible for the Torino area, where, up to 2020, there is a full compensation between forest loss and gain.

In particular, for Roma, Castel Porziano, Castel Volturno, and Appia park, the Fisher–Shannon results were also compared with the results available from other projects/ investigations and field surveys. Concerning Roma, analyses based on Sentinel 2 (obtained by S23 project, https://www.digimat.it/case_study/space-to-tree/, accessed on 6 June 2023) confirmed the ongoing degradation of urban and peri-urban green areas. Detailed investigations based on satellite, close range survey (drone equipped with visible and Thermocamera), and ancillary information were made on the Pinus trees in the Coloseum, as shown in Figure 9. In particular, Figure 9 shows the location of the investigated healthy and unhealthy Pinus trees and clearly highlights the different spectral behaviour of the healthy and unhealthy trees.

For the Appia Park, a visual analysis and comparison of Figure 2d and Figure 8d clearly highlighted two important aspects: (i) there is a prevailing agricultural use, and from the updated Corine Land Cover maps (https://land.copernicus.eu/pan-european/corine-land-cover, accessed on 6 June 2023), no significant changes occurred (as expected, also considering that it is a park); and (ii) the proposed method is capable of suitably removing the NDVI significant seasonal variations, while identifying subtle multi-year trends and changes, that in the case of Appia Park can be mainly related to the anthropogenic activities linked to the prevailing agricultural use of this site.

In the FSIP, Torino and Milano are the sites characterised by higher FIM and lower SEP; a factor indicating that these two sites are characterised by a high level of order and organisation, compared to the other three sites Roma, Castel Porziano and Castel Volturno. The data freely available from forest watch (https://www.globalforestwatch.org/, accessed on 6 June 2023) highlighted that in the analysed time period (2001 to 2020), Torino experienced a net gain of 277 hectares in tree cover gain (315 ha forest gain equal to 2.4% of its total extent, and 38 ha in forest loss equal to 1.7% of its total extent), and in the same time period Milano experienced a net gain of 654 hectares in tree cover gain (664 ha forest gain equal to 3.7% of its total extent, and 10 ha in forest loss equal to 1.7% of its total extent). This clearly fits well with the results obtained from the MODIS NDVI time series analyses based on Singular Spectrum Analysis and the Fisher–Shannon statistical method, from which Torino is relatively more stable, compared to Milano. It is important to remember that, for each study area, the behaviour, and therefore, according to forest watch analyses, the net loss and gain, are related to the averaged time series of NDVIs extracted for each polygon representing each study area (Figure 1). Therefore, the final information is related to a spatial average; namely, an aggregate view (see Section 3) of each site.

As a whole, the proposed approach enabled the identification of subtle signals totally veiled by the seasonality and annual cycle behaviour and could be a tool for early diagnosis of these phenomena.

The NDVI, as any satellite vegetation index, is a superposition of different signal behaviours (seasonal, gradual, and/or abrupt) that are, in part, induced by external driving meteo-climatic mechanisms, and, in part, due to the inner variability of the vegetation time dynamics. To focus on the last, it is necessary to remove the seasonalities (or any periodical cycles). To take into account the great variability exhibited by seasonal variations in NDVI, while identifying subtle multi-year trends and changes, a procedure based on two phases has been devised: (i) the first step is based on SSA applied to the NDVI time series, to detect and remove the annual/seasonal/sub-seasonal (or other potential) periodical cycles; (ii) the second step is based on the processing of the detrended series by using the Fisher–Shannon statistical method to characterise the complexity of the series in terms of their informational properties.

Among the investigated sites, Castel Volturno, Castel Porziano, and Roma were characterised by an inverse proportion between SE and FIM. The pattern shown by these three sites indicates a low level of order and organisation of the NDVI, which could be correlated with ongoing degradation trends. In fact, independent analyses and field surveys highlighted that all these three sites are strongly affected by a vegetation degradation process, due to the parasitic beetle Toumeyella parvicornis that has attacked pine trees in these areas in the last years, dramatically damaging the Pinus trees, as confirmed by the multiple reports in Italian newspapers, by the funds allocated to buffer the problem, and by scientific publications on the subject [51,89].

The capability of this statistical method to characterise the informational status of the NDVI series, and to correlate it to the health status of vegetation, is reinforced by the Fisher–Shannon analysis performed on the other three sites (Torino, Milano, and the Appia park) that show higher FIM and lower Shannon entropy power, in comparison with the first three sites. Torino, Milano, and Appia Park are characterised by a more stable behaviour characterised (up to 2020) by a full compensation of the forest loss and gain, and, thus, by the absence of “pathologic” trends of vegetation deterioration.

Data freely available from forest watch (https://www.globalforestwatch.org/, accessed on 6 June 2023) highlighted that in the analysed time period (2001 to 2020), Milano experienced a net change of 2% compared to its total surface area, and in the same time period Torino experienced a net change of 0.7% compared to its total surface area. This clearly fits well with the results obtained from the MODIS NDVI time series process based on Singular Spectrum Analysis and the Fisher–Shannon statistical method, from which Torino is relatively more stable compared to Milano.

As a whole, the investigations herein conducted highlighted the capability of the method to:

• (i). Identify significant ongoing degradation trends regardless the land cover and land use as evident by the fact that the six investigated pilot areas are characterised by different land use (Figure 2). This is thanks to the first step of data processing which enabled the identification and removal of seasonal and periodical cycles.

• (ii). Provide, for each study area, a net evaluation of the forest loss and gain. This is due to the fact that the area-averaged time series of NDVIs is extracted and processed for each polygon (Figure 1) representing each study area (using GEE).

Therefore, the effective removal of both phenological and seasonal/sub-seasonal/cyclical behaviours enabled the identification of the anomalies linked to degradation, which is the prevailing phenomenon in the case of Castel Volturno, Castel Porziano, and Roma. Whereas, Torino, Milano, and the Appia park show a more stable behaviour characterised (up to 2020) by a full compensation of the forest loss and gain.

4. Conclusions

In this study, a methodological approach, based on the combination of a decompositional method (SSA) and an information method (Fisher–Shannon analysis) for monitoring urban and peri-urban forests using MODIS NDVI satellite products, selected for their long global coverage, was proposed and discussed.

The 2001–2020 temporal variation of the MODIS NDVI for six study areas located from the northern to the southern part of the Italian peninsula was analysed. The pilot sites are urban parks in Milano, Torino, Roma, Appia Ancient park, Castel Porziano, and Castel Volturno. The study areas were specifically selected where there is a combination of urban areas with natural parks. For each site, the focus was on the internal temporal variability of the vegetation given by the detrend series that is not influenced by seasonal or external factors.

The results clearly pointed out by the approach based on the SSA and Fisher Shannon statistical method suitably can (i) provide, for each study area, a net evaluation of the forest loss and gain; and (ii) identify significant ongoing prevailing degradation trends regardless of the land cover and land use.

The results could contribute to the definition of methodologies capable of developing metric and diagnostic tools for monitoring vegetation, in order to assess its status using a synthetic indicator. The importance of this objective laid in the fact that NDVI is one of the most important vegetation key parameters which is extremely important in climate change investigations (hydrological modelling, carbon balance), environmental monitoring (pollution detection, nutrient flows), risk estimation (wild-fire, desert locust), land management, agricultural practices, food security, monitoring of vegetation stress, and disease.

The proposed methodological approach can provide a contribution to the creation of useful tools for:

(i) Early detection of deterioration trends;

(ii) Defining metrics for estimating the recovery/restoration capacity and the effectiveness of the contrast measures adopted to mitigate the phenomena of degradation;

(iii) Creating operational tools for multi-scale, multi-sensor, multi-temporal monitoring of bio-physical parameters relating to the state of vegetation;

(iv) Contributing to the definition of effective indicators for monitoring the natural capital (CN).

Author Contributions

Conceptualization, L.T., R.L. and M.L.; software, L.T., M.L. and N.A.; validation, R.L., N.M. and N.A., investigation, L.T., R.L., M.L., N.M., N.A., A.L., M.P., A.A. and G.C.; resources, R.L.; data curation, A.L., R.L., M.L., N.M., N.A., A.L., M.P., A.A. and G.C., original draft preparation, R.L., L.T., N.A. and N.M.; writing—review and editing, L.T., N.A. and R.L.; visualization, L.T., R.L., M.L., N.M., N.A., A.L., M.P., A.A. and G.C.; project administration, L.T., R.L., M.L., N.M., N.A., A.L., M.P., A.A. and G.C.; funding acquisition, R.L. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Study areas with the combination of urban area and natural parks in northern, central, and southern Italy.

Figure 2. Corine land cover of the study areas: (a) Torino; (b) Milano; (c) Roma; (d) Appia Park; (e) Castel Porziano; (f) Castel Volturno.

Figure 3. Flowchart.

Figure 4. SSA to Torino area: (a) eigenvalue spectrum; (b) reconstructed components.

Figure 5. MDL versus k for Torino NDVI time series. The minimum of the MDL (at k = 10) is used to separate the trend and the detrended series.

Figure 6. Torino detrended NDVI time series.

Figure 7. Fisher–Shannon information plane for all the investigated sites.

Figure 8. Data extracted from Google Earth Engine using the Hansen et al. [83] and Potapov et al. [84] Global Forest Change dataset for the years 2000 to 2021: (a) Torino; (b) Milano; (c) Roma; (d) Appia Park; (e) Castel Porziano; (f) Castel Volturno.

Figure 9. Location and temporal spectral behaviour of the Pinus trees investigated inside the Colosseum Park; in particular: (top) show the location of the sample trees; (bottom) the different temporal spectral behaviour of the healthy and unhealthy Pinus trees investigated.

Table 1

Study areas: area, annual precipitation (average), annual temperature (average), and climate system.

Study Area	Area (Km2)	Annual Precipitation (mm)	Annual Mean Temp. (°C)	Climate System by the Köppen-Geiger [68,69]
Torino	130	1002	12	Csa
Milano	190	1162	13	Csa
Roma	1187	878	15.8	Csa
Appia Park	47	878	15.8	Csa
Castel Porziano	69	878	15.8	Csa
Castel Volturno	12	1078	15.5	Csa

Table 2

Details of the MODIS MOD13Q1.006 data used, from [70].

Dataset	MOD13Q1.006
Platform	Terra
Sensor	Moderate Resolution Imaging Spectroradiometer (MODIS)
Spatial Resolution	250 m
Temporal Resolution	16 days
Coverage	Global
Bands	Band 1: Red (620–670 nm); Band 2: NIR (841–876 nm); Band 3: Blue (459–479 nm); Band 4: Green (545–565 nm); Band 5: MIR (1230–1250 nm); Band 6: MIR (1628–1652 nm); Band 7: SWIR (2105–2155 nm)
Data Period	2000–Present

Table 3

Characteristics of the study areas.

	Roma	Torino	Castel Volturno	Castel Porziano	Appia Park	Milano
kmin	8	10	5	9	10	8

References

1. Chami, R.; Cosimano, T.; Fullenkamp, C.; Nieburg, D. Toward a Nature-Based Economy. Front. Clim.; 2022; 4, 855803. [DOI: https://dx.doi.org/10.3389/fclim.2022.855803]

2. Monfreda, C.; Wackernagel, M.; Deumling, D. Establishing national natural capital accounts based on detailed Ecological Footprint and biological capacity assessments. Land Use Policy; 2004; 21, pp. 231-246. [DOI: https://dx.doi.org/10.1016/j.landusepol.2003.10.009]

3. Jabbar, M.; Yusoff, M.M.; Shafie, A. Assessing the role of urban green spaces for human well-being: A systematic review. GeoJournal; 2022; 87, pp. 4405-4423. [DOI: https://dx.doi.org/10.1007/s10708-021-10474-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34305268]

4. Pimentel, D.; Wilson, C.; McCullum, C.; Huang, R.; Dwen, P.; Flack, J.; Tran, Q.; Saltman, T.; Cliff, B. Economic and Environmental Benefits of Biodiversity. BioScience; 1997; 47, pp. 747-757. [DOI: https://dx.doi.org/10.2307/1313097]

5. Protecting Forests on the Front Line of the Climate-Change Battle|Research and Innovation. Available online: https://ec.europa.eu/research-and-innovation/en/horizon-magazine/protecting-forests-front-line-climate-change-battle (accessed on 27 June 2023).

6. Mercer, C.; Comeau, V.M.; Daniels, L.D.; Carrer, M. Contrasting Impacts of Climate Warming on Coastal Old-Growth Tree Species Reveal an Early Warning of Forest Decline. Front. For. Glob. Chang.; 2022; 4, 775301. [DOI: https://dx.doi.org/10.3389/ffgc.2021.775301]

7. Verkerk, P.J.; Delacote, P.; Hurmekoski, E.; Kunttu, J.; Matthews, R.; Mäkipää, R.; Mosley, F.; Perugini, L.; Reyer, C.P.O.; Roe, S. et al. Forest-Based Climate Change Mitigation and Adaptation in Europe; From Science to Policy European Forest Institute: Joensuu, Finland, 2022.

8. Nunes, L.J.; Meireles, C.I.; Gomes, C.J.P.; Ribeiro, N.M.A. Forest Contribution to Climate Change Mitigation: Management Oriented to Carbon Capture and Storage. Climate; 2020; 8, 21. [DOI: https://dx.doi.org/10.3390/cli8020021]

9. Hansen, A.J.; Neilson, R.P.; Dale, V.H.; Flather, C.H.; Iverson, L.R.; Currie, D.J.; Shafer, S.; Cook, R.; Bartlein, P.J. Global Change in Forests: Responses of Species, Communities, and Biomes. BioScience; 2001; 51, pp. 765-779. [DOI: https://dx.doi.org/10.1641/0006-3568(2001)051[0765:GCIFRO]2.0.CO;2]

10. Nowak, D.J.; Hirabayashi, S.; Bodine, A.; Greenfield, E. Tree and forest effects on air quality and human health in the United States. Environ. Pollut.; 2014; 193, pp. 119-129. [DOI: https://dx.doi.org/10.1016/j.envpol.2014.05.028]

11. Livesley, S.J.; McPherson, E.G.; Calfapietra, C. The Urban Forest and Ecosystem Services: Impacts on Urban Water, Heat, and Pollution Cycles at the Tree, Street, and City Scale. J. Environ. Qual.; 2016; 45, pp. 119-124. [DOI: https://dx.doi.org/10.2134/jeq2015.11.0567]

12. Irga, P.J.; Burchett, M.D.; Torpy, F.R. Does urban forestry have a quantitative effect on ambient air quality in an urban environment?. Atmos. Environ.; 2015; 120, pp. 173-181. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2015.08.050]

13. Kauppi, P.E.; Stål, G.; Arnesson-Ceder, L.; Sramek, I.H.; Hoen, H.F.; Svensson, A.; Wernick, I.K.; Högberg, P.; Lundmark, T.; Nordin, A. Managing existing forests can mitigate climate change. For. Ecol. Manag.; 2022; 513, 120186. [DOI: https://dx.doi.org/10.1016/j.foreco.2022.120186]

14. Kondo, M.C.; Fluehr, J.M.; McKeon, T.; Branas, C.C. Urban Green Space and Its Impact on Human Health. Int. J. Environ. Res. Public Health; 2018; 15, 445. [DOI: https://dx.doi.org/10.3390/ijerph15030445]

15. Lee, A.C.K.; Maheswaran, R. The health benefits of urban green spaces: A review of the evidence. J. Public Health; 2011; 33, pp. 212-222. [DOI: https://dx.doi.org/10.1093/pubmed/fdq068]

16. Zhang, J.; Huang, Y.; Pu, R.; Gonzalez-Moreno, P.; Yuan, L.; Wu, K.; Huang, W. Monitoring plant diseases and pests through remote sensing technology: A review. Comput. Electron. Agric.; 2019; 165, 104943. [DOI: https://dx.doi.org/10.1016/j.compag.2019.104943]

17. Mahlein, A.-K. Plant Disease Detection by Imaging Sensors—Parallels and Specific Demands for Precision Agriculture and Plant Phenotyping. Plant Dis.; 2016; 100, pp. 241-251. [DOI: https://dx.doi.org/10.1094/PDIS-03-15-0340-FE]

18. Riley, J.R. Remote Sensing in Entomology. Annu. Rev. Èntomol.; 1989; 34, pp. 247-271. [DOI: https://dx.doi.org/10.1146/annurev.en.34.010189.001335]

19. Nilsson, H. Remote Sensing and Image Analysis in Plant Pathology. Annu. Rev. Phytopathol.; 1995; 33, pp. 489-528. [DOI: https://dx.doi.org/10.1146/annurev.py.33.090195.002421]

20. Zhang, J.; Huang, Y.; Yuan, L.; Yang, G.; Chen, L.; Zhao, C. Using satellite multispectral imagery for damage mapping of armyworm (Spodoptera frugiperda) in maize at a regional scale. Pest Manag. Sci.; 2016; 72, pp. 335-348. [DOI: https://dx.doi.org/10.1002/ps.4003]

21. Pautasso, M.; Döring, T.F.; Garbelotto, M.; Pellis, L.; Jeger, M.J. Impacts of climate change on plant diseases—Opinions and trends. Eur. J. Plant Pathol.; 2012; 133, pp. 295-313. [DOI: https://dx.doi.org/10.1007/s10658-012-9936-1]

22. Telesca, L.; Aromando, A.; Faridani, F.; Lovallo, M.; Cardettini, G.; Abate, N.; Papitto, G.; Lasaponara, R. Exploring Long-Term Anomalies in the Vegetation Cover of Peri-Urban Parks Using the Fisher-Shannon Method. Entropy; 2022; 24, 1784. [DOI: https://dx.doi.org/10.3390/e24121784]

23. Eklundh, L.; Johansson, T.; Solberg, S. Mapping insect defoliation in Scots pine with MODIS time-series data. Remote Sens. Environ.; 2009; 113, pp. 1566-1573. [DOI: https://dx.doi.org/10.1016/j.rse.2009.03.008]

24. Steddom, K.; Heidel, G.; Jones, D.C.; Rush, C.M.; Garnica, V.C.; Giesler, L.J.; Mahlein, A.-K.; Reynolds, G.J.; Windels, C.E.; MacRae, I.V. et al. Remote Detection of Rhizomania in Sugar Beets. Phytopathology; 2003; 93, pp. 720-726. [DOI: https://dx.doi.org/10.1094/PHYTO.2003.93.6.720] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18943059]

25. Gröll, K.; Graeff, S.; Claupein, W. Use of Vegetation Indices to Detect Plant Diseases; Agrarinformatik im Spannungsfeld Zwischen Regionalisierung und Globalen Wertschöpfungsketten: Stuttgart, Germany, 2007; 98.

26. Zhao, H.; Yang, C.; Guo, W.; Zhang, L.; Zhang, D. Automatic Estimation of Crop Disease Severity Levels Based on Vegetation Index Normalization. Remote Sens.; 2020; 12, 1930. [DOI: https://dx.doi.org/10.3390/rs12121930]

27. Ojdanič, N.; Zelnik, I.; Holcar, M.; Gaberščik, A.; Golob, A. Contrasting Dynamics of Littoral and Riparian Reed Stands within a Wetland Complex of Lake Cerknica. Plants; 2023; 12, 1006. [DOI: https://dx.doi.org/10.3390/plants12051006] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36903867]

28. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring Vegetation Systems in the Great Plains with ERTS. Proceedings of the 3rd Earth Resource Technology Satellite (ERTS) Symposium; Washington, DC, USA, 10–14 December 1973; Volume 1, pp. 48-62.

29. Pervez, S.; McNally, A.; Arsenault, K.; Budde, M.; Rowland, J. Vegetation Monitoring Optimization with Normalized Difference Vegetation Index and Evapotranspiration Using Remote Sensing Measurements and Land Surface Models Over East Africa. Front. Clim.; 2021; 3, 589981. [DOI: https://dx.doi.org/10.3389/fclim.2021.589981]

30. Kumari, N.; Srivastava, A.; Dumka, U.C. A Long-Term Spatiotemporal Analysis of Vegetation Greenness over the Himalayan Region Using Google Earth Engine. Climate; 2021; 9, 109. [DOI: https://dx.doi.org/10.3390/cli9070109]

31. Arekhi, M.; Goksel, C.; Sanli, F.B.; Senel, G. Comparative Evaluation of the Spectral and Spatial Consistency of Sentinel-2 and Landsat-8 OLI Data for Igneada Longos Forest. ISPRS Int. J. Geo-Inf.; 2019; 8, 56. [DOI: https://dx.doi.org/10.3390/ijgi8020056]

32. Lasaponara, R.; Tucci, B.; Ghermandi, L. On the Use of Satellite Sentinel 2 Data for Automatic Mapping of Burnt Areas and Burn Severity. Sustainability; 2018; 10, 3889. [DOI: https://dx.doi.org/10.3390/su10113889]

33. Forkel, M.; Carvalhais, N.; Verbesselt, J.; Mahecha, M.D.; Neigh, C.S.R.; Reichstein, M. Trend Change Detection in NDVI Time Series: Effects of Inter-Annual Variability and Methodology. Remote Sens.; 2013; 5, pp. 2113-2144. [DOI: https://dx.doi.org/10.3390/rs5052113]

34. Schmid, J.N. Using Google Earth Engine for Landsat NDVI Time Series Analysis to Indicate the Present Status of Forest Stands; Georg-August-Universität Göttingen: Basel, Switzerland, 2017; [DOI: https://dx.doi.org/10.13140/RG.2.2.34134.14402/6]

35. Mbatha, N.; Xulu, S. Time Series Analysis of MODIS-Derived NDVI for the Hluhluwe-Imfolozi Park, South Africa: Impact of Recent Intense Drought. Climate; 2018; 6, 95. [DOI: https://dx.doi.org/10.3390/cli6040095]

36. Inghilesi, A.F.; Mazza, G.; Cervo, R.; Gherardi, F.; Sposimo, P.; Tricarico, E.; Zapparoli, M. Alien Insects in Italy: Comparing Patterns from the Regional to European Level. J. Insect Sci.; 2013; 13, 73. [DOI: https://dx.doi.org/10.1673/031.013.7301]

37. De Ros, G.; Conci, S.; Pantezzi, T.; Savini, G. The economic impact of invasive pest Drosophila suzukii on berry production in the Province of Trento, Italy. J. Berry Res.; 2015; 5, pp. 89-96. [DOI: https://dx.doi.org/10.3233/JBR-150092]

38. EFSA Panel on Plant Health (PLH) Bragard, C.; Baptista, P.; Chatzivassiliou, E.; Di Serio, F.; Gonthier, P.; Miret, J.A.J.; Justesen, A.F.; Magnusson, C.S.; Milonas, P. et al. Pest categorisation of Toumeyella parvicornis. EFSA J.; 2022; 20, e07146. [DOI: https://dx.doi.org/10.2903/j.efsa.2022.7146]

39. Hughes, M.J.; Kaylor, S.D.; Hayes, D.J. Patch-Based Forest Change Detection from Landsat Time Series. Forests; 2017; 8, 166. [DOI: https://dx.doi.org/10.3390/f8050166]

40. Huang, C.; Goward, S.N.; Masek, J.G.; Thomas, N.; Zhu, Z.; Vogelmann, J.E. An automated approach for reconstructing recent forest disturbance history using dense Landsat time series stacks. Remote Sens. Environ.; 2010; 114, pp. 183-198. [DOI: https://dx.doi.org/10.1016/j.rse.2009.08.017]

41. Tao, X.; Huang, C.; Zhao, F.; Schleeweis, K.; Masek, J.; Liang, S. Mapping forest disturbance intensity in North and South Carolina using annual Landsat observations and field inventory data. Remote Sens. Environ.; 2019; 221, pp. 351-362. [DOI: https://dx.doi.org/10.1016/j.rse.2018.11.029]

42. Chirici, G.; Giannetti, F.; Mazza, E.; Francini, S.; Travaglini, D.; Pegna, R.; White, J.C. Monitoring clearcutting and subsequent rapid recovery in Mediterranean coppice forests with Landsat time series. Ann. For. Sci.; 2020; 77, 40. [DOI: https://dx.doi.org/10.1007/s13595-020-00936-2]

43. Kennedy, R.E.; Yang, Z.; Cohen, W.B. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sens. Environ.; 2010; 114, pp. 2897-2910. [DOI: https://dx.doi.org/10.1016/j.rse.2010.07.008]

44. Mugiraneza, T.; Nascetti, A.; Ban, Y. Continuous Monitoring of Urban Land Cover Change Trajectories with Landsat Time Series and LandTrendr-Google Earth Engine Cloud Computing. Remote Sens.; 2020; 12, 2883. [DOI: https://dx.doi.org/10.3390/rs12182883]

45. Grabska, E.; Hawryło, P.; Socha, J. Continuous Detection of Small-Scale Changes in Scots Pine Dominated Stands Using Dense Sentinel-2 Time Series. Remote Sens.; 2020; 12, 1298. [DOI: https://dx.doi.org/10.3390/rs12081298]

46. Zhu, Z.; Woodcock, C.E. Continuous change detection and classification of land cover using all available Landsat data. Remote Sens. Environ.; 2014; 144, pp. 152-171. [DOI: https://dx.doi.org/10.1016/j.rse.2014.01.011]

47. Arévalo, P.; Bullock, E.L.; Woodcock, C.E.; Olofsson, P. A Suite of Tools for Continuous Land Change Monitoring in Google Earth Engine. Front. Clim.; 2020; 2, 576740. [DOI: https://dx.doi.org/10.3389/fclim.2020.576740]

48. Jahanifar, K.; Amirnejad, H.; Mojaverian, M.; Azadi, H. Land change detection and effective factors on forest land use changes: Application of land change modeler and multiple linear regression. J. Appl. Sci. Environ. Manag.; 2018; 22, 1269. [DOI: https://dx.doi.org/10.4314/jasem.v22i8.20]

49. Millington, J.D.A.; Perry, G.L.W.; Romero-Calcerrada, R. Regression Techniques for Examining Land Use/Cover Change: A Case Study of a Mediterranean Landscape. Ecosystems; 2007; 10, pp. 562-578. [DOI: https://dx.doi.org/10.1007/s10021-007-9020-4]

50. Telesca, L.; Lasaponara, R.; Faridani, F.; Abate, N.; Lovallo, M. Informational Analysis of MODIS Satellite Evapotranspiration Data of Vegetation Cover: A Method to Reveal the Presence of Plant Diseases. Proceedings of the Eighth International Conference on Advances in Signal, Image and Video Processing; Barcelona, Spain, 13–17 March 2023.

51. Rome, W. In Rome’s Pine Trees Are Dying as Time Runs Out to Save City Skyline. Available online: https://www.wantedinrome.com/news/romes-pine-trees-are-dying-as-time-runs-out-to-save-citys-skyline.html (accessed on 29 June 2023).

52. Tsatsaris, A.; Kalogeropoulos, K.; Stathopoulos, N.; Louka, P.; Tsanakas, K.; Tsesmelis, D.E.; Krassanakis, V.; Petropoulos, G.P.; Pappas, V.; Chalkias, C. Geoinformation Technologies in Support of Environmental Hazards Monitoring under Climate Change: An Extensive Review. ISPRS Int. J. Geo-Inf.; 2021; 10, 94. [DOI: https://dx.doi.org/10.3390/ijgi10020094]

53. Mutanga, O.; Kumar, L. Google Earth Engine Applications. Remote Sens.; 2019; 11, 591. [DOI: https://dx.doi.org/10.3390/rs11050591]

54. Kumar, L.; Mutanga, O. Google Earth Engine Applications Since Inception: Usage, Trends, and Potential. Remote Sens.; 2018; 10, 1509. [DOI: https://dx.doi.org/10.3390/rs10101509]

55. Gorelick, N. Google Earth Engine; American Geophysical Union: Vienna, Austria, 2013; 11997.

56. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ.; 2017; 202, pp. 18-27. [DOI: https://dx.doi.org/10.1016/j.rse.2017.06.031]

57. Mahdianpari, M.; Salehi, B.; Mohammadimanesh, F.; Homayouni, S.; Gill, E. The First Wetland Inventory Map of Newfoundland at a Spatial Resolution of 10 m Using Sentinel-1 and Sentinel-2 Data on the Google Earth Engine Cloud Computing Platform. Remote Sens.; 2019; 11, 43. [DOI: https://dx.doi.org/10.3390/rs11010043]

58. Campos-Taberner, M.; Moreno-Martínez, Á.; García-Haro, F.J.; Camps-Valls, G.; Robinson, N.P.; Kattge, J.; Running, S.W. Global Estimation of Biophysical Variables from Google Earth Engine Platform. Remote Sens.; 2018; 10, 1167. [DOI: https://dx.doi.org/10.3390/rs10081167]

59. Amini, S.; Saber, M.; Rabiei-Dastjerdi, H.; Homayouni, S. Urban Land Use and Land Cover Change Analysis Using Random Forest Classification of Landsat Time Series. Remote Sens.; 2022; 14, 2654. [DOI: https://dx.doi.org/10.3390/rs14112654]

60. Parente, L.; Mesquita, V.; Miziara, F.; Baumann, L.; Ferreira, L. Assessing the pasturelands and livestock dynamics in Brazil, from 1985 to 2017: A novel approach based on high spatial resolution imagery and Google Earth Engine cloud computing. Remote Sens. Environ.; 2019; 232, 111301. [DOI: https://dx.doi.org/10.1016/j.rse.2019.111301]

61. Hansen, C.H. Google Earth Engine as a Platform for Making Remote Sensing of Water Resources a Reality for Monitoring Inland Waters; Department of Civil and Environmental Engineering: Suite, UT, USA, 2015; [DOI: https://dx.doi.org/10.13140/RG.2.1.3688.1047]

62. Deines, J.M.; Kendall, A.D.; Crowley, M.A.; Rapp, J.; Cardille, J.A.; Hyndman, D.W. Mapping three decades of annual irrigation across the US High Plains Aquifer using Landsat and Google Earth Engine. Remote Sens. Environ.; 2019; 233, 111400. [DOI: https://dx.doi.org/10.1016/j.rse.2019.111400]

63. DeVries, B.; Huang, C.; Armston, J.; Huang, W.; Jones, J.W.; Lang, M.W. Rapid and robust monitoring of flood events using Sentinel-1 and Landsat data on the Google Earth Engine. Remote Sens. Environ.; 2020; 240, 111664. [DOI: https://dx.doi.org/10.1016/j.rse.2020.111664]

64. Fattore, C.; Abate, N.; Faridani, F.; Masini, N.; Lasaponara, R. Google Earth Engine as Multi-Sensor Open-Source Tool for Supporting the Preservation of Archaeological Areas: The Case Study of Flood and Fire Mapping in Metaponto, Italy. Sensors; 2021; 21, 1791. [DOI: https://dx.doi.org/10.3390/s21051791]

65. Orengo, H.A.; Conesa, F.C.; Garcia-Molsosa, A.; Lobo, A.; Green, A.S.; Madella, M.; Petrie, C.A. Automated detection of archaeological mounds using machine-learning classification of multisensor and multitemporal satellite data. Proc. Natl. Acad. Sci. USA; 2020; 117, pp. 18240-18250. [DOI: https://dx.doi.org/10.1073/pnas.2005583117]

66. Lasaponara, R.; Abate, N.; Masini, N. On the Use of Google Earth Engine and Sentinel Data to Detect “Lost” Sections of Ancient Roads. The Case of Via Appia. IEEE Geosci. Remote Sens. Lett.; 2021; 19, [DOI: https://dx.doi.org/10.1109/LGRS.2021.3054168]

67. Kumari, N.; Saco, P.M.; Rodriguez, J.F.; Johnstone, S.A.; Srivastava, A.; Chun, K.P.; Yetemen, O. The Grass Is Not Always Greener on the Other Side: Seasonal Reversal of Vegetation Greenness in Aspect-Driven Semiarid Ecosystems. Geophys. Res. Lett.; 2020; 47, e2020GL088918. [DOI: https://dx.doi.org/10.1029/2020GL088918]

68. Beck, H.E.; Zimmermann, N.E.; McVicar, T.R.; Vergopolan, N.; Berg, A.; Wood, E.F. Present and future Köppen-Geiger climate classification maps at 1-km resolution. Sci. Data; 2020; 7, 274. [DOI: https://dx.doi.org/10.1038/s41597-020-00616-w]

69. Arnfield, A.J. Koppen Climate Classification|Definition, System, & Map|Britannica. Available online: https://www.britannica.com/science/Koppen-climate-classification (accessed on 6 June 2023).

70. Didan, K. MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid V061. Distributed Active Archive Center. 2021; Available online: https://lpdaac.usgs.gov/products/mod13q1v061/ (accessed on 11 July 2023).

71. Huete, A.; Justice, C.; van Leeuwen, W. MODIS Vegetation Index (MOD 13). Algorithm Theor. Basis Doc.; 1999; 3, pp. 295-309.

72. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ.; 2002; 83, pp. 195-213. [DOI: https://dx.doi.org/10.1016/S0034-4257(02)00096-2]

73. Li, X.; Lanorte, A.; Lasaponara, R.; Lovallo, M.; Song, W.; Telesca, L. Fisher–Shannon and detrended fluctuation analysis of MODIS normalized difference vegetation index (NDVI) time series of fire-affected and fire-unaffected pixels. Geomat. Nat. Hazards Risk; 2017; 8, pp. 1342-1357. [DOI: https://dx.doi.org/10.1080/19475705.2017.1337652]

74. Meroni, M.; Fasbender, D.; Rembold, F.; Atzberger, C.; Klisch, A. Near real-time vegetation anomaly detection with MODIS NDVI: Timeliness vs. accuracy and effect of anomaly computation options. Remote Sens. Environ.; 2019; 221, pp. 508-521. [DOI: https://dx.doi.org/10.1016/j.rse.2018.11.041] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30774156]

75. Didan, K. MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid V006. NASA Eosdis Land Process. Daac; 2015; 10, 415.

76. Didan, K.; Munoz, A.B.; Solano, R.; Huete, A. MODIS Vegetation Index User’s Guide (MOD13 Series). Available online: https://vip.arizona.edu/documents/MODIS/MODIS_VI_UsersGuide_June_2015_C6.pdf (accessed on 27 June 2023).

77. Vermote, E. MOD09A1 MODIS/Terra Surface Reflectance 8-Day L3 Global 500 m SIN Grid V006. Distributed Active Archive Center. 2015; Available online: https://lpdaac.usgs.gov/products/mod09a1v006/ (accessed on 11 July 2023).

78. Khan, M.A.R.; Poskitt, D.S. Description Length Based Signal Detection in Singular Spectrum Analysis. Monash Econometrics and Business Statistics Working Papers. 2010; Available online: https://www.researchgate.net/profile/D-Poskitt-2/publication/46479471_Description_Length_Based_Signal_Detection_in_singular_Spectrum_Analysis/links/56bc33fc08ae47fa3956cc26/Description-Length-Based-Signal-Detection-in-singular-Spectrum-Analysis.pdf (accessed on 11 July 2023).

79. Vautard, R.; Ghil, M. Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. Phys. D Nonlinear Phenom.; 1989; 35, pp. 395-424. [DOI: https://dx.doi.org/10.1016/0167-2789(89)90077-8]

80. Hassani, H. Singular Spectrum Analysis: Methodology and Comparison. J. Data Sci.; 2021; 5, pp. 239-257. [DOI: https://dx.doi.org/10.6339/JDS.2007.05(2).396]

81. Frieden, B.R. Fisher information, disorder, and the equilibrium distributions of physics. Phys. Rev. A; 1990; 41, pp. 4265-4276. [DOI: https://dx.doi.org/10.1103/PhysRevA.41.4265]

82. Shannon, C.E. A Mathematical Theory of Communication. Bell Syst. Tech. J.; 1948; 27, pp. 379-423. [DOI: https://dx.doi.org/10.1002/j.1538-7305.1948.tb01338.x]

83. Hansen, M.C.; Potapov, P.V.; Moore, R.; Hancher, M.; Turubanova, S.A.; Tyukavina, A.; Thau, D.; Stehman, S.V.; Goetz, S.J.; Loveland, T.R. et al. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science; 2013; 342, pp. 850-853. [DOI: https://dx.doi.org/10.1126/science.1244693]

84. Potapov, P.; Hansen, M.C.; Pickens, A.; Hernandez-Serna, A.; Tyukavina, A.; Turubanova, S.; Zalles, V.; Li, X.; Khan, A.; Stolle, F. et al. The Global 2000-2020 Land Cover and Land Use Change Dataset Derived from the Landsat Archive: First Results. Front. Remote Sens.; 2022; 3, 856903. [DOI: https://dx.doi.org/10.3389/frsen.2022.856903]

85. Schoellhamer, D.H. Singular spectrum analysis for time series with missing data. Geophys. Res. Lett.; 2001; 28, pp. 3187-3190. [DOI: https://dx.doi.org/10.1029/2000GL012698]

86. Telesca, L.; Lovallo, M. On the performance of Fisher Information Measure and Shannon entropy estimators. Phys. A Stat. Mech. Appl.; 2017; 484, pp. 569-576. [DOI: https://dx.doi.org/10.1016/j.physa.2017.04.184]

87. Troudi, M.; Alimi, A.M.; Saoudi, S. Analytical Plug-In Method for Kernel Density Estimator Applied to Genetic Neutrality Study. EURASIP J. Adv. Signal Process.; 2008; 2008, 739082. [DOI: https://dx.doi.org/10.1155/2008/739082]

88. Vignat, C.; Bercher, J.-F. Analysis of signals in the Fisher–Shannon information plane. Phys. Lett. A; 2003; 312, pp. 27-33. [DOI: https://dx.doi.org/10.1016/S0375-9601(03)00570-X]

89. Bertin, S.; Ilardi, F.; Scapini, C.; Simoni, S.; Roversi, P.F. Alien pest Toumeyella parvicornis (cockerell) (hemiptera: Coccidae) on Pinus pinea L.: Short time evaluation of endotherapic treatment. Redia; 2022; 105, pp. 11-16. [DOI: https://dx.doi.org/10.19263/REDIA-105.22.02]