1 Detecting trend and seasonal

2 changes in satellite image time series

3 Jan Verbesselt∗,a , Rob Hyndmanb , Glenn Newnhama , Darius Culvenora

a
4 Remote sensing team, CSIRO Sustainable Ecosystems,
5 Private Bag 10, Melbourne VIC 3169, Australia
b
6 Department of Econometrics and Business Statistics,
7 Monash University, Melbourne VIC 3800, Australia

8 Abstract
9 A wealth of remotely sensed image time series covering large areas is now available to the
10 earth science community. Change detection methods are often not capable of detecting land
11 cover changes within time series that are heavily influenced by seasonal climatic variations.
12 Detecting change within the trend and seasonal components of time series enables the
13 classification of different types of changes. Changes occurring in the trend component often
14 indicate disturbances (e.g. fires, insect attacks), while changes occurring in the seasonal
15 component indicate phenological changes (e.g. change in land cover type). A generic
16 change detection approach is proposed for time series by detecting and characterizing
17 Breaks For Additive Seasonal and Trend (BFAST). BFAST integrates the decomposition
18 of time series into trend, seasonal, and remainder components with methods for detecting
19 change within time series. BFAST iteratively estimates the time and number of changes,
20 and characterizes change by its magnitude and direction. We tested BFAST by simulating
21 16-day Normalized Difference Vegetation Index (NDVI) time series with varying amounts
22 of seasonality and noise, and by adding abrupt changes at different times and magnitudes.
23 This revealed that BFAST can robustly detect change with different magnitudes (> 0.1
24 NDVI) within time series with different noise levels (0.01–0.07 𝜎) and seasonal amplitudes
25 (0.1–0.5 NDVI). Additionally, BFAST was applied to 16-day NDVI Moderate Resolution
26 Imaging Spectroradiometer (MODIS) composites for a forested study area in south eastern
27 Australia. This showed that BFAST is able to detect and characterize spatial and temporal
28 changes in a forested landscape. BFAST is not specific to a particular data type and can be
29 applied to time series without the need to normalize for land cover types, select a reference
30 period, or change trajectory. The method can be integrated within monitoring frameworks
31 and used as an alarm system to flag when and where changes occur.
32 Key words: Change detection, NDVI, time series, trend analysis, MODIS, piecewise
33 linear regression, vegetation dynamics, phenology

∗
Corresponding author. Ph: +61395452265 ; Fax : +61395452448
Email addresses: Jan.Verbesselt@csiro.au (Jan Verbesselt), Rob.Hyndman@buseco.monash.edu.au
(Rob Hyndman), Glenn.Newnham@csiro.au (Glenn Newnham), Darius.Culvenor@csiro.au (Darius
Culvenor)

Preprint submitted to Remote Sensing of Environment September 7, 2009

34 1. Introduction

35 Natural resource managers, policy makers and researchers demand knowledge of land
36 cover changes over increasingly large spatial and temporal extents for addressing many
37 pressing issues such as global climate change, carbon budgets, and biodiversity (DeFries
38 et al., 1999; Dixon et al., 1994). Detecting and characterizing change over time is the
39 natural first step toward identifying the driver of the change and understanding the change
40 mechanism. Satellite remote sensing has long been used as a means of detecting and
41 classifying changes in the condition of the land surface over time (Coppin et al., 2004; Lu
42 et al., 2004). Satellite sensors are well-suited to this task because they provide consistent
43 and repeatable measurements at a spatial scale which is appropriate for capturing the
44 effects of many processes that cause change, including natural (e.g. fires, insect attacks)
45 and anthropogenic (e.g. deforestation, urbanization, farming) disturbances (Jin and Sader,
46 2005).
47 The ability of any system to detect change depends on its capacity to account for
48 variability at one scale (e.g. seasonal variations), while identifying change at another
49 (e.g. multi-year trends). As such, change in ecosystems can be divided into three classes:
50 (1) seasonal change, driven by annual temperature and rainfall interactions impacting
51 plant phenology or proportional cover of land cover types with different plant phenology;
52 (2) gradual change such as interannual climate variability (e.g. trends in mean annual
53 rainfall) or gradual change in land management or land degradation; and (3) abrupt change,
54 caused by disturbances such as deforestation, urbanization, floods, and fires.
55 Although the value of remotely sensed long term data sets for change detection has
56 been firmly established (de Beurs and Henebry, 2005), only a limited number of time series
57 change detection methods have been developed. Two major challenges stand out. First,
58 methods must allow for the detection of changes within complete long term data sets while
59 accounting for seasonal variation. Estimating change from remotely sensed data is not
60 straightforward, since time series contain a combination of seasonal, gradual and abrupt
61 changes, in addition to noise that originates from remnant geometric errors, atmospheric
62 scatter and cloud effects (Roy et al., 2002). Thorough reviews of existing change detection
63 methods by Coppin et al. (2004) and Lu et al. (2004) have shown, however, that most
64 methods focus on short image time series (only 2–5 images). The risk of confounding
65 variability with change is high with infrequent images, since disturbances can occur in

2
66 between image acquisitions (de Beurs and Henebry, 2005). Several approaches have been
67 proposed for analyzing image time series, such as Principal Component Analysis (PCA)
68 (Crist and Cicone, 1984), wavelet decomposition (Anyamba and Eastman, 1996), Fourier
69 analysis (Azzali and Menenti, 2000) and Change Vector Analysis (CVA) (Lambin and
70 Strahler, 1994). These time series analysis approaches discriminate noise from the signal
71 by its temporal characteristics but involve some type of transformation designed to isolate
72 dominant components of the variation across years of imagery through the multi-temporal
73 spectral space. The challenge of these methods is the labeling of the change components,
74 because each analysis depends entirely on the specific image series analyzed. Compared to
75 PCA, Fourier analysis, and wavelet decomposition, CVA allows the interpretation of change
76 processes, but can still only be performed between two periods of time (e.g. between years
77 or growing seasons) (Lambin and Strahler, 1994), which makes the analysis dependent
78 on the selection of these periods. Furthermore, change in time series if often masked by
79 seasonality driven by yearly temperature and rainfall variation. Existing change detection
80 techniques minimize seasonal variation by focussing on specific periods within a year (e.g.
81 growing season) (Coppin et al., 2004), temporally summarizing time series data (Bontemps
82 et al., 2008; Fensholt et al., 2009) or normalizing reflectance values per land cover type
83 (Healey et al., 2005) instead of explicitly accounting for seasonality.
84 Second, change detection techniques need to be independent of specific thresholds or
85 change trajectories. Change detection methods that require determination of thresholds
86 often produce misleading results due to different spectral and phenological characteristics
87 of land cover types (Lu et al., 2004). The determination of thresholds adds significant cost
88 to efforts to expand change detection to large areas. Trajectory based change detection has
89 been proposed to move towards a threshold independent change detection by characterizing
90 change by its temporal signature (Hayes and Cohen, 2007; Kennedy et al., 2007). This
91 approach requires definition of the change trajectory specific for the type of change to be
92 detected and spectral data to be analyzed (e.g. short-wave infrared or near-infrared based
93 indices). Furthermore, the method will only function if the observed spectral trajectory
94 matches one of the hypothesized trajectories. Trajectory based change detection can
95 be interpreted as a supervised change detection method while there is a need for an
96 unsupervised, more generic, change detection approach independent of the data type and
97 change trajectory.

