Identification of Years with Extreme Vegetation State in Central Europe Based on Remote Sensing and Meteorological Data

(1) Eötvös Loránd University, Department of Geophysics and Space Science, Pázmány P. st. 1/A, H-1117 Budapest, Hungary; (2) Croatian Forest Research Institute, Cvjetno naselje 41, HR10450 Jastrebarsko, Croatia; (3) Czech University of Life Sciences, Faculty of Forestry and Wood Sciences, Department of Forest Protection and Entomology, Kamýcká 129, 165 21 Prague 6, Czech Republic; (4) Eötvös Loránd University, Department of Meteorology, Pázmány P. st. 1/A, H-1117 Budapest, Hungary; (5) Eötvös Loránd University, Faculty of Sciences, Excellence Center, H-2462 Martonvásár, Brunszvik u. 2.

Meteorological parameters exhibit strong intra-annual variability as well with extreme events like heatwaves, drought spells, flash floods, extreme cold periods etc. In meteorology, quantification of the extreme events has long tradition [1, 4,5]. Climate indices are calculated typically from daily meteorological records to estimate e.g. strength of heatwaves, number of days with extreme precipitation, length of dry periods etc. [1, 6,7].
Plant processes and productivity have strong economic impacts due to its direct relationship with crop yield, wood production, animal welfare, fodder quantity and quality, ecosystem services, etc. [8]. Thus, quantifying plant production and understanding the processes behind variability of plant status on annual scale is of high importance [9,10].
Studies focusing on ecosystem processes have problems with the existing climate indices and annual means due to the complex interactions between meteorological parameters and the vegetation state [11]. Intra-annual variability of the meteorological parameters might exert strong impact on plant growth and productivity [12]. Actual vegetation state is the result of beneficial and negative effects from the past. In other words, vegetation status is the integrator of all effects from soil water status, temperature, variability of radiation, etc. Therefore extreme indices and annual means are not directly applicable to study plant growth and explain its interannual variability.
One logical way to study the impact of climate variability on plant processes is the effect-based approach, where the plant response is characterized first, then the meteorological cause of the effect is sought (see e.g. [13,14]). Due to economic and political reasons these studies should practically focus on the country scale (see e.g. [15]). However, quantification of plant status on large spatial scales is not straightforward due to the small spatial representativeness and amount of observation based, in situ data (e.g. biomass data, phenology observations, leaf area index measurements, eddy covariance data, crop census data, etc.).
Remote sensing (RS) provides a very convenient solution to this problem. Data obtained from sensors on board Earth Observing satellites like the Moderate resolution Imaging Spectroradiometer (MODIS) provides unique information [16] which can be linked to plant status via vegetation indices like Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI) [17], Leaf Area Index (LAI), Fraction of Photosynthetically Active Radiation (FPAR) [18] or modelled plant production like Gross Primary Production (GPP), and the Net Primary Production (NPP) [19]. They are produced at convenient spatial scale (with finest 250 m -1 km), and today already comprise an exceptional 17-yearslong data record.
The spatial resolution of the MODIS Collection 6 products for NDVI, EVI, LAI, FPAR, GPP, and NPP is 500 m or finer (for NDVI and LAI with 250 m resolution) which is sufficiently small to isolate large number of pixels that contain a single land cover (LC) type. This enables that the analysis of plant behaviour can be studied for each land cover type separately, which is needed due to the substantially different dependence of plant functional types on the environmental conditions [20]. However, the classification of pixels into LC categories is not trivial due to land use change and land cover classification errors. For best results, other sources of land cover type are needed, like national/regional forest management maps, agriculture land use monitoring systems for payment of subsidies (e.g. in EU countries), or CORINE Land Cover database [21]. Here we propose a novel method that would enable selection and identification of the anomalous years based on RS data and RS-derived products (NFVI, EVI, LAI, FPAR, GPP, and NPP) according to different land cover types and spatial coverage. This method overcomes the issues that are related with the study of the meteorological anomalies alone. Based on the area for which the method is implemented, it could help to assess differences at spatial and temporal scale for a given LC within the area of interest. It can be applied at a national level to study the behaviour of each LC type as well as to assess the magnitude of the anomaly that has occurred in terms of "repeat time".
In the present paper country averages are studied to get robust results and to support future research. At smaller spatial scales climate is not necessarily the dominant driver of plant development and vigour. For example forest production can be adversely affected due to pathogens or insects [11]. Cropland management clearly affects plant status at smaller scales; for example winter wheat is typically harvested in the region by June-July which affects overall greenness. Management practices like fertilization, harvest and irrigation, and also soil type clearly modulate the effect of meteorological conditions in grasslands and croplands. Country averages by land cover types are expected to be driven by climate fluctuations and provide robust and clear signal about the overall productivity in the given country (see e.g. [15,22] for such approaches).
The aim of the paper is to present a method to select years that can be characterized as anomalous based on observed plant status and greenness. Using multiple vegetation indices we also test the similarity/dissimilarity of the different vegetation metrics in terms of their usability to detect anomalous plant status. To get robust and easily interpretable results we used country-means. Bosnia and Herzegovina, Croatia and Hungary (and in Supplement Czech Republic, Slovakia and Slovenia) were selected to represent Mediterranean, continental and alpine climates. The reason behind such choice is the relatively small geographical extent and the fact that extreme weather events typically affect large areas and sometimes the whole country [23]. Countries with transitional economies are particularly vulnerable to climate extremes [10]. It means that understanding the cause-effect relationships might support the prevention of adverse effects on plant state and this can support economic growth and human welfare.

