Recent NDVI Trends in Mainland Spain : LandCover and Phytoclimatic-Type Implications

The temporal evolution of vegetation is one of the best indicators of climate change, and many earth system models are dependent on an accurate understanding of this process. However, the effect of climate change is expected to vary from one land-cover type to another, due to the change in vegetation and environmental conditions. Therefore, it is pertinent to understand the effect of climate change by land-cover type to understand the regions that are most vulnerable to climate change. Hence, in this study we analyzed the temporal statistical trends (2001–2016) of the MODIS13Q1 normalized difference vegetation index (NDVI) to explore whether there are differences, by land-cover class and phytoclimatic type, in mainland Spain and the Balearic Islands. We found 7.6% significant negative NDVI trends and 11.8% significant positive NDVI trends. Spatial patterns showed a non-random distribution. The Atlantic biogeographical region showed an unexpected 21% significant negative NDVI trends, and the Alpine region showed only 3.1% significant negative NDVI trends. We also found statistical differences between NDVI trends by land cover and phytoclimatic type. Variance explained by these variables was up to 35%. Positive trends were explained, above all, by land occupations, and negative trends were explained by phytoclimates. Warmer phytoclimatic classes of every general type and forest, as well as some agriculture land covers, showed negative trends.


Introduction
Vegetation and the climate system interact through feedback [1].It is known that plants control water, energy, carbon cycles, the exchange of carbon dioxide with the atmosphere (i.e., photosynthesis, respiration), and the influence on wind circulation [2,3].This explains the importance of plants on Earth.On the other hand, vegetation dynamics are driven by the climate, in particular, incoming radiation, temperature, precipitation, and atmospheric humidity [4].In addition, nutrient availability and short-term natural and anthropogenic disturbances such as fires, volcanic eruptions, logging or insect epidemics, agriculture practices, and the global change affect the health of plants and must be studied in some different spatial-temporal scales [1,[5][6][7].For example, changes at the beginning and end of the growing season interact with the carbon cycle [8,9].In this sense, it is not surprising that, for some years, climate researchers focused on carbon balance and its relationship with the evolution of plants (for instance, References [10][11][12]).Thereby, it seems relevant to highlight that vegetation holds around 42% of the terrestrial carbon storage and assimilates about 20% of the annual anthropogenic emissions of that carbon [7,13].
All of the mentioned studies demonstrated the sensitivity of the terrestrial biosphere to climate fluctuations and the complexity of the relationship between the environment, climate change, and plants functioning.In that way, in recent years, the number of scientific studies exploring the time series of the vegetation index based on satellite data exponentially increased.Studies that focused on evaluating remote-sensing time series were, overall, related to the land-cover change [14], the forest perturbances [15], or the phenological patterns [16] and vegetation index trends over shortand long-term periods [17,18].
The most used vegetation quantity index derived from remote-sensing resources is the normalized difference vegetation index (NDVI) [19].NDVI is an estimate of the fraction of photosynthetically active radiation absorbed by Earth's vegetation [20], which was successfully used to monitor global productivity over the past three decades [21][22][23].The NDVI is strongly related to the biophysical parameter of plants, and, for that reason, this index is used to assess long-term trends of vegetation [24,25].The NDVI index was also used to detect vegetation and climate interactions [26,27], modeling the net primary production [28], and the carbon balance of terrestrial ecosystems [29].Other studies investigated NDVI correlations with droughts, vegetation functional types [30,31], ecosystems [32], and land covers [33].Despite the limitations of this index due to atmospheric and angular effects, saturation, and background reflectance contamination, the use of NDVI increased in recent years because of the availability of large temporal series.
The first applications in remote sensing to evaluate trends and phenological vegetation studies used Landsat sensor data.Then, satellite sensors were used with a more frequent repeat cycle.The longest-running series of high-repeat-frequency sensors is the National Oceanic and Atmospheric Administration's (NOAA) advanced very-high-resolution radiometer (AVHRR).The AVHRR vegetation index datasets are available in a database from 1982 in an 8-km resampling grid covering the globe [19].These sensor studies include References [34,35].Then, with the launch of the Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI 250-m data, time-series data analysis could have a moderate (regional) spatiotemporal resolution.For that reason, since MODIS was launched in December 1999, it was frequently used for vegetation studies [22,[36][37][38][39][40][41][42].Some of the advantages to these works were the MODIS improved geometry, radiometry, and quality combined with its free charge policy distribution.MODIS imagery captured by two sensors (Aqua and Terra satellites) is distributed at various pre-processing levels and, regarding the temporal distribution, data are released as both daily and composite products.These MODIS composited products are also generated at different compositing steps (eight-day, 16-day, monthly) and have some advantages compared to daily data, because the compositing process reduces the effect of clouds, snow, and noise [43].
Regarding the comparability of MODIS and other sensor time series of vegetation index results, it seems interesting to highlight that some previous studies revealed similar spatial patterns in the NDVI between AVHRR and MODIS products throughout the world.We found some examples in which this coherency could be seen for northeastern Brazil [44], China [45,46], and Europe [33].However, this latest work showed lower coherency results for the Iberian Peninsula, especially with phenology parameters.In this sense, in the present paper, we only worked with MODIS because we considered that the coarse spatial and high temporal resolution of that sensor's vegetation products permit us to understand the dynamics of the land surface at a regional scale [47], as applied to mainland Spain and the Balearic Islands.Moreover, MODIS vegetation indices were found to be better correlated to in situ measured vegetation indices, in many cases, as compared to NOAA/AVHRR NDVI [48] or SPOT VGT [49].
Several studies showed the relationship between time series of NDVI values and the climatic factors such as temperature, precipitation, incoming radiation, water vapor content, and droughts [27,47,50,51].On the other hand, there are works in which the authors found a relationship between the trends of the NDVI over the last decades and the land cover and its management [30,[51][52][53].
Furthermore, many works showed an increase in temperature and a regional decrease in precipitation, and models projected an increase in droughts in several parts of the world [54][55][56].In the Iberian Peninsula, an increase in temperature was observed, especially during the summer [57].Changes in precipitation regime, with a mean decrease [58], and an increase in potential evapotranspiration were observed in recent decades [59,60].All of these data are accompanied by several works that warn of the effects on vegetation: increasing insect damage [61,62], reduced growth [53,63], and overall decoupling between plants and territories, which leads to changing climate conditions [64,65].
The climate transition zones are considered as being potentially at risk because of climate change [66].In the Iberian Peninsula, there are three biogeographical regions: Atlantic, Mediterranean, and Alpine (see Figure 1), with a large temperature and precipitation gradient.Phytoclimates try to correlate natural or near-natural habitats or plant associations and climate, using precipitation and temperature thresholds [67,68].In the Iberian Peninsula, the change in these variables led to an increase in drier and warmer types, and a decrease in wetter and fresher ones, that is, a "Mediterranization" [65,69,70].The Co-Ordination of Information on the Environment (CORINE) Land Cover (CLC) is a set of maps of the European environmental landscape which are based on the photo-interpretation of satellite images (SPOT, Landsat TM, and Landsat MSS) and ancillary data (aerial photographs, topographic or vegetation maps, statistics, and local knowledge).The maps present a scale of 1:100.000.The smallest surfaces mapped were 25 hectares, and linear features less than 100 m in width were not mapped [71].The high level of disaggregation and its global coverage make it an interesting coverage as a reference.We used the occupancy map of the CORINE Land Cover 2006 from version 18.5 and three levels of depth (which includes 44 land covers).We used the 2006 version because our In this sense, our objective here was to evaluate the annual NDVI trends in mainland Spain and the Balearic Islands-a zone with high climate anomalies observed in this century-in the relationship with the different land cover and phytoclimatic types present on these territories.

Study Site and Co-Ordination of Information on the Environment (CORINE) Land Cover 2006
The study area in this research was mainland Spain and the Balearic Islands.We can see in Figure 1 that three biogeographical regions are present: Atlantic, Mediterranean, and Alpine.The existence of mountains increases climate variability.For general land-cover types, see Figure 2.

Phytoclimatic Allué Classification
The phytoclimatic regions used in the present study are those delimited in the Allué [68] work for the Iberian Peninsula and the Balearic Islands.The Allué classification was built upon meteorological data, collected from more than 1000 stations, from the Spanish Meteorological Agency and the potential vegetation series elaborated by Reference [72].These phytoclimatic classes enclose a large gradient of bioclimatic conditions, ranging from the Atlantic zone of influence to the most arid areas of the southeast of Spain.The Atlantic area of influence is characterized by high precipitation and mild temperatures throughout the year, and its potential vegetation is the broad-leaved deciduous forest class.In the most arid areas and the Ebro depression, the rainfall is low, the dry season continues beyond the end of summer, and the temperature is higher throughout the year.In the arid areas, the natural vegetation potential corresponds to sparse formations of spiny shrubs.
The resulting Allué phytoclimatic classification consists of 19 different subtypes of vegetation which are associated with different characteristic climatic conditions.Table 2 shows the Allué phytoclimatic zones with a high level of detail, and Figure 3 shows the spatial pattern of the classes that build up to it.
In Table 2, it is necessary to take into account that T´m is the average minimum temperature from the coldest month (normally January), while tf is the lowest average monthly temperature.The time lapse, where the temperature curve is above the rainfall, is represented with a vowel; it is expressed in months and fractions of months, ranging between 11 and zero months.Four types are identified: 3 ≤ a < 11; 1.25 ≤ a < 3; 0 ≤ a < 1.25; a = 0.A value greater than 11 indicates a dry climate, while a low value indicates a wet climate.P is the total annual precipitation value; Hs is the most probable period of frost (in months); Hp is the period (in months) of secure frost; Tc is the greatest monthly average temperature.The precipitation is given in mm, and the temperature is given in °C.The Co-Ordination of Information on the Environment (CORINE) Land Cover (CLC) is a set of maps of the European environmental landscape which are based on the photo-interpretation of satellite images (SPOT, Landsat TM, and Landsat MSS) and ancillary data (aerial photographs, topographic or vegetation maps, statistics, and local knowledge).The maps present a scale of 1:100.000.The smallest surfaces mapped were 25 hectares, and linear features less than 100 m in width were not mapped [71].The high level of disaggregation and its global coverage make it an interesting coverage as a reference.We used the occupancy map of the CORINE Land Cover 2006 from version 18.5 and three levels of depth (which includes 44 land covers).We used the 2006 version because our temporal range was 2001-2016, and changes between 2006 and 2012 were less than 1.5%.Table 1 and Figure 2 show the codes of CORINE Land Cover classes.However, we did not take into account the level 1 artificial surfaces, and from the wetland ones, we only studied the inland marshes.We converted the vector shapefile to raster with the MODIS NDVI resolution.All of the CLC classes used here had a considerable number of pixels.Those with fewer pixels are related to beaches, dunes, and sand with more than 3000 pixels.

Phytoclimatic Allué Classification
The phytoclimatic regions used in the present study are those delimited in the Allué [68] work for the Iberian Peninsula and the Balearic Islands.The Allué classification was built upon meteorological data, collected from more than 1000 stations, from the Spanish Meteorological Agency and the potential vegetation series elaborated by Reference [72].These phytoclimatic classes enclose a large gradient of bioclimatic conditions, ranging from the Atlantic zone of influence to the most arid areas of the southeast of Spain.The Atlantic area of influence is characterized by high precipitation and mild temperatures throughout the year, and its potential vegetation is the broad-leaved deciduous forest class.In the most arid areas and the Ebro depression, the rainfall is low, the dry season continues beyond the end of summer, and the temperature is higher throughout the year.In the arid areas, the natural vegetation potential corresponds to sparse formations of spiny shrubs.
The resulting Allué phytoclimatic classification consists of 19 different subtypes of vegetation which are associated with different characteristic climatic conditions.Table 2 shows the Allué phytoclimatic zones with a high level of detail, and Figure 3 shows the spatial pattern of the classes that build up to it.Allué (1990) regions used in this study.The spatial distribution of the phytoclimatic types is displayed in Figure 3.

The Applied MODIS Dataset
We used the 16-day composite MODIS VI product with 250-m spatial resolution from Terra (MOD13Q1).We worked with 368 images, i.e., all the available images located in our study area that were captured by the MODIS sensor (TERRA) from 2001 to 2016.These images were downloaded from the Earth Data National Aeronautics and Space Administration (NASA) website (https://earthdata.nasa.gov/).
The MOD13Q1 elaborated product (vegetation indices 16-day L3 global 250 m) had one NDVI value (calculated with bands 1 and 2 of the MODIS sensor) every 16 days and allowed the composition of gap-free temporal series [73].This product also brings quality information from the layer reliability (PR).It is important to highlight that, every 16 days, the MODIS team named the image with a nominal date related to the first date of each composition.However, the real acquisition date of each pixel (CDOY) usually differs from the nominal date, although the acquisition may have occurred on any one of the 16 days [73].

Smoothing NDVI Pixel Curves in Spain with TIMESAT Software (Filtering Noise NDVI Profiles).
In the TIMESAT software, Jonsson and Eklundh [16] implemented three different methods.These methods are based on least-squares fitting, as we previously mentioned in Section 1.These methods use NDVI captured by sensor values to determine the shape of the current temporal curve belonging to every pixel.We used the Gaussian filtering TIMESAT algorithm which performs better in heterogeneous environments such as mainland Spain and the Balearic Islands [74].TIMESAT provides a weighting mechanism, although some values in the NDVI time series can be more influential than others.Thus, we assigned high weights for higher-quality retrievals of MODIS data and low weights for lower-quality ones.Therefore, our retrievals with lower quality were those that less influenced the final smoothed curve.After obtaining the smoothing curves of NDVI per pixel, In Table 2, it is necessary to take into account that T'm is the average minimum temperature from the coldest month (normally January), while tf is the lowest average monthly temperature.The time lapse, where the temperature curve is above the rainfall, is represented with a vowel; it is expressed in months and fractions of months, ranging between 11 and zero months.Four types are identified: 3 ≤ a < 11; 1.25 ≤ a < 3; 0 ≤ a < 1.25; a = 0.A value greater than 11 indicates a dry climate, while a low value indicates a wet climate.P is the total annual precipitation value; Hs is the most probable period of frost (in months); Hp is the period (in months) of secure frost; Tc is the greatest monthly average temperature.The precipitation is given in mm, and the temperature is given in • C.
As with the CLC06 classes for the Allué codes, there was a high number of pixels in our study area.The area with fewer pixels is that related to type III (IV) with 2987.

The Applied MODIS Dataset
We used the 16-day composite MODIS VI product with 250-m spatial resolution from Terra (MOD13Q1).We worked with 368 images, i.e., all the available images located in our study area that were captured by the MODIS sensor (TERRA) from 2001 to 2016.These images were downloaded from the Earth Data National Aeronautics and Space Administration (NASA) website (https://earthdata.nasa.gov/).
The MOD13Q1 elaborated product (vegetation indices 16-day L3 global 250 m) had one NDVI value (calculated with bands 1 and 2 of the MODIS sensor) every 16 days and allowed the composition of gap-free temporal series [73].This product also brings quality information from the layer reliability (PR).It is important to highlight that, every 16 days, the MODIS team named the image with a nominal date related to the first date of each composition.However, the real acquisition date of each pixel (CDOY) usually differs from the nominal date, although the acquisition may have occurred on any one of the 16 days [73].In the TIMESAT software, Jonsson and Eklundh [16] implemented three different methods.These methods are based on least-squares fitting, as we previously mentioned in Section 1.These methods use NDVI captured by sensor values to determine the shape of the current temporal curve belonging to every pixel.We used the Gaussian filtering TIMESAT algorithm which performs better in heterogeneous environments such as mainland Spain and the Balearic Islands [74].TIMESAT provides a weighting mechanism, although some values in the NDVI time series can be more influential than others.Thus, we assigned high weights for higher-quality retrievals of MODIS data and low weights for lower-quality ones.Therefore, our retrievals with lower quality were those that less influenced the final smoothed curve.After obtaining the smoothing curves of NDVI per pixel, we calculated the average of every pixel that made the study area (n = 368) per year (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016).Finally, we eliminated those pixels in which the NDVI value was less than zero, because we considered them as errors or bare soil.

Temporal Trend Analyses
From the data of the mean annual values of NDVI obtained, as explained in Section 2.5, we started the yearly trend analyses.Trend analysis employed non-parametric approaches, namely Theil-Sen regressions [75,76].The slopes of the regression approach were tested for their statistical significance using the z-statistic and the p-value of the Mann-Kendall test for slopes.The Mann-Kendall test and the Theil-Sen regression were performed using the Clark Labs TerrSet Earth Trends Modeler (ETM), which is an integrated collection of tools for the analysis of image time-series data.
Mann-Kendall is a temporal trend estimator that is more robust than the least-squares slope because it is much less sensitive to outliers and skewed data (https://clarklabs.org/terrset/).We used this estimator to identify pixels, for which the p-value (Mann-Kendall) was less than 0.1 (90% confidence interval (CI)).After obtaining those pixels that we considered to have a significant short-term trend, we applied the Theil-Shen regression to obtain the Theil-Sen positive or negative slopes of all significant NDVI trend pixels from 2001 to 2016.
Then, with the help of geographic information system (GIS) software, we crossed the Theil-Sen-filtered raster layer with the CLC 06 and Allué phytoclimatic classes.
Figure 4 summarizes the methodological process performed to obtain the Theil-Sen pixel slopes from the 2001-2016 period, belonging to every CLC 06 and Allué phytoclimatic class in our study zone.

Inferential Analyses
From the pixels with a significant trend (90% CI) within the 2001-2016 period, crossed with the CLC 06 and Allué phytoclimatic classes, we could calculate the percentage of pixels with significant positive and significant negative trends by biogeographical region, land cover, and phytoclimatic region within the study area.We carried out analyses of the variance in which the numeric variables were always the percentage of pixels with a significant trend (from the 2001-2016 yearly period) and the categorical variables were the CLC 06 classes and the Allué phytoclimatic regions.We performed ANOVA tests and the Kruskal-Wallis test (a non-parametrical test because our data did not meet the heteroscedasticity assumption).These variance analysis tests were carried out to study, first of all, whether the percentage of significant trend pixels (total, positive, and negative) differed by CLC 06 and by phytoclimatic region in our study area.We also made multiple-regression models (with the same input variables as for the variance tests); thereby, we determined the independent contribution of each explanatory variable, to the response variable of the joint contribution resulting from the correlation with other variables, using the hierarchical partitioning algorithm [77].The hierarchical partitioning tests were performed using R software with the help of the "hier.part" package (version 1.0-4).These multiple-regression models are described below.
In Model 1, the variable used was the percentage of significant pixels.In Model 2, the variable used was the percentage of significant negative pixels, whereas, in model 3 this was the percentage of significant positive pixels.For the three models, the explanatory variables were the phytoclimatic region and the CLC 06 class.

Inferential Analyses
From the pixels with a significant trend (90% CI) within the 2001-2016 period, crossed with the CLC 06 and Allué phytoclimatic classes, we could calculate the percentage of pixels with significant positive and significant negative trends by biogeographical region, land cover, and phytoclimatic region within the study area.We carried out analyses of the variance in which the numeric variables were always the percentage of pixels with a significant trend (from the 2001-2016 yearly period) and the categorical variables were the CLC 06 classes and the Allué phytoclimatic regions.We performed ANOVA tests and the Kruskal-Wallis test (a non-parametrical test because our data did not meet the heteroscedasticity assumption).These variance analysis tests were carried out to study, first of all, whether the percentage of significant trend pixels (total, positive, and negative) differed by CLC 06 Furthermore, we created a statistic contingency table, a type of table in matrix format that displays the CLC and trend frequency distribution of the variables.We made another cross table, but with the phytoclimatic types and the trend frequency distribution of these variables.From this table, we calculated the difference in percentage between the observed and the expected NDVI positive and negative pixels with significant short-term NDVI trends, and made graphs with those differences.Classes with fewer than 30 pixels were excluded from this calculation.

Temporal Short-Term Trend General Results
Figure 5 shows how the positive NDVI trends of the pixels were more common than the negative trends from 2001 to 2016.The total number of pixels used was more than 9,000,000.The percentage of pixels with significant negative trends reached 7.6%, whereas the percentage of pixels with significant positive trends reached 11.8%.Table 3 shows the percentage of pixels, by biogeographical region, that exhibited significant positive and negative NDVI trends.The Atlantic region had 21.1% of pixels with a negative trend and only 6.9% of pixels with a positive trend.The Alpine region showed the opposite result, with only 3.1% of pixels with a negative trend and 16.4% with a positive trend.The spatial pattern of significant trend distribution, shown in Figure 5, indicates how the significant negative trends appear overall on the Spanish Atlantic and Mediterranean coasts, to the south of the Iberian mountain range, in the Sierra Morena mountain range, and in the Balearic Islands.
By contrast, the positive NDVI trends from 2001 to 2016 were mainly located in the northeast of Spain, Northern Plateau, Central System, and Guadalquivir Valley.The results of the parametrical (ANOVA) and non-parametrical (Kruskal-Wallis) variance analyses, performed using the percentage of the total significant pixels by each CLC and Allué phytoclimatic class, differed significantly (p < 0.05) by phytoclimatic type and land occupancy.We shows the pixels with negative Theil-Sen slopes and a Mann-Kendall z-statistic below 0.1 (significant trend).The green color shows the pixels with positive Theil-Sen slopes and a Mann-Kendall z-statistic below 0.1 (significant trend).Finally, the black color shows the pixels with a Mann-Kendall z-statistic higher than 0.1 (non-significant trend).

Inferential Analysis Results of the Pixel Percentage with a Significant Short-Term Trend by Allué Phytoclimatic and CORINE Land-Cover Classes
The results of the parametrical (ANOVA) and non-parametrical (Kruskal-Wallis) variance analyses, performed using the percentage of the total significant pixels by each CLC and Allué phytoclimatic class, differed significantly (p < 0.05) by phytoclimatic type and land occupancy.We could also notice statistical differences in the average between the groups of these two classes (CLC and Allué phytoclimatic types) when we performed the analyses focused first on the significant negative trend sample and then on the significant positive trend sample.The positive pixel percentage and the negative pixel percentage also differed significantly by CLC and Allué phytoclimatic regions.
Figure 6 shows the interval mean plots (95% CI) from the variance analyses in which the percentage of significant trend pixels was taken into account.The Allué phytoclimatic type average percentages ranged from 10% to 25%; however, the CLC land-cover percentages ranged from 11% to approximately 37%.The two Allué phytoclimatic types with lower significant pixel percentages were, on average, types X(IX)1 and X(IX)2, i.e., the Oromarticoids or those less influenced by the Mediterranean climate and located in high mountainous areas.By contrast, the Allué phytoclimatic types with a higher mean of significant trend pixels were those related to types III(IV) and IV(VI)2, VI(IV)3, VI(IV)4, VI(V), and VI(VII).Types III(IV) and IV(VI)2 are the most arid types from our study zone; types VI(IV)3 and VI(IV)4 are between the Mediterranean and nemoral type, i.e., between the Mediterranean and a continental climate.Finally, types VI(V) and the VI(VII) are within the nemoral category, similar to a continental climate (see Figure 6a).By CLC class, it can be seen that the fruit trees and berry plantations (222), beaches, dunes, and sand (331), and inland marshes (411) were those with a higher percentage of significant trends from 2001 to 2016.By contrast, the pastures (231), agro-forestry (244), natural grasslands (321), bare rocks (332), and burnt areas (334) were those that showed the lowest percentage of significant trends.
Figure 6 shows the interval mean plots (95% CI) from the variance analyses in which the percentage of significant negative and positive pixels was taken into account.At a glance, it can be noted that, by Allué phytoclimatic regions, the highest rate of negative trends was found in type III(IV), i.e., arid, types IV(III) and IV(VI)2 in the Mediterranean class, type VI (IV)3 in the nemo-Mediterranean type and VI(V) in the nemoral class (rounding the 15% on average).By contrast, the highest percentage of positive trend pixels was found in types VI(IV)4 and VI (VII), belonging to the nemo-Mediterranean and the Mediterranean types (more than 20%) (see Figure 6a).On the other hand, by CLC class, the greatest percentage of negative trends was found in rice fields (213), fruits (222), and inland marshes (411) (>13%) while, the highest positive trend percentage was found in olive groves (223), beaches and dunes (331), and inland marshes (411) (> 20%) (See Figure 6b).
As shown in Figure 7, the greatest difference between the positive pixel percentage and the negative pixel percentage was found in type IV(III), i.e., arid, types VI(IV)1, VI(IV)2, and VI(IV)4 in the nemo-Mediterranean type, type VI(VII) in the nemoral class, and type X(VIII) in the Oroboreal category.Moreover, all of these classes show how the positive percentage was higher than the negative one.The same occurred with the CLC classes; the greatest differences in the percentage of positive and negative pixels were found in non-irrigated arable land (211), vineyards (221), olive groves (223), agro-forestry areas (244), beaches and dunes (331), and bare rocks (332).In all of these CLC categories, the percentage of positive pixels was higher than the percentage of negative pixels.Figure 6 shows the interval mean plots (95% CI) from the variance analyses in which the percentage of significant negative and positive pixels was taken into account.At a glance, it can be noted that, by Allué phytoclimatic regions, the highest rate of negative trends was found in type III(IV), i.e., arid, types IV(III) and IV(VI)2 in the Mediterranean class, type VI (IV)3 in the nemo-Mediterranean type and VI(V) in the nemoral class (rounding the 15% on average).By contrast, the  Figure 7b shows that the classes of fruit trees (222), pastures (231), complex cultivation patterns (242), and broad-leaved (311) and coniferous (312) forests had a higher negative NDVI percentage from 2001-2016 than positive NDVI percentage.This negative percentage was higher than one-third the positive percentage in the case of the broad-leaved and coniferous masses.Although it is not shown in Figure 7, for the positive trend pixels, there were some combinations of CLC/Allué regions in which the percentage was particularly high.The positive statistical trends were evident in moors and heathland (322) with type IV2, which is from the Mediterranean phytoclimatic type; and in beaches, dunes, and sand (331) with types VI and X(VIII) belonging to the nemoral and Oroboreal classes, respectively.For the significant negative trends, the highest Although it is not shown in Figure 7, for the positive trend pixels, there were some combinations of CLC/Allué regions in which the percentage was particularly high.The positive statistical trends were evident in moors and heathland (322) with type IV2, which is from the Mediterranean phytoclimatic type; and in beaches, dunes, and sand (331) with types VI and X(VIII) belonging to the nemoral and Oroboreal classes, respectively.For the significant negative trends, the highest percentages were found in fruit trees with type VIII(IV) (Oroboreal class) and in beaches, dunes, and sand with type IV(VI)1 (within the Mediterranean class).
Figure 8 shows the percentage of significant trends by CLC class (a) and phytoclimatic types (b), and the difference charts (y-axis) with bars that allow us to identify CLC classes (c) and phytoclimatic types (d) with the highest percentage of the difference between observed and expected significant negative and positive NDVI trends.We changed the sign for negative trends to show categories with the harmful vegetation dynamic (for negative trends, percentages observed higher than expected percentages and vice versa for positive trends) in the lower part of the figure.These results showed some differences by CLC class and also by phytoclimatic type.By land-cover class, the highest percentage of negative trends was found in rice fields (213), fruits (222), pastures (231), complex cultivation patterns (242), all the forests (311, 312, and 313), and moors and heathland (322).The highest percentage of positive trends was found in permanently irrigated areas ( 212 Figure 8 shows the percentage of significant trends by CLC class (a) and phytoclimatic types (b), and the difference charts (y-axis) with bars that allow us to identify CLC classes (c) and phytoclimatic types (d) with the highest percentage of the difference between observed and expected significant negative and positive NDVI trends.We changed the sign for negative trends to show categories with the harmful vegetation dynamic (for negative trends, percentages observed higher than expected percentages and vice versa for positive trends) in the lower part of the figure.These results showed some differences by CLC class and also by phytoclimatic type.By land-cover class, the highest percentage of negative trends was found in rice fields (213), fruits (222), pastures (231), complex cultivation patterns (242), all the forests (311, 312, and 313), and moors and heathland (322).The highest percentage of positive trends was found in permanently irrigated areas ( 212  In Figure 9, it can be observed that the percentages of the difference between significant positive and negative NDVI short-term trends, shown by combinations of CLC classes and phytoclimatic types, were also different.For instance, in the general graph from Figure 8a, the 311, 312, and 313 CLC classes presented a percentage of pixels with a negative trend greater than the mean and up to 50% of the difference between observed and expected NDVI pixels with significant NDVI short-term trends.Thus, the evolution of the NDVI values from 2001 to 2016 was weak in the aforementioned In Figure 9, it can be observed that the percentages of the difference between significant positive and negative NDVI short-term trends, shown by combinations of CLC classes and phytoclimatic types, were also different.For instance, in the general graph from Figure 8a, the 311, 312, and 313 CLC classes presented a percentage of pixels with a negative trend greater than the mean and up to 50% of the difference between observed and expected NDVI pixels with significant NDVI short-term trends.Thus, the evolution of the NDVI values from 2001 to 2016 was weak in the aforementioned CLC classes.However, if we focus on Figure 9a-c, we can observe that these values occurred in some phytoclimatic types, but not in others.For instance, the behavior of classes 311 and 312 was similar, except for phytoclimates IV(III), with a good dynamic for 311 and a bad dynamic for 312; type VIII(VI) contrasts this finding.Negative trends in IV(VI)2 and IV(VI)2 did not occur for 313, 311, and 312.In all of the phytoclimates, we found a higher percentage of significant negative NDVI trends except in type IV (VI)1 (more than 1200% between observed and expected pixels) and in type IV3 (100% of the difference between observed and expected pixels).
ISPRS Int.J. Geo-Inf.2018, 7, x; doi: FOR PEER REVIEW www.mdpi.com/journal/ijgiCLC classes.However, if we focus on Figures 9a-c, we can observe that these values occurred in some phytoclimatic types, but not in others.For instance, the behavior of classes 311 and 312 was similar, except for phytoclimates IV(III), with a good dynamic for 311 and a bad dynamic for 312; type VIII(VI) contrasts this finding.Negative trends in IV(VI)2 and IV(VI)2 did not occur for 313, 311, and 312.In all of the phytoclimates, we found a higher percentage of significant negative NDVI trends except in type IV (VI)1 (more than 1200% between observed and expected pixels) and in type IV3 (100% of the difference between observed and expected pixels).The three hierarchical partitioning (HP) tests carried out-with the percentage of the total significant trend pixels, by CLC and Allué phytoclimatic combinations, as the y variable-showed that the total variance explained by these CLC and Allué phytoclimatic variables reached 35; this also occurred for the HP of the negative trend percentage numerical variable (see Table 4).However, the CLC and Allué phytoclimatic classes could explain 46% of the positive trend percentages.Regarding the individual variability explained in each of the three HP models (I), it is important to highlight the following: in the HP made to explain the percentage of trend pixels, the CLC explained three times more variability than the phytoclimatic types.In the HP test results related to the negative trends, the phytoclimatic types explained a higher percentage of variability than the CLC classes; on the contrary, in the HP test performed to explain the positive trend percentages, the CLC explained more variability than the phytoclimatic types.The three hierarchical partitioning (HP) tests carried out-with the percentage of the total significant trend pixels, by CLC and Allué phytoclimatic combinations, as the y variable-showed that the total variance explained by these CLC and Allué phytoclimatic variables reached 35; this also occurred for the HP of the negative trend percentage numerical variable (see Table 4).However, the CLC and Allué phytoclimatic classes could explain 46% of the positive trend percentages.Regarding the individual variability explained in each of the three HP models (I), it is important to highlight the following: in the HP made to explain the percentage of trend pixels, the CLC explained three times more variability than the phytoclimatic types.In the HP test results related to the negative trends, the phytoclimatic types explained a higher percentage of variability than the CLC classes; on the contrary, in the HP test performed to explain the positive trend percentages, the CLC explained more variability than the phytoclimatic types.

Discussion
In this study, we utilized MODIS NDVI data to examine long-term changes in the vegetation greenness in mainland Spain and the Balearic Islands from 2001 to 2016.We evaluated the percentage of significant negative and positive trends by land-cover occupancies and phytoclimatic types to find patterns in our study area.
In our study, we found that, in mainland Spain and the Balearic Islands, the annual average NDVI dataset showed that the area with significant positive trends was more important (11.8%)than the territory with negative trends (7.6%).This increasing photosynthetic activity agrees with other analyses made by satellite time series [21,23,[78][79][80][81][82][83] and analyses of time series of phenological field observations (e.g., References [84][85][86][87]), which provide evidence of the start of the growing season in the 1980s, mainly related to increased spring temperatures.Moreover, these shifts correspond to the increased annual accumulation of growing degree days and are most highly concentrated in eastern United States and Europe.The general growth of the photosynthetic activity in Europe was also explained by the increase in the forested area, the CO 2 fertilization, the elevated atmospheric nitrogen deposits, the juvenile age structure, and climate change [88].
It is important to remark that our significant negative NDVI trend percentage (7.6%)was higher than the negative NDVI trend percentage found in other studies for the same site in the latter part of the 20th century [57,[89][90][91].The result of the Atlantic biogeographic region is quite unexpected, reaching a percentage of 21.6% for the pixels with a negative trend.For instance, in Alcaraz-Segura et al. [89], the authors evaluated the consistency of the 1982-1999 NDVI trends in the Iberian Peninsula across time series derived from the AVHRR sensor, also using a similar methodology as that used here, but for an earlier period.They also found that the percentage of positive trends in mainland Spain was higher than the percentage of positive trends.However, the percentage of significant positive trends they found was around 23%, while the percentage of significant negative trends was lower than 1%.Furthermore, Fesholt et al. [92] calculated vegetation trends in drylands and semi-arid regions across the world, through the MOD13Q1 product, finding an overall greening for semi-arid areas across the globe over 27 years.In the Fesholt et al. [92] results related to mainland Spain, it can be seen how the majority of significant NDVI trends were positive from 1982 to 2007.García-Ruíz et al. [58] found a generalized positive trend but significant negative trends in the SW.The Atlantic region showed negative trends in none of these works.
The analyses of NDVI significant trends from 2001-2016 in mainland Spain and the Balearic Islands showed that the phytoclimatic and CLC classes had different percentages of significant negative and positive trends.The land-cover and phytoclimatic types with a higher number of NDVI positive or negative trends are, in part, explained by climatic factors and, in part, by anthropological factors, as other authors reported in their articles in this field.In more detail, other authors reported that worldwide vegetation dynamics are generally driven by climate, in particular, by precipitation, incoming radiation, and temperature [1,4,81,[93][94][95][96][97].Additionally, nutrient availability and short-term natural and anthropogenic disturbances (e.g., fires, logging, insect epidemics, land-use change, and agricultural management) can be crucial at various spatiotemporal scales [5][6][7]98].Hill et al. [52] and Serra et al. [99] reported that the rural exodus that occurred in the last 30 years of the 20th century provoked land-cover changes from crops to transitional wood-shrub zones and from burnt areas to transitional wood-shrub sites.These changes to shrubs resulted in an increase in the NDVI values in the last 30 years of the 20th century, and these results coincide, in part, with those of our significant positive NDVI study area.However, they did not explain our negative NDVI trends.
The negative trends found in our study area could be caused by the state of the North Atlantic Jet.It is known that there were anomalies in the trajectory of the North Atlantic Jet that provoked more severe droughts in some critical months of the year for the growth of the vegetation [100,101].Moreover, some variability in the precipitation cycle was observed in Spain.In some areas, the precipitation decreased in Spain during the period that we took into account [102][103][104].However, it is unclear whether there was a general negative trend in the precipitation of the last 20 years in the Mediterranean areas [105].What is clearer than the decrease in precipitation is the global increase in temperature since the last century [1,57,[102][103][104].In Khorchani et al. [57] or Papagiannopoulou et al. [1], it was reported that, in a general pattern in the northern Peninsular of Spain, the limiting factors on vegetation growth are the temperature and the irradiation of the sun, while, in the center and south of this zone, the limiting factor is the precipitation.In other studies, we found that, in general, when the limiting factor was the precipitation, an increase in temperature translated into a decrease in NDVI value.Meanwhile, when the limiting factor on vegetation growth was the temperature, an increase in temperature resulted in a decrease in NDVI value or a positive NDVI trend in our study area [57].This was not entirely in accordance with our results.One example was found in our Atlantic region, where we did not expect to find a major decrease in NDVI values, as we saw in the period from 2001 to 2016.
The percentage of significant NDVI trends observed by phytoclimatic types revealed that there are statistical differences and that, in some of these types, the positive trends predominated and, in others, the significant negative trends predominated (see Figure 8).In our study area, the vegetation tendencies of the Mediterranean types showed different behaviors.It is interesting to remark that the warmer Mediterranean phytoclimatic types (tf > 9.5 • : IV(III) and IV2), the warmer types between Mediterranean and nemo-Mediterranean (tf > 7.5 • : IV(VI)2 and VI (IV)3), and the warmer nemoral types showed negative trends.This pattern found in our study must be explained by the fact that an increase in temperature in the summer [57], in addition to some droughts reported in those areas in which these phytoclimatic types are located, resulted in a decrease in photosynthetic activity [106,107].These results coincide with results from other authors who present models that revealed how sub-Mediterranean zones would suffer most with an increase in temperature in addition to the increased frequency of droughts in the summer [65].
The nemo-Mediterranean types presented higher pixels with more significant positive trends than negative trends (see Figure 7a).The differences in the mean NDVI nemo-Mediterranean types between the negative and the positive trends were high, as can be seen in Figure 7a.Therefore, this could be explained by the fact that the precipitation is not a limiting factor on these types (see Table 2) and an increase in the temperature could be related to the increases in vegetation activity, as other authors mentioned for the temperate zones of the world (e.g., References [21, 23,57,[78][79][80]82,108]. The same occurred for the nemoral phytoclimatic types that are more similar to the nemo-Mediterranean types, i.e., the VI(V) and the VI phytoclimatic class-the explanation seems to be the same.However, we could not note this pattern in type VI(V).The last mentioned phytoclimate is geo-located in the northwest of mainland Spain (see Figure 3) and coincides with this North Atlantic zone in which we did not expect to find a predominance of negative NDVI trends.Those trends might be explained by the anomalies in the trajectory of the North Atlantic Jet [101], which could be related to a decrease of the cloudy atmosphere or which could provoke more severe droughts in some critical months of the year for the growth of vegetation.Other factors could be land degradation and insect plagues, in addition to the high rates of land abandonment and fragmentation, severe overgrazing, and the increase in wildfire occurrence in recent decades [109][110][111].Furthermore, in this area, other studies suggest that the Fagus sylvatica species in the northern zones of Spain are less resistant to climate change than those located in the southern regions of Spain [112].Consequently, it should be noted that we did not observe a long period, and future research on that area might be necessary to evaluate and monitor vegetation in that Atlantic zone, which is less studied than drylands or semi-arid zones.
The Oromaticoid phytoclimatic types are those that match with the highest mountain areas of our study site.In these zones, we found a predominance of significant positive NDVI trends.These results match with those obtained by other authors for mountains areas [88,89].
By land cover, we also found statistical differences in the percentage of positive and negative trends of NDVI in the period studied.Three agriculture covers-permanently irrigated areas (212), rice fields (213), and fruit trees (222)-showed a large percentage of both negative and positive trends.These land covers depend greatly on the availability of water, and different restrictions were imposed regionally.However, when water is not a problem, irrigated vegetation can take advantage of higher temperatures, as occurred for the broad-leaved (311) and coniferous masses (312) (see Figure 7b).Regarding these two CLCs, to explain the decrease in the vegetation pattern, it is essential to mention that the percentages of negative trends shown in Figure 7b were negative overall because of the local negative trend percentages of the two classes.Additionally, the phytoclimate IV(VI)2 belongs to the Mediterranean class, and type VI (V) was in the nemoral class (higher than 21% in both cases).Therefore, the broad-leaved and coniferous forests from the Mediterranean north coastline and the Atlantic north coastline could be explained by the aforementioned effect of diseases and defoliation of trees in those areas.
The cover types that had a higher percentage of positive trends than negative trends were those belonging to vineyards (221), olive groves (223), and beaches, dunes, and sand (331).Vineyards and olive groves in mainland Spain and the Balearic islands had permanently increasing irrigated areas; thus, the increase in temperature could be translated into an increase in vegetation activity.
Finally, it is crucial to remark that our HP model results show that the variation in the negative percentage was higher, as explained by the Allué phytoclimatic types, while the positive trends were mostly explained by the land-cover types.As in the general statistical model, we had more positive than negative trends.The results show that, without making the distinction between negative and positive trends, the land-cover types explained most of the variability in the percentage of significant trend pixels.As we noted with our inferential test and the explanations found in the present paper, the land-cover types and the phytoclimatic classes showed marked differences in their percentages of significant trends that could be explained by different climate and human causes.We determined that the remote-sensing datasets are among the best options to assess the trends in the vegetation dynamics.Thus, it is crucial to further research this topic to understand and effectively manage our territory in the face of future challenges.

Conclusions
In this study, we evaluated the NDVI trends from 2001 to 2016 in mainland Spain and the Balearic Islands, a region with three biogeographic zones and with great climatic and vegetation variability.The most important findings of this study are as follows:

•
Significant trends showed a non-random distribution.We found great differences between biogeographical regions, and phytoclimatic and land-cover types.

•
The unexpected negative trend of the Atlantic region could be a response to climate variability, with a global temperature increase and more precipitation variability.

•
Land cover and phytoclimates explained approximately 35% of the variance.Land cover explained most of the positive trends and phytoclimates explained most of the negative trends.

•
Warmer types of each general phytoclimate showed a negative vegetation dynamic.For instance, from the nemo-Mediterranean type, the warmest type showed negative trends, while the other types showed positive NDVI trends.

•
Forest land cover had negative trends, especially for phytoclimates IV2, IV(VI)2, IV(VI)3.and VI(V).For grasslands, only type IV(VI)1 showed clear negative values.It is interesting to highlight that, for some phytoclimates, most land cover had the same behavior as NDVI trends, but, for others, the behavior was the opposite.

28 Figure 2 .
Figure 2. Map of Co-Ordination of Information on the Environment (CORINE) Land Cover (CLC) 2006 classes at level 1 of disaggregation from our study area.

Figure 2 .
Figure 2. Map of Co-Ordination of Information on the Environment (CORINE) Land Cover (CLC) 2006 classes at level 1 of disaggregation from our study area.

Figure 3 .
Figure 3. Map of Allué phytoclimatic types or zones in mainland Spain and the Balearic Islands.

Figure 3 .
Figure 3. Map of Allué phytoclimatic types or zones in mainland Spain and the Balearic Islands.

Figure 4 .
Figure 4. Data processing scheme.TS refers to Theil-Sen slope data and MK refers to the Mann-Kendall test.

Figure 4 .
Figure 4. Data processing scheme.TS refers to Theil-Sen slope data and MK refers to the Mann-Kendall test.

28 Figure 5 .
Figure 5. Spatial pattern that resulted from intersecting the Theil-Sen slopes with the Mann-Kendall test (z-statistic), testing annual averages from 2001-2016 Moderate Resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) values.The red colorshows the pixels with negative Theil-Sen slopes and a Mann-Kendall z-statistic below 0.1 (significant trend).The green color shows the pixels with positive Theil-Sen slopes and a Mann-Kendall z-statistic below 0.1 (significant trend).Finally, the black color shows the pixels with a Mann-Kendall z-statistic higher than 0.1 (non-significant trend).

Figure 5 .
Figure 5.Spatial pattern that resulted from intersecting the Theil-Sen slopes with the Mann-Kendall test (z-statistic), testing annual averages from 2001-2016 Moderate Resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) values.The red color shows the pixels with negative Theil-Sen slopes and a Mann-Kendall z-statistic below 0.1 (significant trend).The green color shows the pixels with positive Theil-Sen slopes and a Mann-Kendall z-statistic below 0.1 (significant trend).Finally, the black color shows the pixels with a Mann-Kendall z-statistic higher than 0.1 (non-significant trend).

Figure 6 .
Figure 6.Interval plots (95% confidence interval (CI)) from the variance analyses.The x-axis labels are the Allué phytoclimatic regions (a) and the CORINE Land-Cover code; (b) the y-axis shows the percentage of significant (p < 0.1) trend pixels.

Figure 6 .
Figure 6.Interval plots (95% confidence interval (CI)) from the variance analyses.The x-axis labels are the Allué phytoclimatic regions (a) and the CORINE Land-Cover code; (b) the y-axis shows the percentage of significant (p < 0.1) trend pixels.

Figure
Figure7ashows slightly higher negative trend percentages than positive trend percentages in type IV(III), the most arid from the Mediterranean type; types IV(VI)1 and 2 (also from the Mediterranean type, but the least arid in this category); type VI(V) from the nemoral class; and type VIII(VI) from the Oroboreal type.Figure7bshows that the classes of fruit trees (222), pastures (231), complex cultivation patterns (242), and broad-leaved (311) and coniferous (312) forests had a higher negative NDVI percentage from 2001-2016 than positive NDVI percentage.This negative percentage was higher than one-third the positive percentage in the case of the broad-leaved and coniferous masses.

Figure 7 .
Figure 7. Interval plots (95% CI) from the variance analyses.The x-axis labels are the Allué phytoclimatic regions (a) and the CLC codes; (b) the y-axis shows the percentage of significant (p < 0.1) positive (in green color) and negative (in red color) trend pixels.

Figure 7 .
Figure 7. Interval plots (95% CI) from the variance analyses.The x-axis labels are the Allué phytoclimatic regions (a) and the CLC codes; (b) the y-axis shows the percentage of significant (p < 0.1) positive (in green color) and negative (in red color) trend pixels.

Figure 8 .
Figure 8. Percentage of pixels with a significant trend by CORINE Land-Cover class (a) and by phytoclimatic type (b).Negative is represented by negative (in red color) and positive is green in color.Lines are mean values.Difference in percentage (y-axis) between the observed and expected negative NDVI, with negative (in red color) and positive (in green color) short-term trend pixels from our study area by CORINE Land-Cover class (c) and by phytoclimatic type (d).

Figure 8 .
Figure 8. Percentage of pixels with a significant trend by CORINE Land-Cover class (a) and by phytoclimatic type (b).Negative is represented by negative (in red color) and positive is green in color.Lines are mean values.Difference in percentage (y-axis) between the observed and expected negative NDVI, with negative (in red color) and positive (in green color) short-term trend pixels from our study area by CORINE Land-Cover class (c) and by phytoclimatic type (d).

Figure 9 .
Figure 9. Difference in percentage (y-axis) between the observed and expected NDVI negative (in red color) and positive (in green color) short-term trend pixels from our study area, by phytoclimatic type, within the 311 (a), the 312 (b), the 313 (c), and the 321 (d) CORINE Land-Cover classes.

Figure 9 .
Figure 9. Difference in percentage (y-axis) between the observed and expected NDVI negative (in red color) and positive (in green color) short-term trend pixels from our study area, by phytoclimatic type, within the 311 (a), the 312 (b), the 313 (c), and the 321 (d) CORINE Land-Cover classes.

Table 1 .
Codes of Co-Ordination of Information on the Environment (CORINE) Land Cover (CLC) 2006 classes: levels 1 and 3 of disaggregation.

Table 2 .
Summary of the phytoclimatic

Table 3 .
Percentage of pixels, by biogeographical region, that exhibit positive and negative significant normalized difference vegetation index (NDVI) trends.