3
 98 The purpose of this research was to develop a generic change detection approach for
99 time series, involving the detection and characterization of Breaks For Additive Seasonal
100 and Trend (BFAST). BFAST integrates the iterative decomposition of time series into
101 trend, seasonal and noise components with methods for detecting changes, without the
102 need to select a reference period, set a threshold, or define a change trajectory. The main
103 objectives are:

104 (1) the detection of multiple abrupt changes in the seasonal and trend components of the
105 time series; and
106 (2) the characterization of gradual and abrupt ecosystem change by deriving the time,
107 magnitude, and direction of change within the trend component of the time series.

108 We assessed BFAST for a large range of ecosystems by simulating Normalized Difference
109 Vegetation Index (NDVI) time series with varying amounts of seasonal variation and noise,
110 and by adding abrupt changes with different magnitudes. We applied the approach on
111 MODIS 16-day image composites (hereafter called 16-day time series) to detect major
112 changes in a forested area in south eastern Australia. The approach is not specific to
113 a particular data type and could be applied to detect and characterize changes within
114 other remotely sensed image time series (e.g. Landsat) or be integrated within monitoring
115 frameworks and used as an alarm system to provide information on when and where
116 changes occur.

117 2. Iterative change detection

118 We propose a method that integrates the iterative decomposition of time series into
119 trend, seasonal and noise components with methods for detecting and characterizing
120 changes (i.e. breakpoints) within time series. Standard time series decomposition methods
121 assume that trend and seasonal components are smooth and slowly changing, and so
122 these are not directly applicable to the problem of identifying change. For example, the
123 Seasonal-Trend decomposition procedure (STL) is capable of flexibly decomposing a series
124 into trend, seasonal and remainder components based on a LOcally wEighted regreSsion
125 Smoother (LOESS) (Cleveland et al., 1990). This smoothing prevents the detection of
126 changes within time series.

4
127 2.1. Decomposition model

128 We propose an additive decomposition model that iteratively fits a piecewise linear
129 trend and seasonal model. The general model is of the form 𝑌𝑡 = 𝑇𝑡 + 𝑆𝑡 + 𝑒𝑡 , 𝑡 = 1, . . . , 𝑛,
130 where 𝑌𝑡 is the observed data at time 𝑡, 𝑇𝑡 is the trend component, 𝑆𝑡 is the seasonal
131 component, and 𝑒𝑡 is the remainder component. The remainder component is the remaining
132 variation in the data beyond that in the seasonal and trend components (Cleveland et al.,
133 1990). It is assumed that 𝑇𝑡 is piecewise linear, with break points 𝑡∗1 , . . . , 𝑡∗𝑚 and define
134 𝑡∗0 = 0, so that
𝑇𝑡 = 𝛼𝑗 + 𝛽𝑗 𝑡 (1)

135 for 𝑡∗𝑗−1 < 𝑡 ≤ 𝑡∗𝑗 and where 𝑗 = 1, . . . , 𝑚. The intercept and slope of consecutive linear
136 models, 𝛼𝑗 and 𝛽𝑗 , can be used to derive the magnitude and direction of the abrupt change
137 (hereafter referred to as magnitude) and slope of the gradual change between detected break
138 points. The magnitude of an abrupt change at a breakpoint is derived by the difference
139 between 𝑇𝑡 at 𝑡∗𝑗−1 and 𝑡∗𝑗 , so that

Magnitude = (𝛼𝑗−1 − 𝛼𝑗 ) + (𝛽𝑗−1 − 𝛽𝑗 )𝑡 (2)

140 and the slopes of the gradual change before and after a break point are 𝛽𝑗−1 and 𝛽𝑗 .
141 This technique represents a simple yet robust way to characterize changes in time series.
142 Piecewise linear models, as a special case of non-linear regression (Venables and Ripley,
143 2002), are often used as approximations to complex phenomena to extract basic features of
144 the data (Zeileis et al., 2003).
145 Similarly, the seasonal component is fixed between break points, but can vary across
146 break points. Furthermore, the seasonal break points may occur at different times from
147 the break points detected in the trend component. Let the seasonal break points be given
148 by 𝑡# # # # #
1 , . . . , 𝑡𝑝 , and define 𝑡0 = 0. Then for 𝑡𝑗−1 < 𝑡 ≤ 𝑡𝑗 , we assume that

⎧
⎨ 𝛾𝑖,𝑗 if time 𝑡 is in season 𝑖, 𝑖 = 1, . . . , 𝑠 − 1;
𝑆𝑡 = (3)
⎩ − ∑𝑠−1 𝛾 if time 𝑡 is in season 0,
𝑖=1 𝑖,𝑗

149 where 𝑠 is the period of seasonality (e.g. number of observations per year) and 𝛾𝑖,𝑗 denotes
150 the effect of season 𝑖. Thus, the sum of the seasonal component, 𝑆𝑡 across 𝑠 successive
151 times is exactly zero for 𝑡# #
𝑗−1 < 𝑡 ≤ 𝑡𝑗 . This prevents apparent changes in trend being

5
152 induced by seasonal breaks happening in the middle of a seasonal cycle. The seasonal term
153 can be re-expressed as
𝑠−1
∑
𝑆𝑡 = 𝛾𝑖,𝑗 (𝑑𝑡,𝑖 − 𝑑𝑡,0 ) (4)
𝑖=1