Study area
The target area of this study is the broader region of the Carpathian Basin, determined by the coverage of the applied meteorological dataset (see below). Bosnia and Herzegovina, Croatia and Hungary were selected to represent Mediterranean and continental climate as well, with high biodiversity and variability in the meteorological conditions. Results for Czech Republic, Slovakia and Slovenia are presented in the Supplemental Material. The reason behind chosen countries is the spatial domain of the applied meteorological database, covering fully only the selected countries.

Meteorological Database
To investigate the effects of weather variability on the vegetation activity and greenness, we used the FORESEE database [24]. This freely available meteorological database contains observed and projected daily maximum/minimum temperature and precipitation fields for Central Europe on a regular grid with a spatial resolution of 1/6° × 1/6°, covering the years between 1951 and 2100. For the 1951-2014 period, FORESEE provides observation-based interpolated meteorological fields for the wide region of the Carpathian Basin.
In order to match the temporal and spatial resolution of the applied remote sensing datasets and of the FORESEE, 8-day mean temperature values and precipitation sums were calculated on the finer grid of the MODIS products based on the methodology of Kern et al. [25]. Country averaged mean and anomaly values were calculated for the above mentioned six countries within the study area.
Mean temperature and precipitation fields (respectively) calculated for the entire domain area of the FORESEE database are presented in Figure 1 and 2 for the study period of 2000-2014. Yearly anomaly fields, shown in Figure 3-6, are illustrations of years which might be considered as extreme from the meteorological perspective. Figure 3 and 4 show the anomaly maps of temperature and precipitation (respectively) for 2011, while Figure 5 and 6 shows the anomaly maps of temperature and precipitation (respectively) for 2014 (note that the temperature and precipitation anomaly maps have different legends).

