Evaluation of the potential of MODIS satellite data to predict vegetation phenology in
different biomes: An investigation using ground-based NDVI measurements
G. Hmimina a
, E. Dufrêne a
, J.-Y. Pontailler a
, N. Delpierre a
, M. Aubinet b
, B. Caquet d
, A. de Grandcourt d
,
B. Burban e
, C. Flechard f
, A. Granier c
, P. Gross c
, B. Heinesch b
, B. Longdoz c
, C. Moureaux b
, J.-M. Ourcival g
,
S. Rambal g
, L. Saint André h
, K. Soudani a,
⁎
a
University of Paris-Sud, CNRS, AgroParisTech, Laboratoire Ecologie Systematique et Evolution, Faculty of Sciences of Orsay, France
b
University of Liège — Gembloux Agro-Bio Tech (GxABT), Passage des Déportés 2, Gembloux, Belgium
c
INRA, UMR EEF 1137, INRA/University of Nancy, Champenoux, France
d
CIRAD/CRDPI, France
e
INRA, UMR Ecofog, Kourou, Guyane Française, French Guiana
f
INRA, Agrocampus Ouest, UMR 1069 SAS, Rennes, France
g
CNRS, Centre d'Ecologie Fonctionnelle et Evolutive, Montpellier, France
h
INRA, Unité Biogéochimie des Ecosystèmes Forestiers, Champenoux, France
a b s t r a c ta r t i c l e i n f o
Article history:
Received 5 June 2012
Received in revised form 17 December 2012
Accepted 12 January 2013
Available online 11 February 2013
Keywords:
Ground-based NDVI
MODIS
Deciduous forests
Evergreen forests
Crops
Phenology
Vegetation phenology is the study of the timing of seasonal events that are considered to be the result of adaptive
responses to climate variations on short and long time scales. In the field of remote sensing of vegetation phenology,
phenological metrics are derived from time series of optical data. For that purpose, considerable effort has been specifically
focused on developing noise reduction and cloud-contaminated data removal techniques to improve the
quality of remotely-sensed time series. Comparative studies between time series composed of satellite data acquired
under clear and cloudy conditions and from radiometric data obtained with high accuracy from ground-based measurements
constitute a direct and effective way to assess the operational use and limitations of remote sensing for
predicting the main plant phenological events. In the present paper, we sought to explicitly evaluate the potential
use of MODerate resolution Imaging Spectroradiometer (MODIS) remote sensing data for monitoring the seasonal
dynamics of different types of vegetation cover that are representative of the major terrestrial biomes, including
temperate deciduous forests, evergreen forests, African savannah, and crops. After cloud screening and filtering,
we compared the temporal patterns and phenological metrics derived from in situ NDVI time series and from
MODIS daily and 16-composite products. We also evaluated the effects of residual noise and the influence of data
gaps in MODIS NDVI time series on the identification of the most relevant metrics for vegetation phenology monitoring.
The results show that the inflexion points of a model fitted to a MODIS NDVI time series allow accurate
estimates of the onset of greenness in the spring and the onset of yellowing in the autumn in deciduous forests
(RMSE≤one week). Phenological metrics identical to those provided with the MODIS Global Vegetation Phenology
product (MDC12Q2) are less robust to data gaps, and they can be subject to large biases of approximately two weeks
or more during the autumn phenological transitions. In the evergreen forests, in situ NDVI time series describe the
phenology with high fidelity despite small temporal changes in the canopy foliage. However, MODIS is unable to
provide consistent phenological patterns. In crops and savannah, MODIS NDVI time series reproduce the general
temporal patterns of phenology, but significant discrepancies appear between MODIS and ground-based NDVI
time series during very localized periods of time depending on the weather conditions and spatial heterogeneity
within the MODIS pixel. In the rainforest, the temporal pattern exhibited by a MODIS 16-day composite NDVI
time series is more likely due to a pattern of noise in the NDVI data structure according to both rainy and dry seasons
rather than to phenological changes. More investigations are needed, but in all cases, this result leads us to conclude
that MODIS time series in tropical rainforests should be interpreted with great caution.
© 2013 Elsevier Inc. All rights reserved.
1. Introduction
Vegetation phenology is the study of the timing of seasonal events,
such as leaf budburst and leaf senescence, that are considered to be the
result of adaptive responses to climatic constraints. As such, an understanding
of phenology brings important insights into both climate and
Remote Sensing of Environment 132 (2013) 145–158
⁎ Corresponding author. Tel.: +33 1 69 15 64 27.
E-mail address: kamel.soudani@u-psud.fr (K. Soudani).
0034-4257/$ – see front matter © 2013 Elsevier Inc. All rights reserved.
http://dx.doi.org/10.1016/j.rse.2013.01.010
Contents lists available at SciVerse ScienceDirect
Remote Sensing of Environment
journal homepage: www.elsevier.com/locate/rse
vegetation interactions and their impacts on matter and energy exchange
processes at local, regional and global scales. Because field phenological
observations are work intensive and cannot be easily generalized,
remote-sensing tools were developed to track Earth surface changes.
The use of satellite-derived vegetation indices is now frequent in the literature
and has been closely linked to canopy foliage biomass (Soudani et
al., 2006), the onset of leaf greenness in the spring and the onset of leaf
coloring in the autumn (Soudani et al., 2008; Zhang & Goldberg, 2011;
Zhang et al., 2003). Remote sensing-based phenology began with the Advanced
Very High Resolution Radiometer (AVHRR) (Reed et al., 1994) and
has been significantly improved with the Moderate Resolution Imaging
Spectroradiometer (MODIS) onboard Terra and Aqua satellites (Zhang
et al., 2003). Data are acquired daily by AVHRR and MODIS sensors, but
MODIS represents a significant improvement in terms of spatial resolution
(250 m to 1 km vs. 1 km), spectral resolution (36 spectral bands
vs. 6), geolocation accuracy [50 m at nadir (Wolfe et al., 2002) vs. 1 to
2 km (Box et al., 2006)], the atmospheric correction scheme and cloud
screening (Heidinger et al., 2001) and sensor calibration (Justice et al.,
1998). MODIS data are now used routinely for building the MODIS global
vegetation phenology product that provides estimates of the timing of
main vegetation seasonal cycles events at global scales. The first version
of this product (MOD12Q2) was already evaluated, particularly in the
studies of Zhang et al. (2003) and Soudani et al. (2008). Since 2009, a
new version of the global vegetation phenology product (MCD12Q2)
has been available that covers the period from 2001 through 2006.
Compared to the first version, MCD12Q2 uses MODIS with both Aqua
and Terra platforms at higher spatial and temporal resolutions (500 m
vs. 1 km and 8 days vs. 16 days). The first validation studies of this product
are underway (Ganguly et al., 2010).
In the field of remote vegetation phenology sensing, considerable
effort has been focused on developing noise reduction and cloudcontaminated
data removal techniques [e.g., Best Index Slope Extraction
(BISE) (Viovy et al., 1992), a CVA-MVC compositing algorithm
used to produce MODIS-based global vegetation phenology products
(Huete et al., 2002), an adaptive Savitzky–Golay filter (Chen et al.,
2004) and a mean value iteration filter (Ma & Veroustraete, 2006)]. Different
phenological markers were then derived from remotely-sensed
time series data after filtering and noise reduction pre-processing.
These phenological markers may be categorized as follows (Soudani
et al., 2008): (1) user-defined thresholds separating growing and dormancy
seasons (Chen et al., 2004; Delbart et al., 2006; Schwartz et al.,
2002; Studer et al., 2007; Suzuki et al., 2003; White & Nemani, 2006;
White et al., 1997, 2002); (2) markers based on significant and rapid increases
in remotely-sensed signals (Kaduk & Heimann, 1996; Moulin et
al., 1997; Schwartz et al., 2002) and (3) parameters directly determined
from functions fitted to remotely-sensed time series data (Beck et al.,
2006; Fisher et al., 2006; Jönsson & Eklundh, 2002; Soudani et al.,
2008; Zhang et al., 2003). These phenological markers are related to
the vegetation cover types characterized by strong and rapid changes
in leaf density that are sufficient to be detected by remote sensing
sensors. These phenological markers focus on the beginning and end
of the vegetation season, that is, the beginning and end of the period
of canopy photosynthesis, respectively. These events are characteristic
of the phenology of deciduous species. The timing of the beginning of
the photosynthetically active period is associated with the emergence
of buds and the first leaves. The timing of the end of this period is characterized
by depigmentation, leaf yellowing and then leaf fall under the
control of abscission processes. For evergreen species that show less
seasonal change in foliage biomass, the noise inherent to satellitebased
radiance measurements may completely mask the seasonal
variations (Moulin et al., 1997). This interference may explain the fact
that few studies have been devoted to the evergreen vegetation and
that the potential use of remote sensing to monitor the seasonal dynamic
of these biomes has not been sufficiently assessed.
Despite the technological maturity and significant progress achieved
over the last 10 years, there remains a strong need for an effective and
unbiased assessment of the potential and practical use of remotelysensed
data to monitor vegetation phenology. Indeed, the consequences
of applying pre-processing techniques (atmospheric corrections,
noise filtering, and compositing methods) on the performance
of remotely-sensed time series for detecting phenological events have
been evaluated under specific conditions through limited comparisons
of one method against others without referring to field observations
(Chen et al., 2004) or through comparisons with field observations
that are themselves subject to multiple sources of uncertainty (operator
bias, sampling density, temporal frequency, data compilation process,
etc.). However, the multitude of remote sensing-based phenological
metrics used can also make an accurate evaluation of the applicability
of remote sensing for the detection of key vegetation phenological
events much more difficult (White et al., 2009). Finally, from a practical
point of view, phenological metrics provide many estimates that correspond
to different phenological situations, making their practical use in
other studies problematic. Therefore, comparative studies between
time series composed of satellite data in clear and cloudy conditions
and of high-accuracy radiometric data obtained from ground-based
measurements constitute a direct and effective way to assess the operational
use and limitations of remote sensing in predicting the main plant
phenological events. In this study, we sought to explicitly evaluate the
potential use of MODIS remote sensing data for monitoring the seasonal
dynamics of vegetation cover from in situ NDVI measurements in different
vegetation cover types that are representative of the major terrestrial
biomes, including temperate deciduous forests of oak and beech, an evergreen
forest, a tropical rainforest, an African savannah, and a succession
of crops. This assessment relies on tower-based measurements of
NDVI at a half-hourly time step. After cloud screening and filtering, we
will 1) compare temporal patterns and phenological metrics derived
from in situ NDVI time series and from a MODIS daily and 16-day composite
product and 2) evaluate the effects of residual noise and the influence
of data gaps in the MODIS NDVI time series to identify the most
relevant metrics for vegetation phenology monitoring.
2. Materials and methods
2.1. Study sites
This study was undertaken at seven experimental sites that are members
of FLUXNET, the global network of eddy covariance flux towers
measuring carbon, water and energy fluxes between the vegetation
and atmosphere. These study sites cover three main bioclimatic regions
(temperate, Mediterranean, and tropical) and the major plant functional
types encountered: deciduous and evergreen forests, tropical moist evergreen
forest, African savannah, and crops (Table 1). More details about
these sites are provided in Soudani et al. (2012).
Briefly, the temperate deciduous forests are situated in the Fontainebleau,
Hesse, and Fougeres regions in the Northern France. The two main
overstory species are: sessile and pedunculate Oaks [Quercus petraea
(Matt.) Liebl and Quercus robur L.] in the Fontainebleau forest, and
beech (Fagus Sylvatica L.) in the Hesse and Fougeres forests. The Hesse
site is described in more details in Granier et al. (2008).
The evergreen forests are situated in the Puechabon region in Southern
France. The Puechabon site is a holm oak (Quercus Ilex L.) evergreen
broadleaf forest. It is located on the northern Mediterranean coast, and is
representative of the whole region. Holm oak is emblematic of Mediterranean
sclerophyllous vegetation and is encountered in Southern Europe
and the Arab Maghreb region in Northern Africa.
The tropical rainforest site is located in Paracou in French Guiana
(Guyaflux experimental site). It is a mature forest with unknown
human disturbance over the past centuries. This forest is characterized
by a high diversity of plant species. The major species are in
the families of Caesalpiniaceae, Lecythidaceae, Chrysobalanacae, and
Sapotaceae. More details about this experimental site are provided
in Bonal et al. (2008).
146 G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
The African savannah is located in the Tchizalamou study site (North–
North-East of Pointe Noire — Congo). This site is part of the CARBOAFRICA
network. It is composed of grassland dominated by the grass Loudetia
simplex (Nees) Hubb., one of the most common species in this region
of West Africa, with sparse shrubs of Annona arenaria (Schumach. &
Thonn.). This area is burnt every year almost at the same date at the
end of June. More details about this experimental site are provided in
Castaldi et al. (2010).
The succession of crops is located in Lonzee in the Belgian province
of Namur. It is composed of a succession of annual crops over three
years at the same location as wheat (2007), sugar beet (2008), and
wheat and mustard (2009). This site is part of CarboEurope-IP, and
more details about this site and the farming operations are provided
in Aubinet et al. (2009) and Dufranne et al. (2011).
For each site except the Guyaflux site, spatial representativity of
the main species present in the MODIS pixel and “seen” by the in
situ sensor is shown in Table 1. It is calculated as the ratio of basal
areas (or biomass in Tchizalamou site) of species present in the
field of view of in situ sensor to total basal area within MODIS pixel
from field inventories. In forest stands field inventories were done
on the whole MODIS pixel in Hesse and Fontainebleau sites (for
Hesse forest, more details are given in Granier et al. (2008) and for
Fontainebleau forest, data are unpublished but details about the
forest stand are given in Delpierre et al. (2009)).In Fougere forest,
field inventories were done in December 2010 on an area of approximately
10,000 m2
surrounding the flux tower. Both in and outside
the sampled area, over the MODIS pixel, the stand is almost monospecific,
composed of beech (unpublished data). In the Puechabon
forest, inventory was made on 12 plots of 100 m2
each around the
tower. As shown in Table 1, the holm oak occupies about 96% of
total basal area (31 m2
/ha versus 32.3 m2
/ha). This composition is
considered as homogeneous on the 50 ha of the forest management
unit and beyond, because the management over the past century was
exactly the same across the whole forest region as shown in Goerner
et al. (2009). In the tropical rainforest, as underlined above, there is a
high specific diversity, including up to 180 different trees/ha, but no
species is dominant. In the grassland savanna, Loudetia simplex,
followed by Ctenium newtonii Hack. (Poaceae) accounted for more than
75% of the aerial total phytomass regardless of the season. The Lonzee
herbaceous site, as shown in Dufranne et al. (2011), is a homogeneous
agricultural area of about 12 ha.
2.2. In situ NDVI measurements and pre-processing
In situ measurements of NDVI were achieved using a laboratory-made
sensor. More details about this NDVI sensor and the in situ measurement
protocol are provided in Soudani et al. (2012). Briefly, NDVI sensors are
made according to the design described in Pontailler et al. (2003) and
Pontailler and Genty (1996). The body of the sensor is made of Teflon®
installed into a stainless steel cylinder having a diameter of 3 cm and a
height of 9 cm and is equipped with two photodiodes having spectral
sensitivity in red and near-infrared bands. The two photodiodes are
covered with two filters resulting in bandwidths of 640–660 nm and
780–920 nm for red and near infrared, respectively. The technical specifications
of the components of the sensor are provided in Soudani et al.
(2012). The sensor was calibrated against a spectroradiometer (LI-1800,
LI-COR, Inc.). NDVI sensors are installed on towers above the vegetation,
directed downward at a height of several meters above the top of the
canopy. The sensor is inclined at approximately 30° from the vertical
and oriented towards the south to avoid a hotspot effect. The field of
view is 100°, but it is often collimated to account for viewing constraints
encountered at each site. The area viewed is approximately tens to hundreds
of square meters, depending on the site. The data were recorded
at half-hour time steps in a data-storage central unit. The NDVI is computed
for measured radiances reflected by the canopy.
The processing of in situ measurements of NDVI was achieved
according to Soudani et al. (2012). Only daily radiance measurements
acquired under clear sky conditions during the time of the MODIS overpass
are considered in this study.
2.3. MODIS NDVI data and pre-processing
For this study, we used two MODIS products: MODIS NDVI MOD13Q1
V005 16-day composites and MOD09GQK daily surface reflectance at a
250 m resolution. The MOD13Q1 V005 16-day 250 m NDVI for years
2000 to 2009 were obtained through the MODIS Subset gateway for the
pixels centered on the study sites. The MOD13Q1 NDVI values were
built using the Constrained View Angle Maximum Value Composite
(CV-AMVC) algorithm on a 16-day compositing period described in
Huete et al. (2002). The day of the year for each retained NDVI value is
provided. The MODIS/TERRA surface reflectance daily L2G global 250 m
V005 product (MOD09GQK) was obtained through the Earth Observing
Table 1
Main characteristics of the study sites. Percent coverage quantifies the spatial representativity of the species “seen” by the in situ sensor over the MODIS pixel. This representativity
is calculated on the basis of basal areas (or biomass in savanna) of the species present in the MODIS pixel from field inventories.
Site name Type of
biome
Lat/long Altitude
(m)
Average
temperature
(°C)
Average
precipitation
(mm)
Main vegetation species Age
(years)
Maximum
Leaf
Area Index
(m2
/m2
)
Percent
coverage
of main species
Fontainebleau Deciduous
broadleaf
48°28′35″ N
2°46′48″E
120 10.2 720 Sessile and pedunculate Oaks
(Quercus petraea (Matt.) Liebl
and Quercus robur L.)
145 5 80%
Hesse Deciduous
broadleaf
48°40′27″ N
7°03′56″E
300 9.2 820 European beech (Fagus sylvatica L.) 44 5.6 95%
Fougeres Deciduous
broadleaf
48°22′59″ N
1°11′05″ W
140 11.2 900 European beech 40 – 95%
Puechabon Evergreen
broadleaf
forest
43°44′29″ N
3°35′45″ E
270 13.4 907 Holm Oak (Quercus ilex L.) 70 2.9 96%
French
Guiana
Tropical
rain
forest
5°16′54″ N
52°54′44″
W
29 25.7 3136 150 species (DBH>10 cm)/ha – 7 –
Tchizalamou
(Congo)
Herbaceous
Savanna
4°17′210″ S
11°39′23″E
82 25.7 1150 Loudetia simplex – 1.6 75%
Lonzee
(Belgium)
Succession
of
crops
50° 33′ 8″ N
04° 44′ 42″
E
165 10 800 Succession Wheat/Sugar beet/Wheat and
mustard
– – 100%
147G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
System Data Gateway. We used the MODIS Reprojection Tool to extract
the red and near-infrared bands and MODIS per-pixel quality assurance.
For each site, we used the pixel centered on the study site for each
band and year, and we compiled them using MATLAB software (MATLAB
2008a, MathWorks, Natick, Massachusetts, USA). The NDVI was computed
using MOD09GQK MODIS/Terra bands 1 (red: 630–690 nm) and 2
(near infrared: 780–900 nm) data produced at “ideal global quality”.
NDVI ¼
NIR−RED
NIR þ RED
ð1Þ
Where NIR stands for the measured reflectance in the near infrared
band and RED stands for the measured reflectance in the red band.
Unlike the MOD13Q1 products, which were already filtered using
the 16-day CVA-MVC algorithm, daily MOD09GQK NDVI measurements
are highly noisy despite the removal of pixels not flagged as “produced
at ideal quality” in the pixel-level quality assurance (QA) image that
provides QA descriptions about each pixel. Therefore, in this study, we
developed the following scheme to improve the NDVI time series
data. This scheme includes two steps:
- The removal of all values that can be considered unlikely considering
what is known about the annual cycle of vegetation greenness using
a Gaussian mixture model (GMM, McLachlan & Peel., 2000). For all
studied biomes, we considered that the NDVI distribution is bimodal.
The two modes correspond to low and high NDVI values. Low
values coincide with the winter dormancy period in deciduous forests
and periods of bare soils and low vegetation cover in crops
and savannah. In evergreen forests, low NDVI values coincide with
the short periods of leaf renewal that generally occur in the spring.
- The reduction of random noise using a moving-window mean filter
based on the one described in Soudani et al. (2008).
First, to remove unlikely NDVI values contaminated by clouds and
snow, a GMM was fitted on the distribution of NDVI values for each
site and each year.
The two resulting Gaussian densities N(μ1, σ1) and N(μ2, σ2) with
μ1 bμ2 provide three classes of NDVI values defined for each year, as
follows:
- High NDVI values found in the [μ2 −2∗σ2; μ2 +2∗σ2] interval.
- Intermediate values found in the ]μ1 +2*σ1;μ2 −2*σ2[ interval.
- Low NDVI values found in the [μ1 −2∗σ1; μ1 +2∗σ1] interval.
All NDVI values that were not within these bounds were discarded.
Second, we checked the distribution of NDVI acquisition dates in
each of the three NDVI classes defined above. We assumed that the
NDVI time series should exhibit coherent temporal patterns that predominate
over noise and that most noisy observations correspond to
low NDVI values. The few low NDVI observations occurring during periods
that exhibit mostly high NDVI values are hereby considered as
noise and vice versa. Formally, the NDVI filtering process was applied
according to the following rules:
- High NDVI values acquired on the dates that were closer to the modal
date of the low NDVI value class than to the modal date of the high
NDVI value class were discarded and vice versa.
- Intermediate values acquired during the dates that fell in the period
of the high NDVI class were discarded.
- Low NDVI values acquired during the dates that fell in the period of
the intermediate NDVI class were discarded.
After this first step of processing, the retained NDVI values were filtered
and smoothed according to the algorithm presented in Soudani et
al. (2008) using an 11-day moving window and excluding NDVI values
lower than the average value minus the standard deviation of the NDVI
values within the moving window.
2.4. Deriving phenological metrics from NDVI time series
For deciduous canopies, an asymmetric double-sigmoid function
(ADS) was fitted to the NDVI time series independently for each year
using the following equation:
NDVI tð Þ ¼ p à t þ a þ cð Þ þ
1
2
à a−cð Þ Ã tanh b à t−μð Þð Þð Þ−
1
2
à a−eð Þ Ã tanh d à t−vð Þð Þ ð2Þ
tanh is the hyperbolic tangent, t is the time (day of year) and a, b, c, d, e,
u, v, and p are the fitting parameters, where (a+c) is the winter NDVI
value and (a−e) is the amplitude of the NDVI variation. u and v are
the dates corresponding to the highest rates of change of NDVI(t) (maximum
and minimum peaks of the first derivative of NDVI(t)) (Fig. 1).
They are the dates of the two inflexion points when NDVI(t) increases
during leaf expansion (u) and decreases during leaf senescence (v). u
and v correspond to S2 and A2 shown in Fig. 1, respectively.
Note that Eq. (2) is based on the equation of Zhang et al. (2003),
rewritten as in Soudani et al. (2008) but modified by including two
new parameters (e and p). The additional parameter e allows Eq. (2)
to fit two different winter NDVI minima for the start and end of the
year, and the parameter p accounts for the slight linear decrease observed
in the NDVI time series during the winter and summer seasons
(Fig. 1).
For the deciduous species and based on the fitted NDVI time series,
we derived 6 metrics for both spring (S1 to S3) and autumn (A1 to A3)
seasons (Fig. 1). S2 and A2 are the dates of maximum increase and
decrease of the NDVI (inflexion points) during the leaf expansion and
leaf senescence phases, and they are directly given by u and v fitting
parameters, respectively. The metrics S1 and S3 are the days of the
year delimiting the leaf expansion phase in the spring. A1 and A3 are
the days of year delimiting the leaf senescence phase in the autumn.
S1, S3, A1, and A3 are determined from the local extrema of the third
derivative of the fitted NDVI time series (Soudani et al., 2008). Note
that S1, S3 and A1, A3 correspond to the onset of the greenness increase,
the onset of the greenness maximum, the onset of the greenness
decrease and the onset of the greenness minimum, respectively, given
in the MODIS Vegetation phenology product (MCD12Q2). S2 and A2
are not given in the MCD12Q2 product.
For evergreen species whose phenological characteristics exhibit
little change through time, and for crops and savannahs whose phenological
characteristics exhibit irregular temporal patterns that cannot
be modeled by the ADS function (Eq. 2), we used cubic splines to
Day of Year
0 30 60 90 120 150 180 210 240 270 300 330 360
InsituNDVI
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
FirstderivativeoffittedNDVI
-0.8
-0.6
-0.4
-0.2
0.0
0.2
S1 S2 S3 A1 A2 A3
Fig. 1. In situ measured (gray square) and fitted NDVI time series (bold curve) over an
oak forest stand in Fontainebleau during the year 2006. The rate of change given by the
first derivative of the fitted NDVI time series (continuous curve — right axis). The
vertical lines show six phenological markers derived from the fitted NDVI time-series
(S and A refer to spring and autumn phenological events, respectively).
148 G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
fit the NDVI time series because of their great flexibility for fitting
data and because they are differentiable, allowing for the implementation
of the following phenological date detection procedure:
- In a first step, the NDVI transition phases are localized by determining
all local extrema. The duration of each phase is given by
the duration between two successive extrema. Thus, each NDVI
phase corresponds to a period of time where the NDVI varies in
a monotonous manner.
- In a second step, all dates for which the ratio of the noise-to-signal
metric described below (Section 2.5.1) is at a minimum and lower
than the duration of the NDVI phase are considered as having a
phenological significance. In other words, metrics that exhibit a
noise-to-signal ratio higher than the duration of the NDVI transition
phase were considered non-significant.
Fig. 2 illustrates an in situ NDVI time series over an evergreen
forest.
For every site, the dates of phenological events detected using
MODIS daily and MODIS 16-days were compared to those detected
using in situ NDVI time-series.
2.5. Theoretical assessment of the predictive power of vegetation phenology
from in situ and satellite-based NDVI time series
2.5.1. Predictive power of vegetation phenological markers derived from
in situ and satellite-based NDVI time series
Uncertainties on the dates of phenological events derived from NDVI
time series are determined as follows. The NDVI time series may be
written as:
NDVI tð Þ ¼ fNDVI tð Þ þ ε ð3Þ
fNDVI(t) is the signal given by the fitted curve (ADS or cubic spline), and
ε (fitting error) is a random noise with zero mean. This last assumption
is discussed below.
The total variance of NDVI around time t describes the total information
(signal and noise) contained in NDVI(t) and is given by the following
equation (assuming independence between NDVI and ε):
Var NDVIð Þ ¼ Var fNDVI tð Þð Þ þ Var εð Þ ð4Þ
Var(fNDVI(t)) may be locally approximated by:
Var fNDVI tð Þð Þ≈
dfNDVI tð Þ
dt
 2
à Var tð Þ ð5Þ
dfNDVI tð Þ
dt
 2
gives the NDVI variance per unit of time and may be used as a
measure of the information in the NDVI signal observed at time t. Var(t)
is the variance around time (t).
Var(ε) corresponds to the variance of residuals between the predicted
NDVI from the fit fNDVI(t) and the observed NDVI values.
Var(ε) is equal to RMSE(t)2
— i.e., the mean square error.
We define the predictive power of the NDVI time series locally at
time (t) using the following expression:
pp tð Þ ¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
RMSE2
dfNDVI tð Þ
dt
 2
v
u
u
u
t ð6Þ
pp is expressed in days (time unit) and can be read as the number of
days during which a monotonous variation in NDVI(t) should continue
until it exceeds the noise. Therefore, pp is the theoretical uncertainty estimation
of a given date (t) based on the NDVI change from the NDVI fit
and taking into account the noise in the data around the fit. More generally,
it corresponds to the noise-to-signal ratio and expresses the
number of days necessary to obtain an NDVI temporal change higher
than the NDVI noise.
The denominator is the first derivative of the NDVI fit because the
estimation of the phenological metrics is based on the detection of an
NDVI increase or decrease. The numerator term (RMSE) measures the
NDVI noise. The noise is assumed to be constant over the year. For
deciduous forests, this assumption is not strong after modifying the
ADS function according to Eq. (2). Indeed, the examination of the
NDVI time series shows a slight monotonic decrease in the NDVI
throughout the seasons of winter and summer. The two parameters p
and e were introduced to homogenize the distribution of residuals
around the fit. However, over evergreen forests, the distribution of the
residuals around the fit is not always homogeneous, as the flexibility
of the spline that is used is not sufficient to account for all of the NDVI
signal variations, which are particularly fast and short.
Eq. (6) was used to assess the potential error around an estimated
phenological metric over both deciduous and evergreen forests.
2.5.2. Sensitivity analysis of MODIS-derived phenological metrics to data
gaps in NDVI time series
To estimate the influence of data gaps due to clouds and snow in
remotely-sensed NDVI time series on the ability to predict phenological
events, we artificially introduced data gaps with several different
lengths in in situ NDVI time series.
In a first step, the probability density function (pdf) of clear sky in
the year at the time of the MODIS overpass was established using
measurements of the sunshine duration acquired using a BF3 sunshine
sensor (Delta-T devices) at a half-hour time step in the Fontainebleau
tower-flux site during the year 2008. The year 2008 may
be considered as a “normal year” representative of the regional
cloud cover regime. We defined a clear sky as that in which the sunshine
duration is at least 10 min per 30 min of measurements using
the BF3 sensor. The pdf of clear sky days was determined using a
non-parametric kernel smoothing density estimation fitted to the
empirical histogram of clear days (Fig. 3). The pdf determines the
conditional probability that the sky will be clear at the time of the
MODIS overpass on a given day of the year. For example, if it is assumed
that half of the year is covered by clouds during the MODIS
overpass and that all days are equally likely to be cloudy, then the
probability of having a clear sky for a given day is: (1/365)×0.5.
In a second step, an in situ NDVI time series is used to randomly
generate an NDVI time series with gaps. We generate an NDVI time
series with between 30 and 189 observations according to the probability
density function for clear sky determined above. For each length
of the NDVI time series, this operation is repeated 100 times. As a result,
a total of 16,000 NDVI time series were created. Then, the ADS
Day of year
30 60 90 120 150 180 210 240 270 300 330 360
InsituNDVI
0.0
0.2
0.4
0.6
0.8
1.0
FirstderivativeoffittedNDVItime-series
-0.010
-0.005
0.000
0.005
0.010
Fig. 2. In situ measured (gray square) and fitted NDVI time series using cubic splines
(continuous curve) over an evergreen broadleaf forest in Puechabon during the year
2009. First derivative of the fitted NDVI time series (continuous curve — right axis).
149G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
(Eq. 2) was fitted to each sample of the simulated NDVI time series,
and the six phenological metrics defined above were estimated.
To quantify the strength of the NDVI signal of each simulated NDVI
time series, we calculate the first derivative of the fitted NDVI at the
date of each NDVI observation of the dataset. The absolute values of
the first derivatives averaged over the series are used as a simple
measure of the NDVI signal strength. For example, for a time series
composed of very close NDVI values, it is expected that the curve
fitted to the NDVI time series will be flat, the first derivative calculated
for each NDVI observation will be zero or close to zero and the average
signal strength over all NDVI observations will also be zero or
close to zero.
Finally, according to the strength of the NDVI signal, we grouped
all samples into 100 classes using an unweighted pair group method
with the arithmetic mean (Sokal & Michener, 1958). For each class,
we determined the average length of the simulated NDVI time series
within the class and the average RMSE between the estimates of the
phenological marker determined from the simulated time series and
full NDVI time series.
3. Results
3.1. Comparison between ground- and MODIS-based NDVI time series
Fig. 4 shows the overall comparison between the in situ, daily and
16-day MODIS NDVI data.
For all vegetation types, the MODIS NDVI values are generally
higher than those measured by in situ sensors. For the deciduous forests,
the overall agreement between in situ and MODIS NDVI data is
very good (R2
=0.91, pb0.00). In the holm oak forest at Puechabon,
the relationship is highly scattered due to small temporal changes
in the canopy foliage area compared to the magnitude of the noise affecting
the NDVI signal, even after filtering (R2
=0.13, pb0.00). In the
tropical forest, no significant relationship could be found (R2
=0,
p>0.75). In the succession of crops in Lonzee (R2
=0.56, pb0.00)
and in the savannah site (R2
=0.45, pb0.00), the relationships between
in situ and MODIS NDVI measurements are significant but
noisy despite the high in situ NDVI temporal changes.
In the deciduous forests, the linear relationships between the in
situ and daily MODIS NDVI and between the in situ NDVI and
MODIS 16-day composite data are similar for all sites (Fig. 4). While
both filtering methods – CVA-MVC and GMM – appear to perform
equally well, the filtering method based on a GMM applied to the daily
MODIS NDVI time series provides a better temporal resolution and, as
expected, preserves the intermediate values, which generally correspond
to key phenological phases (leaf expansion and leaf senescence), as
shown in Fig. 4.
While the filtered MODIS NDVI observations are linearly correlated
to the in situ NDVI measurements for deciduous forests, the evergreen
forests exhibit a high level of noise in regard to the low in situ NDVI
variability. The noise levels as estimated using RMSE between the
observed and fitted NDVI time series in the evergreen forest in the
Puechabon and deciduous forests are similar (0.03 versus 0.04), and
the lack of correlation is mainly due to a low NDVI temporal variability
(low seasonal variation of phenology). In the tropical forest of French
Guyana, no significant relationship could be found between the
MODIS and in situ NDVI observations (p>0.75). The RMSE is larger
than in the deciduous forests (0.12 versus 0.04), possibly due to the effects
of sub-pixel cloud contamination not detected by a MODIS cloud
mask algorithm and a failure in the filtering process.
Despite the noise in the data, the in situ and MODIS data reproduce
similar temporal patterns in deciduous forests (Fig. 5). The NDVI time
series show with high temporal resolution the phenological seasonality
of these species, which is characterized by two phases: the growing season
from mid-spring to summer and the dormancy season during late
autumn and winter. These main phases are separated by two short transition
phases delimited by two main phenological events: the onset of
greenness when budburst starts in the spring and the onset of senescence
in the autumn.
For the evergreen broadleaf forest of Puechabon (Fig. 6), the in situ
NDVI time series show clear phases of NDVI decrease during the spring
despite small NDVI variations. This pattern is consistent with the phenology
of the holm oak characterized by the partial foliage renewal
each year from March to the middle of June. Note that in the Puechabon
forest, an unexplained sensor dysfunction coinciding with strong rains
explains the gap in NDVI measurements during the autumn of the
year 2008. The MODIS daily and MODIS 16-day composite NDVI time
series show small signal variations that do not coincide with the in
situ NDVI measurements.
In the tropical forest (Fig. 7), in situ NDVI time series show two
periods characterized by declines in the NDVI of variable magnitudes
occurring around the middle of March for the first period and around
days 300–320 (October) for the second period. For the first period, the
decline in NDVI is clearly visible in 2007 and 2008. For the second
period, the decline in NDVI appears only in 2008 and 2009. The first
period of NDVI decline was much shorter than the second one. The
second decline was more pronounced in 2009. Contrary to the in situ
NDVI measurements, MODIS NDVI time series (Fig. 7) include so
much noise that none of the tested filtering methods could provide a usable
signal. The temporal pattern from the MODIS 16-day composite
NDVI data is inconsistent with the in situ NDVI time series. This pattern
is mainly characterized by two periods that coincide with the main
rainy season from December to July (sometimes interrupted by a
short dry season in March called the little summer of March) and the
main dry season (July–November). The rainy season is characterized
by abnormally low values of NDVI, and the second period is characterized
by NDVI values at the same level as the daily MODIS NDVI data.
This temporal pattern may arise from variations in noise intensity. As
shown in Fig. 7, none of the tested filters could provide a good agreement
between the in situ and MODIS observations.
During the succession of crops at the Lonzee site (Fig. 8), the in
situ NDVI measurements started in 2007 at the end of March, during
the growth of winter wheat. For this crop, the in situ NDVI increases
during the spring, reaches a peak at the end of April and then decreases
during June and July to reach a minimum value a few weeks
before the harvest at the beginning of August. After the harvest, the
NDVI peaks again in the first week of September due to a re-growth
of wheat and weeds. In 2008, during the growth of the sugar beet
crop, the NDVI increases, reaches its maximum at the end of June
and remains almost constant during the summer until the harvest
at the beginning of November. In 2009, the in situ NDVI time series
Day of Year (2008 - Fontainebleau)
0 30 60 90 120 150 180 210 240 270 300 330 360
ClearSkydaysProbabilityDensityFunction
0.0005
0.0010
0.0015
0.0020
0.0025
0.0030
0.0035
0.0040
InsitumeasuredNDVI
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Fig. 3. In situ measured NDVI in Fontainebleau forest for the year 2008 (empty square,
left axis). The filled area under the continuous curve is the probability density of clear
sky for each day of the year (the sky is considered clear when the sunshine duration is
at least 10 min per period of 30 min at the time of the MODIS overpass (right axis)).
150 G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
is bimodal, reproducing the phenology of a succession of two crops of
winter wheat and mustard.
In the succession of crops, the MODIS daily and 16-day composite
signals (Fig. 8) exhibit strong noise, but the main temporal patterns
associated with the phenology of this vegetation could be identified.
In the grass savannah at the Tchizalamou site (Fig. 8), the temporal
patterns of the in situ NDVI measurements in 2008 and 2009 are
similar, with the exception that the NDVI remains high during February
and March 2009. The NDVI is at its maximum during the main wet
season from October to May and at its minimum during the main dry
0.4 0.5 0.6 0.7 0.8 0.9
0.4
0.5
0.6
0.7
0.8
0.9
Deciduous forest
R²=0.91
MODISdailyNDVI
0.4 0.5 0.6 0.7 0.8
0.4
0.5
0.6
0.7
0.8
0.9
R²=0.92
MODIS16daysNDVI
Ground-based NDVI
0.65 0.7 0.75 0.8 0.85
0.65
0.7
0.75
0.8
0.85
Evergreen forest (Puechabon)
R²=0.13
0.6 0.65 0.7 0.75 0.8
0.6
0.65
0.7
0.75
0.8
R²=0.14
Ground-based NDVI
0 0.2 0.4 0.6 0.8
0
0.2
0.4
0.6
0.8
Crops
R²=0.56
0 0.2 0.4 0.6 0.8
0
0.2
0.4
0.6
0.8
R²=0.55
Ground-based NDVI
-0.2 0 0.2 0.4 0.6
-0.2
0
0.2
0.4
0.6
Savannah
R²=0.45
-0.2 0 0.2 0.4 0.6
-0.2
0
0.2
0.4
0.6
R²=0.49
Ground-based NDVI
0.5 0.6 0.7 0.8 0.9
0.5
0.6
0.7
0.8
0.9
Tropical forest
R²=0
0.4 0.5 0.6 0.7 0.8 0.9
0.4
0.5
0.6
0.7
0.8
0.9
R²=0
Ground-based NDVI
Fig. 4. In situ NDVI (x-axis) versus daily (upper figures) and 16-day composite MODIS NDVI data for the different biomes. R2
: coefficient of determination.
Fig. 5. Time series of in situ NDVI (gray squares), daily MODIS NDVI (empty circles) and 16-day composite MODIS NDVI (filled circles) over deciduous forests.
151G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
season from May to October. The first decrease of the NDVI in 2008 is
due to a short dry season, which may occur in February–March. For
the two years, the sudden drop of NDVI at the end of the wet season
is due to human-induced fire. Nevertheless, it is important to note
that the mismatch between the in situ and MODIS observations in
the savannah site occur during both rainy and dry seasons, and it is
most likely due to contamination of the data by clouds that are not
detected by the filtering process.
3.2. Comparison of phenological metrics estimates derived from in situ
and MODIS daily NDVI time series
For the deciduous species, the comparative analysis for the 6 considered
phenological metrics are shown in Fig. 9 and Table 2.
The best agreement between the predictions of the phenological
dates based on the MODIS time series and in situ NDVI measurements
is found for the two inflexion points S2 and A2 during the leaf expansion
and the leaf senescence phases, respectively. The bias between the
MODIS predictions and those based on in situ NDVI measurements is
positive for S2, S3, A2, and A3 (MODIS-based phenological markers
occur later). It is less important for S2 and A2 than for S3 and A3, which
delimit the end of the two phases, i.e., the end of the leaf expansion
phase in the spring (S3) and the end of the leaf senescence phase in
the autumn (A3). In comparison with the daily MODIS series, the bias
for S2 is positive at approximately 4 days, and it is also positive for
A2 at 2.5 days. On both sides of the two inflexion points, the bias is
negative for S1 in early spring (MODIS-based phenology estimates are
earlier) and positive for S3 at the end of the leaf expansion phase.
Furthermore, it is negative for A1 in early autumn and positive for A3
at the end of the leaf senescence phase. The relationships between the
ground-based and MODIS daily spring metrics (S1, S2, S3) are statistically
significant (pb0.01), while no statistically significant relationship
could be found for the autumn metrics (A1, A2, A3). The metrics derived
from the MODIS 16-day series exhibit lower R2
, and only the S2 and S3
metrics could be significantly related to the ground-based metrics
(pb0.02 and pb0.04, respectively).
3.3. Theoretical analysis of the predictive power of NDVI time series for
phenology detection
Because the phenological metrics are based on the derivatives of the
fitted NDVI time series, the RMSE of the first derivative ratio was used to
quantify the expected theoretical uncertainty of any particular phenological
metric regardless of the signal noise (Eq. 6). As described above, we
Fig. 6. Time series of in situ NDVI (gray squares), daily MODIS NDVI (empty circles) and 16-day composite MODIS NDVI (filled circles) over an evergreen broadleaf forest in
Puechabon. The continuous curves represent the series of cubic splines fitted to the NDVI data.
Fig. 7. Time series of in situ NDVI (gray squares), daily MODIS NDVI (empty circles) and 16-day composite MODIS NDVI (filled circles) over the tropical forest in French Guyana.
Continuous curves represent the series of cubic splines fitted to the NDVI data.
152 G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
recall that this ratio corresponds to the number of days needed to obtain
an NDVI temporal change higher than the NDVI noise. Fig. 10A and B
illustrates the application of this method to track features and to estimate
the uncertainties (in days) of main phenological metrics in a
deciduous forest in Fontainebleau and in an evergreen forest in
Puechabon (Fig. 10).
Summary statistics of uncertainty in the main phenological metrics
over all deciduous forests and for all years are provided in Table 3.
The uncertainties of phenological dates determined from the MODIS
daily and 16-day NDVI products are similar. Nevertheless, the standard
deviation is larger using the MODIS 16-day NDVI time series. In comparison
with the in situ NDVI-derived phenological metrics, the theoretical
uncertainty based on the daily and MODIS 16-day composite NDVI
observations is higher, particularly for phenological markers associated
with autumn. The phenological dates given by the inflexion points during
the leaf expansion phase in the spring and the leaf senescence phase
in the autumn are significantly more accurate than the other metrics.
This result is of great importance because it demonstrates that the inflexion
point metric is more robust for tracking the phenology from
the NDVI time series in temperate broadleaf deciduous forests.
For the other biomes, summary statistics of uncertainty assessed
over evergreen forests, savannah and crops are provided in Table 4.
For crops and the herbaceous savannah, the uncertainties obtained
for the most significant transition dates derived from the in situ NDVI
time series are comparable to those obtained for deciduous forests. For
the evergreen forest of holm oak and the rainforest, the uncertainties
are higher, differing from in situ NDVI measurements by approximately
10 days. In contrast, the daily and 16-day MODIS NDVI time series are
not able to describe the phenology of these ecosystems with sufficient
accuracy.
3.4. Influence of data gaps in the MODIS NDVI time series on the prediction
accuracy of phenological metrics in deciduous forests
Based on the Fontainebleau 2008 in situ NDVI time series considered
as a reference and by inserting artificial gaps into the actual data using
the method described in Section 2.5.2, Fig. 11 shows the differences between
the phenological estimates from the full NDVI time series and the
simulated (with data gaps) NDVI time series for the different phenological
markers defined in Fig. 1.
In Fig. 11, the abscissa (x-axis) corresponds to the mean of the absolute
values of the first derivatives of ADS fitted to the simulated NDVI
time series. The first derivative is numerically calculated for every day
of NDVI observation. Low values of the x-axis represent a loss of information
(loss of NDVI signal) due to a bad compositing of the NDVI
time series (decrease of the proportion of informative observations
that are acquired during leaf expansion and leaf senescence phases),
while the increase represents a relative gain of information due to the
removal of uninformative observations (i.e., NDVI values during the
winter and summer seasons). The two ordinate axes (left y-axis and
right y-axis) correspond to the average RMSE between the estimates
of the phenological marker determined from the simulated and full
NDVI time series (left y-axis) and the average length of the simulated
NDVI time series (right y-axis).
Fig. 8. Time series of in situ NDVI (gray squares), daily MODIS NDVI (empty circles) and 16-day composite MODIS NDVI (filled circles) over herbaceous species in Lonzee (succession
of crops) and at the Tchizalamou site (African savannah). Continuous curves represent the series of cubic splines fitted to the NDVI data.
Predicted dates from in situ NDVI
70 80 90 100 110 120
PredicteddatesfromMODISdata
70
80
90
100
110
120
MODIS daily
MODIS 16-days
100 105 110 115 120 125 130
100
105
110
115
120
125
130
110 120 130 140 150 160
110
120
130
140
150
160
Predicted dates from in situ NDVIPredicted dates from in situ NDVI
Spring minimum Spring maximumSpring inflection point
Fig. 9. In situ NDVI-derived metrics (day of year) versus MODIS NDVI-derived metrics (day of year) for deciduous forests.
153G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
The general form of the relationships between the RMSE and the
mean of the absolute values of the first derivatives of the fitted NDVI
time series is concave up. The RMSE is minimal around the middle of
the x-axis and then increases rapidly on both sides for low and high
values. In contrast, the general form of the relationships between the
average length of simulated NDVI time series and the x-axis is concave
down, indicating that the prediction errors of phenological metrics are
lowest when the length of the simulated NDVI time series is high.
In Fig. 11, both the level and the extent of the central part of the
RMSE curve, characterized by stable and low values, vary according to
the phenological metric considered. The RMSE is lower and the region
of error stability is wider for phenological metrics based on the inflexion
points (S2 and A2) during the spring and autumn, respectively. For the
spring phenological metrics, the point of inflexion S2 is subject to
the smallest error (b1 day) for a length of NDVI time series varying
between 35 and 140 observations. For the autumn phenological metrics,
the prediction error at the point of inflexion (A2) is higher and
less stable at approximately 3 days for a range of sample sizes varying
from 60 to 140 observations.
The extent of the portion of the curve where the RMSE is stable
indicates the sensitivity of the phenological metrics to the length of
the NDVI time series and to signal degradation due to cloud cover.
Table 5 gives the range of the length of the simulated NDVI time
series that bounds the region of RMSE stability (defined as the region
having the lowest RMSE±1 day). This range gives the minimum number
of NDVI observations retained without significantly degrading the
predictive quality of the NDVI time series.
The results in Table 5 show that the average length of the simulated
NDVI time series is considerably lower for S2 and A2 than for the
other metrics, indicating that the phenological metrics based on the
inflexion points are less sensitive to data gaps and signal degradation
due to cloud cover. The prediction error remains at its minimum by
keeping only approximately 27% and 40% of the total number of NDVI
observations in the best case (gaps or cloudy observations concentrated
during the summer and winter) or approximately 40% and 61% of the
total number of observations in the worst case (gaps or cloudy observations
concentrated during the spring or autumn).
4. Discussion
In situ NDVI measurements are made only a few meters above the
canopy, and because NDVI is a normalized index, the effects of the sky
conditions produce little noise. In situ NDVI measurements can thus be
carried out under diffuse sky conditions, allowing for the monitoring of
vegetation phenology at high temporal frequency in deciduous and evergreen
forests for which the phenological variations are less pronounced.
These data may be considered as a reference offering adequate empirical
and theoretical frameworks for directly assessing the potential use of
satellite data to predict vegetation phenology in different biomes and
under different sky conditions. Nevertheless, when comparing coarse
satellite data to ground measurements, spatial heterogeneity can become
an important source of uncertainty in predicting phenology if certain
precautions are not taken. In this study, the ground-based NDVI
measurements benefited from an existing network of seven eddy
covariance flux towers that were installed on flat terrain with relatively
homogeneous vegetation cover, specifically chosen to satisfy the assumptions
of the eddy covariance method and to avoid scaling issues
and plant species heterogeneity (Chen et al., 2009; Metzger et al.,
2012). In addition, for each study site, the homogeneity of the vegetation
composition within the MODIS 250 m pixel was checked by visual photointerpretation
complemented by the use of Landsat TM/ETM+ based
NDVI subsets and ancillary ground-based observations. Nevertheless, it
is important to note that these precautions do not completely remove
the residual uncertainty due to the spatial heterogeneity in the MODIS
pixel. We emphasize that ground-based NDVI measurements are acquired
at a constant viewing angle, while MODIS data are acquired
with different viewing geometries. BRDF effects lead to uncertainties
of variable magnitude in the seasonal course of surface reflectance,
Table 2
Comparison between the phenological metrics derived from the in situ NDVI measurements (considered as reference) and MODIS data. (S1, S2, S3) for the spring and (A1, A2, A3) for
the autumn. MAE (days): mean absolute error, bias (days): (+) MODIS overestimation, (−) underestimation, RMSE (days): root mean square error. R2
: coefficient of determination
of the regression between the phenological markers based on the in situ and MODIS NDVI time series. Significant values are noted in bold (pb0.05).
S1 S2 S3 A1 A2 A3
MODIS
daily
MODIS
16-day
MODIS
daily
MODIS
16-day
MODIS
daily
MODIS
16-day
MODIS
daily
MODIS
16-day
MODIS
daily
MODIS
16-day
MODIS
daily
MODIS
16-day
MAE 2.5 9 4 3 10 7 11 11 4 5.5 15 12
Bias −2 −3.5 4 1.5 10 7 −4 7.5 2.5 4 9 0.5
RMSE 3. 11 4.5 4 10 10.5 14.5 14 6 8 22 14
R2
0.8 0.2 0.9 0.7 0.97 0.59 0.42 0.76 0.56 0.47 0.04 0.01
p 0.007 0.312 0.001 0.019 0.000 0.043 0.115 0.102 0.052 0.09 0.668 0.825
Day of Year
0 30 60 90 120 150 180 210 240 270 300 330 360
InsituNDVI
0.0
0.2
0.4
0.6
0.8
1.0
RMSEtofirstderivative(days)
0
2
4
6
8
10
12
14
16
Deciduous forest
(Fontainebleau)
Day of year
30 60 90 120 150 180 210 240 270 300 330 360
InsituNDVI
0.0
0.2
0.4
0.6
0.8
1.0
RMSEtofirstderivative(days)
0
2
4
6
8
10
12
14
16
Evergreen forest
(Puechabon)
Fig. 10. In situ measured (gray squares) and fitted NDVI time series over a deciduous forest in Fontainebleau for year 2006 (left) and over an evergreen broadleaf forest of holm oak
in Puechabon forest (right) during the year 2009. The continuous curve is the ratio between the RMSE and first derivative (right axis).
154 G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
and they may cause bias in the identification of vegetation phenological
events (Fensholt et al., 2010; Hird & McDermid, 2009; Sims et al.,
2011; Tan et al., 2006). In this study, we used the MODIS 16-day
product derived from the CVA-MVC compositing methodology that
preferentially select highest NDVI values with zenith view angle
closest to nadir view (Huete et al., 2002). No specific constraint on
the viewing angle has been applied to daily MODIS data previously.
Nevertheless, daily 250 m MODIS NDVI data were previously filtered
using GMM and a moving-window mean filter to minimize the total
noise due to variations in the atmospheric conditions, mismatch in
the spatial scales between the MODIS data sets, and differences in
radiometric data acquisition geometry.
In this study, a new filtering method based on mixture Gaussian
models has been developed to remove spurious MODIS NDVI data in
deciduous forests without altering their phenological patterns. This
method showed good performance in terms of the similarity between
the in situ and daily MODIS NDVI time series (Fig. 5). However, in evergreen
forests, this method has a limited efficiency for filtering daily and
16-day MODIS NDVI time series because the magnitude of noise is of the
same order as the phenological signal (Figs. 6 & 7).
After the removal of spurious NDVI observations, both MODIS and in
situ NDVI time series allow us to predict with good accuracy the two
main phenological events in temperate deciduous forests: the date of
the onset of greenness in the spring and the date of the onset of leaf
senescence in the autumn (Fig. 9, Table 2). The MODIS daily NDVI
derived onset of greenness metrics was shown to be well correlated
to the in situ NDVI metrics, while the related senescence metrics may
still be challenging. The use of the MODIS 16-day NDVI series yielded
more variable results, which is probably due to the loss of intermediate
NDVI values.
The inflexion points during NDVI increase and NDVI decrease phases
in the spring and autumn constitute the best predictors in terms of
robustness to data gaps and prediction accuracy. These results are in
agreement with Fisher et al. (2006), Soudani et al. (2008), and Busetto
et al. (2010). For inflexion-point-based phenological metrics, the biases
between in situ and MODIS-based NDVI time series estimates are positive
and vary between 2 and 4 days for the daily and 16-day composite
MODIS NDVI time-series, respectively. For the spring minimum (S1)
and autumn maximum (A1) metrics derived from the daily MODIS
time series (Table 2), the biases are negative, meaning that MODIS
tends to detect the onset of greenness earlier. Negative bias in S1 was
also reported by Fisher et al. (2006) and Soudani et al. (2008). This
result can be explained by the greater sensitivity of (S1, S3) and (A1,
A3) to a lack of NDVI observations during short periods at the beginning
and end of the leaf expansion and leaf yellowing phases, respectively. A
lack of NDVI measurements has the effect of shifting the start of the sigmoid
to the left at S1 and A1 and to the right at S3 and A3. This point will
be discussed in detail below.
The use of the method based on the noise-to-signal ratio developed
in this study (Eq. 6) provided a quantitative insight about the uncertainty
of each phenological metric and its reliability, accounting for the initial
NDVI signal, filtering, and fitting techniques in both deciduous and
evergreen forests. In deciduous forests (Table 3), whether the data
were obtained from the in situ or MODIS NDVI time series, the theoretical
uncertainty yielded significantly lower values for the phenological
transition dates based on the inflexion points for both spring and autumn
phases. This result may be explained by two reasons. First, the rate of
change of the NDVI during these two periods is higher, and thus, the
noise-to-signal ratio is lower, allowing a better fit of the data. Second,
the inflexion point of the NDVI curve is more stable due to the constraint
of symmetry around this position. During the leaf expansion phase, this
date is constrained by two NDVI plateaus in the winter and summer.
During leaf senescence, it is constrained by two other plateaus in the
summer and autumn/winter. For these two reasons, a relatively small
number of NDVI measurements that are of good quality and are well
distributed over the seasons (winter, leaf expansion and senescence
phases, summer and autumn) may be sufficient to obtain good estimates.
The dates of the NDVI minimum increase and the NDVI maximum
are not constrained, and it is necessary to have high-quality NDVI data
during these periods to obtain accurate estimates.
The conclusions underlined above are confirmed by the results of
the sensitivity analysis of the phenological markers to the lack of
NDVI data conducted on the 2008 year in Fontainebleau site, which
exhibit strong phenological pattern and low noise, as summarized in
Fig. 11 and Table 5. The dates of spring and autumnal phenological transitions
derived from inflexion points are the most accurate and most
robust. The use of the inflexion point may even be necessary to estimate
the date of leaf senescence with sufficient accuracy because of the
strong instability of the other two indices (A1 and A2), as shown in
Fig. 11 and Table 5. However, during the autumnal phase, the NDVI
decline is generally slower and less pronounced than during leaf
expansion in the spring because it depends on biological and physical
mechanisms (leaf yellowing, browning, leaf fall, marcescence, and the
mechanical influence of wind and precipitation) that may vary from
year to year. Marcescence, which means that leaves die but do not fall
off of trees in the autumn, is frequent in temperate deciduous forests
and may influence the NDVI decline. In addition, it is highly likely that
the contribution of the soil covered with newly fallen leaves may also
significantly affect the NDVI signal and may explain (at least partially)
the slow decline of the NDVI during the autumn and throughout the
winter, as shown in previous studies (Nagler et al., 2000; Van
Leeuwen & Huete, 1996). These factors may shift the position of the
inflexion point to the right and cause an overestimation of the date of
leaf yellowing and senescence in the autumn.
The results of the sensitivity analysis (Fig. 11) also show that the
inflexion point is the most robust remote sensing-based phenological
metric to gaps in the NDVI time series. The spring phenological transition
prediction error remains less than one week when the number of
Table 3
Average and standard deviation of the theoretical uncertainties (from Eq. 6) calculated for six phenological metrics (Fig. 1) derived from the fitted NDVI time series over deciduous
species and for all years (sample size=10). (S1, S2, S3) for the spring and (A1, A2, A3) for the autumn.
Theoretical uncertainty
(days)
S1 S2 S3 A1 A2 A3
In situ NDVI Average 2 0.7 2 6.5 2.5 6.0
Standard deviation 2.5 0.7 2.5 6.5 2.5 6.5
MODIS daily NDVI Average 7.5 2 7.5 13 5 13
Standard deviation 8.5 2.5 8.5 16.5 6 16.5
MODIS 16-day NDVI Average 8 2.5 8 14 5 14
Standard deviation 12 3.5 12 22 7 22
Table 4
Average theoretical uncertainty (from Eq. 6) for the significant phenological transitions
derived from the fitted NDVI time series over crops, herbaceous savannah and evergreen
forests.
Theoretical uncertainty
(days)
Savannah Crops Evergreen
forests
In situ NDVI Average 3.5 2.5 9.1
MODIS daily NDVI Average 6.6 10.5 19.5
MODIS 16-days NDVI Average 13.5 24.2 162.5
155G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
NDVI observations is less than 30, corresponding to one observation
every two weeks, which is equivalent to the multi-temporal MODIS
NDVI 16-day composite product. This leads us to the conclusion that
the MODIS 16-day composite NDVI data may allow accurate predictions
of spring phenology using the inflexion point of the NDVI curve provided
that the NDVI observations are not contaminated by clouds and that
they are well distributed over the main transition phases. When those
two conditions are not met, the use of MODIS 16-day composite NDVI
may lead to hazardous prediction, as shown in the left part of the graphs
(Fig. 11). However, it is difficult to predict the timing of autumnal
phenological transition with one week of accuracy using the MODIS
16-day NDVI data, and it is still highly unlikely that such accuracy can
be achieved by using the criteria A1 and A3, even when the MODIS
daily NDVI data are used (Table 3).
In the evergreen forest in Puechabon and the tropical rainforest
(Figs. 6, 7 & Table 4), the in situ NDVI time series show low NDVI variations.
In the evergreen Mediterranean forest of holm oak, the NDVI
variations are consistent with the phenology of this species, which is
mainly characterized by two major events: the sprouting of leaves and
shoots in the spring and the shedding of leaves, which is particularly important
during the phase of leaf sprouting in the spring and occasionally
autumn (La Mantia et al., 2003; Soudani et al., 2012). In the tropical
rainforest in French Guiana, the interpretation of the NDVI temporal
patterns is more complex because of the high species diversity in such
forests. Nevertheless, the two periods of NDVI decline, which are
observed occasionally during the first short dry season in February–
March and during the second (main) dry season from the end of August
to the end of October, are concomitant with two periods of lower rainfall
and higher solar radiation. The NDVI decline during the second dry
season coincided with a peak of litterfall, as shown by measurements
of the litterfall regime based on the use of litter traps placed beneath
the canopy in this forest (Soudani et al., 2012).
The seasonal phenological features derived from the daily MODIS
NDVI time series measured over the evergreen forests are quite poor
because only a slight decrease of the NDVI in the spring in the
Puechabon forest was detected. However, 16-day composite MODIS
NDVI time series could not provide a sufficient certainty to precisely
detect any of the phenological features of the evergreen forests.
In the tropical forest (Fig. 7), the MODIS NDVI time series exhibited
strong noise, so none of the temporal features detected in ground-based
NDVI time series could be found in the MODIS NDVI data. In contrast,
the wave shapes observed in the 16-day composite and daily (with a
lower magnitude for the latter) NDVI time series are mainly driven by
seasonal variations in noise intensity. We note that during the main dry
season, this pattern is opposite to that observed in the in situ NDVI time
series. It is also important to note that the seasonal patterns of the
MODIS 16-day composite NDVI shown in Fig. 7 are similar to those
obtained in previous works (Huete et al., 2006; Saleska et al., 2007).
These studies concluded that there is an increase in the canopy greenness
during dry periods. Our results suggest that this pattern may not reflect a
phenological signal but a variation of the noise intensity in the NDVI
observations. The use of the MODIS daily or 16-day composite data without
any ground-based reference may therefore be misleading.
The results from the savannah and crop sites (Fig. 8, Table 4) pinpoint
an important limitation of the MODIS NDVI time series in the detection
phenological features: while the actual features were detected, strong
Fig. 11. Relationships of phenology prediction error (days) (blue line, left y-axis) and the length of the simulated NDVI time-series (black line — right y-axis) versus the average of the
absolute values of the first derivatives of fits of simulated NDVI time-series (x-axis). The red lines are the confidence limits (95%) of the phenology prediction error.
Table 5
Sample sizes defining the stability region of the root mean square errors (lower RMSE±1 day) of the phenological estimates in deciduous forests due to the introduction of artificial data gaps in
the NDVI time series (first line). (Second line): Proportions of the remaining NDVI data (percentage of the total number of observations of the entire in situ NDVI time series, n=189). (S1, S2, S3)
for the spring and (A1, A2, A3) for the autumn. * [left–right] corresponds to both sides of the stability region for each phenological metric (see Fig. 11).
Phenological metrics: S1 S2 S3 A1 A2 A3
Average length of simulated NDVI time series (n)
[left–right*]
83–73 75–51 94–97 124–118 115–73 139–98
(In % of full NDVI dataset)
[left–right]
44%–39% 40%–27% 49%–51% 65%–62% 61%–40% 73%–51%
156 G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
errors occurred due to mixed pixels or bad sky conditions coinciding
with phenological events. These errors could not be addressed by the
tested filters or by the noise-to-signal ratio. At the savannah site, which
may be heterogeneous at the MODIS pixel scale, some comparable features
could be found between the in situ and MODIS series such as the
NDVI drop due to a short dry period in March 2008, and the green-up
at the end of 2009, while there were localized mismatches, notably
around June, when the area is burnt. This observation may indicate a possible
scale mismatch between the in situ and MODIS observations.
5. Conclusion
In this study, in situ NDVI time series allowed us to directly assess
the accuracy of MODIS-derived phenological estimates. In deciduous
forests, inflexion points of a double sigmoid model fitted to NDVI
data allow for the most accurate estimates of the onset of greenness
in the spring and the onset of yellowing in the autumn (RMSE≤one
week). Phenological metrics delimiting the leaf expansion phase in
the spring and the leaf senescence phase in the autumn, which are
identical to those provided in MODIS Global Vegetation Phenology
product (MDC12Q2), are less robust to data gaps, and they can be
subject to large biases of approximately two weeks or more during
the autumn phenological transitions. The inflexion point detection
was shown to be more precise and less sensitive to data gaps than
these metrics. However, the use of the date at the beginning of the
NDVI decline in the autumn (identical to the onset of the greenness decrease
in MDC12Q2) instead of the date at the inflexion point can be
justified because of the slow, monotonic decline of the NDVI during
the autumn and winter, which could be due to the contribution of freshly
fallen leaf litter and because the phenomenon of marcescence can
cause a shift to the right of the inflexion point that could lead to
overestimation of the onset of leaf yellowing. In the evergreen forests,
in situ NDVI time series describe the phenology with high fidelity despite
small temporal changes in the canopy foliage. However, MODIS
is unable to provide consistent phenological patterns. In savannah and
crops, the detection of phenological patterns could be achieved but
was hampered by a seasonal variation of noise amplitudes. Similarly,
in the tropical rainforest, the temporal pattern exhibited in the MODIS
16-day composite NDVI time series is more likely due to a pattern of
noise in the NDVI data, structured according to both rainy and dry seasons
rather than to phenological changes. More investigations are needed,
but in all cases, this result leads us to conclude that the MODIS time
series in tropical rainforests should be interpreted with great caution.
Acknowledgments
The authors thank GIP ECOFOR and F-ORE-T «Observatoires de
Recherche en Environnement (ORE) sur le Fonctionnement des
Écosystèmes Forestiers» ECOFOR, INSU, Ministère de l'Enseignement
Supérieur et de la Recherche for funding the project that enabled
manufacturing and deployment of the NDVI sensors in different vegetation
types. The last author would also like to thank the University of
Paris Sud for the PhD grant given to Mr Gabriel Hmimina. We would
like to express our profound gratitude to Laurent Vanbostal and Daniel
Berveiller for their help in manufacturing the NDVI sensors and measurements.
We also thank all those involved in the data collection process
at all study sites. We are very grateful for the thorough and helpful
comments from the reviewers of the manuscript.
Appendix A. Supplementary data
Supplementary data associated with this article can be found in the
online version, at http://dx.doi.org/10.1016/j.rse.2013.01.010. These
data include Google map of the most important areas described in this
article.
References
Aubinet, M., Moureaux, C., Bodson, B., Dufranne, D., Heinesch, B., Suleau, M., et al.
(2009). Carbon sequestration by a crop over a 4-year sugar beet/winter wheat/seed
potato/winter wheat rotation cycle. Agricultural and Forest Meteorology, 149, 407–418.
Beck, P. S. A., Atzberger, C., Høgda, K. A., Johansen, B., & Skidmore, A. K. (2006). Improved
monitoring of vegetation dynamics at very high latitudes: A new method using
MODIS NDVI. Remote Sensing of Environment, 100, 321–334.
Bonal, D., Bosc, A., Goret, J. Y., Burban, B., Gross, P., Bonnefond, J. M., et al. (2008). The
impact of severe dry season on net ecosystem exchange in the Neotropical rainforest
of French Guiana. Global Change Biology, 14, 1917–1933.
Box, J. E., Bromwich, D. H., Veenhuis, B. A., Bai, L. S., Stroeve, J. C., Rogers, J. C., et al.
(2006). Green land Ice Sheet Surface Mass Balance Variability (1988–2004) from
calibrated Polar MM5 Output. Journal of Climate, 19, 2783–2800.
Busetto, L., Colombo, R., Migliavacca, M., Cremonese, E., Meroni, M., Galvagno, M., et al.
(2010). Remote sensing of larch phenological cycle and analysis of relationships
with climate in the Alpine region. Global Change Biology, 16, 2504–2517.
Castaldi, S., de Grandcourt, A., Rasile, A., Skiba, A., & Valentini, R. (2010). Fluxes of CO2, CH4
and N2O from soil of burned grassland savannah of central Africa. Biogeosciences
Discussions, 7, 4089–4126.
Chen, B., Black, T. A., Coops, N. C., Hilker, T., Trofymow, (Tony) J. A., & Morgenstern, K.
(2009). Assessing tower flux footprint climatology and scaling between remotely
sensed and Eddy covariance measurements. Boundary-Layer Meteorology, 130,
137–167.
Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsushita, B., & Eklundh, L. (2004). A simple
method for reconstructing a high-quality NDVI time series data set based on the
Savitzky–Golay filter. Remote Sensing of Environment, 91, 332–344.
Delbart, N., Le Toan, T., Kergoat, L., & Fedotova, V. (2006). Remote sensing of spring
phenology in boreal regions: A free of snow-effect method using NOAA-AVHRR
and SPOT-VGT data (1982–2004). Remote Sensing of Environment, 101, 52–62.
Delpierre, N., Soudani, K., François, F., Köstner, B., Pontailler, J. Y., Aubinet, M., et al.
(2009). Exceptional carbon uptake in European forests during the warm spring
of 2007: A data-model analysis. Global Change Biology, 15, 1455–1474.
Dufranne, D., Moureaux, C., Vancutsem, F., Bodson, B., & Aubinet, M. (2011). Comparison
of carbon fluxes, growth and productivity of a winter wheat crop in three contrasting
growing seasons. Agriculture, Ecosystems & Environment, 141, 133–142.
Fensholt, R., Sandholt, I., Proud, S. R., Stisen, S., & Rasmussen, M. O. (2010). Assessment
of MODIS sun-sensor geometry variations effect on observed NDVI using MSG
SEVIRI geostationary data. International Journal of Remote Sensing, 31, 6163–6187.
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,
265–279.
Ganguly, S., Friedl, M. A., Tan, B., Zhang, X., & 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.
Goerner, A., Reichstein, M., & Rambal, S. (2009). Tracking seasonal drought effects on
ecosystem light use efficiency with satellite-based PRI in a Mediterranean forest.
Remote Sensing of Environment, 113, 1101–1111.
Granier, A., Bréda, N., Longdoz, B., Gross, P., & Nao, J. (2008). Ten years of fluxes and stand
growth in a young beech forest at Hesse, North-eastern France. Annals of Forest Science,
65(7), 704.
Heidinger, A. K., Anne, V. R., & Dean, C. (2001). Using MODIS to estimate cloud contamination
of the AVHRR data record. Journal of Atmospheric and Oceanic Technology, 19,
586–601.
Hird, J., & McDermid, G. J. (2009). Noise reduction of NDVI time-series: An empirical
comparison of selected techniques. Remote Sensing of Environment, 113, 248–258.
Huete, A. R., 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, 95–213.
Huete, A. R., Didan, K., Shimabukuro, Y. E., Ratana, P., Saleska, S. R., Hutyra, L. R., et al.
(2006). Amazon rain forests green-up with sunlight in dry season. Geophysical
Research Letters, 33, L06405, http://dx.doi.org/10.1029/2005GL025583.
Jönsson, P., & Eklundh, L. (2002). Seasonality extraction by function fitting to timeseries
of satellite sensor data. IEEE Transactions on Geoscience and Remote Sensing,
40, 1824–1832.
Justice, C. O., Vermote, E., Townshend, J. R. G., Defries, R., Roy, D. P., Hall, D. K., et al.
(1998). The Moderate Resolution Imaging Spectroradiometer (MODIS): Land remote
sensing for global change research. IEEE Transactions on Geoscience and Remote Sensing,
36, 1228–1249.
Kaduk, J., & Heimann, M. (1996). A prognostic phenology scheme for global terrestrial
carbon cycle models. Climate Research, 6, 1–19.
La Mantia, T., Cullotta, S., & Garfi, G. (2003). Phenology and growth of Quercus ilex (L.) in
different environmental conditions in Sicily (Italy). Ecologia Mediterranea, 29, 15–25.
Ma, M., & Veroustraete, F. (2006). Reconstructing pathfinder AVHRR land NDVI timeseries
data for the Northwest of China. Advances in Space Research, 37, 835–840.
McLachlan, G., & Peel, D. (2000). Finite mixture models. New York: John Wiley and Sons.
Metzger, S., Junkermann, W., Mauder, M., Beyrich, F., Butterbach-Bahl, K., Schmid, H. P.,
et al. (2012). Eddy Covariance flux measurements with a weight-shift microlight
aircraft, Atmos. Atmospheric Measurement Techniques, 5, 2591–2643.
Moulin, S., Kergoat, L., Viovy, N., & Dedieu, G. (1997). Global-scale assessment of vegetation
phenology using NOAA/AVHRR satellite measurements. Journal of Climate,
10, 1154–1170.
Nagler, P. L., Daughtry, C. S. T., & Goward, S. N. (2000). Plant litter and soil reflectance.
Remote Sensing of Environment, 71, 207–215.
Pontailler, J. -Y., & Genty, B. (1996). A simple red:far-red sensor using gallium arsenide
phosphide detectors. Functional Ecology, 10, 535–540.
157G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158
Pontailler, J. -Y., Hymus, G. J., & Drake, B. G. (2003). Estimation of leaf area index using
ground-based remote sensed NDVI measurements: Validation and comparison
with two indirect techniques. Canadian Journal of Remote Sensing, 29, 381–387.
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.
Saleska, S. R., Didan, K., Huete, A. R., & da Rocha, H. R. (2007). Amazon forests green-up
during 2005 drought. Science, http://dx.doi.org/10.1126/science.1146663.
Schwartz, M. D., Reed, B. C., & White, M. A. (2002). Assessing satellite-derived start-of-season
measures in the conterminous USA. International Journal of Climatology, 22, 1793–1805.
Sims, D. A., Rahman, A. F., Vermote, E. F., & Jiang, Z. (2011). Seasonal and inter-annual
variation in view angle effects on MODIS vegetation indices at three forest sites.
Remote Sensing of Environment, 115, 3112–3120.
Sokal, R. R., & Michener, C. D. (1958). A statistical method for evaluating systematic
relationships. University of Kansas Scientific Bulletin, 28, 1409–1438.
Soudani, K., François, C., le Maire, G., Le Dantec, V., & Dufrêne, E. (2006). Comparative analysis
of IKONOS, SPOT, and ETM+ data for leaf area index estimation in temperate coniferous
and deciduous forest stands. Remote Sensing of Environment, 102, 161–175.
Soudani, K., Hmimina, G., Delpierre, N., Pontailler, J. -Y., Aubinet, M., Bonal, D., et al.
(2012). Ground-based Network of NDVI measurements for tracking temporal dynamics
of canopy structure and vegetation phenology in different biomes. Remote
Sensing of Environment, 123, 234–245.
Soudani, K., le Maire, G., Dufrêne, E., François, C., Delpierre, N., Ulrich, E., et al. (2008).
Evaluation of the onset of green-up in temperate deciduous broadleaf forests derived
from Moderate Resolution Imaging Spectroradiometer (MODIS) data. Remote Sensing
of Environment, 112, 2643–2655.
Studer, S., Stöckli, R., Appenzeller, C., & Vidale, P. (2007). A comparative study of satellite
and ground-based phenology. International Journal of Biometeorology, 51, 405–414.
Suzuki, R., Nomaki, T., & Yasunari, T. (2003). West–east contrast of phenology and climate
in northern Asia revealed using a remotely sensed vegetation index. International
Journal of Biometeorology, 47, 126–138.
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, 98–114.
Van Leeuwen, W. J. D., & Huete, A. R. (1996). Effects of standing litter on the biophysical
interpretation of plant canopies with spectral indices. Remote Sensing of Environment,
55, 123–138.
Viovy, N., Arino, O., & Belward, A. S. (1992). The Best Index Slope Extraction (BISE): A
method for reducing noise in NDVI time-series. International Journal of Remote Sensing,
13, 1585–1590.
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., & Nemani, R. R. (2006). Real-time monitoring and short-term forecasting
of land surface phenology. Remote Sensing of Environment, 104, 43–49.
White, M. A., Nemani, R. R., Thornton, P. E., & Running, S. W. (2002). Satellite evidence of
phenological differences between urbanized and rural areas of the Eastern United
States deciduous broadleaf forest. Ecosystems, 5, 260–273.
White, M. A., Thorntorn, P. E., & Running, S. W. (1997). A continental phenology model
for monitoring vegetation responses to interannual climatic variability. Global
Biogeochemical Cycles, 11, 217–234.
Wolfe, R. E., Nishihama, M., Fleig, A. J., Kuyper, J. A., Roy, D. P., Storey, J. C., et al. (2002).
Achieving sub-pixel geolocation accuracy in support of MODIS land science. Remote
Sensing of Environment, 83, 31–49.
Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C. F., Gao, F., et al. (2003).
Monitoring vegetation phenology using MODIS. Remote Sensing of Environment,
84, 471–475.
Zhang, X., & Goldberg, M. D. (2011). Monitoring fall foliage coloration dynamics using
time-series satellite data. Remote Sensing of Environment, 115, 382–391.
158 G. Hmimina et al. / Remote Sensing of Environment 132 (2013) 145–158