154 where 𝑑𝑡,𝑖 = 1 when 𝑡 is in season 𝑖 and 0 otherwise. Therefore, if 𝑡 is in season 0, then
155 𝑑𝑡,𝑖 − 𝑑𝑡,0 = −1. For all other seasons, 𝑑𝑡,𝑖 − 𝑑𝑡,0 = 1 when 𝑡 is in season 𝑖 =
∕ 0. 𝑑𝑡,𝑖 is
156 often referred to as a seasonal dummy variable (Makridakis et al., 1998, pp.269-274); it
157 has two allowable values (0 and 1) to account for the seasons in a regression model. The
158 regression model expressed by Eq. 4 can also be interpreted as a model without intercept
159 that contains 𝑠 − 1 seasonal dummy variables.

160 2.2. Iterative algorithm to detect break points

161 Our method is similar to that proposed by Haywood and Randal (2008) for use with
162 monthly tourism data. Following Haywood and Randal (2008), we estimate the trend and
163 seasonal components iteratively. However, we differ from their method by: (1) using STL to
164 estimate the initial seasonal component (𝑆ˆ𝑡 ); (2) using a robust procedure when estimating
165 the coefficients 𝛼𝑗 , 𝛽𝑗 and 𝛾𝑖,𝑗 ; (3) using a preliminary structural change test; and (4) forcing
166 the seasonal coefficients to always sum to 0 (rather than adjusting them afterward). An
167 alternative approach proposed by Shao and Campbell (2002) combines the seasonal and
168 trend term in a piecewise linear regression model without iterative decomposition. This
169 approach does not allow for an individual estimation of breakpoints in the seasonal and
170 trend component.
171 Sequential test methods for detecting break points (i.e. abrupt changes) in a time series
172 have been developed, particularly within econometrics (Bai and Perron, 2003; Zeileis et al.,
173 2003). These methods also allow linear models to be fitted to sections of a time series, with
174 break points at the times where the changes occur. The optimal position of these breaks
175 can be determined by minimizing the residual sum of squares, and the optimal number of
176 breaks can be determined by minimizing an information criterion. Bai and Perron (2003)
177 argue that the Akaike Information Criterion usually overestimates the number of breaks,
178 but that the Bayesian Information Criterion (BIC) is a suitable selection procedure in
179 many situations (Zeileis et al., 2002, 2003; Zeileis and Kleiber, 2005). Before fitting the
180 piecewise linear models and estimating the breakpoints it is recommended to test whether
181 breakpoints are occurring in the time series (Bai and Perron, 2003). The ordinary least

6
182 squares (OLS) residuals-based MOving SUM (MOSUM) test, is selected to test for whether
183 one or more breakpoints are occurring (Zeileis, 2005). If the test indicates significant
184 change (𝑝 < 0.05), the break points are estimated using the method of Bai and Perron
185 (2003), as implemented by Zeileis et al. (2002), where the number of breaks is determined
186 by the BIC, and the date and confidence interval of the date for each break are estimated.
187 The iterative procedure begins with an estimate of 𝑆ˆ𝑡 by using the STL method, where
188 𝑆ˆ𝑡 is estimated by taking the mean of all seasonal sub-series (e.g. for a monthly time series
189 the first subseries contains the January values). Then it follows these steps.

190 Step 1 If the OLS-MOSUM test indicates that breakpoints are occurring in the trend
191 component, the number and position of the trend break points (𝑡∗1 , . . . , 𝑡∗𝑚 ) are
192 estimated from the seasonally adjusted data, 𝑌𝑡 − 𝑆ˆ𝑡 .

193 Step 2 The trend coefficients, 𝛼𝑗 and 𝛽𝑗 for 𝑗 = 1, . . . , 𝑚, are then computed using robust
194 regression of Eq. 1 based on M-estimation (Venables and Ripley, 2002). The trend
195 ˆ 𝑗 + 𝛽ˆ𝑗 𝑡 for 𝑡 = 𝑡∗𝑗−1 + 1, . . . , 𝑡∗𝑗 .
estimate is then set to 𝑇ˆ𝑡 = 𝛼

196 Step 3 If the OLS-MOSUM test indicates that breakpoints are occurring in the seasonal
197 component, the number and position of the seasonal break points (𝑡# #
1 , . . . , 𝑡𝑝 ) are

198 estimated from the detrended data, 𝑌𝑡 − 𝑇ˆ𝑡 .

199 Step 4 The seasonal coefficients, 𝛾𝑖,𝑗 for 𝑗 = 1, . . . , 𝑚 and 𝑖 = 1, . . . , 𝑠 − 1, are then
200 computed using a robust regression of Eq. 4 based on M-estimation. The seasonal
estimate is then set to 𝑆ˆ𝑡 = 𝑠−1 ˆ𝑖,𝑗 (𝑑𝑡,𝑖 − 𝑑𝑡,0 ) for 𝑡 = 𝑡# #
∑
201
𝑖=1 𝛾 𝑗−1 + 1, . . . , 𝑡𝑗 .

202 These steps are iterated until the number and position of the breakpoints are unchanged.
203 We have followed the recommendations of Bai and Perron (2003) and Zeileis et al. (2003)
204 concerning the fraction of data needed between the breaks. For 16-day time series, we used
205 a minimum of one year of data (i.e. 23 observations) between successive change detections,
206 corresponding to 12% of a 9 year data span (2000–2008). This means that if two changes
207 occur within a year, only the most significant change will be detected.

208 3. Validation

209 The proposed approach can be applied to a variety of time series, and is not restricted
210 to remotely sensed vegetation indices. However, validation has been conducted using

7
211 Normalized Difference Vegetation Index (NDVI) time series, the most widely used vegetation
212 index in medium to coarse scale studies. The NDVI is a relative and indirect measure of
213 the amount of photosynthetic biomass, and is correlated with biophysical parameters such
214 as green leaf biomass and the fraction of green vegetation cover, whose behavior follows
215 annual cycles of vegetation growth (Myneni et al., 1995; Tucker, 1979).
216 We validated BFAST by (1) simulating 16-day NDVI time series, and (2) applying
217 the method to 16-day MODIS satellite NDVI time series (2000–2008). Validation of
218 multi-temporal change-detection methods is often not straightforward, since independent
219 reference sources for a broad range of potential changes must be available during the change
220 interval. Field validated single-date maps are unable to represent the type and number
221 of changes detected (Kennedy et al., 2007). We simulated 16-day NDVI time series with
222 different noise, seasonality, and change magnitudes in order to robustly test BFAST in a
223 controlled environment. However, it is challenging to create simulated time series that
224 approximate remotely sensed time series which contain combined information on vegetation
225 phenology, interannual climate variability, disturbance events, and signal contamination
226 (e.g. clouds) (Zhang et al., 2009). Therefore, applying the method to remotely sensed data
227 and performing comparisons with in-situ data remains necessary. In the next two sections,
228 we apply BFAST to simulated and MODIS NDVI time series.