Vegetation Related MODIS products
In the present study we used the latest version (Collection 6) MODIS NDVI, EVI, FPAR, LAI, GPP, and NPP with 500 × 500 m spatial resolution, as part of the MOD13A1, MOD15A2H, MOD17A2H and MOD17A3H official products [26,27] derived from the measurements of MODIS sensor on board satellite Terra. The longest possible datasets (covering the period of 2000-2014) were chosen to match the availability of the MODIS products and the temporal coverage of the applied meteorological dataset.
Quality filtering of each dataset was performed using the quality flag information included in the datasets, based on the method described by Kern et al. [25]. Besides this, to filter out unrealistic sudden increases and decreases in the state of the vegetation the so called Best Index Slope Extraction (BISE) method [28] was applied afterwards on pixel level for NDVI and EVI data [25]. While FPAR, LAI and GPP has the same 8-day temporal resolution, NDVI and EVI are 16-day composite products, therefore temporal resampling based on the Julian date information included in the MOD13 datasets was also necessary to create NDVI and EVI dataset with the same 8-day temporal resolution [25].
It has to be noted that although most of the MODIS products are state-of-the-art, standardized, well documented data sets, the quality of the model based products is expected to be lower and that affects their applicability [19]. Most notably, MODIS annual NPP products suffers from a major error (with large unexpected positive bias for the first four years (2000)(2001)(2002)(2003) as it is presented in Figure 7 and should not be used in the FORESEE domain. Therefore, the anomaly values from the yearly NPP products were calculated relative to the period of 2004-2014.

Methods and Metrics for Defining anomalous Years
From the quality (and in the case of MOD13A1 data BISE-filtered) 8-day temporal resolution datasets country averaged multiannual means and yearly anomaly values were calculated for the whole vegetated area of the selected countries within the Carpathian Basin. Land cover specific country averaged values were also derived for broadleaf and coniferous forests, croplands and grasslands. The yearly anomalies were calculated without the first 40 and last 40 days of the year to avoid any misleading result originating from the effect of snow cover. Using the derived yearly anomaly values, relative anomaly values were calculated by dividing the yearly anomaly values by the maximum of the absolute anomaly values during the investigated 2000-2014 period.
To analyse the spatial distribution of the yearly anomaly fields, we calculated the 0.5, 2, 9, 25, 50, 75, 91, 98 and 99.5 percentiles values of all yearly mean anomaly values for all land cover specific pixels separately within a given country. The selected percentiles define thresholds for the classification of anomalies into the following 10 categories: The derived percentile values are appropriate to des-cribe the total distribution of the anomaly values of a given land cover type (within a given country) during the study period of 2000-2014. Based on the proposed categorization we calculated for every year separately the number of the pixels (as the percentage for a given LC type) within all the created percentile ranges. Note that there might be situations when e.g. strong positive anomaly was not occurring in the entire time period for a given LC type at all. In this study the maximum anomaly is used independent of its magnitude.
The advantages of the proposed categorization is the following: (1) the method we applied is the same for all LC types, but it yields LC specific results; (2) it enables the identification of areas within the selected geographical area (e.g. country), where mean anomaly is more positive or more negative with respect to the country average, indicating that <0.5p most extreme negative anomaly, 0.5p -2p extremely negative anomaly, 2p -9p strong negative anomaly, 9p -25p negative anomaly, 25p -50p } normal range 50p -75p 75p -91p positive anomaly, 91p -98p strong positive anomaly, 98p -99.5p extremely positive anomaly, >99.5p most extreme positive anomaly. some areas are more/less productive, more/less prone or more/less sensitive to meteorological anomalies (where a special attention should be paid to croplands due to the possible yearly crop type change); (3) it is appropriate to be done for any other vegetation related characteristics.
In this paper we propose the following classification of the years, which is based on the percentage of the area showing positive and negative anomalies of a given vegetation characteristics (NDVI, EVI, FPAR, LAI, GPP, and NPP). Using the yearly anomaly fields of given vegetation 12 [30] were used to identify pixels with high probability of being broadleaf forests, coniferous forests, croplands and grasslands. From the five types of land cover classification contained in the MCD12Q1 product [31] we used the Type-1 (International Geosphere Biosphere Programme -IGBP) to identify croplands and grasslands and the Type-3 (MODIS-derived LAI/fPAR scheme) classifications for broadleaf and coniferous forests. The categories of the land cover classification for the used MODIS pixels are given in the translation matrix between the different LC schemes ( Table 2). The reason behind the usage of Type-3 for forests was that the widely used Type-1 IGBP land cover classification suffers from well-documented errors [32,33] with 75% overall accuracy [31] including misclassification of the mixed forest types as well [34]. We found that the misclassification is especially evident for the lowland forests along rivers, for example like the Gemenc forest by the Danube in Hungary, or the river Spačva basin in eastern Croatia, extending also to Vojvodina in Serbia. According to taBlE 1. Thresholds to the selection of the anomalous years based on the percentage of the affected area characteristic, for a specific land cover type, the chosen year can be graded based on the thresholds presented in Table 1. It should be noted that the grade for a given year at a given location depends on the selected reference period (e.g. MODIS era, in our case 2000-2014) and the area representing a domain (e.g. region, country or a continent).

