Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI

Remote Sensing of Environment 100 (2006) 321 – 334 www.elsevier.com/locate/rse

Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI Pieter S.A. Beck a,T, Clement Atzberger b, Kjell Arild Høgda c, Bernt Johansen c, Andrew K. Skidmore b b

a Institute of Biology, University of Tromsø, N-9037 Tromsø, Norway International Institute for Geo-Information Science and Earth Observation (ITC), P.O. Box 6, 7500 AA Enschede, The Netherlands c NORUT IT AS, N-9291 Tromsø, Norway

Received 6 June 2005; received in revised form 27 September 2005; accepted 15 October 2005

Abstract Current models of vegetation dynamics using the normalized vegetation index (NDVI) time series perform poorly for high-latitude environments. This is due partly to specific attributes of these environments, such as short growing season, long periods of darkness in winter, persistence of snow cover, and dominance of evergreen species, but also to the design of the models. We present a new method for monitoring vegetation activity at high latitudes, using Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI. It estimates the NDVI of the vegetation during winter and applies a double logistic function, which is uniquely defined by six parameters that describe the yearly NDVI time series. Using NDVI data from 2000 to 2004, we illustrate the performance of this method for an area in northern Scandinavia (35
162 km2, 68- N 23- E) and compare it to existing methods based on Fourier series and asymmetric Gaussian functions. The double logistic functions describe the NDVI data better than both the Fourier series and the asymmetric Gaussian functions, as quantified by the root mean square errors. Compared with the method based on Fourier series, the new method does not overestimate the duration of the growing season. In addition, it handles outliers effectively and estimates parameters that are related to phenological events, such as the timing of spring and autumn. This makes the method most suitable for both estimating biophysical parameters and monitoring vegetation phenology. D 2005 Elsevier Inc. All rights reserved. Keywords: Boreal forests; Tundra; Climate change; Green up; Fennoscandia

1. Introduction The monitoring of global vegetation activity using time series of satellite data indicates that vegetation at high latitudes is displaying considerable variation both annually and interannually (Goward & Dye, 1987; Lucht et al., 2002; Myneni et al., 1997; Shabanov et al., 2002; Tucker et al., 2001). Zhou et al. (2001), for example, observed an overall greening and longer growing seasons in the period 1982 – 1999. Such trends may be caused by global warming (Braswell et al., 1997; Zhou et al., 2003). For very high latitudes, potential scenarios for the future include a northward shift of the area suitable for boreal forest by 150 –550

* Corresponding author. E-mail address: [email protected] (P.S.A. Beck). 0034-4257/$ - see front matter D 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2005.10.021

km over the next century and an upward shift of the tree line (Kirschbaum et al., 1996; Skre, 2001). 1.1. Characterising vegetation using the NDVI Most monitoring of large-scale vegetation activity is based on normalized difference vegetation index (NDVI) time-series datasets acquired by sensors aboard satellites with short revisit periods, in particular the National Oceanic and Atmospheric Administration (NOAA) Advanced Very High Resolution Radiometer (AVHRR) (Goward & Prince, 1995). This sensor has a spatial resolution of roughly 1 km
1 km and with its large swath width it can sense almost the entire Earth on a daily basis. The NDVI relies on the absorption of red radiation by chlorophyll and other leaf pigments, and the strong scattering of near-infrared radiation by foliage. As a consequence, red reflectance (q Red) tends to decrease as the amount of green

322

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

vegetation in a pixel increases. At the same time the structural properties of a denser canopy (e.g. increasing leaf area index and crown coverage) cause an increase in near-infrared reflectance (q NIR). The NDVI quantifies the q Red – q NIR contrast by calculating the ratio of the red and near-infrared surface reflectance: NDVI ¼ ðqNIR  qRed Þ=ðqNIR þ qRed Þ:

ð1Þ

The interpretation of NDVI in terms of structural properties of vegetation (e.g. leaf area index—LAI) and its photosynthetic activity (i.e. the fraction of Absorbed Photosynthetically Active Radiation—fAPAR) is supported by empirical studies, as well as simulations with physically based radiative transfer models (e.g. Baret & Guyot, 1991; Guyot et al., 1989). However, it has to be pointed out that the relation between optical properties of forests and variables such as LAI and fAPAR depend on a variety of other canopy characteristics as well (Nagler et al., 2004; Running et al., 1989; Schlerf & Atzberger, accepted for publication); e.g. leaf optical properties, foliage clumping, and leaf angle distribution. Unfortunately, the NDVI is also sensitive to external factors not related to the forest canopy itself (Gamon et al., 2004; Peterson, 1991). Depending on the illumination and observation geometry, the spectral bi-directional reflectance of forests may decrease or increase, even without any change in the forest characteristics due to BRDF effects (Deering et al., 1999). Recent work has shown, however, that the BRDF effects of many vegetation types are quite similar (Bacour & Bre´on, 2005). Due to additive path radiance, the atmosphere leads to an increase in red reflectance. At the same time, the near-infrared reflectance decreases due to a reduced atmospheric transmission. The effect is that the top-of-atmosphere NDVI is generally lower than the top-of-canopy NDVI (Guyot et al., 1989). Consequently, atmospheric ‘‘noise’’ in the NDVI is often considered negatively biased. Another fundamental challenge to the retrieval of useful land cover information from remote sensing is the complex stand structure, particularly the two layered canopy (e.g., conifer of birch forests on top of a surface moss or lichen layer and/or snow layer), and the varying degree of canopy closure. For example, Chen and Cihlar (1996) and Hu et al. (2000) showed that the variability in the understory largely determines the variation in the optical signatures of open boreal forests. When combined with changing sun and view angles, these structural complexities yield spatially and temporally variable reflectance patterns that often confound the estimation of land cover information from remote sensing (Gamon et al., 2004). 1.2. Modelling vegetation phenology using NDVI time series Seasonal variations of NDVI are closely related to vegetation phenology, such as green-up, peak and offset of development (McCloy & Lucht, 2004). Specific algorithms use NDVI time series to derive parameters related to vegetation phenology and production that are essential for modelling vegetation distribution and dynamics (reviewed in Pettorelli et al., 2005). These algorithms either model the