229 3.1. Simulation of NDVI time series

230 NDVI time series are simulated by extracting key characteristics from MODIS 16-
231 day NDVI time series. We selected two MODIS NDVI time series (as described in 3.2)
232 representing a grassland and a pine plantation (Fig. 1), expressing the most different
233 phenology in the study area, to extract seasonal amplitude, noise level, and average value.
234 Simulated NDVI time series are generated by summing individually simulated seasonal,
235 noise, and trend components. First, the seasonal component is created using an asymmetric
236 Gaussian function for each season. This Gaussian-type function has been shown to perform
237 well when used to extract seasonality by fitting the function to time series (Jönsson and
238 Eklundh, 2002). The amplitude of the MODIS NDVI time series was estimated using the
239 range of the seasonal component derived with the STL function, as shown in Fig. 2. The
240 estimated seasonal amplitudes of the real forest and grassland MODIS NDVI time series
241 were 0.1 and 0.5 (Fig. 1). Second, the noise component was generated using a random
242 number generator that follows a normal distribution N(𝜇 = 0, 𝜎 = 𝑥), where the estimated

8
243 𝑥 values were 0.04 and 0.02, to approximate the noise within the real grass and forest
244 MODIS NDVI time series (Lhermitte et al., submitted). Vegetation index specific noise was
245 generated by randomly replacing the white noise by noise with a value of −0.1, representing
246 cloud contamination that often remains after atmospheric correction and cloud masking
247 procedures. Third, the real grass and forest MODIS NDVI time series were approximated
248 by selecting constant values 0.6 and 0.8 and summing them with the simulated noise and
249 seasonal component. A comparison between real and simulated NDVI time series is shown
250 in Fig. 1.
251 Based on the parameters required to simulate NDVI time series similar to the real grass
252 and forest MODIS NDVI time series (Fig. 1), we selected a range of amplitude and noise
253 values for the simulation study (Table 1). These values are used to simulate NDVI time
254 series of different quality (i.e. varying signal to noise ratios) representing a large range of
255 land cover types.
Table 1: Parameter values for simulation of 16-day NDVI time series

Parameters Values
Amplitude 0.1, 0.3, 0.5
𝜎 Noise 0.01, 0.02, . . . , 0.07
Magnitude −0.3, −0.2, −0.1, 0

256 The accuracy of the method for estimating the number, timing and magnitude of abrupt
257 changes was assessed by adding disturbances with a specific magnitude to the simulated
258 time series. A simple disturbance was simulated by combining a step function with a
259 specific magnitude (Table 1) and linear recovery phase (Kennedy et al., 2007). As such,
260 the disturbance can be used to simulate, for example, a fire in a grassland or an insect
261 attack on a forest. Three disturbances were added to the sum of simulated seasonal, trend,
262 and noise components using simulation parameters in Table 1. An example of a simulated
263 NDVI time series with three disturbances is shown in Fig. 3. A Root Mean Square Error
264 (RMSE) was derived for 500 iterations of all the combinations of amplitude, noise and
265 magnitude of change levels to quantify the accuracy of estimating: (1) the number of
266 detected changes, (2) the time of change, and (3) the magnitude of change.

267 3.2. Spatial application on MODIS image time series

268 We apply BFAST to real remotely sensed time series, and compare the detected changes
269 with a spatial validation data set. BFAST provides information on the number, time,

9
270 magnitude and direction of changes in the trend and seasonal components of a time series.
271 We focussed on the timing and magnitude of major changes occurring within the trend
272 component.
273 We selected the 16-day MODIS NDVI composites with a 250m spatial resolution
274 (MOD13Q1 collection 5), since this product provides frequent information at the spatial
275 scale at which the majority of human-driven land cover changes occur (Townshend and
276 Justice, 1988). The MOD13Q1 16-day composites were generated using a constrained view
277 angle maximum NDVI value compositing technique (Huete et al., 2002). The MOD13Q1
278 images were acquired from the February 24th of 2000 to the end of 2008 (23 images/year
279 except for the year 2000) for a multi-purpose forested study area (Pinus radiata plantation)
280 in South Eastern Australia (Lat. 35.5∘ S, Lon. 148.0∘ E). The images contain data from
281 the red (620–670nm) and near-infrared (NIR, 841–876nm) spectral wavelengths. We used
282 the binary MODIS Quality Assurance flags to select only cloud-free data of optimal quality.
283 The quality flags, however, do not guarantee cloud-free data for the MODIS 250 m pixels
284 since that algorithms used to screen clouds use bands at coarse resolution. Missing values
285 are replaced by linear interpolation between neighboring values within the NDVI series
286 (Verbesselt et al., 2006).
287 The 16-day MODIS NDVI image series were analyzed, and the changes revealed were
288 compared with spatial forest inventory information on the ’year of planting’ of Pinus
289 radiata. Time of change at a 16-day resolution was summarized to a yearly temporal
290 resolution to facilitate comparison with the validation data. The validation protocol was
291 applied under the assumption that no other major disturbances (e.g. tree mortality) would
292 occur that would cause a change in the NDVI time series bigger than the change caused
293 by harvesting and planting activities.

294 4. Results

295 4.1. Simulated NDVI time series

296 Fig. 3 illustrates how BFAST decomposes and fits different time series components. It
297 can be seen that the fitted and simulated components are similar, and that the magnitude
298 and timing of changes in the trend component are correctly estimated. The accuracies
299 (RMSE) of the number of estimated changes are summarized in Fig. 4. Only results for
300 seasonal amplitude 0.1 and 0.5 are shown but similar results were obtained for 0.3 NDVI

10
301 amplitude. Three properties of the method are illustrated. First, the noise level only has
302 an influence on the estimation of the number of changes when the magnitude of the change
303 is -0.1, and is smaller than the overall noise level. The noise level is expressed as 4 𝜎, i.e.
304 99% of the noise range, to enable a comparison with the magnitude (Fig. 4). Second, the
305 noise level does not influence the RMSE when no changes are simulated (magnitude =
306 0), indicating a low commission error independent of the noise level. Third, the seasonal
307 amplitude does not have an influence on the accuracy of change detection. In Fig. 5 only
308 simulation results for an amplitude 0.1 are shown, since similar results were obtained for
309 other amplitudes (0.3 and 0.5). Overall, Fig. 5 illustrates that the RMSE of estimating the
310 time and magnitude of change estimation is small and increases slowly for increasing noise
311 levels. Only when the magnitude of change is small (−0.1) compared to the noise level
312 (> 0.15), the RMSE increases rapidly for increasing noise levels.