Synergy of the MODIS and CORInE land Cover Databases
In order to study the response of the various vegetation types to meteorological anomalies we distinguished the main land cover types based on the synergistic use of two land cover datasets. The so called MCD12 land cover products based on MODIS observations [29] and the MODIS Type-1 land cover classification both forests have been categorized as mixed forest, while they are broadleaf forests in reality. In fact Spačva forest with its area of 40 kHa size is the largest complex of pedunculate oak forest in Europe [35]. In contrast to this, the Type-3 classification doesn't have mixed forests category, resulting in all forested pixels categorized as either broadleaf or coniferous category.
We have to note that in Type-3 classification croplands and grasslands are categorized together as "Grasses/Cereal crops", therefore Type-3 classification cannot be used to discriminate croplands and grasslands. From the yearly MCD12Q1 land cover datasets with 500 m × 500 m spatial resolution covering the period of 2001-2013 we selected the pixels (hereafter called stable LC pixels) which had no land cover change during the entire period. The usage of the temporal stability selection criteria excluded pixels which underwent any kind of land cover change during the investigated period, either actual, or resulting from the error in classification. Implementation of such a strict rule resulted with dropping significant number of the pixels. The percentage of stable LC pixels was only ~56% of the whole study area based on Type-1 classification. Specifically, 32%, 34%, 53%, 69%, and 28% of the pixels remained for the categories of broadleaf forests, coniferous forests, mixed forests, croplands, and grasslands, respectively. Based on Type-3 classification the percentages of stable LC pixels in the whole study were 48% and 42% for broadleaf and coniferous forests, respectively. The percentages were calculated relatively to the mean number of the pixels during the 13 years of the total dataset.
In order to increase the reliability of the applied LC categorization we used the CORINE land cover dataset as a reference LC dataset for the year 2012. The accuracy of CORINE is 87.0 ± 0.7% [21], which is significantly better than the 75% accuracy of MODIS [31]. Using GIS software [QGIS 2.16], we intersected the vector-based CORINE layer with the grid of the MODIS pixels at 500 m × 500 m spatial resolution and obtained the share of every of the 44 CORINE land cover types present within each of the MODIS grid-cells (i.e. pixel). Due to differences in land cover categorization between MODIS and CORINE, we used translation matrix (Table 2) to unify the classification. Based on the area share information of CORINE we selected so called CLC2012 pure pixels for deciduous forests, coniferous forests, croplands and grasslands, which contained at least 90% area share from the given LC type (Table 2, column 1). In the case of croplands and grasslands the 90% threshold was applied for the sum of the area shares of the different CLC categories (separately for croplands and grasslands) listed in Table 2 (column 3).
Using the derived set of MODIS stable pixels and CORINE pure pixels we selected "reliable" pixels which were classified by both land cover database to the same vegetation type of broadleaf forests, coniferous forests, croplands and grasslands. Figure 8 shows the location of the remaining reliable pixels of the main land cover types based on MODIS IGBP (MCD12 Type-1), LAI/FPAR scheme (MCD12 Type-3) classifications and CLC2012 database. This land cover map shows the location of the pixels which were finally used in the present study as reliable pixels, being constant during the study period and having probably right classification. Based on the information of the CORINE database regarding the percentage values of the presence of a given land cover type the mean percentage for the selected broadleaf forests, coniferous forests, croplands and grasslands pixel was 99.0%, 99.0%, 99.3% and 98.2%, respectively. The problem of using the MODIS Type-1 classification for forests is illustrated in Figure 9. The presented map shows the location of the pixels which are classified as broadleaf forests and coniferous based on the MODIS Type-3 and CORINE classifications, but not based on MODIS Type-1 classification. The number of the forested pixels (in the study area) which are not represented well in Type-1 classification is considerable: 59.5% and 66.3% for broadleaf and coniferous forests, respectively, relatively to the number of the derived reliable broadleaf and  coniferous pixels, which were selected finally to our further investigations.
After the synergetic use of the two MODIS datasets with CORINE dataset the number of the selected pixels was 6.3% for broadleaf, 3.1% for coniferous forest, 17.1% for cropland and 0.1% for grassland pixels relative to the all vegetated pixels in the whole domain. It means, with the final pixel selection we were taking into account 26.6% of the whole domain in a total value. The country specific pixel numbers are presented in Table 3.