entire time series using a mathematical function (Jo¨nsson & Eklundh, 2002; Los et al., 1994; Running & Nemani, 1988; Sto¨ckli & Vidale, 2004) or estimate the timing of single phenological dates and characteristics, such as the start of the growing season (Badeck et al., 2004; Chen et al., 2001; Duchemin et al., 1999; Kaduk & Heimann, 1996; Karlsen et al., accepted for publication; Lloyd, 1990; Reed et al., 1994; Xin et al., 2002). Methods for analysing entire time series include principle component analysis (Hirosawa et al., 1996), Fourier analysis (Azzali & Menenti, 2000; Sto¨ckli & Vidale, 2004), harmonic analysis (Jakubauskas et al., 2001), wavelet decomposition (Anyambe & Eastman, 1996; Sakamoto et al., 2005), and curve fitting (Jo¨nsson & Eklundh, 2004; Zhang et al., 2003). Modelling NDVI time series as such has the advantage of conserving a maximum amount of information in the NDVI data, while reducing the dimensionality of the data (Jo¨nsson & Eklundh, 2004). Therefore, in addition to the phenological dates, other parameters can be estimated from the models output. Due to the close link between NDVI and fAPAR, for a variety of cover types, net primary production and CO2 fluxes can be estimated from time integrated NDVI (Running & Nemani, 1988; Tucker et al., 1983; Wylie et al., 2003). NDVI values are also associated with actual and potential evapotransporation (Choudhury et al., 1994; Cihlar et al., 1991). These parameters can in turn serve as biosphere-related input for biochemical and climate models (Cramer, 1997). 1.3. Compositing and describing NDVI time series As mentioned, NDVI images are often degraded by atmospheric conditions such as cloud cover, dust and aerosols (Forster, 1984), as well as by sensor-related artefacts such as line dropout, band-to-band calibration, and so on (Richards, 1986). Several methods exist to minimize these perturbations. For example, Cihlar (1996) developed a method for cloud estimation from composite imagery based on albedo and NDVI trends. Roerink et al. (2000) focussed on harmonic analysis of time series, whereas Viovy et al. (1992) developed the best index slope extraction (BISE). One often applied technique is maximum value compositing (MVC, Holben, 1986). It selects for each pixel the highest observation from a predefined compositing period to represent this period. This is based on the logic that low-value observations either are erroneous or have less vegetation vigour for the period under consideration. Applying techniques such as maximum value compositing has resulted in a number of corrected NDVI datasets (e.g. Global Inventory Monitoring and Modelling Studies NDVI (GIMMS, 2005) and Pathfinder NDVI (James & Kalluri, 1994)). For these in turn, algorithms have been created that increase the spatial and temporal consistency in the data and reduce data dimensionality (e.g. Chen et al., 2000; Jo¨nsson & Eklundh, 2002). The Fourier Adjustment, Solar zenith angle, correction, Interpolation and Reconstruction NDVI (FASIR-NDVI, Sellers et al., 1996) and the European Fourier-Adjusted and Interpolated NDVI (EFAINDVI, Sto¨ckli & Vidale, 2004) datasets were produced by such algorithms and originate from GIMMS NDVI data for the

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

323

Fig. 1. Map of Northern Europe showing the location of the 35
162 km2 study area, as well as a digital elevation model and a simplified vegetation map of the study area.

entire globe and Pathfinder data for Europe, respectively. They have been used for vegetation monitoring (Sto¨ckli & Vidale, 2004) and climate modelling (Sellers et al., 1997). The algorithms used to create these datasets typically consist of operations to eliminate or replace spurious values, and the description of the NDVI time series using mathematical functions. These functions need to result in smooth curves consistent with the seasonal behaviour of vegetation. For example, in humid temperate areas, the NDVI of natural vegetation increases in spring, stays relatively constant during summer, and decreases again in autumn. The FASIR and the related EFAI algorithms both rely on second-order Fourier series, whereas Jo¨ nsson and Eklundh (2002) applied a combination of Gaussian-type functions. 1.4. NDVI time series analysis in high-latitude environments High-latitude environments pose particular challenges to modelling the annual behaviour of vegetation using vegetation indices. They are characterized by a short growing season, often shorter than four months, due to long cold winters and the persistence of snow (Tuhkanen, 1980; Wielgolaski & Inouye, 2003). This limits the opportunities to acquire reliable NDVI measurements over active vegetation (Justice 1986). Evergreen coniferous forests, which dominate boreal areas, maintain their foliage during winter but show higher needle-to-shoot ratios in summer, which is reflected in vegetation indices (Chen & Cihlar, 1996; Huemmrich & Goward, 1997). Hence, one needs to estimate the NDVI outside the growing season (later referred to as winter NDVI). The estimated winter NDVI should represent the NDVI of coniferous forests when the needle-toshoot ratio is at its lowest and be close to zero for deciduous forests. During winter, the polar night at very high latitudes makes it impossible to measure the NDVI. Snow cover causes a

negative bias on NDVI values, further limiting the possibilities of observing the winter NDVI. This bias causes NDVI values to rise during snowmelt, even if this is not related to increased vegetation activity. Estimating winter NDVI is therefore crucial to distinguishing the effective greening up of the vegetation during spring from snowmelt. Many of the algorithms used for vegetation monitoring are ill-suited for boreal regions (e.g. Andres et al., 1994; Menenti et al., 1993; Moody & Johnson, 2001; Olsson & Eklundh, 1994; Sellers et al., 1996; Sto¨ckli & Vidale, 2004). These algorithms rely on second-order Fourier series which cannot represent short growing seasons well, especially when the increase and decrease in NDVI in spring and autumn occur abruptly. In this case, the fitted Fourier series tends to overestimate the duration of the growing season (Sto¨ckli & Vidale, 2004). Other algorithms have not yet been applied in high latitude environments. This is striking, given the susceptibility of these areas to environmental change (Cramer, 1997; Foley et al., 1994). For more southern biomes some algorithms have been developed that can model short growing seasons as well as abrupt changes in NDVI during spring and autumn. Jo¨nsson and Eklundh (2002), for example, used asymmetric Gaussian functions with Pathfinder NDVI data of Africa and Zhang et al. (2003) used piecewise logistic functions to model phenological transition dates in the northeastern United States. Agricultural vegetation can display a short growing season and rapid changes in vigour, related to the emergence, greening up, wilting, or harvesting of the crop. To describe seasonal variation in the NDVI of agricultural crops and of shrubland and grassland, double logistic functions of time have been used (Fischer, 1994a, 1994b; Paruelo & Lauenroth, 1998). The potential suitability of double logistic functions to describe vegetation dynamics of subarctic

324

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

vegetation was originally suggested by Shibayama et al. (1999). 1.5. Objectives Here, we test whether a method using a double logistic function is more suitable for modelling the yearly NDVI time series of boreal and tundra vegetation than approaches based on Fourier series or asymmetric Gaussian functions, and compare the method against current state-of-the-art alternatives. We present a method for modelling growing season dynamics of high-latitude environments using yearly NDVI time series acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS). Compared with the AVHRR sensors, the MODIS sensors aboard the TERRA and AQUA satellites launched in 1999 and 2002, respectively, have improved calibration and atmospheric correction. Furthermore, MODIS NDVI MVC products exist for the entire globe at spatial resolutions as high as 250 and 500 m (compared with 8 km for Pathfinder and GIMMS). This makes studies of regional vegetation dynamics possible, but increases the need for algorithms that reduce the size of the datasets. 2. Developing the method 2.1. Test site The test site is a 35
162 km2 area located on the borders of Finland, Norway and Sweden (geographical corners of the area: 68-57V30WN, 22-34V17WE; 67-30V1WN, 21-10V39WE; 68-57V30WN, 23-26V28WE; 67-30V1WN, 21-59V36.22WE; Fig. 1). The area is covered by diverse boreal vegetation, including deciduous and evergreen forest dominated by Betula pubescens and Pinus sylvestris respectively, bogs, minerotrophic mires, and alpine tundra vegetation dominated by shrub heaths. 2.2. Data preprocessing We used the 16-day composite MODIS NDVI dataset at a 250-m resolution, spanning the period from 2000 to 2004 (the MOD13Q1 product, Huete et al., 2002). It is produced from surface reflectance data corrected for molecular scattering, ozone absorption, and aerosols (Vermote et al., 2002). Prior to compositing the NDVI data, a filter selects cloud-free nearnadir observations. Vegetation indices (VI) often increase as they are measured more obliquely. To buffer against this effect, the compositing algorithm used for the MODIS VI constrains the variation in viewing angle in the data. It compares the two highest NDVI values after filtering and selects the observation closest to nadir view to represent the 16-day composite cycle. The MVC dataset contains 23 NDVI observations per pixel per year. These observations are accompanied by auxiliary information on different parameters, such as aerosol quantity and likelihood of snow cover, as well as an overall indication of the usefulness of each data point (scaled from 15 to 0, with 0 indicating Fperfect quality_).