313 4.2. Spatial application on MODIS image time series

314 The application of BFAST to MODIS NDVI time series of a Pinus radiata plantation
315 produced estimates of the time and magnitude of major changes. These results are shown
316 spatially in Figs. 6 and 7. The time of change estimated by BFAST is summarized
317 each year to facilitate comparison. Only areas for which we had validation data available
318 were visualized in Figs. 6 and 7. The overall similarity between the time of planting
319 and time of detected change illustrates how BFAST can be used to detect change in a
320 forest plantation (Fig. 6). However, differences in the estimated time of change can be
321 interpreted using differences in the magnitude of change estimated by BFAST. Fig. 7
322 shows that detected changes can have either a positive or a negative magnitude of change.
323 This can be explained by the fact that planting in pine plantations in the study area
324 corresponds with a harvesting operation in the preceding year (personal communication
325 with C. Stone). Harvesting operations cause a significant decrease in the NDVI times series,
326 whereas planting causes a more gradual increase in NDVI. Firstly, if planting occurred
327 before 2002, the NDVI time series did not contain any significant decrease in NDVI caused
328 by the harvesting operations, since the MODIS NDVI time series only started in early
329 2000. BFAST therefore detected change with a positive magnitude, indicating regrowth
330 (Fig. 7), corresponding to a time of change during or later than the plant date (Fig. 6).
331 Fig. 8 (top) illustrates detected changes within a NDVI time series extracted from a single
332 MODIS pixel within a pine plantation with a planting activity during 2001. Secondly,

11
333 if planting occurred after 2003, the time series contained a significant decrease in NDVI
334 caused by the harvesting operations. Major change detected as a consequence are changes
335 corresponding to harvesting preceding the planting operation, and are therefore detected
336 before the planting date (Fig. 6) and have a negative magnitude (Fig. 7). Fig. 8 (middle)
337 illustrates detected changes within a NDVI time series with harvesting operation activity
338 during 2004. These points illustrate BFAST’s capacity to detect and characterize change,
339 but also confirm the importance of simulating time series in a controlled environment, since
340 it is very difficult to find validation data to account for all types of change occurring in
341 ecosystems.
342 Fig. 8 (bottom) shows an example of changes detected by BFAST in an area where
343 harvesting and thinning activities were absent. Fig. 9 illustrates how BFAST decomposed
344 the NDVI time series and fitted seasonal, trend and remainder components. In 2002 and
345 2006 the study area experienced a severe drought, which caused the pine plantations to
346 be stressed and the NDVI to decrease significantly. Severe tree mortality occurred in
347 2006, since trees were drought-stressed and not able to defend themselves against insect
348 attacks (Verbesselt et al., in press). This explains why the change detected in 2006 is
349 bigger (magnitude of the abrupt change) and the recovery (slope of the gradual change) is
350 slower than the change detected in 2003, as shown in (Fig. 9). This example illustrates
351 how the method could be used to detect and characterize changes related to forest health.

352 5. Discussion and further work

353 The main characteristics of BFAST are revealed by testing the approach using simulated
354 time series and by comparing detected changes in 16-day MODIS NDVI time series with
355 spatial forest inventory data. Simulation of NDVI time series illustrated that the iterative
356 decomposition of time series into a seasonal and trend component was not influenced by
357 the seasonal amplitude and by noise levels smaller than the simulated change magnitude.
358 This enabled the robust detection of abrupt and gradual changes in the trend component.
359 As such, full time series can be analyzed without having to select only data during a
360 specific period (e.g. growing season), or can avoid the normalization of reflectance values
361 for each land cover type to minimize seasonal variability (Healey et al., 2005). Seasonal
362 adjustment by decomposing time series, as implemented in the BFAST approach, facilitates
363 the detection of change in the trend component independent of seasonal amplitude or land

12
364 cover type information. Considerations for further research fall into four main categories:

365 (1) Further research is necessary to study BFAST’s sensitivity to detecting phenological
366 change in the seasonal component. This research has focussed on the detection and
367 characterization of changes within the trend component of 16-day NDVI time series.
368 Changes in the seasonal component were not simulated, and BFAST’s sensitivity to
369 detecting seasonal changes using simulated data was not assessed. However, changes
370 occurring in the seasonal component can be detected using BFAST. The application of
371 BFAST to 16-day MODIS NDVI time series on a forested area (40000ha) revealed that
372 seasonal breaks were detected in only 5% of the area. The small number of seasonal
373 breaks occurring in the study area could be explained by the fact that a seasonal
374 change is only detected when a change between land cover types with a significantly
375 different phenology occurs. Time series with a higher temporal resolution (e.g. daily or
376 8-day) could increase the accuracy of detecting seasonal changes but might also impact
377 the ability to detect subtle changes due to higher noise levels. Zhang et al. (2009)
378 illustrated that vegetation phenology can be estimated with high accuracy (absolute
379 error of less than 3 days) in time series with a temporal resolution of 6–16 days, but
380 that accuracy depends on the occurrence of missing values. It is therefore necessary to
381 study BFAST’s capacity to detect phenological change caused by climate variations or
382 land use change in relation to the temporal resolution of remotely sensed time series.
383 (2) Future algorithm improvements may include the capacity to add functionality to identify
384 the type of change with information on the parameters of the fitted piecewise linear
385 models (e.g. intercept and slope). In this study we have focussed on the magnitude of
386 change, but the spatial application on MODIS NDVI time series illustrated that change
387 needs to be interpreted by combining the time and magnitude of change. Alternatively,
388 different change types can be identified based on whether seasonal and trend breaks
389 occur at the same time or not and whether a discontinuity occurs (i.e. magnitude
390 > 0) (Shao and Campbell, 2002). Parameters of the fitted piecewise linear models can
391 also be used to compare long term vegetation trends provided by different satellite
392 sensors. Fensholt et al. (2009), for example, used linear models to analyze trends in
393 annually integrated NDVI time series derived from Advanced Very High Resolution
394 Radiometer (AVHRR), SPOT VEGETATION, and MODIS data. BFAST enables the
395 analysis of long NDVI time series and avoids the need to summarize data annually