RESUltS anD DISCUSSIOn problems with the Selection of Extreme Years Based Solely on Meteorological Data
In order to study the dependence of the strength of meteorological anomalies on the selected integration time period, first we present country averaged mean temperature and precipitation anomalies in Figure 10 (left and right side images, respectively) during 2000-2014 for Bosnia and Herzegovina, Croatia and Hungary (and in Figure S1 in the Supplementary Material for Czech Republic, Slovakia and Slovenia).
While the columns show whole year anomalies in Figure  10, the curves indicate periods with various length (starting at different date and all ending with 7 th of October within the year), representing gradually only the predominant part of the growing season. 7 th of October was used as the end of the integration period as we can hypothesize that meteorology in the dormant season affects plant state to a lesser extent. The beginning was varied with ~2 weeks periods and ranging from 1 st January to 26 th June.
The figure shows that even if there is an extreme temperature or precipitation anomaly for a given year, it does not necessarily mean that the vegetation was exposed to similar meteorological anomalies during the growing season. A good example is 2003, when the countryaveraged annual temperature anomaly was negative in the presented countries (especially for Hungary), but a strong, mostly positive anomaly was detectable for the growing season. On the contrary, year 2004 or 2005 (with similar negative annual temperature anomalies) showed only negative temperature anomalies during the shorter time periods until the end of the growing season. In terms of annual temperature anomaly year 2007 was very similar to year 2000 for all three selected countries, but with much higher intra-annual variability, on average having cooler than usual period in summer and early autumn. On the contrary, in 2012 (with annual temperature anomaly similar to the one in 2007) the summer-early autumn period was characterized with the highest positive anomaly during the year, when the mean temperature of the countries was continuously higher than the average. Finally, year 2014 is worth mentioning as well. Though it had the largest annual anomaly, its temporal evolution was not consistent. Precipitation showed similar features, however the intra-annual variability was much lower. It is important to mention the two consecutive years of 2010 and 2011, which was characterized by very strong positive (2010) and then negative precipitation anomaly (2011) all over the Carpathian-Basin.
These results illustrate the problems in defining temperature and precipitation anomalies as they strongly depend on the selected integration period. In other words, we cannot unambiguously select extreme years based on the yearly (or monthly/seasonal) mean meteorological conditions that can be related with e.g. vegetation state, crop yield, forest productivity, outbreak of insects or other phenomena. To find the extreme years in plant greenness and describe the response of the vegetation greenness, it is more straightforward to study anomalous behaviour of the plant state (or other phenomena that is related with ongoing climate anomalies).