Using the Fquality information_ provided by the MODIS NDVI MVC product, we first eliminated lower-quality data (usefulness index 7 Faverage quality_ or worse). Next, we removed outliers (NDVI > 0.98). We then estimated the winter NDVI, which was assumed to be pixel-specific but constant across years. At the end of winter in northern boreal areas, the disappearance of snow from the landscape and bud break in deciduous trees can be less than a few days apart (Wielgolaski & Inouye, 2003). Hence, the only time window available for observing the winter NDVI is between the end of the growing season and the arrival of snow or the polar night (e.g. at 68- N the latter occurs around 10 November). The first snow fall mostly occurs several weeks after deciduous trees have shed their leaves. Ground observations indicate that leaf fall in northern Norway starts no later than mid-October (NML, 2005). Therefore, we considered a window of one month, between 16 October and 17 November, as suitable for observing the winter NDVI, although some snow fall might occur during this time. The MODIS quality data indicate that in the second half of this period pixels are likely to contain snow. This is problematic for the estimation of the winter NDVI as it could introduce a negative bias. As the MODIS time series used in this study consists of only 5 years, we decided to nevertheless use all data points within the time window in order to base the winter NDVI on enough data points. As the MODIS time series grows, this element of the algorithm could be adjusted to resolve this issue. Each year, two MVC NDVI values represent this time window, namely the images 19 and 20 in the time series. Per pixel, the 5-year dataset provides two times five NDVI values from which to estimate the winter NDVI. We selected the highest values observed in each of these two sets, giving the highest value observed in image 19 and in image 20 across years. If any of these maxima were negative, they were assigned a value of 0, because negative NDVI values indicate the absence of green vegetation. We then calculated the geometric mean of the two values and used it as a pixelspecific estimate of winter NDVI (Fig. 2): pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi winter NDVI ¼ max ðimage19Þ
max ðimage20Þ: ð2Þ This method of calculation assigns a value of 0 to the winter NDVI if either of the values is 0. Applying the geometric mean rather than the arithmetic mean thus ensures that only pixels displaying a prolonged period of greenness after the growing season are assigned positive winter NDVI values. Next, we replaced all values lower than the winter NDVI with the winter NDVI (Fig. 2b). No useful NDVI recordings were ever obtained in December because of the polar night. As a consequence, very few NDVI recordings were available after the growing season—which complicates the curve fitting for this period. Hence, we set the missing recordings for December to the winter NDVI (Fig. 2b). For each year, pixels where less than six NDVI observations exceeded the winter NDVI were excluded from further analysis. After pre-processing the data as outlined, a double logistic curve was fitted to the data (Fig. 2c). This procedure is described in the following section.

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

325

a)

NDVI

1

0

2000

2001

2002

b)

NDVI

0

0 j

f m a m j

j

a s o n d

2004

c)

1

NDVI

1

2003

j

time (months)

f m a m j

j

a s o n d

time (months)

Fig. 2. Filtering and modelling of yearly NDVI series. (a) Five-year NDVI series for a single pixel. In each year, the period between 16 October and 17 November is represented by two, or fewer, reliable NDVI observations (O). The winter NDVI (dotted line) is calculated as the geometric mean of the highest NDVI values within these two sets of images (bold circles), after setting negative values to zero. The probable presence of snow is indicated by (
). Note that while one of the maxima is likely to be affected by snow in the pixel, it is nevertheless included for the calculation of the winter NDVI (see text for details). (b) For each year (e.g. 2004), values lower than the winter NDVI were adjusted, as were the missing data points for December. (c) A double logistic curve was fitted to the adjusted data (dashed curve), negative outliers during the growing season (o) were eliminated, and the curve was refitted (solid line).

We propose a double logistic function to describe NDVI time series spanning one year. The double logistic function models the NDVI as a function of time (t) using six parameters (Fig. 3): the winter NDVI (wNDVI); the maximum NDVI (mNDVI) during the growing year; two inflection points, one as the curve rises (S) and one as it drops (A); and the rate of increase or decrease (mS and mA) of the curve at the inflection points. NDVIðt Þ ¼ wNDVI þ ðmNDVI  wNDVIÞ  1
1 þ expð  mS
ðt  S ÞÞ  1 1 : ð3Þ þ 1 þ expðmA
ðt  AÞÞ For each pixel and year, we used the earlier derived winter NDVI and estimated the remaining five parameters using an iterative non-linear least squares procedure. This resulted in a double logistic curve fitted to the NDVI time series. Based on the residuals produced by this fit, we set up a weighting scheme for the observations and in a second iteration reestimated all parameters except the winter NDVI. The weighting scheme considers observed NDVI values that are initially overestimated by the fitted curve and assigns lower weights to them for the second iteration (Sellers et al., 1994). This type of weighting is commonly applied in NDVI modelling algorithms and rests on the assumption that anomalies in the NDVI dataset are negatively biased (Sellers et al., 1994). We only applied weighting to the observations during the estimated growing season (delimited by the initial estimates of the inflection points S and A in Fig. 3). Visual

inspection of the data revealed that exceptional variation in the NDVI occurred early during the growing season in some years (e.g. Fig. 4b). These variations are likely to be erroneous, as the NDVI is expected to evolve in a relatively smooth manner. As the quality information did not allow us to determine which of the values were most likely to be wrong, we limited the weighting early in the growing season. The weighting scheme assigned a weight W to each observation, based on the initial residuals r and on tr, the relative timing of the observation during the growing season (Fig. 5). The calculated timing (tr) varied between 0 (at the first inflection point S) and 100 (at the second inflection point A).

mNDVI

~mA

NDVI

2.3. Curve fitting

~mS wNDVI

S

A

time Fig. 3. Example of the double logistic function (Eq. 1) used to model the yearly NDVI time series. It is defined by six parameters: the winter and maximum NDVI (wNDVI and mNDVI), the spring and autumn inflection points (S and A), and the rate of increase or decrease in NDVI at S and A, respectively (mS and mA).

326

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

time (months)

a)

j

m

a

m

j

j

a

s

o

n

d

time (months)

b) 1.0

0.8

0.8

0.6

0.6

NDVI

NDVI

1.0

f

0.4

0.4

0.2

0.2

0.0

j

f

m

a

m

j

j

a

s

o

n

d

0.0 1 2 3 4 5 6 7 8 9 10 111213141516171819 202122 23

1 2 3 4 5 6 7 8 9 1011 121314151617 1819 20 2122 23

image number

image number

Fig. 4. One-year NDVI time series for two pixels in the study area after data filtering and correction (dots), with fitted curves using a second-order Fourier function (thin line), a double logistic function (thick line), and asymmetric Gaussian functions (dashed line) all with weighting schemes applied.

1

W ¼

2

ðtr
r þ 1Þ W ¼ 1;

if r > 0:01

;