13
396 (i.e. loss of information) by accounting for the seasonal and trend variation within
397 time series. This illustrates that further work is needed to extend the method from
398 detecting change to classifying the type of change detected.
399 (3) Evaluating BFAST’s behavior for different change types (e.g. fires versus desertification)
400 in a wide variety of ecosystems remains important. BFAST is tested by combining
401 different magnitudes of an abrupt change with a large range of simulated noise and
402 seasonal variations representing a wide range of land cover types. BFAST is able to
403 detect different change types, however, it remains important to understand how these
404 change types (e.g. woody encroachment) will be detected in ecosystems with drastic
405 seasonal changes (e.g. strong and variable tropical dry seasons) and severe noise in the
406 spectral signal (e.g. sun angle and cloud cover in mountainous regions).
407 (4) The primary challenge of MODIS data, despite its high temporal resolution, is to
408 extract useful information on land cover changes when the processes of interest operate
409 at a scale below the spatial resolution of the sensor (Hayes and Cohen, 2007). Landsat
410 data have been successfully applied to detect changes at a 30m spatial resolution.
411 However, the temporal resolution of Landsat, i.e. 16-day, which is often extended by
412 cloud cover, can be a major obstacle. The fusion of MODIS with Landsat images to
413 combine high spatial and temporal resolutions has helped to improve the mapping of
414 disturbances (Hilker et al., 2009). It is our intention to use BFAST in this integrated
415 manner to analyze time series of multi-sensor satellite images, and to be integrated
416 with data fusion techniques.

417 This research fits within an Australian forest health monitoring framework, where
418 MODIS data is used as a ‘first pass’ filter to identify the regions and timing of major
419 change activity (Stone et al., 2008). These regions would be targeted for more detailed
420 investigation using ground and aerial surveys, and finer spatial and spectral resolution
421 imagery.

422 6. Conclusion

423 We have presented a generic approach for detection and characterization of change in
424 time series. ‘Breaks For Additive Seasonal and Trend’ (BFAST) enables the detection of
425 different types of changes occurring in time series. BFAST integrates the decomposition of
426 time series into trend, seasonal, and remainder components with methods for detecting

14
427 multiple changes in time series. BFAST iteratively estimates the dates and number
428 of changes occurring within seasonal and trend components, and characterizes changes
429 by extracting the magnitude and direction of change. Changes occurring in the trend
430 component indicate gradual and abrupt change, while changes occurring in the seasonal
431 component indicate phenological changes. The approach can be applied to other time
432 series data without the need to select specific land cover types, select a reference period,
433 set a threshold, or define a change trajectory.
434 Simulating time series with varying amounts of seasonality and noise, and by adding
435 abrupt changes at different times and magnitudes, revealed that BFAST is robust against
436 noise, and is not influenced by changes in amplitude of the seasonal component. This
437 confirmed that BFAST can be applied to a large range of time series with varying noise
438 levels and seasonal amplitudes, representing a wide variety of ecosystems. BFAST was
439 applied to 16-day MODIS NDVI image time series (2000–2008) for a forested study area
440 in south eastern Australia. This showed that BFAST is able to detect and characterize
441 changes by estimating time and magnitude of changes occurring in a forested landscape.
442 The algorithm can be extended to label changes with information on the parameters
443 of the fitted piecewise linear models. BFAST can be used to analyze different types of
444 remotely sensed time series (AVHRR, MODIS, Landsat) and can be applied to other
445 disciplines dealing with seasonal or non-seasonal time series, such as hydrology, climatology,
446 and econometrics. The R code (R Development Core Team, 2008) developed in this paper
447 is available by contacting the authors.

448 7. Acknowledgements

449 This work was undertaken within the Cooperative Research Center for Forestry Pro-
450 gram 1.1: Monitoring and Measuring (www.crcforestry.com.au). Thanks to Dr. Achim
451 Zeileis for support with the ‘strucchange’ package in R, to professor Nicholas Coops, Dr.
452 Geoff Laslett, and the four anonymous reviewers whose comments greatly improved this
453 paper.

454 References

455 Anyamba, A., Eastman, J. R., 1996. Interannual variability of NDVI over Africa and
456 its relation to El Nino Southern Oscillation. International Journal of Remote Sensing

15
457 17 (13), 2533–2548.

458 Azzali, S., Menenti, M., 2000. Mapping vegetation-soil-climate complexes in southern
459 Africa using temporal Fourier analysis of NOAA-AVHRR NDVI data. International
460 Journal of Remote Sensing 21 (5), 973–996.

461 Bai, J., Perron, P., 2003. Computation and analysis of multiple structural change models.
462 Journal of Applied Econometrics 18 (1), 1–22.

463 Bontemps, S., Bogaert, P., Titeux, N., Defourny, P., 2008. An object-based change detection
464 method accounting for temporal dependences in time series with medium to coarse spatial
465 resolution. Remote Sensing of Environment 112 (6), 3181–3191.

466 Cleveland, R. B., Cleveland, W. S., McRae, J. E., Terpenning, I., 1990. STL: A seasonal-
467 trend decomposition procedure based on loess. Journal of Official Statistics 6, 3–73.

468 Coppin, P., Jonckheere, I., Nackaerts, K., Muys, B., Lambin, E., 2004. Digital change
469 detection methods in ecosystem monitoring: a review. International Journal of Remote
470 Sensing 25 (9), 1565–1596.

471 Crist, E. P., Cicone, R. C., 1984. A physically-based transformation of thematic mapper
472 data – The TM tasseled cap. IEEE Transactions on Geoscience and Remote Sensing
473 22 (3), 256–263.

474 de Beurs, K. M., Henebry, G. M., 2005. A statistical framework for the analysis of long
475 image time series. International Journal of Remote Sensing 26 (8), 1551–1573.

476 DeFries, R. S., Field, C. B., Fung, I., Collatz, G. J., Bounoua, L., 1999. Combining satellite
477 data and biogeochemical models to estimate global effects of human-induced land cover
478 change on carbon emissions and primary productivity. Global Biogeochemical Cycles
479 13 (3), 803–815.

480 Dixon, R. K., Solomon, A. M., Brown, S., Houghton, R. A., Trexier, M. C., Wisniewski, J.,
481 1994. Carbon pools and flux of global forest ecosystems. Science 263 (5144), 185–190.