Selection of anomalous Years Based on the Overall Impact on Vegetation
Yearly relative anomalies of the vegetation related characteristics (such as NDVI, EVI, FPAR, LAI, GPP, and NPP) during 2000-2014 are presented in Figure 11 for croplands and grasslands. Figure 12 shows the same for broadleaf and coniferous forests for Bosnia and Herzegovina, Croatia and Hungary ( Figure S2 and S3 in the Supplementary Material illustrate results for the Czech Republic, Slovakia and Slovenia). The yearly, land cover specific anomalies of the meteorological variables (temperature and precipitation) are also presented in Figure 11 and 12.
Based on the results, the relative anomalies for the different vegetation characteristics could be described as quite consistent for croplands and grasslands, showing the same direction of relative anomaly in most of the cases. It is also notable that GPP and NPP show opposite character for taBlE 3. Number of the pixels at 500 m × 500 m spatial resolution which complied with the selection criteria separately to the main LC types and countries (selection criteria: MODIS LC type for a given pixel did not change during 2001-2013 and its LC classification with 90% share corresponds to that of CORINE 2012 -see Table 2  For a given land cover type in some area (country) within the reference period 2000-2014, the calculated relative anomalies facilitate straightforward comparison between years, as well as between different land cover types within the same year. For example, we can select years when croplands and grasslands were both anomalous (like 2000, 2003 and 2012 with negative anomalies, or 2010 and 2014 with positive anomalies), but in some years (such as 2001 and 2004) croplands and grasslands showed different response to the environmental conditions. These differences might be related to land cover type specific features, species distribution and human disturbances (harvest, fertilization), and maybe other factors like soil water holding capacity, etc. Additionally, croplands in some countries cover C4 plants with higher drought resistance, which might affect the results. Also, croplands and grasslands might exhibit different critical time periods when the ongoing weather affects their productivity [36]. Clearly, it is not only the yearly anomalies of the meteorological parameters that affect the plant status, but the within-year temporal pattern of the meteorological elements. Note that the number of pixels (see Table 3) was low for grasslands in all countries except Bosnia and Herzegovina and Croatia, which can affect the reliability of the results for grasslands in those countries.
While in the case of croplands and grasslands the vegetation related characteristics are mostly consistent, forests reveal different behaviour ( Figure 12). As it is obvious from the figure, in the case of forests the different characteristics are not showing similarity for the anomalies, as it was found for croplands and grasslands. It is especially true for GPP and NPP, which are in fact models and not indices that are derived directly from the reflectance data [19]. It means that GPP and NPP results must be treated with caution. For example, in 2012 the reason behind the large disagreement between the mean anomaly of NDVI, EVI, LAI and FPAR, and the mean behaviour of GPP and NPP might be related with the extreme meteorological conditions during the growing season, causing significant drought in the Carpathian-Basin and a consequently simulated lower productivity by the GPP-model. In addition, forests with typically deeper root zone do not react in the same way to changing meteorological conditions as the shallow-rooted croplands and grasslands [16]. Forest productivity depends on the available soil water content in the topsoil to a lesser extent, because deeper groundwater supplies can be accessible to trees. Trees store a relatively large amount of non-structural carbohydrate -which is the product of the previous year -to mitigate negative effect of shortage in nutrients and photosynthetic carbon uptake. Therefore they are likely to be less affected by the shorter term weather anomalies, but showing stronger exposure to the longer term changes. Based on the observed results, we could not select the same years as anomalous years for herbaceous vegetation, which were highlighted in the case of croplands and grassland over the studied countries. The role of the carbohydrate reserves seems to be important which might impose lagged effect [13]. This can be recognised in 2013 as the effect of the previous, extreme year (evident in strong negative anomaly for croplands, see Figure 11) which has contributed to a negative anomaly in broadleaf forests (see Figure 12).
Co-variation between the meteorological and vegetation anomalies does exist and it can give us some basic information at the annual scale. However, the main aim of the present study was not to execute in-depth analysis of cause-effect of environmental variables on the detectable plant state anomalies, but to select anomalous years. In order to quantify their co-variance Pearson's r values (correlation coefficients) were calculated and compared. The correlation coefficients between the yearly anomalies of the vegetation characteristics and the yearly meteorological anomalies were calculated for the different plant types separately. Table 4 present the correlation coefficients calculated between the NDVI anomalies and the temperature and precipitation anomalies, where the statistically significant (P<0.05) values are indicated. In the case of temperature, the calculated r values show the strongest correlation for forests (especially with broadleaf forests), while in the case of precipitation the largest correlation is present for croplands and grasslands. Both the direction and the strength of the correlation depend strongly on the land cover type, indicating complex relationship between the vegetation and meteorological conditions. It also implies that large scale and land cover specific studies should not be made together, but separately and at finer temporal and spatial scale.

Selection of anomalous Years Based on the Magnitude and Spatial Extent of nDVI anomalies
The investigation of the country averaged mean yearly anomaly does not provide information about its spatial distribution. Therefore, pixel-level investigations were performed using percentiles analysis for the area of a given land cover type affected with severe anomaly (see Section 2.4 in the Materials and Methods). Here we present the results based only on the NDVI, while the results for all other metrics are provided in the Supplementary material (S6, S8, S10, S12).
Different land cover types react with varying intensity to meteorological anomalies as it can be inferred by careful inspection of Figure 13. It shows that, for example, broadleaf forests (top left) maintain rather stable mean annual NDVI and the corresponding anomalies are in the range of +/-0.05 of the mean in 99% of the cases. Coniferous forests ( Figure 13, top right), on the other hand, show higher variability of the mean NDVI compared to broadleaf forests, but the magnitude is still lower than for grasslands ( Figure  13, bottom right), and considerably lower than in the case of croplands ( Figure 13, bottom left).
Differences between countries exist, but they are rather small, with increasing discrepancies toward the extremes (2p and less or 98p and more). However, this could be the result of relatively small number of pixels with extreme values and should be interpreted with caution. This is particularly the case for grasslands, where the number of pixels was large enough only for Bosnia and Herzegovina and Croatia (see Table 3). We would also like to point out the differences among countries in the case of the croplands, and in particular to Slovenia. It can be seen that the high values for the anomalies in NDVI are occurring in Slovenia less frequently than e.g. in Slovakia or in Czech Republic. The reason behind this could be climate (relatively large amount of precipitation in Slovenia -see Figure 2) but also it could be due to the differences in crop type or cultivation method. Somewhat similar behaviour can be observed also for the coniferous forests in Hungary, where the observed magnitude of their positive anomaly is apparently lower than for coniferous forests in the other countries. However, here the caution is needed, because it could be an artefact due to the modest number of available pixels and this would require further investigation. In any case, these results corroborate the logic of country level analysis, where the considered area is large enough to reduce the effects of random noise (caused e.g. by species or management differences), but small enough to retain the information on the exiting differences among the countries. The relatively uniform distance between the values of the mean anomalies for the shown percentiles ( Figure 13) indicates that the choice of the percentiles, representing the border between the classes of anomalies (see Section 2.4 in the Materials and Methods for the list) has been made properly.
In an effort to identify specific years that can be characterized as extremes according to a given RS characteristic, we performed a land cover and country specific analysis. The results for Bosnia and Herzegovina, Croatia and Hungary based on NDVI are given in Figures  14-16 and Table 5   Here the case of Slovenia is somewhat easier to interpret. Namely, in January-February 2014 large part of Slovenia and part of Croatia was experiencing severe case of freezing rain that caused extensive damage to forests of the affected region [37,38].
Regarding the negative extreme in the case of forests, the situation is more complex. Unlike for croplands and grasslands, year 2003 was not the worst in all cases except for broadleaf forests in Hungary and coniferous forests in Hungary and Croatia. Interestingly, year 2004 and 2005 are indicated as worst or second worst years for forests in different countries. This could imply that the negative effects of year 2003 (which was in fact an extremely warm and dry year all over Europe [39] on forest ecosystems were not immediately visible, but the consequences of the unfavourable weather was "remembered" by the ecosystem and manifested itself in the following years. This "dampening" effect of the negative anomalies in the case of forests is very important. It shows how difficult is to make straightforward conclusions based only on the current observation without taking into account the events of the past. Table 6 shows the summary of country-specific extreme years for the different land cover types and different vegetation indices (based on the approach given in Table 1

COnClUDIng REMaRKS
In this study we do not attempt to discriminate the different RS-based indices, which means that the indices are not ranked or qualified in any way. NDVI, EVI, LAI and FAPAR are all related to plant greenness and leaf development status, which is in turn related to photosynthetic capacity and plant productivity. GPP and NPP are based on a series of model assumptions and they provide measures that are directly related to productivity. Therefore, we assume that they all serve as proxies of plant processes, and as such we can expect from them to respond to the climate fluctuations. The selection of the relatively large number of indices was done due to their different complexity in their algorithm and assumptions which were used during their calculations.
For example, NDVI and EVI are indices obtained directly from reflectances provided by the polar orbiting MODIS instrument, which can already be related to the potential plant productivity. On the other hand, FPAR and LAI (which are also obtained from remotely sensed data) are also related to plant productivity and plant development status, but unlike the former two, FPAR and LAI are calculated taking into account the biome type derived from the land cover information [18]. Thus they embedded additional information but possibly additional errors as well (either random due to e.g. misclassification of the pixel, or bias due to the possible bias in biome taBlE 4. Correlation coefficients (r) between the yearly NDVI and meteorological (temperature and precipitation) anomalies for the main land cover type in the case of the different countries (note, for Czech Republic and Slovenia the number of the selected grassland pixels were less than 5 (  0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% parameters used in FPAR and LAI calculation). Finally, MOD17 products, namely GPP and NPP, are likewise not independent from FPAR [19] but they also carry additional information (and probably additional errors) due to the meteorology information used for their calculation.
Considering the above, it is clear that all of the indices or metrics potentially have advantages and setbacks and it is difficult to select the best one. Therefore, in our work we decided to investigate all and treat them as an ensemble, similarly to the ensemble technique used in numerical weather prediction.
The selection of land cover specific anomalous years highlighted that plant response to unusual climate conditions strongly depend on the land cover type. Intraspecies differences also exist but in this study this was not addressed. The advantage of the country-mean studies is the robustness which was demonstrated earlier in other studies [15,22]. Using country means the spatial differences are most likely diminished to some extent, and in fact regions with considerably higher anomalous behaviour might exist within the countries that might require further studies.
We need to mention that in the present approach anomalies were defined for the entire growing season. It is clear that the growing season (from start of season to cessation) might be split into different time periods that needs further investigation. These anomalies associated with smaller temporal scale might provide additional information on the nature of anomalies and their meteorological driver.
The presented results might provide invaluable information for researchers associated with plant production (ecologists, agronomists, foresters, etc.). As ecosystem services are closely tied to plant productivity [8], stakeholders might also find the presented information useful. With the availability of additional years, the time series should be extended and thus the trend in the strength of the anomalies might be estimated.
The present study highlighted the difficulties related with the selection of appropriate time periods and climate indices to define extreme weather from the point of view of ecosystems. The effect-based approach relying on RS data resolves this difficulty as the observed anomalies act as integrators of the past environmental conditions thus they are good indicators of unusual conditions. Explicit selection of anomalous years indicated that the unusual plant state is not independent of the plant functional type. It is also clear that large differences exist between herbaceous and woody vegetation, where the latter is also associated with legacy effect from the previous year [40]. This legacy effect is one major issue which needs further research. In any case, legacy effect (which is well documented) alone questions the pure meteorological approach for the selection of anomalous years, as climate anomalies in the last year might exert stronger impact on the current plant growth than the present year. We propose to further refine the effect-based approach separately for the different land cover types to better understand the main drivers of plant growth in Central Europe and also worldwide.