if 0V r V0:01

W ¼ 4;

ð4Þ

if r < 0

where r = predicted NDVI  measured NDVI, tr = relative time during the growing season. Using this scheme, NDVI values that were lower than the initially fitted curve (r > 0) were ignored, particularly if they occurred later in the growing season. Values that were slightly overestimated (0  r  0.01) or underestimated (r < 0) were given extra weight. The result is a smooth curve that describes the yearly NDVI data, allowing for large increases and decreases in spring and autumn, respectively, yet ignoring abrupt drops in the NDVI during the growing season.

2 1

0 ≤ r < 0.01

0

weight

3

4

r<0

r = 0.01 r = 0.1 0

20

40

60

80

100

timing during the growing season Fig. 5. Weighting applied to the NDVI values during the growing season for the double logistic curve fitting, whereby the growing season is delimited by the first inflection point (timing = 0) and the second inflection point (timing = 100) in the unweighted double logistic fit. To minimize the influence of large NDVI variations early in the growing season, the weights are adjusted according to when the NDVI values were recorded. Further, the weights are a function of r, the residuals produced by the unweighted fit, calculated as measured NDVI  predicted NDVI.

Fourier series consist of a summation of sine and cosine functions. Second-order Fourier series require five parameter estimates (a 0, a 1, a 2, b 1, b 2) and are suitable for describing seasonal NDVI trajectories of many different vegetation types (Los, 1998). They model the NDVI at time t as: 2 X NDVIðt Þ ¼ aj
cosð j
Ut Þ þ bj
sinð j
Ut Þ ð5Þ j¼0

where U ¼ 2p
ðt  1Þ=n. For comparison purposes, we fitted second-order Fourier series to the NDVI time series after data filtering and replacement with winter NDVI (n = 23, t ranges from 1 to 23). We then applied the weighting scheme of the EFAI algorithm (Sto¨ckli & Vidale, 2004), with slight modifications to optimize it for our data. Here too, weighting was only applied to the growing season. Estimates of the start and end of the growing season cannot be derived from the parameters of the fitted curves. They were therefore based on the days when the first derivative of the curve first exceeded 0.003 or  0.003/ day, as these are the days when the NDVI, and thus vegetation activity, starts to increase/decrease. These dates were then shifted 16 days toward the beginning/end of the year instead of 10 days as in the original EFAI algorithm, according to the compositing period in the dataset. We also compared our method with a new algorithm based on fitting an asymmetric Gaussian function to filtered NDVI data, as implemented in the TIMESAT software (Jo¨nsson & Eklundh, 2004): NDVIðt Þ ¼ wNDVI þ ðmNDVI  wNDVIÞ
g ðt Þ where     t  a1 a3 gðt Þ ¼ exp  ; if t > a1 a2     a1  t a5 gðt Þ ¼ exp  ; if t < a1 a4

ð6Þ

a 1 determines the position of the maximum NDVI (mNDVI), with respect to t. The other five parameters are similar to the parameters in the double logistic function: a 2 and a 4 are namely related to the

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

A estimated in the unweighted double logistic function were used to delimit the growing season. The root mean square error (RMSE) calculated from observed and estimated NDVI values was then used as a measure of model performance. A Wilcoxon rank sum test for paired observations was used to test whether different methods generated different RMSEs, corresponding to a better or worse fit. To further investigate the capacity of the methods to describe the NDVI trajectories, we used the residual NDVI values. The residuals produced by a model were checked for the presence of autocorrelation in time, as this can indicate that the model is misspecified (Gujarati, 2003). This was done using the autocorrelation function (ACF). The ACF quantifies the similarity between residuals that are a certain lag apart in the time series. An ACF close to zero indicates nearly no

position of the spring and autumn inflection points relative to the position of the maximum, while a 3 and a 5 are related to the slope of the curve at these points. We fitted unweighted asymmetric Gaussian functions to the corrected NDVI time series to test whether Gaussian functions describe the yearly time series better than double logistic functions. Afterwards, we applied the default weighting scheme implemented in the TIMESAT software and refit the asymmetric Gaussian curves. 2.4. Evaluation of the models Outside the growing season, the double logistic function should fit the observed NDVI values well, as it includes the winter NDVI as a parameter. Hence, we only quantified the fit of the curves during the growing season. The parameters S and

0

0.5

Fourier Fit

0.5

ACF of residual NDVI

327

1

2

3

4

5

6

7

8

9

10

11

12

13

10

11

12

13

11

12

13

0

0.5

Double Logistic Fit

0.5

ACF of residual NDVI

lag (images)

1

2

3

4

5

6

7

8

9

0

0.5

Asymmetric Gaussian Fit

0.5

ACF of residual NDVI

lag (images)

1

2

3

4

5

6

7

8

9

10

lag (images) Fig. 6. Autocorrelation plots for the unweighted Fourier series, double logistic function and asymmetric Gaussian function based models. Bars and whiskers indicate the mean and standard deviation, respectively, of the estimates of autocorrelation produced by 500 pixels uniformly distributed over the entire study area and uniformly drawn from the years 2000 to 2004.

328

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

time (months) j

f

m

a

m

j

j

a

s

o

n

d

time (months)

b)

j

f

m

a

m

j

j

a

s

o

n

d

time (months)

c)

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4 0.2

0.2

1 2 3 4 5 6 7 8 9 10 11 12 1314 1516 171819 20 2122 23

f

m

a

m

j

j

a

s

o

n

d

0.4

0.0 1 2 3 4 5 6 7 8 9 10 11 12 1314 15 16 1718 19 20 21 22 23

image number

j

0.2

0.0

0.0

1.0

NDVI

1.0

1.0

NDVI

NDVI

a)

1 2 3 4 5 6 7 8 9 1011 12 13 1415 161718 19 20 21 22 23

image number

image number

Fig. 7. Fit of unweighted and weighted (a) Fourier series, (b) double logistic functions and (c) asymmetric Gaussian functions to NDVI time series after data filtering and correction. Unweighted functions are shown as thin lines and weighted functions as thick lines. The weighting algorithm increases the positive bias in NDVI estimates outside the growing season, particularly for the Fourier series.

correlation, and as ACF values approach T 1 they indicate increasing correlation/anticorrelation. To reveal bias and error in the estimates of the NDVI values throughout the year, we used the mean and variation in the residuals for each of the 23 NDVI observations per year. For the evaluation based on the RMSE and the residuals, we used a subset of the data, consisting of 500 one-year time series. These were extracted from 100 pixels uniformly distributed over the entire study area and drawn from the years 2000 to 2004. Note, however, that some of these time series consisted of too few observations to be successfully modelled. They were therefore excluded from the evaluation. 3. Results and discussion The double logistic function describes the NDVI data better than the Fourier series and slightly better than the asymmetric

0.5 0.55

0

residual NDVI

0

residual NDVI

0.55

0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

image number

image number

Fourier Fit

Double Logistic Fit

0.55

0

residual NDVI

0.5

Asymmetric Gaussian Fit

0.5

residual NDVI

0.55

0

0

residual NDVI

0.55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

image number

0.5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

0.55

residual NDVI

TIMESAT weighted Gaussian Fit

0.5

Weighted Double Logistic Fit

0.5

Weighted Fourier Fit

