Despite its small area, Slovenia is exceptionally diverse in terms of its physical geography.1 This is confirmed by the fact that its territory is categorized into various landscape types and regions at the European level.23.–4 The area of Slovenia is on the border of northern and southern climatic division (Fig. 1),5 what is another proof of high natural diversity. In classifying Slovenia’s small area with very diverse landscape (that is evident also in a small scale), the question arises whether it is also necessary to take into account natural elements that are commonly known as appropriate for landscape classification at continental or worldwide level (e.g., latitude, climate). One such element is the delayed appearance of phenological phases with regard to vertical and, especially, horizontal distance from the sea (e.g., onset of leafing out and onset of leaf loss). These phenomena are also known as the green wave and brown wave.6 Climate (e.g., temperature) is decisive for the development of plants and affects phenology,7,8 and climate is also an important factor in classification at the smallest scale.9 Regarding the significance of other natural elements for spatial classification, climate and the related climatic vegetation zones connected9,10 are more appropriate for analysis at the continental level, and not for small countries. Judging by climate data—specifically, the ratio between average April and October temperatures, which can show the influence of a continental location or distance from the sea—Slovenia is quite diverse despite its small area (). The ratio between average April and October temperatures in Portorož, where the cool sea prevents rapid warming of the air in early spring, is 0.81, whereas the ratio for Murska Sobota, which lies in the northeast part of the country and where the air warms up more quickly in the spring but also cools off more quickly in the fall, is 1.04.11 Similar differences can be observed between Lisbon, Portugal (0.82) and Kyiv, Ukraine (1.08), for example.12
This diversity of Slovenia raises the question whether its territory has an observable delay in phenological events as distance from the sea increases, that is, whether certain natural factors (that are usually variable at greater distance) retain high variation in a small area as well and are also observable in a small area. At the same time, phenological events can be compared with elevation, namely, temperatures change more drastically with elevation than with latitude or distance from the sea. Specifically, change in the vertical direction usually has a greater influence than in the horizontal direction.10 Because of this it can be expected that changes in phenological phenomena are more distinct in the vertical direction than in the horizontal direction (i.e., with distance from the sea). An important argument for carrying out phenological research with satellite images in Slovenia is that Slovenia is one of the most heavily forested countries in Europe.13 Along with this, some farmed areas are being abandoned and becoming forested (again).14,15
The main goal of this study is to use phenological data (especially on the onset of greenness and greenness decrease) obtained from satellite images to analyze phenological phenomena with regard to distance from the sea and elevation and to examine whether green wave is noticed in vertical and horizontal direction.
Remote Sensing and Phenology
The spectral ratio between infrared and red bands has great importance for observing the state of vegetation because vegetation reflects infrared and absorbs red part of the electromagnetic spectrum.16 The ratio enables observation of changes in structural, phenological, and biophysical parameters between seasons and years.17 It should also be added that the spectral vegetation indices are more stable across time than individual spectral bands.18
The main advantage of using remote sensing in phenological studies is the continuous capture of phenological samples in the larger areas and the ability to analyze past periods from archival images.19 The seasonal sample of variation in areas with vegetation acquired using remote sensing data is defined as landsurface phenology. This kind of phenology differs from the traditional phenology of plants, which deals with life-cycle phenomena (such as flowering, leafing out, and leaf loss) with the help of in situ observation of individual plant species.20 Specifically, remote sensing makes possible the shift from observing individual plants to broader phenological samples of the area.19 The phenology of entire ecosystems was less studied until the introduction of remote sensing.6 Satellite sensors are capable of detecting and measuring changes in the landscape at a broader scale (Fig. 2), which may not coincide with the phenological incidence of individual plants because they describe the state of the entire ecosystem.6 Due to this fact landsurface phenology is indicating a mixture of different land covers and is distinct from traditional species-centric phenology.21 Such analysis can also be carried out in areas with deficient phenological data.22 Different approaches were used to define the key dates (events) of phenological phenomena in terrain phenology.6,18,20,23
In addition to the method of defining the key dates, there are also several different phenological measures. Phenological measures are divided into three groups: temporal measures, vegetation index measures, and derived measures (Table 1). These measures usually do not directly correspond to an actual phenological phenomenon, but they are indicators of ecosystem dynamics. These phenomena can be detected with satellite observations as a sudden increase in the vegetation index and they signal a change in the ecosystem.6
Measures derived from vegetation index data (Ref. 6).
|Group of metrics||Metric||Phenological interpretation|
|Temporal vegetation index metrics||Time of onset of greenness||Beginning of measureable photosynthesis|
|Time of end of greenness||Cessation of measureable photosynthesis|
|Duration of greenness||Duration of photosynthetic activity|
|Time of maximum vegetation index||Time of maximum measureable photosynthesis|
|Vegetation index value metrics||Value of onset of greenness||Level of photosynthetic activity at the beginning of growing season|
|Value of end of greenness||Level of photosynthetic activity at the end of growing season|
|Value of maximum vegetation index||Maximum measureable level of photosynthetic activity|
|Range of vegetation index||Range of measureable photosynthetic activity|
|Derived metrics||Time-integrated vegetation index||Net primary production|
|Rate of greenup||Acceleration of photosynthesis|
|Rate of senescence||Deceleration of photosynthesis|
|Modality||Periodicity of photosynthetic activity|
This article uses measures that could be categorized into the group of temporal measures: we studied the onset of greenness increase, onset of greenness decrease, and duration of the greenness (i.e., vegetation period).
Phenology studies using satellite data go back relatively far. Rouse et al. carried out one of the first studies to determine phenological events using Landsat multispectral scanner (MSS) images.25 This was one of the first studies that used the vegetation index obtained from the red and near-infrared spectrum to define the green wave and brown wave.6,25 Later, these types of phenomena were comprehensively studied by using data from the advanced very high resolution radiometer (AVHRR)6,22,26 and the moderate resolution imaging spectroradiometer (MODIS) sensors.18,27
Data for Analysis
The information useful for regional phenology studies is captured by several satellite systems (Table 2), among which those with higher spatial resolution are especially appropriate for use.19 Some sensors have more spectral bands than others (e.g., MODIS has more spectral bands in comparison to AVHRR).28
Satellite sensors and datasets useful for studying phenological terrain features (Refs. 19 and 29).
|Landsat||MSS||1973 to 1985||79 m||18 days|
|Landsat||TM||1984 to 2012||30 m||16 days|
|Landsat||ETM+||1999 to present||30 m||16 days|
|Landsat||OLI||2013 to present||30 m||16 days|
|SPOT||Vegetation||1999 to present||1 km||1 to 2 days|
|NOAA||AVHRR||1982 to present||8 km||2×per month|
|NOAA||AVHRR||1989 to present||1 km||2×per week|
|Terra||MODIS||2000 to present||250 m, 500 m, 1 km||1 to 2 days|
|Aqua||MODIS||2002 to present||250 m, 500 m, 1 km||1 to 2 days|
|Envisat||MERIS||2002 to 2012||300 m||1 to 3 days|
|Proba-V||Vegetation||2013 to present||330 m||1 day|
For our analysis, we used MODIS satellite products on phenological phenomena (onset of greenness increase and onset of greenness decrease). We also included forest area, and digital elevation model and horizontal distance from the sea (Table 3). All of the layers were masked to exclude nonforest areas. For validation we used phenological data recorded by observers from the network of the Slovenian Environment Agency (code ARSO).30
List of data for carrying out the analysis.
|Data layer||Measurement unit||Original data layer|
|Onset of greenness increase||Sequential day of the year||MODIS data set MCD12Q2 (land cover dynamics yearly L3 global 500 m SIN grid)31|
|Onset of greenness decrease||Sequential day of the year||MODIS data set MCD12Q2 (land cover dynamics yearly L3 global 500 m SIN grid)31|
|Forest area||Forest area; nonforest area||Gozdno-vegetacijska karta Slovenije (Forest vegetation map of Slovenia)32|
|Elevation||Meters||Digital elevation model 25 m (Ref. 33)|
|Distance from sea||Meters||National general map of the Republic of Slovenia34|
|Leafing out of beeches (Fagus sylvatica)||Sequential day of the year||Meteorološki letopis Slovenije (Meteorological annals of Slovenia)30|
|General yellowing of beeches (Fagus sylvatica)||Sequential day of the year||Meteorološki letopis Slovenije (Meteorological annals of Slovenia)30|
Ground Phenological Data and Relief Data
Since 1980, there have been 61 phenological stations in Slovenia operated by the Slovenian Environment Agency, distributed according to a regional climate key. In the stations, for selected plants, the day when observed phenological phases occur (e.g., beginning of leaf folding, beginning of flowering, beginning of ripening) are recorded.35 Parallel to this monitoring, phenological observations are also carried out by the Forestry Institute of Slovenia as part of intense monitoring of the state of forest ecosystems.36
Data on elevation were obtained from a digital elevation model with a 25 m resolution,33 and data on distance from the sea were calculated based on the National general map of the Republic of Slovenia.34
MODIS Satellite Data and Products
Because of its excellent temporal and good spatial resolution, we used data and products from the MODIS system.37 MODIS is an instrument with high radiometric (12 bit) and spectral resolution: 36 spectral bands in a wavelength range from 0.4 to 14.4 μm. Its spatial resolution is somewhat smaller, at 250 to 1000 m. The MODIS instrument is onboard of two satellites: Terra, which was placed in orbit in December 1999, and Aqua, which was launched in May 2002. Both of them have near-polar orbits, are sun-synchronous, and together they image the entire Earth every one to two days.38
The applicability of products from MODIS satellite images in phenology was extensively studied for North and South America, especially in comparison with aerial photos and terrestrial measurements. It was determined that the vegetation indices calculated from the MODIS images are sensitive to seasonable vegetation variability, land cover, and biophysical parameters.17 However, it should be noted that vegetation indexes can be different, depending on sensor.39
The product range of the land cover dynamics yearly L3 global 500 m SIN grid, also known as the MODIS global vegetation phenology product, coded MCD12Q2 version 005, shows the course of phenological development of vegetation through the year at the global level. It contains data on the range and value of the enhanced vegetation index and it offers information on vegetation growth, maturity, and aging. Data on the enhanced vegetation index (EVI), calculated from the MODIS product nadir bidirectional reflectance distribution function adjusted reflectance (NBAR), were used to produce the MCD12Q2 data set. NBAR values were calculated from eight daily periods with a spatial resolution of 500 m.31 It was proved that product MCD12Q2 can be generally reliable when analyzing crop phenology.40
To define the dates of phenological transitions (or phenomena), the level of change in the curve of the model was used.41 The dates correspond to the time when the degree of change in the EVI value was lowest or greatest. Four transitional phenological events are defined: onset of greenness increase, onset of greenness maximum, onset of greenness decrease, and onset of greenness minimum (Fig. 3).41 The data layers containing the daily phenological phenomena define the approximate onset and end of the growth period, the length of the growth period, and the time of the greatest value of greenness in the year.
For all of the MCD12Q2 data sets used, we applied a uniform processing procedure. We first converted the data from hierarchical data format (HDF) format into Idrisi image data file (IMG) and transferred them to a suitable coordinate system (geographic coordinate system D48 with transverse mercator projection and Bessel spheroid 1841). Then we excluded pixels with a NoData value and masked pixels in forest areas based on the forest vegetation map. We then compared the pixel values with elevation (code DMV25) and distance from the sea (code M_DIST), and calculated the correlation coefficient.
We analyzed the dependence of phenological phenomena on elevation and distance from the sea for the following phenological data for the year 2006: onset of greenness increase [Fig. 4(a)], onset of greenness decrease [Fig. 4(b)], and length of the vegetation period [i.e., difference between onset of greenness decrease and onset of greenness increase; Fig. 4(c)].
We prepared the data and calculated the correlations using the Idrisi Taiga (version 16.05) and SPSS (version 16.0).
First, we calculated the average value for every phenological data layer and looked for the lowest and highest values. In this way, we obtained a basic overview of the phenomenon (Table 4).
Basic statistical features of phenological data from MODIS satellite images (MCD12Q2 data set). Onset of greenness decrease is more variable than onset of greenness increase.
|Day in the year|
|Average value||95% value range||Actual value range|
|Onset of greenness increase||97||77 to 143||29 to 250|
|Onset of greenness decrease||218||163 to 266||62 to 337|
|Length of growth period (onset of decrease–onset of increase)||121||60 to 179||11 to 244|
With regard to the average value of all pixels in the data layer, it can be seen that the average onset of greenness increase is the 97th day of the year, or sometime in the beginning of April. Onset of greenness decrease is defined on average on the 218th day of the year, or sometime in August, but its variability is greater than in the previous case. Length of the growing period is defined as the length of the period between leafing out and yellowing of foliage.42 In our case, we subtracted the day of the onset of greenness decrease from the day of the onset of greenness increase in order to calculate the length of the vegetative period, a consequence of which is a significantly shorter vegetative period than the actual one. The average length is, thus, 121 days, which is approximately four months.
We analyzed correlations with elevation and distance from the sea for each phenological data layer (onset of greenness increase, onset of greenness decrease, and length of the vegetation period). The correlation indicator was selected after we checked whether the variables had a normal distribution (by calculating the coefficients of skewness and kurtosis). Variables that have a normal distribution normally have coefficients of skewness and kurtosis between and 1. The basic features of certain variables in the MODIS data layers show (Table 5) that not all of them have a normal distribution, which prevents the use of parametric methods (e.g., calculating the Pearson correlation coefficient and the coefficient of partial correlation) for all variables. Therefore, we calculated Spearman’s correlation coefficient (Table 6) because this does not require a normal distribution and allows a comparison of correlation among all variables. Because the number of units (pixels) was too large, we randomly divided pixels in 10 sets. In each one of them there were of all pixels. We then calculated Spearman coefficient for each of the 10 sets; all the pixels in our study area were, thus, included in the calculation.
Basic statistical features of MODIS data (standardized values). The number of pixels is 961,570. There are no missing values. Non-normal data distribution, which is determined from the coefficients of skewness and kurtosis, prevents the use of parametric statistical methods.
|Elevation||Distance from sea||Onset of greenness increase||Onset of greenness decrease||Vegetation period|
Maximum and minimum values of Spearman’s correlation coefficients based on 10 samples with 10% of all pixels in each sample being statistically significant at p<0.01. This is a two-tailed test. Onset of greenness increase and elevation correlate most strongly (>0.7. The negative correlation coefficient between onset of greenness increase and distance from the sea is partially a consequence of the fact that in western Slovenia (which is relatively close to the sea), the land is considerably higher than in the east and northeast. Onset of greenness decrease correlates very weakly with elevation and also with distance from the sea. As a result, the correlation of the vegetation period with these data layers is also low.
|Distance from sea||Onset of greenness increase||Onset of greenness decrease||Vegetation period|
|Distance from sea||Maximum value||1.000||−0.398||0.253||0.329|
As can be seen in Table 6, there is a high positive correlation between onset of greenness increase and elevation () and, to a considerably less extent, also the correlation between elevation and vegetation period (). The latter is lower than the former primarily due to the small correlation between elevation and onset of greenness decrease, which amounts from up to and is properly negative because the vegetation period is shorter with increasing elevation. On the basis of these results, relative strong vertical green wave can be confirmed, but not also vertical brown wave.
Because elevation rises sharply across the Karst rim with distance from the sea (toward the northeast), and then gradually decreases from the Dinaric karst plateaus onward (Fig. 5), it is necessary to exclude the influence of elevation for suitably calculating the correlation of phenological phenomena with distance from the sea. Spearman’s correlation coefficient points to the decrease in elevation with distance from the sea; its values are from to , which is statistically significant. Unfortunately, the non-normal distribution of certain data makes it impossible to calculate a partial correlation, which would provide more reliable results. The fact that the onset of greenness increase occurs earlier with distance from the sea and the onset of greenness decrease occurs later is, therefore, probably a result of higher elevation in western Slovenia. Because of this, the vegetation period becomes longer with distance from the sea in our data.
In order to exclude the influence of altitude, correlation between phenological phenomena and distance from the sea only for pixels with similar elevation (Table 7) was calculated. It is evident that the correlations for most of the elevation classes are and only for three of them the correlations are just . These relatively weak connections at only certain areas do not support our hypothesis that horizontal green wave is visible in small area. In comparison to horizontal, the vertical green wave is much more evident.
Correlations between different phenological phenomena and distance from the sea for different altitude classes.
|Altitude class||Phenological phenomena|
|From (m)||To (m)||Pixels||Percentage||Onset of greenness increase||Onset of greenness decrease||Vegetation period|
Not significant at level 0.01.
For validation, it is necessary to have a good network of national and global in situ phenological and climatic observations or measures. In reality, these are usually sparse or incomplete and many studies have been carried out without validation.19,20 In Slovenia, we can use data collected by the Slovenian Environment Agency. We can also use the results of these measurements, which are available in the literature or in the form of phenological maps.7
Our results were confirmed by the facts from the literature. The onset of the growth period becomes later with increasing elevation, and the growth period is also shorter with elevation.7 Because of the low level of correlation between distance from the sea and phenological phenomena, we are unable to confirm the dependence of the onset and length of the growth period on distance from the sea.
The phenological map of the leafing out of beeches, for example, agrees with elevation, and the northeastern (continental) and southwestern (sub-Mediterranean) parts of Slovenia are classified into the same category (Fig. 6), which, thus, shows that the significance of distance from the sea is much smaller than elevation.
We used data on onset of leafing out of beeches and onset of yellowing of beeches30 to calculate the correlation in 2006 between the elevation and distance from the sea of phenological stations and the day of onset of leafing out or the onset of yellowing of leaves. Beech is one of the most widespread tree species in Slovenia, and therefore, we have data for this species for most of Slovenia. We were, thus, able to additionally assess the results obtained using data derived from MODIS satellite images.
Basic statistical features of ARSO data (standardized values). The units total 46.
|Elevation (m)||Distance (m)||Onset of leafing out (day in year)||Onset of yellowing (day in year)|
Analysis of ARSO data with Spearman’s correlation coefficients. There is a visible correlation between elevation and leafing out: the higher the elevation, the later in the year leafing out begins. The correlation between elevation and the onset of yellowing is the opposite. This is a two-tailed test. The statistical significance is given in parentheses.
|Onset of leafing out||Onset of yellowing|
|Elevation||0.66 (p<0.001)||−0.44 (p=0.002)|
|Distance from sea||Statistically insignificant||Statistically insignificant|
The calculated values show that the data from the phenological stations make it statistically possible to prove the correlation between elevation and onset of leafing out of beeches, and also between elevation and onset of yellowing of beeches. The latter correlation is relatively low in comparison to the former one. Vertical green wave proved to be more distinctive than the vertical brown wave. Similar relationships were demonstrated by the satellite image analysis. At higher elevations, leafing out occurs later in the year and yellowing begins earlier. With regard to the calculated correlation, distance from the sea does not have a statistically significant role in leafing out and yellowing. In this manner, the correlation of MODIS data on the onset of the growth period, onset of decrease in greenness, and elevation are again confirmed. We can also confirm the fact that distance from the sea does not have as much (or as strong a) significance as elevation. Additionally, we also calculated Spearman’s coefficient between distance from sea and onset of leafing out and onset of yellowing for each 100 m elevation class. None of the correlations was statistically significant at .
Direct Correlation of ARSO and MODIS Data
We also checked how the ARSO observation data and satellite image data correlate and also how individual variables differ. We compared the onset of leafing out of beeches and onset of yellowing of beeches from ARSO data and onset of greening increase and onset of greening decrease from MODIS data. In total, 44 stations were included (two stations were excluded because we did not have satellite data for them).
Because the variables did not have a normal distribution, we calculated Spearman’s correlation coefficient. This one was 0.42 for onset of leafing out of beeches (ARSO) and onset of greening increase (MODIS) and is statistically significant at the level of . The correlation between onset of yellowing of beeches and onset of greening decrease is not statistically significant, and therefore, no reliable conclusions can be drawn about it.
When comparing the features of variables from both data sources, there is a greater difference in recording the end of the vegetation period (Table 10). From this it can be concluded that the data on the onset of the vegetation period (onset of leafing out of beeches and onset of greening increase) are more comparable than are the data on the end of the vegetation period (onset of yellowing of beeches and onset of greening decrease). At the same time, this also confirms the fact that onset of greening increase is better correlated with elevation than is onset of greening decrease.
Basic statistical differences (in days) between ARSO and MODIS data for the location of phenological stations.
|Minimum value||Maximum value||Average||Standard deviation||Variance||Coefficient of variation|
|Difference between onset of leafing out of beeches (ARSO) and onset of greening increase (MODIS) (days)||5||46||24.70||10.405||108.260||42.1|
|Difference between onset of yellowing of beeches (ARSO) and onset of greening decrease (MODIS) (days)||20||161||83.11||44.743||2001.964||53.8|
The differences between ARSO and MODIS data are expected because the first involves the recording of phenological phenomena of one plant species (beeches), and the second involves the increase or decrease in greenness for the larger area, where various plant species grow—both deciduous trees as well as coniferous.
We have compared MODIS MCD12Q2 products, which contain information on various phenological phenomena, with elevation and distance from the sea for Slovenia. We confirmed the general knowledge that increasing elevation delays onset of greening and shortens the growth period. We additionally confirmed this with ARSO phenological data. The degree of correlation between elevation and onset of greening increase that was recorded with satellite images is 0.712 to 0.720, and the degree of correlation between elevation and onset of leafing out of beeches recorded by ARSO is 0.655. Onset of greening decrease recorded with satellite images and onset of yellowing recorded by ARSO are less correlated with elevation. This is primarily true for satellite image data, for which the correlation is exceptionally low (only to ), and not so much for the onset of yellowing from ARSO data, for which the correlation is . Regarding the onset and end of the vegetation period, the onset is, therefore, more comparable. We can conclude that vertical green wave is more distinct than vertical brown wave.
On the basis of the sample correlations (Table 6) between onset of greening increase and distance from the sea, it is observable that the onset of the vegetation period begins earlier with greater distance from the sea, but this is more a result of the fact that in Slovenia general elevation decreases from the southwest to the northeast. In addition, the correlations between distance from the sea and ARSO phenological phenomena are not statistically significant. Because of methodological difficulties (the non-normal distributions of variables do not permit the calculation of a partial correlation), we evaluated the correlation between MODIS phenological data and distance from the sea excluding the influence of elevation on the basis of elevation classes. After the calculation of correlations for each altitude class, it was evident that for the majority of the classes weak correlation was shown, and only for three of them correlation between distance from the sea and the onset of greenness increase was just above 0.4, which is relatively weak in comparison to the correlation between onset of greenness increase and altitude.
Despite high correlation between onset of greening increase (data recorded with satellite images) and onset of leafing out of beeches (data recorded in the field), there are differences in the days. The resolution of the satellite images used is 500 m,31 which means that the satellite images cover a greater area with various vegetation types—including many coniferous forests—and do not reflect changes for individual plants, which are recorded by ARSO observations.
This article has determined that elevation is a very important factor for the sequence of natural phenomena and the vertical green wave is clearly seen in Slovenia. We did not confirm strong brown wave in autumn when leaves fall off. We also observed that the area of Slovenia is apparently too small to observe such strong trend in the horizontal direction (horizontal green wave), particularly in the direction away from the sea toward the interior (from southwest to northeast). One of the possibilities for determining these trends would be to study the changes by days, but the problem is that monitoring of daily changes is greatly affected by weather conditions, especially cloud cover.
The research was performed in the frame of the research program Anthropological and Spatial Studies (P6-0079) and the research project Determination of natural landscape types of Slovenia using geographic information system (L6-3643) financed by Slovenian Research Agency.
Rok Ciglič received his PhD in geography from the University of Ljubljana. He works as a research fellow at Research Centre of the Slovenian Academy of Sciences and Arts. His research fields are landscape classification, natural hazards, and environment protection. He is especially interested in usage of GIS, remote sensing, machine learning, and other quantitative methods in landscape classification. He is a member of the editorial board of scientific reviews Geographical Bulletin (Geografski vestnik) and Acta Geographica Slovenica.
Krištof Oštir received his PhD in remote sensing from the University of Ljubljana. His main research fields are geometric and radiometric image correction, object-based land use classification, and change detection. He is active in the development of a small satellite for Earth observation. He is employed at the Research Centre of the Slovenian Academy of Sciences and Arts, Centre of Excellence Space-SI and as an associate professor at the University of Ljubljana.