Detecting interannual variation in deciduous broadleaf forest phenology using
Landsat TM/ETM+ data
Eli K. Melaas ⁎, Mark A. Friedl, Zhe Zhu
Department of Earth and Environment, Boston University, 675 Commonwealth Avenue, Boston, MA 02215, United States
a b s t r a c ta r t i c l e i n f o
Article history:
Received 6 July 2012
Received in revised form 31 December 2012
Accepted 12 January 2013
Available online 13 February 2013
Keywords:
Budburst
Landsat
Phenology
Deciduous forest
New England
Observations of vegetation phenology provide valuable information regarding ecosystem responses to climate
variability and change. Phenology is also a first-order control on terrestrial carbon and energy budgets,
and remotely sensed observations of phenology are often used to parameterize seasonal vegetation dynamics
in ecosystem models. Current land surface phenology products are only available at moderate spatial resolution
and possess considerable uncertainty. Higher resolution products that resolve finer spatial detail are
therefore needed. A need also exists for data sets and methods that link ground-based observations of phenology
to moderate resolution land surface phenology products. Data from the Landsat TM and ETM+ sensors
have the potential to meet these needs, but have been largely unexplored by the phenology research
community. In this paper we present a method for characterizing both long-term average and interannual
dynamics in the phenology of temperate deciduous broadleaf forests using multi-decadal time series of
Landsat TM/ETM+ images. Results show that spring and autumn phenological transition dates estimated
from Landsat data agree closely with in-situ measurements of phenology collected at the Harvard Forest in
central Massachusetts, and that Landsat-derived estimates for the start and end of the growing season in
Southern New England varied by as much as 4 weeks over the 30-year record of Landsat images. Application
of this method over larger scales has the potential to provide valuable information related to landscape-scale
patterns and long term dynamics in phenology, and for bridging the gap between in-situ phenological measurements
collected at local scales and land surface phenology metrics derived from moderate spatial resolution
of instruments such as MODIS and AVHRR.
© 2013 Elsevier Inc. All rights reserved.
1. Introduction
In most temperate ecosystems, forest canopy processes related to
leaf development and senescence are strongly controlled by air temperature
(Chuine et al., 2010; Delpierre et al., 2009). As a result,
spatio-temporal patterns in phenology are widely viewed to be important
indicators of climate change (Cleland et al., 2007). Growing
season dynamics are also tightly coupled with biosphere–atmosphere
interactions (Churkina et al., 2005; Richardson et al., 2009), and accurate
representation of phenology is therefore important in land surface
models (Richardson et al., 2012). To support studies of how
ecosystems are responding to climate change and reduce uncertainty
in terrestrial energy and carbon budgets, better datasets characterizing
the response of vegetation phenology to variations in climate
are required (Morisette et al., 2009).
Phenological observations are traditionally collected using two
main approaches: (1) surface observation networks (Schwartz et al.,
2012), and (2) high frequency, coarse spatial resolution satellite remote
sensing (e.g., Jonsson & Eklundh, 2002). Surface observations
provide detailed information related to the timing of leaf development
and flowering phenology for individual plants. However, the
utility of such data for large-scale monitoring and model development
is constrained by data availability, the limited spatial extent of
available samples, and biases inherent to methods used to characterize
the phenology of plants (Cleland et al., 2007).
Remotely sensed data collected by instruments such as the Advanced
Very High Resolution Radiometer (AVHRR) and the Moderate
Resolution Imaging Spectroradiometer (MODIS) provide near-daily
global observations of vegetation dynamics at 250 m to 8 km spatial
resolution. Exploiting this high frequency acquisition strategy, numerous
remote sensing products and algorithms have been developed
over the past two decades that use time series of vegetation
indices such as the normalized difference vegetation index (NDVI)
to track seasonal plant activity (e.g., de Beurs & Henebry, 2010;
Ganguly et al., 2010; Jonsson & Eklundh, 2002; Moulin et al., 1997;
Reed et al., 1994; White et al., 1997). These products and algorithms
have been shown to successfully capture regional-to-global phenological
patterns. However, they are less reliable at local scales and in areas
with heterogeneous land cover where remotely sensed measurements
Remote Sensing of Environment 132 (2013) 176–185
⁎ Corresponding author at: 675 Commonwealth Ave. Room 132, Boston University,
Boston, MA 02215, United States. Tel.: +1 248 343 9278; fax: +1 617 353 8399.
E-mail addresses: emelaas@bu.edu (E.K. Melaas), friedl@bu.edu (M.A. Friedl),
zhuzhe@bu.edu (Z. Zhu).
0034-4257/$ – see front matter © 2013 Elsevier Inc. All rights reserved.
http://dx.doi.org/10.1016/j.rse.2013.01.011
Contents lists available at SciVerse ScienceDirect
Remote Sensing of Environment
journal homepage: www.elsevier.com/locate/rse
reflect mixtures of land cover types, plant species, and plant functional
types (Badeck et al., 2004; Liang et al., 2011).
Bridging the gap between ground observations and moderate resolution
remote sensing approaches, Fisher et al. (2006) used a time
series of 57 Landsat images to generate a map of long-term average
spring phenology across Southeastern New England. At 30 m spatial
resolution, this map yielded substantial information related to spatial
patterns and sources of local variability in leaf phenology that are not
observable from lower spatial resolution instruments such as MODIS.
Landsat spatial resolution also supports analyses at regional scales
(1 scene ~31,000 km2
) and allows factors that control phenology
over relatively short spatial scales to be explored (e.g., topography,
land use and urban heat islands, and coastal effects). However, because
Landsat image acquisitions are relatively infrequent compared
to moderate resolution sensors such as MODIS, Landsat-based characterization
of phenology at annual time scale using conventional
methodologies such as those described by Jonsson and Eklundh
(2002) or Zhang et al. (2003) is challenging. Some efforts have
addressed this using data fusion algorithms that blend MODIS with
Landsat data (e.g., Walker et al., 2012). However, retrieval of phenological
information from such fused data sets is complicated by land
cover heterogeneity below the spatial resolution of MODIS and uncertainty
introduced by the data fusion algorithm (Tan et al., 2006;
Zhu et al., 2010).
In this paper we extend the work of Fisher et al. (2006) to develop
an algorithm that uses time series of Landsat images to characterize
both long-term average and interannual variability in vegetation phenology.
Our analysis was made possible by the opening of the Landsat
archive (Woodcock et al., 2008), which provides new opportunities
for higher spatial resolution analyses of phenology. We exploit this
potential using a 30-year time series of Landsat images to develop
and test a simple algorithm that efficiently and accurately estimates
both long-term average phenology and annual spring and autumn
transition dates in temperate deciduous forests.
2. Materials and methods
2.1. Landsat data
We used all available TM/ETM+ L1T images from 1982 to 2011 with
cloud cover less than 80% (a total of 541 images; Fig. 1) for a Landsat
scene centered over southern New England (Path 12 Row 31; Fig. 2).
The land area in this scene encompasses a mix of land cover types characteristic
of Southern New England, including substantial areas of deciduous
broadleaf forest. The northeast quadrant of this scene contains the
Harvard Forest Long Term Ecological Reserve site (http://harvardforest.
fas.harvard.edu/), where more than 20 years of in-situ vegetation phenology
measurements have been collected and are available for comparison
with time series of Landsat measurements. To reduce atmospheric
effects, DN values from each image were converted to units of surface reflectance
using the Landsat Ecosystem Disturbance Adaptive Processing
System (LEDAPS) atmosphere correction tool (Vermote et al., 1997;
Masek et al., 2008), which uses the MODIS/6S radiative transfer model.
Clouds, cloud shadows, and snow were screened using an algorithm
that has recently been developed at Boston University called Fmask
(Zhu & Woodcock, 2012).
2.2. Phenology algorithm
Most satellite-based algorithms for monitoring phenology use
either the Normalized Difference Vegetation Index (NDVI) or the
Enhanced Vegetation Index (EVI). While both of these indices capture
seasonal changes in vegetation properties, the EVI provides a
larger dynamic range than the NDVI over vegetation with high leaf
area (Huete et al., 2002). Here we used the EVI to develop an algorithm
that detects spatio-temporal patterns in phenology using two main
steps.
The first step in our algorithm estimates a model of long-term average
annual phenology as a function of day of year (DOY) using all cloud-free
1985 1990 1995 2000 2005 2010
50
100
150
200
250
300
350
Year
DOY
Fig. 1. Schematic showing the number and timing of Landsat TM/ETM+ images used in the study. Each point represents a single Landsat image. Note how the density of observations
increases after the launch of Landsat 7.
177E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185
Fig. 2. Study region located in southern New England. Green pixels were classified as deciduous forest based on the method described in Section 2.3, and pixels labeled as “Other”
were classified as non-deciduous forest. Pixels labeled as “No Observations” either fall outside the Landsat scene's spatial extent or have ≤100 cloud-free observations. Water areas
were mapped using to the 2006 National Land Cover Dataset (Xian et al., 2009).
0 50 100 150 200 250 300 350
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
DOY
EVI
Removed Obs.
Other Obs.
Spring Obs.
Autumn Obs.
S(t)
A(t)
Long−term Mean
Transitions
Fig. 3. An example of multi-year phenology for one Landsat pixel located at Harvard Forest, MA. Green and red dots identify observations used to estimate interannual anomalies in
transition dates for spring and autumn, respectively. Spring and autumn anomalies are calculated as the difference between the date of each observation and the date where the
logistic curve (S(t) for spring; A(t) for autumn) reaches the corresponding EVI value. Transition dates for each year are then estimated as the sum of the long-term mean transition
date (yellow dots) and the estimated anomaly. Points identified with x's are outside the date range of the analysis.
178 E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185
EVI values at each pixel (Fig. 3). To do this, we fit separate spring and autumn
logistic functions to EVI time series at each pixel, and use these functions
to estimate long-term average spring onset and autumn offset dates
(ps and pa, respectively). Spring greenup was modeled using a 4 parameter
logistic function (Zhang et al., 2003):
S tð Þ ¼ m1 þ
m2
1 þ e−m3 t−m4ð Þ
ð1Þ
where t is time (expressed as DOY), m1 is the pre-greenup background
EVI, m1+m2 is the maximum EVI, and m3 and m4 control the slope and
phase, respectively. Long-term mean spring onset dates were estimated
as the DOY corresponding to when the first derivative of S(t) is maximum,
which coincides with the most stable and well-constrained portion of
functions estimated using Eq. (1). Uncertainty in logistic function fits is
therefore lower during this period relative to other portions of the curve.
As described by Elmore et al. (2012), a gradual decrease in EVI
(“greendown”) is commonly observed during the summer in many
deciduous broadleaf forests. As a result, Eq. (1) sometimes provides
a sub-optimal representation of autumn dynamics in EVI. To account
for this, we modeled EVI dynamics during the mid-summer and autumn
using the 5 parameter logistic function proposed by Elmore et
al. (2012):
A tð Þ ¼ n1 þ
n2t þ n3
1 þ en4 n5−tð Þ
ð2Þ
where n1 is the leaf-off season background EVI, n2 and n3 control the
trajectory of EVI during the mid-growing season, and n4 and n5 control
the slope and phase, respectively. Following Elmore et al. (2012) autumn
offset dates were estimated using the value of n5 in Eq. (2), which corresponds
to the DOY when the derivative of Eq. (2) is maximum. To optimize
the logistic model fits in both spring and fall, nonlinear regressions
were estimated using the Levenberg–Marquardt algorithm (Levenberg,
1944).
Estimation of coefficients for Eqs. (1) and (2) is sensitive to the
range of dates and density of EVI observations. To determine the
best date for transitioning from S(t) to A(t), we compute the slope
of EVI for moving windows composed of 21 observations starting on
DOY 90. The DOY corresponding to the mid-point for the first window
detected to have a negative slope was identified as the transition
point between the spring and autumn logistic functions. Observations
before DOY 80 and after DOY 340 were excluded during the early
spring and late autumn periods, respectively, and pixels with fewer
than 100 cloud-free observations were excluded from the analysis.
However, because Eqs. (1) and (2) are estimated using a different
number of observations at each pixel, it was important to evaluate
the sensitivity of our approach to variation in the number of EVI observations
at each pixel. To do this, we generated 2500 unique EVI
time series with n ranging from 100 to 288 by randomly sampling
from the pixel with the largest number of cloud free observations in
the scene (n=288). We then estimated spring and autumn logistic
models for each time series and used the results to characterize how
variation in the number of samples affects uncertainty in our results.
In the second step of our algorithm, we use the estimated long-term
mean phenology models at each pixel to quantify interannual variability
in phenology during both spring and autumn. To accomplish this, we examined
all observations acquired within ±20 days of the long-term
mean transition dates whose EVI values met the following criteria for
spring (a and b) and autumn (a and c):
a: EVI≥min[F(t)]+0.2×(max[F(t)]−min[F(t)])
b: EVI≤max[F(t)]+0.2×(max[F(t)]−min[F(t)])
c: EVI≥max[F(t)]−0.4×(max[F(t)]−min[F(t)])
where F(t) is Eq. (1) for spring observations and Eq. (2) for autumn observations.
These criteria are designed to include observations acquired
during the most stable portion of each logistic curve. For each observation
that meets these criteria, we compute the deviation between the
DOY on which the observation was acquired and the DOY corresponding
to the same EVI value for the logistic function(s) fit across all years
(i.e., in Fig. 3, the horizontal deviation between each selected EVI
value and the fitted logistic function). This deviation provides an estimate
of the annual anomaly in the timing of spring onset and autumn
offset dates relative to the long-term average. If more than one observation
from the same year meet the criteria described above, we use the
mean deviation across all qualifying observations. The spring onset
and autumn offset dates for each year are then computed by adding
the estimated annual deviation to the mean transition date; i.e. ps +Δt
and pa+Δt, where ps and pa are the long term average spring onset
and fall offset dates, respectively, and Δt is the annual anomaly in spring
or fall, as defined above.
2.3. Deciduous forest stratification
The method we describe above is designed to characterize both
long-term average and interannual dynamics in the phenology of deciduous
broadleaf forests. To do this, our algorithm exploits the large
amplitude in EVI at each pixel that results from the seasonal emergence
and loss of leaves. Before applying our method, we exploited
this property to distinguish deciduous forest from other land cover
types (e.g., Donohue et al., 2009; Fisher et al., 2006). Specifically,
using the long-term lower (m1) and upper asymptote (m1 +m2) of
EVI derived from the estimated spring logistic function, we classified
pixels as deciduous forest if the lower asymptote of EVI was between
0.1 and 0.25 and the upper asymptote was between 0.6 and 0.9. Both
thresholds were chosen based on comparison of EVI time series with
high-resolution imagery acquired during the leaf-off season (Fig. 4a).
To exclude pixels affected by land cover change, we computed the
95th percentile of EVI values at each pixel identified as deciduous forest
for two periods (1982–2005 and 2006–2011). Pixels for which the
95th percentile in EVI decreased by more than 0.2 between these two
periods were assumed to have experienced some type of land cover
change, and we excluded from further analysis.
2.4. Algorithm assessment
To assess the ability of our algorithm to capture interannual variability
in spring and autumn phenology, we used a 22-year time series of
in-situ observations collected at the Harvard Forest in Petersham, MA
from 1990 to 2011. Tree species composition at Harvard Forest is representative
of temperate forests in southern New England and is dominated
by red oak (Quercus rubra, 36% of basal area), red maple (Acer rubra,
22% of basal area) and yellow birch (Betula allaghaniensis, 14% of basal
area), with smaller quantities of other hardwoods and some conifer, primarily
eastern hemlock (Tsuga canadensis) (Richardson et al., 2009).
Stem density in the forest is approximately 660 stems ha−1
, which corresponds
to roughly 60 stems per Landsat pixel (Moore et al., 1996),
and green leaf area index in deciduous stands reaches a seasonal maximum
of ~5.5 (Urbanski et al., 2007). The in-situ measurements we use
here consist of visual observations of leaf length and coloration collected
during spring and autumn (respectively) for at least three individuals of
eleven tree species at 3–7 day intervals (O'Keefe, 2000).
Spring and autumn transition dates estimated by our Landsatbased
method were compared against ground observations using a
60×60 pixel window (~3.2 km2
) centered over the area on the
ground where in-situ measurements were collected. To perform this
comparison, we used the basal area proportions for each of the three
dominant species identified above to compute site-wide weighted averages
for leaf length in spring and leaf coloration in fall on each ground
observation date. We then compared spring onset and autumn offset
dates estimated from Landsat data with the resulting time series of
leaf length and leaf color measurements at 5% intervals from 5% to
179E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185
95%. This analysis indicated that spring onset and fall offset dates were
most strongly related to the DOY when leaves reached 25% of their maximum
length and when 90% of leaves reached peak coloration, respectively
(Fig. 5). In the results below, these two metrics are a basis for
assessing results from Landsat.
3. Results
Fig. 3 plots Landsat EVI values versus DOY along with estimated
spring and fall logistic fits for a representative deciduous forest pixel
at Harvard Forest. Both spring and autumn logistic curves realistically
capture the average timing of increase and decrease in EVI using the algorithm
described in Section 2.2. The yellow dots in Fig. 3 identify the
long-term mean spring onset (DOY 138=May 18) and autumn offset
dates (DOY 285=Oct 12). The green and red dots identify EVI observations
that meet the criteria defined in Section 2.2 and were used to calculate
annual onset and offset dates at this pixel. Observations located
to the left of the fitted logistic functions indicate that spring greenup
or autumn offset occurred earlier than average, while data points
located to the right of the fitted functions indicate later development
or senescence of leaves in the spring and autumn, respectively.
Fig. 4a presents an aerial photograph of Harvard Forest (obtained
from Bing Maps: www.bing.com/maps) that was acquired during the
leaf-off season when evergreen and deciduous forest areas are distinguishable
from each another. Fig. 4b shows the long-term average amplitude
in EVI estimated from the fitted spring logistic curve (i.e., m2 in
Eq. (1) at each Landsat pixel in a 60×60 window co-located with the
image in Fig. 4a. Red and orange colored areas in this figure identify deciduous
broadleaf stands (larger amplitude), while blue colored areas
have low EVI amplitude and correspond to evergreen needleleaf stands
or other non-deciduous land cover types. Yellow regions have moderate
amplitudes that are characteristic of mixed forests. Based on the
stratification procedure described in Section 2.3, 37% of pixels in
the 60×60 pixel window are deciduous broadleaf forest.
Fig. 6 shows EVI time series for 3 pixels that span the range of
cloud-free observations in the scene, along with a histogram showing
the distribution of cloud free observations between DOY 80 and 340
at each pixel across the entire scene. The median number of cloud
Fig. 4. (a) Aerial photo of Harvard Forest during the leaf-off season. Light brown areas are mainly deciduous forest and green areas are evergreen forests; (b) corresponding annual
EVI amplitude across a 60×60 pixel window (3.24 km2
) covering the same area. Note that pixels with high amplitude (>0.45) correspond to areas with deciduous forest.
125 130 135 140 145 150
125
130
135
140
145
150
Mean 25% Leaf Length (DOY)
LandsatSpringOnset(DOY)
275 280 285 290 295 300 305
275
280
285
290
295
300
305
Mean 90% Leaf Coloring (DOY)
LandsatAutumnOffset(DOY)
a b
R
2
= 0.80
RMSE = 3.0
n = 9
y = 1.16x−21.1
R
2
= 0.80
RMSE = 2.0
n = 7
y = 0.95x+15.3
Fig. 5. (a) Landsat-derived spring onset date versus the mean date when leaf length of common deciduous tree species reached 25% of maximum for one pixel at Harvard Forest;
(b) Landsat-derived autumn offset versus the mean date when leaves of deciduous species reached 90% coloring. Vertical error bars represent 95% confidence intervals for Landsat
estimated transition dates. Horizontal bars show ±1 standard deviation for spring and autumn in-situ phenology metrics. Dashed lines are 1:1.
180 E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185
free EVI values was 240, the maximum was 288, and the vast majority
of pixels had more than 200 cloud free observations available. Fig. 7
presents results showing the sensitivity of logistic function coefficients
to variation in the number of samples at each pixel. These results
show that the slope and phase coefficients for both spring and
autumn logistic functions are relatively insensitive to the number of
EVI time series observations at each pixel, at least for the range of
time series considered here. Specifically, 95% confidence intervals
for the slope coefficient in spring and fall were ±0.029 and ±0.023
(respectively), which translate into uncertainties in spring onset
and autumn offset of ±1.0 days. Similarly, 95% confidence intervals
for the phase coefficient were ±2.0 days for spring and ±2.3 days
for autumn. Uncertainty in model coefficients gradually increases as
the number of observations decreases below 200. However, as we describe
above, the vast majority of pixels in the Landsat scene considered
here have at least 200 cloud-free observations.
Landsat-based estimates of annual spring and autumn transition
dates correspond closely to the timing of in-situ observations at
Harvard Forest. Fig. 8 plots the RMSE and R2
between ground observations
and Landsat-derived estimates of spring phenology at each
deciduous forest pixel in Fig. 4b, which corresponds to the area within
Harvard Forest where the in-situ measurements were collected.
Agreement between ground observations and Landsat retrievals is
uniformly high over the window: average RMSE and R2
values in
spring were 3.4±0.8 days and 0.79±0.10, respectively. Autumn
transition dates estimated from Landsat (not shown) were also highly
correlated with observed dates (RMSE=4.7±1.8; R2
=0.64±0.21).
Spatial variation in Fig. 8a and b is controlled by several factors
including pixel-to-pixel variation in forest phenology and uncertainty
in the species composition, noise inherent to the Landsat EVI data, and
uncertainty in the algorithm itself. Note that depending on the density
of EVI data in the spring and fall of any year (which varies across pixels
because of cloud cover and missing data introduced by the ETM+ scan
line corrector problem), it is not possible to estimate onset and offset
dates in all years. To illustrate, Fig. 9 shows the number of spring and
autumn retrievals (out of a possible 30 years) at each deciduous forest
pixel. On average, the algorithm was able to estimate spring onset and
autumn offset dates in 10 and 9 years, respectively, out of thirty.
4. Discussion
Observations of phenology have been made using remote sensing
for over three decades. Indeed, Multispectral Scanner System data
was used to monitor crop phenology early in the Landsat program
(e.g., Badhwar, 1980). Recently, however, the vast majority of remote
sensing applications related to phenology have used coarse or
moderate spatial resolution data from sensors such as the AVHRR
or MODIS (White et al., 2009; Zhang et al., 2003). Relative to these
instruments, higher spatial resolution sensors such as the Landsat
TM and ETM+ acquire imagery at much less frequent intervals and
have been largely ignored as a source of data for phenological studies.
The availability of long and dense time series of Landsat data, especially
in regions of the world such as North America where the
Landsat archive is densely populated, provides new opportunities
to exploit Landsat data in regional (and potentially continental)
studies of vegetation phenology.
In this paper, we present and test a new method for detecting
long-term average and interannual variation in spring and autumn
phenology of temperate deciduous broadleaf forests using Landsat
data. The algorithm is simple, computationally efficient, and based
100 150 200 250 300 350
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
DOY
EVI
a
n = 101
100 150 200 250 300 350
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
DOY
EVI
b
n = 198
100 150 200 250 300 350
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
DOY
EVI
c
n = 288
100 150 200 250 300
0
2
4
6
8
10
x 10
5
Cloud−free observations DOY 80−340
PixelCount
d
Fig. 6. Multi-year phenological variability for Landsat pixels with the (a) minimum, (b) median, and (c) maximum numbers of cloud-free observations between DOY 80-340 in the
Landsat scene used in this study. The number of EVI observations (n) is displayed in the upper left corner of each plot; (d) histogram showing the number of cloud-free observations
for all deciduous forest pixels across the entire Landsat scene.
181E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185
on comparison of results against field data collected at Harvard Forest,
provides remarkably accurate results. Relative to coarser spatial resolution
data sources such as MODIS and AVHRR, Landsat is able to resolve
much greater fine-scale geographic variability in phenology. Landsat
data therefore has substantial advantages in regions with heterogeneous
land cover where topography, land use, and local microclimates
can lead to substantial variation in phenology over relatively short spatial
scales (Elmore et al., 2012; Fisher & Mustard, 2007). Equally important,
the Landsat TM/ETM+ time series is now over 30 years long, which
makes it an excellent resource for studies of long-term phenological
responses to climate change.
The 16-day revisit period provided by Landsat combined with missing
data arising from cloud cover and the scan line corrector failure on
the Landsat 7 ETM+ significantly reduces the frequency with which annual
spring onset and autumn offset dates can be estimated — for the
scene we considered here, it was possible to estimate these dates in
only about one of every three years. One way to overcome this limitation
is by spatially aggregating retrievals to provide statistical moments
100 150 200 250 300
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
S(t)slope[m3]
100 150 200 250 300
130
132
134
136
138
140
142
S(t)phase[m4]
MC Simulation
Mean
± 95% CI
100 150 200 250 300
−0.5
−0.4
−0.3
−0.2
−0.1
A(t)slope[n4]
100 150 200 250 300
290
295
300
305
310
Cloud−free observations DOY 80−340
A(t)phase[n5]
ba
c d
Fig. 7. Results from uncertainty analysis for slope and phase coefficients of logistic functions. Plots show the range of the slope (a, c) and phase coefficients (b, d) for spring (a, b) and
autumn (c, d) logistic functions with respect to the number of available cloud-free observations between DOY 80-340. Green horizontal lines show the mean values for each coefficient
and horizontal lines show 95% confidence intervals.
Fig. 8. (a) Root mean square error (days) between Landsat-derived spring onset date and the average date when the leaf length of deciduous tree species reached 25% of maximum
across the 60×60 pixel window from Fig. 4; (b) R2
for the same metrics. Note that only pixels classified as deciduous forest were included in the analysis.
182 E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185
of onset and offset dates at coarser spatial resolution (e.g., the mean,
median, and standard deviation in onset and offset dates at 1–5 km).
However, devising the optimal method for doing this requires further
research. Fig. 10, for example, presents boxplots for spring onset and
autumn offset dates for the 60×60 pixel window centered over
Harvard Forest shown in Figs. 8 and 9; even at this scale (~3 km2
),
significant gaps exist in the time series where no retrievals were
made (e.g., 1982–1990 and 1996–1998 in spring; 1994–1998 in
fall). Despite these gaps, it is striking to note that median spring
onset and autumn offset dates in Fig. 10 vary by over 4 weeks.
Consistent with previous studies, spring onset and autumn offset
dates estimated by our method reflect particular phenological metrics
for the start and end of season (e.g., Hufkens et al., 2012a; Schwartz &
Hanes, 2010). Specifically, the Landsat-based retrievals of spring
onset dates at Harvard Forest provided unbiased and highly correlated
estimates for the DOY when leaves reached 25% of their maximum
size. However, depending on the application, different information
related to phenology may be required (e.g., budburst date, which
tends to occur 1–2 weeks earlier at Harvard Forest), and alternative remote
sensing-based metrics may be more appropriate (e.g., Jonsson &
Eklundh, 2002). Similarly, the timing and rate of senescence in autumn
can vary by several weeks both among and within species (John
O'Keefe, pers. comm.), making comparison of in-situ metrics with
satellite-based measurements difficult, even at 30 m spatial resolution.
Thus, more work is needed to develop flexible remote sensing-based
metrics and methods that address the diverse needs of the ecological
community.
Despite the encouraging results we show here, accurate characterization
of uncertainty in remotely sensed measures of phenology remains
an important challenge. For the method we describe in this
paper, uncertainty in retrieved onset and offset dates is controlled
by two main factors: (1) uncertainty in estimated slope and phase coefficients
in Eqs. (1) and (2), and (2) noise in the EVI data. Our analysis
suggests that uncertainty in estimates of the phase and slope
coefficients translates into relatively modest uncertainties for annual
spring onset and autumn offset dates. Our method also assumes that
undetected clouds, cloud shadows, and other sources of noise that
tend to reduce EVI value (thereby biasing estimates for spring onset
and advancing autumn offset dates) are negligible. However, for the
window shown in Figs. 8 and 9, the average rate of change in spring
for the fitted spring logistic functions is 0.013 EVI/day. A decrease in
EVI of 0.1 caused by clouds, shadows, or other sources of noise
would therefore bias the estimated spring onset by 7 days (late).
Our results do not show evidence of any significant bias, but high
quality atmospheric correction and cloud screening is essential to
the success of this method, and future efforts using this algorithm
will need to assess this carefully.
The algorithm we present in this paper assumes that the rate of
change of increase (spring) and decrease (fall) in EVI is consistent
across years. However, recent studies have demonstrated that the
rate of leaf growth can be much higher during anomalously warm
springs than in more normal years (e.g., Barr et al., 2004; Hufkens et
al., 2012b). To assess this, we fit logistic functions to annual time
series of leaf length measurements for individual oak trees in the
Harvard Forest field data set (not shown), and used the results to estimate
uncertainty in the timing associated with 25% leaf length.
Across all trees and all years, the standard deviation of the model
slopes was ±0.08, which corresponds to an uncertainty of ±3.5 days
for the DOY on which leaves achieve 25% their maximum length. This
result supports our assumption that the rate of change in EVI (or leaf
length) is relatively constant across years.
Finally, the assessment we present in this paper is limited to a
temperate deciduous forest site in Southern New England. While
the results we present are encouraging, future efforts are needed to
more fully assess the method over a wider range of geographic and
Fig. 9. Number of successful annual retrievals for (a) spring onset and (b) autumn offset derived from Landsat across the 60×60 pixel window in Fig. 4.
260
270
280
290
300
310
AutumnOffset(DOY)
120
130
140
150
1985 1990 1995 2000 2005 2010
SpringOnset(DOY)
Fig. 10. Boxplots for autumn offset and spring onset dates across all deciduous forest pixels in
the 60×60 Landsat pixel window in Fig. 4 from 1982 to 2011. The horizontal line in the
center of each box is the median, the edges of each box are the 25th and 75th percentiles,
and the whiskers extend to 1.5 times the interquartile range. Red points outside the whiskers
are potential outliers.
183E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185
ecological conditions. Areas located in the regions where adjacent
Landsat scenes overlap are of particular interest, since twice as many
data points are available for those areas. It may be possible to create
gapless time series of spring onset and fall offset dates at these locations,
supporting more extensive analysis of long-term phenological trends
associated with climate variability and change. Evaluation is also needed
across a broader range of ecosystems where phenology is strongly
governed by air temperature (e.g., boreal forests) (Barr et al., 2004;
Melaas et al., 2013), and more generally, in moisture-controlled ecosystems,
where the coupling between climate and phenology is also strong,
but where interannual variation in phenology can be significantly
higher (>1 month; Williams et al., 1997).
5. Conclusions
Phenology is widely viewed to be an important diagnostic of climate
change and is also a first-order control on biosphere–atmosphere interactions.
As a result, remote sensing-based approaches for mapping and
monitoring phenology have been the focus of substantial research over
the past decade. To date, however, most efforts have focused on moderate
spatial resolution data sources derived from sensors such as AVHRR,
SPOT, MERIS and MODIS. These sensors provide frequent observations
at global scale. However, substantial variation in phenology occurs
below the spatial resolution of such instruments and cannot be retrieved
from these data sources.
In this paper we describe and test a method to detect the timing of
spring and autumn phenology for temperate deciduous forests using
Landsat TM/ETM+ data. Our results show that retrievals of spring and
autumn transition dates correspond closely to phenological observations
collected on the ground at Harvard Forest. Relative to coarser spatial
resolution data sources that are commonly used to monitor
phenology, Landsat is able to resolve much finer spatial detail and
therefore captures geographic variability in phenology that is not observable
from coarser spatial resolution sensors. Landsat-derived measurements
of phenology therefore have advantages in regions with
heterogeneous land cover, where the presence of evergreen forests
can significantly complicate estimation of phenological metrics, or
where other factors such as topography, land use, and microclimates
create substantial variation in phenology over short spatial scales.
By providing information related to phenology at much finer spatial
resolution, the 30+ year archive of Landsat data has the potential
to significantly improve understanding of how climate controls phenology
at interannual-to-decadal time scales and at regional-tocontinental
spatial scales. From a purely empirical perspective, longterm
field observations related to phenology are quite rare, especially in
North America. As a result, progress within the community focused on
phenological theory and model development has been limited by available
data. Recent initiatives prioritizing collection of ground or nearsurface
remote sensing based phenology data such as the PhenoCam
(http://klima.sr.unh.edu/) and National Phenology Networks (www.
usanpn.org/) are rapidly expanding available datasets (Betancourt et
al., 2005; Richardson et al., 2009). In the near future, however, these efforts
will not provide long-term observations that are urgently needed
to improve understanding of how climate change affects the phenology
of terrestrial ecosystems. The method we describe in this paper has the
potential to substantially improve this data and knowledge gap. More
effort is required to test this approach over a broader range of conditions,
but at a minimum, the results presented here suggest that it
should be possible to generate high-quality multi-decadal observations
of phenology at 30 m spatial resolution over large areas of deciduous
forest ecosystems.
Acknowledgments
This work was partially supported support by NASA Cooperative
Agreement number NNX08AE61A and by NSF Macrosystems Biology
award number 1064614. The authors gratefully acknowledge data
provided by John O'Keefe at the Harvard Forest.
References
Badeck, F. W., Bondeau, A., Bottcher, K., Doktor, D., Lucht, W., Schaber, J., et al. (2004).
Responses of spring phenology to climate change. New Phytologist, 162(2),
295–309.
Badhwar, G. D. (1980). Crop emergence data determination from spectral data. Photogrammetric
Engineering and Remote Sensing, 46, 369–377.
Barr, A. G., Black, T. A., Hogg, E. H., Kljun, N., Morgenstern, K., & Nesic, Z. (2004).
Inter-annual variability in the leaf area index of a boreal aspen-hazelnut forest in
relation to net ecosystem production. Agricultural and Forest Meteorology, 126,
237–255.
Betancourt, J. L., Schwartz, M. D., Breshears, D. D., Cayan, D. R., Dettinger, M. D., Inouye, D.
W., et al. (2005). Implementing a U.S. National Phenology Network. EOS Transactions
AGU, 86, 539–541.
Chuine, I., Morin, X., & Bugmann, H. (2010). Warming, photoperiods, and tree phenology.
Science, 329(5989), 277–278.
Churkina, G., Schimel, D., Braswell, B. H., & Xiao, X. M. (2005). Spatial analysis of growing
season length control over net ecosystem exchange. Global Change Biology,
11(10), 1777–1787.
Cleland, E. E., Chuine, I., Menzel, A., Mooney, H. A., & Schwartz, M. D. (2007). Shifting
plant phenology in response to global change. Trends in Ecology & Evolution,
22(7), 357–365.
De Beurs, K. M., & Henebry, G. M. (2010). Spatio-temporal statistical methods for
modelling land surface phenology. In Hudson, & Keatley (Eds.), Phenological
research: Methods for environmental and climate change analysis (pp. 177–208).
Delpierre, N., Dufrene, E., Soudani, K., Ulrich, E., Cecchini, S., Boe, J., et al. (2009). Modelling
interannual and spatial variability of leaf senescence for three deciduous tree
species in France. Agricultural and Forest Meteorology, 149(6–7), 938–948.
Donohue, R. J., McVicar, T. R., & Roderick, M. L. (2009). Climate-related trends in
Australian vegetation cover as inferred from satellite observations, 1981–2006.
Global Change Biology, 15, 1025–1039.
Elmore, A. J., Guinn, S. M., Minsley, B. J., & Richardson, A. D. (2012). Landscape controls
on the timing of spring, autumn, and growing season length in mid-Atlantic
forests. Global Change Biology, 18(2), 656–674.
Fisher, J. I., & Mustard, J. F. (2007). Cross-scalar satellite phenology from ground,
Landsat, and MODIS data. Remote Sensing of Environment, 109(3), 261–273.
Fisher, J. I., Mustard, J. F., & Vadeboncoeur, M. A. (2006). Green leaf phenology at
Landsat resolution: Scaling from the field to the satellite. Remote Sensing of Environment,
100(2), 265–279.
Ganguly, S., Friedl, M. A., Tan, B., Zhang, X. Y., & Verma, M. (2010). Land surface phenology
from MODIS: Characterization of the Collection 5 global land cover dynamics
product. Remote Sensing of Environment, 114, 1805–1816.
Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002). Overview
of the radiometric and biophysical performance of the MODIS vegetation indices.
Remote Sensing of Environment, 83(1–2), 195–213.
Hufkens, K., Friedl, M. A., Keenan, T. F., Sonnentag, O., Bailey, A., O'Keefe, J., et al. (2012).
Ecological impacts of a widespread frost event following early spring leaf-out.
Global Change Biology, 18(7), 2365–2367.
Hufkens, K., Friedl, M., Sonnentag, O., Braswell, B. H., Milliman, T., & Richardson, A. D.
(2012). Linking near-surface and satellite remote sensing measurements of deciduous
broadleaf forest phenology. Remote Sensing of Environment, 117, 307–321.
Jonsson, 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.
Levenberg, K. (1944). A method for the solution of certain problems in least squares.
Quarterly of Applied Mathematics, 2, 164–168.
Liang, L. A., Schwartz, M. D., & Fei, S. L. (2011). Validating satellite phenology through
intensive ground observation and landscape scaling in a mixed seasonal forest. Remote
Sensing of Environment, 115, 143–157.
Masek, J. G., Huang, C. Q., Wolfe, R., Cohen, W., Hall, F., Kutler, J., et al. (2008). North
American forest disturbance mapped from a decadal Landsat record. Remote Sensing
of Environment, 112(6), 2914–2926.
Melaas, E. K., Richardson, A. D., Friedl, M. A., Dragoni, D., Gough, C. M., Herbst, M.,
Montagnani, L., & Moors, E. (2013). Using FLUXNET data to improve models of
springtime vegetation activity onset in forest ecosystems. Agricultural and Forest
Meteorology, 171–172, 46–56.
Moore, K. E., Fitzjarrald, D. R., Sakai, R. K., Goulden, M. L., Munger, J. W., & Wofsy, S. C.
(1996). Seasonal variation in radiative and turbulent exchange at a deciduous forest
in central Massachusetts. Journal of Applied Meteorology, 35(1), 122–134.
Morisette, J. T., Richardson, A. D., Knapp, A. K., Fisher, J. I., Graham, E. A., Abatzoglou, J.,
et al. (2009). Tracking the rhythm of the seasons in the face of global change:
phenological research in the 21st century. Frontiers in Ecology and the Environment,
7, 253–260.
Moulin, S., Kergoat, L., Viovy, N., & Dedieu, G. (1997). Global-scale assessment of vegetation
phenology using NOAA/AVHRR satellite measurements. Journal of Climate,
10(6), 1154–1170.
O'Keefe, J. (2000). Phenology of woody species. Harvard Forest Data Archive: HF003.
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.
Richardson, A. D., Anderson, R. S., Arain, M. A., Barr, A. G., Bohrer, G., Chen, G., et al. (2012).
Land surface models need better representation of vegetation phenology: Results
184 E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185
from the North American Carbon Program Site Synthesis. Global Change Biology, 18,
566–584.
Richardson, A. D., Hollinger, D. Y., Dail, D. B., Lee, J. T., Munger, J. W., & O'Keefe, J. (2009).
Influence of spring phenology on seasonal and annual carbon balance in two contrasting
New England forests. Tree Physiology, 29(3), 321–331.
Schwartz, M. D., Betancourt, J. L., & Weltzin, J. F. (2012). From Caprio's lilacs to the USA
National Phenology Network. Frontiers in Ecology and the Environment, 10,
324–327.
Schwartz, M. D., & Hanes, J. M. (2010). Intercomparing multiple measures of the onset
of spring in eastern North America. International Journal of Climatology, 30,
1614–1626.
Tan, B., Woodcock, C. E., Hu, J., Zhang, P., Ozdogan, M., Huang, D., et al. (2006). The impact
of gridding artifacts on the local spatial properties of MODIS data: Implications for
validation, compositing, and band-to-band registration across resolutions. Remote
Sensing of Environment, 105(2), 98–114.
Urbanski, S., Barford, C., Wofsy, S., Kucharik, C., Pyle, E., Budney, J., et al. (2007). Factors controlling
CO2 exchange on timescales from hourly to decadal at Harvard Forest. Journal
of Geophysical Research, 112, G02020. http://dx.doi.org/10.1029/2006/JG000293.
Vermote, E. F., El-Saleous, N., Justice, C. O., Kaufman, Y. J., Privette, J. L., Remer, L., et al.
(1997). Atmospheric correction of visible to middle-infrared EOS-MODIS data over
land surfaces: Background, operational algorithm and validation. Journal of Geophysical
Research-Atmospheres, 102(D14), 17131–17141.
Walker, J. J., de Beurs, K. M., Wynne, R. H., & Gao, F. (2012). Evaluation of Landsat and
MODIS data fusion products for analysis of dryland forest phenology. Remote Sensing
of Environment, 117, 381–393.
White, M. A., de Beurs, K. M., Didan, K., Inouye, D. W., Richardson, A. D., Jensen, O. P.,
et al. (2009). Intercomparison, interpretation, and assessment of spring phenology
in North America estimated from remote sensing for 1982–2006. Global Change
Biology, 15, 2335–2359.
White, M. A., Thornton, P. E., & Running, S. W. (1997). A continental phenology model
for monitoring vegetation responses to interannual climatic variability. Global
Biogeochemical Cycles, 11, 217–234.
Williams, R. J., Myers, B. A., Muller, W. J., Duff, G. A., & Eamus, D. (1997). Leaf phenology of
woody species in a North Australian tropical savanna. Ecology, 78(8), 2542–2558.
Woodcock, C. E., Allen, R., Anderson, M., Belward, A., Bindschadler, R., Cohen, W., et al.
(2008). Free access to Landsat imagery. Science, 320(5879), 1011.
Xian, G., Homer, C., & Fry, J. (2009). Updating the 2001 National Land Cover Database
land cover classification to 2006 by using Landsat imagery change detection
methods. Remote Sensing of Environment, 113, 1133–1147.
Zhang, X. Y., 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(3), 471–475.
Zhu, X. L., Chen, J., Gao, F., Chen, X. H., & Masek, J. G. (2010). An enhanced spatial and
temporal adaptive reflectance fusion model for complex heterogeneous regions.
Remote Sensing of Environment, 114(11), 2610–2623.
Zhu, Z., & Woodcock, C. E. (2012). Object-based cloud and cloud shadow detection in
Landsat imagery. Remote Sensing of Environment, 118, 83–94.
185E.K. Melaas et al. / Remote Sensing of Environment 132 (2013) 176–185