Gaussian functions, both outside and during the growing season (RMSE (double logistic) < RMSE (Gaussian) < RMSE (Fourier), p < 0.001 in all cases). The Fourier series cannot mimic the shape of the annual NDVI curves as well as the double logistic and asymmetric Gaussian functions can, which is indicated by the temporal autocorrelation in the NDVI residuals (Fig. 6). In particular, the Fourier series is unable to describe the abrupt start and end of the growing season (Fig. 4). This results in estimated values that are too high before the spring increase in NDVI (images 7 and 8 in Fig. 7). On the other hand, the estimated values at the end of spring, as the NDVI ceases to increase, are too low (images 9 and 11). Accordingly, the NDVI is underestimated as autumn starts (images 17 and 18) and overestimated as it ends (images 19 and 20). The double logistic function produces the same trend in the residuals around the spring periods. However, the errors are smaller than those with the Fourier-based method and extend over two compositing

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

image number

image number

image number

Fig. 8. Residuals for Fourier series, double logistic function and asymmetric Gaussian functions based models of yearly NDVI time series. Positive residuals represent overestimated NDVI values; negative residuals represent underestimated values. Whiskers span 4
interquartile range and outliers are not drawn. Data are presented for 500 pixels uniformly distributed over the entire study area and uniformly drawn from the years 2000 to 2004.

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

periods only (images 7 and 8) rather than three or more (images 7 –10). In the case of the double logistic function, the errors are caused by the fact that this function, by definition, displays an equal but opposite curvature at the start of spring and at the end of spring. The observed rise in NDVI at the start of spring is generally more abrupt, however, than the levelling out at the end of spring. As a result, the fitted function slightly overestimates values right before the start of spring, and underestimates the values at the start of spring (Fig. 7). The double logistic and asymmetric Gaussian functions produced similar results (Fig 8), as expected based on their definition. However, the logistic functions mimicked the NDVI trajectories after and at the very end of the growing season slightly better than the asymmetric Gaussian function did. This is indicated by the larger residuals the asymmetric Gaussian functions produced after the growing season (images 18, 19 and 20). The weighted methods are designed to reduce the effect of relatively low NDVI values on the fitted values (Fig. 7). The latter are identified as those observations that have positive residuals after the non-weighted fit. A general effect of the weighting procedure is an overall increase in residuals (Fig. 8). For the Fourier-based method this leads to better estimates of NDVI at the end of spring (images 9 –11) and the start of autumn (images 17 and 18). These values were initially underestimated. The better estimates for spring and autumn NDVI values are at the cost of the accuracy of the NDVI estimates right before and after the growing season (images 7 – 8 and 19– 21, respectively). Although the weighting scheme results in values outside the growing season remaining unchanged, these values are down-weighted relative to the succeeding underestimated values in the case of spring, and the preceding underestimated values in the case of autumn. This causes the weighting algorithm to increase the positive bias in their estimates (Fig. 8). In the case of the weighted double logistic function and asymmetric Gaussian function, similar trends are observed. The added bias at the beginning and end of the growing season is much smaller, however, because the initial fit of the function is better. It indicates that the weighting schemes should be used with care, especially if descriptions of vegetation behaviour during spring and autumn are the main aim of the modelling effort. All methods produce smaller residuals during winter because of the constancy of the winter NDVI. The estimate of the winter NDVI was based on the NDVI time series only. We expected areas with more coniferous forest to display higher winter NDVI, as these forests retain their leaves in winter. A comparison with a vegetation map based on a Landsat 7 Enhanced Thematic Mapper image taken on 27 July 2000 confirms this, as it shows similar patterns in the winter NDVI and the presence of coniferous forest (Fig. 9). The use of a baseline value for the NDVI during winter relies on the knowledge that vegetation in boreal areas goes through a period of stable dormancy during winter. Using data for all years to calculate the winter dormancy ensures that the limited opportunity to measure the winter NDVI is maximally exploited. As it produces a common winter NDVI for all

329

Fig. 9. Agreement between the presence of evergreen coniferous vegetation and winter NDVI. The cover of coniferous forest was mapped using a classified Landsat Enhanced Thematic Mapper image, taken 27 June 2000, and the winter NDVI was derived from the MODIS NDVI time series.

years per pixel, curves fitted to yearly NDVI time series can be joined to describe the seasonal behaviour of vegetation over several years. However, it does assume that the greenness of the vegetation during winter stays constant across years. This might not be the case for growing coniferous forests. At high latitudes, solar zenith angles are typically lower during spring than during autumn. In our test area the maximum solar zenith angle is 82- on the 25th of October and 57- on the 25th of April. This influences reflectance values because of different atmospheric path lengths and nonLambertian properties of vegetation. While the ratioing

330

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

properties of the NDVI reduce these influences (Kimes et al., 1984), the winter NDVI values could be biased, as they can only be based on NDVI observations during autumn. Applying a bidirectional reflectance distribution function (Rahman et al.,

1993) with parameters for high latitude pine forests (Deering et al., 1995) indicated that different sun angles result in 10% higher NDVI at the end of the dormant period (April) compared to the start of the dormant period (October). Using

2003

2004

2003

2004

2003

2004

2003

2004

Fig. 10. Estimates of the timing of spring and autumn in the study area for 2003 and 2004. Values in the histograms and maps represent the parameters S (middle of spring) and A (middle of autumn). These were estimated when fitting the weighted double logistic function to the NDVI time series. To substitute missing values, a 5
5 pixel median filter was applied to the maps.

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

10% lower winter NDVIs, the estimates of the inflection points in spring were on average 1.5 days earlier (standard error = 0.002) and the ones in autumn 0.6 days later (S.E. = 0.01). The effect a 10% bias in the winter NDVI has on the predicted curve parameters depended on the winter NDVI itself. For pixels with a 0.1 higher winter NDVI the inflection point in spring was predicted 10 (S.E. = 0.1) days earlier and the one in autumn 3.7 (S.E = 0.1) days later when using the lower winter NDVI values. Altogether, the deviations are significantly smaller than the 16-day uncertainty inherent to the MODIS 16-day MVC data and are not likely to bias the modelled growing season characteristics. One of the advantages of using the double logistic function for describing a unimodal growing season is that the estimated parameters can be directly interpreted in terms of vegetation phenology (Fischer, 1994a). The parameters S and A, for example, provide the two inflection points of the fitted curve and can serve as estimates of the timing of spring and autumn, respectively. Mapping these parameters for an area or different years enables the monitoring of large-scale vegetation dynamics. A map of the midpoints of spring and autumn for 2003 and 2004 (Fig. 10) reveals considerable differences across the study area and between years. The timing of autumn appears to vary less in space than the timing of spring. This might be because the environmental and physiological cues for leaf emergence in spring and leaf senescence in autumn show different spatial patterns. Spring phenophases differ between species and correlate with winter and spring temperatures, which vary locally following altitude and geographical gradients (Taulavuori et al., 2004). In addition, spring phenophases within single tree species vary up to several weeks owing to genetic variation (Kramer, 1995). The controls over plant phenology in autumn are little understood, as correlations between meteorological factors and leaf senescence are weak (Menzel, 2002). Nevertheless, the lack of geographical variation in the timing of autumn indicates that this timing is driven by large-scale processes such as longer-term temperature regimes and day length, as suggested by Smit-Spinks et al. (1985). In our study, the modelled timing of autumn correlates with altitude, particularly in 2003 (compare DEM in Fig. 1 and autumn 2003 in Fig. 10). This could indeed be the result of mean temperatures, which decrease with altitude, driving autumn phenology. Vegetation communities of high latitude environments typically show single modes of growth and senescence within the short yearly growing season (Wielgolaski & Inouye, 2003). Seasonal patterns become more complex however when viewed at coarser resolutions, due to the presence of different communities and local climatic gradients. Earlier, this posed problems when using NDVI datasets of relatively low spatial resolution such as GIMMS NDVI with 8 km spatial resolution (Karlsen et al., accepted for publication). The high spatial resolution of the MODIS dataset used here, namely 250 m reduces such effects. Hence, it increases the likelihood of observing unimodal NDVI curves, representing the growthsenescence cycle of the canopy, which we found to be well described by our algorithm. Therefore, when used with the