482 Fensholt, R., Rasmussen, K., Nielsen, T. T., Mbow, C., 2009. Evaluation of earth ob-
483 servation based long term vegetation trends – intercomparing NDVI time series trend
484 analysis consistency of Sahel from AVHRR GIMMS, Terra MODIS and SPOT VGT
485 data. Remote Sensing of Environment 113 (9), 1886 – 1898.

16
486 Hayes, D. J., Cohen, W. B., 2007. Spatial, spectral and temporal patterns of tropical forest
487 cover change as observed with multiple scales of optical satellite data. Remote Sensing
488 of Environment 106 (1), 1–16.

489 Haywood, J., Randal, J., 2008. Trending seasonal data with multiple structural breaks.
490 NZ visitor arrivals and the minimal effects of 9/11. Research report 08/10, Victoria
491 University of Wellington, New Zealand.
492 URL http://msor.victoria.ac.nz/twiki/pub/Main/ResearchReportSeries/
493 mscs08-10.pdf

494 Healey, S. P., Cohen, W. B., Zhiqiang, Y., Krankina, O. N., 2005. Comparison of Tasseled
495 Cap-based Landsat data structures for use in forest disturbance detection. Remote
496 Sensing of Environment 97 (3), 301–310.

497 Hilker, T., Wulder, M. A., Coops, N. C., Linke, J., McDermid, G., Masek, J. G., Gao, F.,
498 White, J. C., 2009. A new data fusion model for high spatial- and temporal-resolution
499 mapping of forest disturbance based on Landsat and MODIS. Remote Sensing of Envi-
500 ronment 113 (8), 1613–1627.

501 Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., Ferreira, L. G., 2002. Overview
502 of the radiometric and biophysical performance of the MODIS vegetation indices. Remote
503 Sensing of Environment 83 (1-2), 195–213.

504 Jin, S. M., Sader, S. A., 2005. MODIS time-series imagery for forest disturbance detection
505 and quantification of patch size effects. Remote Sensing of Environment 99 (4), 462–470.

506 Jönsson, P., Eklundh, L., 2002. Seasonality extraction by function fitting to time-series
507 of satellite sensor data. IEEE Transactions on Geoscience and Remote Sensing 40 (8),
508 1824–1832.

509 Kennedy, R. E., Cohen, W. B., Schroeder, T. A., 2007. Trajectory-based change detection
510 for automated characterization of forest disturbance dynamics. Remote Sensing of
511 Environment 110 (3), 370–386.

512 Lambin, E. F., Strahler, A. H., 1994. Change-Vector Analysis in multitemporal space - a
513 tool to detect and categorize land-cover change processes using high temporal-resolution
514 satellite data. Remote Sensing of Environment 48 (2), 231–244.

17
515 Lhermitte, S., Verbesselt, J., Verstraeten, W. W., Coppin, P., submitted. Comparison of
516 time series similarity measures for monitoring ecosystem dynamics: a review of methods
517 for time series clustering and change detection. Remote Sensing of Environment.

518 Lu, D., Mausel, P., Brondizio, E., Moran, E., 2004. Change detection techniques. Interna-
519 tional Journal of Remote Sensing 25 (12), 2365–2407.

520 Makridakis, S., Wheelwright, S. C., Hyndman, R. J., 1998. Forecasting: methods and
521 applications, 3rd Edition. John Wiley & Sons, New York.

522 Myneni, R. B., Hall, F. G., Sellers, P. J., Marshak, A. L., 1995. The interpretation of
523 spectral vegetation indexes. IEEE Transactions on Geoscience and Remote Sensing
524 33 (2), 481–486.

525 R Development Core Team, 2008. R: A Language and Environment for Statistical Com-
526 puting. R Foundation for Statistical Computing, Vienna, Austria.
527 URL www.R-project.org

528 Roy, D. P., Borak, J. S., Devadiga, S., Wolfe, R. E., Zheng, M., Descloitres, J., 2002. The
529 MODIS Land product quality assessment approach. Remote Sensing of Environment
530 83 (1-2), 62–76.

531 Shao, Q., Campbell, N. A., 2002. Modelling trends in groundwater levels by segmented
532 regression with constraints. Australian & New Zealand Journal of Statistics 44 (2),
533 129–141.

534 Stone, C., Turner, R., Verbesselt, J., 2008. Integrating plantation health surveillance
535 and wood resource inventory systems using remote sensing. Australian Forestry 71 (3),
536 245–253.

537 Townshend, J. R. G., Justice, C. O., 1988. Selecting the spatial-resolution of satellite
538 sensors required for global monitoring of land transformations. International Journal of
539 Remote Sensing 9 (2), 187–236.

540 Tucker, C. J., 1979. Red and photographic infrared linear combinations for monitoring
541 vegetation. Remote Sensing of Environment 8 (2), 127–150.

542 Venables, W. N., Ripley, B. D., 2002. Modern applied statistics with S, 4th Edition.
543 Springer-Verlag, pp. 156–163.

18
544 Verbesselt, J., Jönsson, P., Lhermitte, S., van Aardt, J., Coppin, P., 2006. Evaluating
545 satellite and climate data derived indices as fire risk indicators in savanna ecosystems.
546 IEEE Transactions on Geoscience and Remote Sensing 44 (6), 1622–1632.

547 Verbesselt, J., Robinson, A., Stone, C., Culvenor, D., in press. Forecasting tree mortality
548 using change metrics derived from MODIS satellite data. Forest Ecology and Management,
549 doi:10.1016/j.foreco.2009.06.011.

550 Zeileis, A., 2005. A unified approach to structural change tests based on ML scores, F
551 statistics, and OLS residuals. Econometric Reviews 24 (4), 445 – 466.

552 Zeileis, A., Kleiber, C., 2005. Validating multiple structural change models – A case study.
553 Journal of Applied Econometrics 20 (5), 685–690.

554 Zeileis, A., Kleiber, C., Krämer, W., Hornik, K., 2003. Testing and dating of structural
555 changes in practice. Computational Statistics and Data Analysis 44, 109–123.

556 Zeileis, A., Leisch, F., Hornik, K., Kleiber, C., 2002. strucchange: An R package for testing
557 for structural change in linear regression models. Journal of Statistical Software 7 (2),
558 1–38.

559 Zhang, X., Friedl, M., Schaaf, C., 2009. Sensitivity of vegetation phenology detection to
560 the temporal resolution of satellite data. International Journal of Remote Sensing 30 (8),
561 2061–2074.