331

MODIS 16 day MVC product at its highest resolution, the presented method should be suitable for modelling and mapping the growing season for larger high latitude areas, also outside Scandinavia. The method is designed to describe the general characteristics of the growing season, such as its start, end, and peak. It does not allow capturing variations within the growing season that do not correspond with the unimodal shape of the double logistic curve. Such variations in NDVI could occur due to infestation of insects, and exceptional drought or temperatures which can cause diminished photosynthetic activity, leaf discolouration, and leaf fall (Tømmervik et al., 2001). As the method assumes the NDVI time series to follow a predefined pattern, deviating courses will not be modelled well. On the other hand, the overall fit of the modelled functions will decrease in such cases. This, in turn, allows for further inspection of the NDVI data to assess if the unexpected patterns in the NDVI time series are due to errors in the data or land cover related effects. 4. Conclusion The use of a double logistic function or asymmetric Gaussian function is appropriate for describing vegetation dynamics at high latitudes, as quantified by NDVI time series. Existing methods relying on second-order Fourier functions are unsuited for modelling the short growing seasons in these parts of the world. The method as outlined in this paper holds the additional advantage that the estimated parameters are of phenological significance, and can potentially reveal spatial and temporal trends in NDVI time-series data. While measuring the NDVI of boreal vegetation during winter remains problematic, it can apparently be estimated from the MODIS NDVI data alone, without requiring a vegetation map. The method presented here is the first algorithm specifically designed for vegetation monitoring at very high latitudes. It overcomes problems related to algorithms based on secondorder Fourier series and makes use of MODIS NDVI data to monitor vegetation phenology at a higher spatial resolution and with higher accuracy than hitherto possible. The results show exceptionally early onset of autumn in the northern part of the study area. This might be a natural phenomenon or an artefact of undetected measurement errors. As a next step in our research, we will develop ways to assign confidence intervals to the estimated model parameters. These confidence intervals should integrate the effects of data quality and quantity, the weighting scheme applied, as well as the variation in date of image acquisition and unknown stochastic effects.

Acknowledgements We thank Stein Rune Karlsen and two anonymous referees for their comments on the manuscript. Funding for this research was obtained from the Norwegian Research Council through Phenoclim project.

332

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

References Andres, L., Salas, W., & Skole, D. (1994). Fourier analysis of multi-temporal AVHRR data applied to land cover classification. International Journal of Remote Sensing, 15, 1115 – 1121. Anyambe, A., & Eastman, J. R. (1996). Interannual variability of NDVI over Africa and its relation to El Nino˜-southern oscillation. International Journal of Remote Sensing, 17, 2533 – 2548. Azzali, A., & Menenti, M. (2000). Mapping vegetation – soil complexes in southern Africa using temporal Fourier analysis of NOAA AVHRR NDVI data. International Journal of Remote Sensing, 21, 973 – 996. Bacour, C., & Bre´on, F. M. (2005). Variability of biome reflectance directional signatures as seen by POLDER. Remote Sensing of Environment, 89, 80 – 95. Badeck, F. -W., Bondeau, A., Bo¨ttcher, K., Doktor, D., Lucht, W., Schaber, J., et al. (2004). Responses of spring phenology to climate change. New Phytologist, 162, 295 – 309. Baret, F., & Guyot, G. (1991). Potentials and limits of vegetation indices for LAI and APAR assessment. Remote Sensing of Environment, 35, 161 – 173. Braswell, B. H., Schimel, D. S., Linder, E., & Moore III, B. (1997). The response of global terrestrial ecosystems to interannual temperature variability. Science, 278, 870 – 872. Chen, J. M., & Cihlar, J. (1996). Retrieving leaf area index of boreal forests using Landsat TM images. Canadian Journal of Botany, 55, 153 – 162. Chen, X., Tan, Z., Schwartz, M. D., & Xu, C. (2000). Determining the growing season of land vegetation on the basis of plant phenology and satellite data in Northern China. International Journal of Biometeorology, 44, 97 – 101. Chen, X., Xu, C., & Tan, Z. (2001). TAKE DOWN An analysis of relationships among plant community phenology and seasonal metrics of Normalized Difference Vegetation Index in the northern part of the monsoon region of China. International Journal of Biometeorology, 45, 170 – 177. Choudhury, B. J., Ahmed, N. U., Idso, S. B., Reginato, R. J., & Daughtry, C. S. T. (1994). Relations between evaporation coefficients and vegetation indices studied by model simulations. Remote Sensing of Environment, 50, 1 – 17. Cihlar, J. (1996). Identification of contaminated pixels in AVHRR composite images for studies of land biosphere. Remote Sensing of Environment, 56, 149 – 153. Cihlar, J., St.-Laurent, L., & Dyer, J. A. (1991). Relation between the normalized difference vegetation index and ecological variables. Remote Sensing of Environment, 35, 279 – 298. Cramer, W. (1997). Modelling the possible impact of climate change on broad-scale vegetation structure: Examples from northern Europe. In W. Oechel, T. Callaghan, T. Gilmanov, J. I. Holten, B. Maxwell, U. Molau, & B. Sveinbjo¨rnsson (Eds.), Global change and arctic terrestrial ecosystems (pp. 381 – 401). New York’ Springer. Deering, D. W., Ahmad, S. P., Eck, T. F., & Banerjee, B. P. (1995). Temporal attributes of the bidirectional reflectance for three boreal forest canopies. IEEE Proceedings of IGARSS ’95, Florence, 1239 – 1241. Deering, D. W., Eck, T. F., & Banerjee, B. (1999). Characterization of the reflectance anisotropy of three boreal forest canopies in spring – summer. Remote Sensing of Environment, 67, 205 – 229. Duchemin, B., Goubier, J., & Courrier, G. (1999). Monitoring phenological key stages and cycle duration of temperate deciduous forest ecosystems with NOAA/AVHRR data. Remote Sensing of Environment, 67, 68 – 82. Fischer, A. (1994a). A model for the seasonal variations of vegetation indices in coarse resolution data and its inversion to extract crop parameters. Remote Sensing of Environment, 48, 220 – 230. Fischer, A. (1994b). A simple model for the temporal variations of NDVI at regional scale over agricultural countries. Validation with ground radiometric measurements. International Journal of Remote Sensing, 15, 1421 – 1446. Foley, J. A., Kutzbach, J. E., Coe, M. T., & Levis, S. (1994). Feedbacks between climate and boreal forests during the Holocene epoch. Nature, 371, 52 – 54. Forster, B. C. (1984). Derivation of atmospheric correction procedures for Landsat MSS with particular reference to urban data. International Journal of Remote Sensing, 5, 799 – 817.

Gamon, J. A., Huemmrich, K. F., Peddle, D. R., Chen, J., Fuentes, D., Hall, F. G., et al. (2004). Remote sensing in BOREAS: Lessons learned. Remote Sensing of Environment, 89, 139 – 162. GIMMS (Global Inventory Monitoring and Modeling Studies) at the National Aeronautics and Space Administration (2005), GIMMS NDVI. http://ltpwww.gsfc.nasa.gov/gimms/htdocs/ndvi/ndviGIMMS.htm Goward, S. N., & Dye, D. G. (1987). Evaluating North American net primary productivity with satellite observations. Advances in Space Research, 7, 165 – 174. Goward, S. N., & Prince, S. D. (1995). Transient effects of climate on vegetation dynamics: Satellite observations. Journal of Biogeography, 22, 549 – 563. Gujarati, D. N. (2003). Autocorrelation: What happens if the error terms are correlated? Basic econometrics (pp. 441 – 505). New York’ McGraw-Hill. Guyot, G., Gugon, D., & Riom, J. (1989). Factors affecting the spectral response of forest canopies: A review. Geocarta International, 3, 43 – 60. Hirosawa, Y., Marsh, S. E., & Kliman, D. H. (1996). Application of standardized principle component analysis to land-cover characterization using multitemporal AVHRR data. Remote Sensing of Environment, 58, 267 – 281. Holben, B. N. (1986). Characteristics of maximum-value composite images for temporal AVHRR data. International Journal of Remote Sensing, 7, 1435 – 1445. Hu, B., Inannen, K., & Miller, J. R. (2000). Retrieval of leaf area index and canopy closure from CASI data over the BOREAS flux tower sites. Remote Sensing of Environment, 74, 255 – 274. Huemmrich, K. F., & Goward, S. N. (1997). Vegetation canopy PAR absorptance and NDVI: An assessment for ten tree species with the SAIL model. Remote Sensing of Environment, 61, 254 – 269. Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83, 195 – 213. Jakubauskas, M. E., Legates, D. R., & Kastens, J. H. (2001). Harmonic analysis of time-series AVHRR NDVI data. Photogrammetric Engineering and Remote Sensing, 67, 461 – 470. James, M. E., & Kalluri, S. N. V. (1994). The Pathfinder AVHRR land data set: An improved coarse resolution data set for terrestrial monitoring. International Journal of Remote Sensing, 15, 3347 – 3364. Jo¨nsson, P., & Eklundh, L. (2002). Seasonality extraction by function fitting to time-series of satellite sensor data. IEEE Transactions on Geoscience and Remote Sensing, 40, 1824 – 1832. Jo¨nsson, P., & Eklundh, L. (2004). TIMESAT—a program for analyzing time-series of satellite sensor data. Computers and Geosciences, 30, 833 – 845. Kaduk, J., & Heimann, M. (1996). A prognostic classification of phenology model for global terrestrial carbon cycle models. Climate Research, 6, 1 – 19. Karlsen, S. R., Elvebakk, A., Høgda, K. A., & Johansen, B. (accepted for publication), Satellite based mapping of the growing season and bioclimatic zones in Fennoscandia. Global Ecology and Biogeography. Kimes, D. S., Holben, B. N., Tucker, C. J., & Newcomb, W. W. (1984). Optimal directional view angles for remote sensing images. International Journal of Remote Sensing, 5, 887 – 908. Kirschbaum, M. U. F., Fischlin, A., Cannell, M. G. R., Cruz, R. V. O., & Galinski, W. (1996). Climate change impacts on forests. In R. T. Watson, M. C. Zinyowera, & R. H. Moss (Eds.), Climate change 1995: Impacts, adaptations and mitigation of climate change: Scientific-technical analyses. Contribution of working group II to the second assessment report of the intergovernmental panel on climate change (pp. 95 – 129). Cambridge, United Kingdom and New York, NY, USA’ Cambridge University Press. Kramer, K. (1995). Phenotypic plasticity of the phenology of seven European species in relation to climate warming. Plant, Cell and Environment, 18, 93 – 104. Lloyd, D. (1990). A phenological classification of terrestrial vegetation cover using shortwave vegetation index imagery. International Journal of Remote Sensing, 11, 2269 – 2279.

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334 Los, S. O. (1998). Linkages between global vegetation and climate. An analysis based on NOAA Advanced Very High Resolution Radiometer data. Amsterdam, Vrije Universiteit, PhD. dissertation, p. 179. Los, S. O., Justice, C. O., & Tucker, C. J. (1994). A global 1 degree by 1 degree NDVI data set for climate studies derived from the GIMMS continental NDVI. International Journal of Remote Sensing, 15, 3493 – 3518. Lucht, W., Prentice, R. B., Myneni, R. B., Sitch, S., Friedlingstein, P., Cramer, W., et al. (2002). Climate control of the high-latitude vegetation greening trend and Pinatubo effect. Science, 296, 1687 – 1689. McCloy, K. R., & Lucht, W. (2004). Comparative evaluation of seasonal patterns in long time series of satellite image data and simulations of a global vegetation model. IEEE Transactions on Geoscience and Remote Sensing, 42, 140 – 153. Menenti, M., Azzali, S., Verhoef, W., & van Swol, R. (1993). Mapping agroecological zones and time lag in vegetation growth by means of Fourier analysis of time series of NDVI images. Advances in Space Research, 13, 233 – 237. Menzel, A. (2002). Phenology: Its importance to the global change community, an editorial comment. Climatic Change, 54, 379 – 398. Moody, A., & Johnson, D. M. (2001). Land-surface phenologies from AVHRR using the discrete Fourier transform. Remote Sensing of Environment, 75, 305 – 323. Myneni, R. B., Keeling, C. D., Tucker, C. J., Asrar, G., & Nemani, R. R. (1997). Increased plant growth in the northern high latitudes from 1981 – 1991. Nature, 386, 698 – 702. Nagler, P. L., Glenn, E. P., Lewis Thompson, T., & Huete, A. (2004). Leaf area index and normalized difference vegetation index as predictors of canopy characteristics and light interception by riparian species on the Lower Colorado River. Agricultural and Forest Meteorology, 125, 1 – 17. NML (Nettverk for miljølære) (2005), Fenologi- Brukerstyrt visning. Bergen: Skolelaboratoriet i realfag, Universitetet i Bergen: http://nml.uib.no/data/ ut/land/natur/ln14/?visselektert=1 Olsson, L., & Eklundh, L. (1994). Fourier-series for analysis of temporal sequences of satellite sensor imagery. International Journal of Remote Sensing, 15, 3735 – 3741. Paruelo, J. M., & Lauenroth, W. K. (1998). Interannual variability of NDVI and its relationship to climate for North American shrublands and grasslands. Journal of Biogeography, 25, 721 – 733. Peterson, D. L. (1991). Applications in forest science and management. In G. Asrar (Ed.), Theory and applications of optical remote sensing (pp. 429 – 473). New York’ Wiley. Pettorelli, N., Vik, J. O., Mysterud, A., Gaillard, J. -M., Tucker, C. J., & Stenseth, N. C. (2005). Using the satellite-derived NDVI to assess ecological responses to environmental change. Trends in Ecology and Evolution, 20, 503 – 510. Rahman, H., Pinty, B., & Verstraete, M. M. (1993). Coupled SurfaceAtmosphere Reflectance (CSAR) Model 2. Semi-empirical surface model usable with NOAA AVHRR advanced very high resolution radiometer data. Journal of Geophysical Research, 98, 20791 – 20801. Reed, B. C., Brown, J. F., VanderZee, D., Loveland, T. R., Merchant, J. W., & Ohlen, D. O. (1994). Measuring phenological variability from satellite imagery. Journal of Vegetation Science, 5, 703 – 714. Richards, J. A. (1986). Remote sensing and digital image analysis (p. 281). New York’ Springer-Verlag. Roerink, G. J., Menenti, M., & Verhoef, W. (2000). Reconstructing cloudfree NDVI composites using Fourier analysis of time series. International Journal of Remote Sensing, 21, 1911 – 1917. Running, S. W., & Nemani, R. R. (1988). Relating seasonal patterns of the AVHRR vegetation index to simulated photosynthesis and transpiration of forests in different climates. Remote Sensing of Environment, 24, 347 – 367. Running, S. W., Nemani, R. R., Peterson, D. L., Band, L. E., Potts, D. F., Pierce, L. L., et al. (1989). Mapping regional forest evapotranspiration and photosynthesis by coupling satellite data with ecosystem simulation. Ecology, 70, 1090 – 1101. Sakamoto, T., Yokozawa, M., Toritani, H., Shibayama, M., Ishitsuka, N., & Ohno, H. (2005). A crop phenology detection method using time-series MODIS data. Remote Sensing of Environment, 96, 366 – 374.