19
562 Figures

563 For interpretation of the references to color in this figure legend, the reader is referred
564 to the web version of this article.

0.9
0.7
NDVI

0.5
0.3
0.90

Real
Simulated
0.80
NDVI

0.70

2000 2002 2004 2006 2008

Time

Figure 1: Real and simulated 16-day NDVI time series of a grassland (top) and pine
plantation (bottom).

20
 0.90
data

0.80
0.70
seasonal

0.02
−0.04
0.82
trend

0.78
remainder

0.05
−0.05
2000 2002 2004 2006 2008

time

Figure 2: The STL decomposition of a 16-day NDVI time series of a pine plantation
into seasonal, trend, and remainder components. The seasonal component is
estimated by taking the mean of all seasonal sub-series (e.g. for a monthly time
series the first sub-series contains the January values). The sum of the seasonal,
trend, and remainder components equals the data series. The solid bars on the
right hand side of the plot show the same data range, to aid comparisons. The
range of the seasonal amplitude is approximately 0.1 NDVI.

21
 0.9
data

0.6
0.3

0.20
seasonal

0.05
−0.10
0.70
trend

0.55
0.40
remainder

0.00
−0.10
2000 2002 2004 2006 2008

Time

Figure 3: Simulated 16-day MODIS NDVI time series with a seasonal amplitude = 0.3,
𝜎 = 0.02 and change magnitude = -0.3. The simulated data series is the sum
of the simulated seasonal, trend and noise series (- - -), and is used as an input
in BFAST. The estimated seasonal, trend and remainder series are shown in
red. Three break points are detected within the estimated trend component (⋅ ⋅ ⋅ ).
The solid bars on the right hand side of the plot show the same data range, to
aid comparisons.

22
 m = −0.1 m = −0.2 m = −0.3 m=0

a = 0.1 a = 0.5
2.0

1.5
RMSE

1.0

0.5

0.0

0.05 0.10 0.15 0.20 0.25 0.05 0.10 0.15 0.20 0.25

Noise

Figure 4: RMSEs for the estimation of number of abrupt changes within a time series, as
shown in Figure 3 (a = amplitude of the seasonal component, m = magnitude
of change). The units of the 𝑥 and 𝑦-axes are 4𝜎 (noise) and the number
of changes (RMSE). See Table 1 for the values of parameters used for the
simulation of the NDVI time series. Similar results were obtained for a = 0.3

23
 m = −0.1 m = −0.2 m = −0.3 m=0
Time Magnitude
a = 0.1 a = 0.1

0.08
10
8

0.06
6
RMSE

0.04
4

0.02
2
0

0.05 0.10 0.15 0.20 0.25 0.05 0.10 0.15 0.20 0.25

Noise

Figure 5: RMSEs for the estimation of the time and magnitude of abrupt changes within
a time series (a = amplitude of the seasonal component, m = magnitude of
changes). The units of the 𝑥-axis are 4𝜎 NDVI, and 𝑦-axis are relative time
steps between images (e.g. 1 equals a 16-day period) (left) and NDVI (right).
See Table 1 for the values of parameters used for the simulation of NDVI time
series. Similar results were obtained for a = 0.3 and 0.5.

24
 148 0’0"E 148 2’0"E 148 4’0"E 148 6’0"E 148 8’0"E
35 35’0"S 35 34’0"S 35 33’0"S 35 32’0"S 35 31’0"S 35 30’0"S 35 29’0"S 35 28’0"S 35 27’0"S 35 26’0"S 35 25’0"S

35 35’0"S 35 34’0"S 35 33’0"S 35 32’0"S 35 31’0"S 35 30’0"S 35 29’0"S 35 28’0"S 35 27’0"S 35 26’0"S 35 25’0"S
Year of planting
2001
2002
2003
2004
2005
2006
Year of change
2001
2002
2003
2004
2005
2006

0 1 2 4 6

Kilometers

148 0’0"E 148 2’0"E 148 4’0"E 148 6’0"E 148 8’0"E

Figure 6: Comparison between the year of Pinus radiata planting derived from spatial
forest inventory data and the BFAST estimate of the year of major change
occurring in MODIS NDVI image time series (2000–2008) for a forested area
in south eastern Australia.

25
 148 0’0"E 148 2’0"E 148 4’0"E 148 6’0"E 148 8’0"E
35 35’0"S 35 34’0"S 35 33’0"S 35 32’0"S 35 31’0"S 35 30’0"S 35 29’0"S 35 28’0"S 35 27’0"S 35 26’0"S 35 25’0"S

35 35’0"S 35 34’0"S 35 33’0"S 35 32’0"S 35 31’0"S 35 30’0"S 35 29’0"S 35 28’0"S 35 27’0"S 35 26’0"S 35 25’0"S
NDVI
-0.48 - -0.37
-0.36 - -0.29
-0.28 - -0.2
-0.19 - -0.1
-0.09 - 0
0.01 - 0.12
0.13 - 0.23

0 1 2 4 6

Kilometers

148 0’0"E 148 2’0"E 148 4’0"E 148 6’0"E 148 8’0"E

Figure 7: BFAST estimated magnitudes of major changes occurring in MODIS NDVI

image time series (2000–2008) for a forested area in south eastern Australia.
Negative values generally indicate harvesting, while positive values indicate
forest growth.

26
 0.9
0.7
NDVI

0.5
0.3

0.9
0.7
NDVI

0.5
0.3
0.90
NDVI

0.80
0.70

2000 2002 2004 2006 2008

Figure 8: Detected changes in the trend component (red) of 16-day NDVI time series
(black) extracted from a single MODIS pixel within a pine plantation, that
was planted in 2001 (top), harvested in 2004 (middle), and with tree mortality
occurring in 2007 (bottom). The time of change (- - -), together with its
confidence intervals (red) are also shown.

27
 0.90
data

0.80
0.70

0.04
seasonal

0.00
−0.04
0.82
trend

0.76
remainder

0.05
−0.05
2000 2002 2004 2006 2008

Time

Figure 9: Fitted seasonal, trend and remainder (i.e. estimated noise) components for
a 16-day MODIS NDVI time series (data series) of a pine plantation in the
northern part of the study area. Three abrupt changes are detected in the trend
component of the time series. Time (- - -), corresponding confidence interval
(red), direction and magnitude of abrupt change and slope of the gradual change
are shown in the estimated trend component. The solid bars on the right hand
side of the plot show the same data range, to aid comparisons.

28