333

Schlerf, M. & Atzberger, C. (2006). Inversion of a forest reflectance model to estimate structural canopy variables from hyperspectral remote sensing data. Remote Sensing of Environment. 100, 281-294 Sellers, P. J., Dickinson, R. E., Randall, D. A., Betts, A. K., Hall, F. G., Berry, J. A., et al. (1997). Modeling the exchanges of energy and water, and carbon between continents and atmosphere. Science, 275, 502 – 509. Sellers, P. J., Los, S. O., Tucker, C. J., Justice, C. O., Dazlich, D. A., Collatz, G. J., et al. (1996). A revised land surface parameterization (SiB2) for atmospheric GCMS: Part II. The generation of global fields of terrestrial biophysical parameters from satellite data. Journal of Climate, 9, 706 – 737. Sellers, P. J., Tucker, C. J., Collatz, G. J., Los, S. O., Justice, C. O., Dazlich, D. A., et al. (1994). A global 1- by 1- NDVI data set for climate studies: 2. The generation of global fields of terrestrial biophysical parameters from the NDVI. International Journal of Remote Sensing, 151, 3519 – 3545. Shabanov, N. V., Zhou, L., Knyazikhin, Y., Myneni, R. B., & Tucker, C. J. (2002). Analysis of interannual changes in northern vegetation activity observed in AVHRR data from 1981 to 1994. IEEE Transactions on Geoscience and Remote Sensing, 40, 115 – 130. Shibayama, M., Salli, A., Ha¨me, T., Iso-Iivari, L., Heino, S., Alanen, M., et al. (1999). Detecting phenophases of subarctic shrub canopies by using automated reflectance measurements. Remote Sensing of Environment, 67, 160 – 180. Skre, O. (2001). Climate change on mountain birch ecosystems. In F. E. Wielgolaski (Ed.), Nordic mountain birch ecosystems (pp. 343 – 357). Paris, New York and London’ UNESCO, Parthenon Publ. Group. Smit-Spinks, B., Swanson, B. T., & Markhart, A. H. I. (1985). The effect of photoperiod and thermoperiod on cold acclimation and growth of Pinus sylvestris. Canadian Journal of Botany, 15, 453 – 460. Sto¨ckli, R., & Vidale, P. L. (2004). European plant phenology and climate as seen in a 20-year AVHRR land-surface parameter dataset. International Journal of Remote Sensing, 10, 3303 – 3330. Taulavuori, K. M., Taulavuori, E. B., Skrer, O., Nilsen, J., Igeland, B., & Laine, K. M. (2004). Dehardning of mountain birch (Betula pubescens ssp. czerepanovii) ecotypes at elevated winter temperatures. New Phytologist, 162, 427. Tømmervik, H., Høgda, K. A., & Karlsen, S. R. (2001). Using remote sensing for detection of caterpillar outbreaks in mountain birch forests, a new approach. In F. -E. Wielgolaski (Ed.), Nordic mountain birch ecosystems (pp. 241 – 249). Paris and New York, London’ UNESCO and Parthenon Publishing Group. Tucker, C. J., Slayback, D. A., Pinzon, J. E., Los, S. O., Myneni, R. B., & Taylor, M. G. (2001). Higher northern latitude normalized difference vegetation index and growing season trends from 1982 to 1999. International Journal of Biometeorology, 45, 184 – 190. Tucker, C. J., Vanpraet, C., Boerwinkel, E., & Gaston, A. (1983). Satellite remote sensing of total dry matter production in the Senegalese Sahel. Remote Sensing of Environment, 13, 461 – 474. Tuhkanen, S. (1980). Climatic parameters and indices in plant geography. Acta Phytogeographica Suecica, 67, 1 – 105. Vermote, E. F., El Saleous, N. Z., & Justice, C. O. (2002). Atmospheric correction of MODIS data in the visible to middle infrared: first results. Remote Sensing of Environment, 83, 97 – 111. Viovy, N., Arino, O., & Belward, A.S. (1992). The best index slope extraction (BISE): a method for reducing noise in NDVI time-series. International Journal of Remote Sensing, 13, 1585 – 1590. Wielgolaski, F. -E., & Inouye, D. W. (2003). High latitude climates. In M. D. Schwartz (Ed.), Phenology: An integrative environmental science (pp. 175 – 194). Dordrecht, Boston, London’ Kluwer Academic Publishers. Wylie, B. K., Johnson, D. A., Laca, E., Saliendra, N. Z., Gilmanov, T. G., Reed, B. C., et al. (2003). Calibration of remotely sensed, coarse resolution NDVI to CO2 fluxes in a sagebrush-steppe ecosystem. Remote Sensing of Environment, 85, 243 – 255. Xin, J., Yu, Z., van Leeuwen, L., & Driessen, P. M. (2002). Mapping crop key phenological stages in the North China Plain using NOAA time series

334

P.S.A. Beck et al. / Remote Sensing of Environment 100 (2006) 321 – 334

images. International Journal of Applied Earth Observation and Geoinformation, 4, 109 – 117. Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C. F., Gao, F., et al. (2003). Monitoring vegetation phenology using MODIS. Remote Sensing of Environment, 84, 471 – 475. Zhou, L., Kaufmann, R. K., Tian, Y., Myneni, R. B., & Tucker, C. J. (2003). Relation between interannual variations in satellite measures of northern

forest greenness and climate between 1982 and 1999. Journal of Geophysical Research, 108. doi:10.1029/2002JD002510. Zhou, L., Tucker, C. J., Kaufmann, R. K., Slayback, D., Shabanov, N. V., & Myneni, R. B. (2001). Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981 to 1999. Journal of Geophysical Research, 106, 20069 – 20083.