# Evaluation of Temporal Resolution Effect in Remote Sensing Based Crop Phenology Detection Studies

Hu Zhao, Zhengwei Yang, Liping Di, Zhiyuan Pei

<figure><img src="image_1.png" /></figure>

► To cite this version:

Hu Zhao, Zhengwei Yang, Liping Di, Zhiyuan Pei. Evaluation of Temporal Resolution Effect in Remote Sensing Based Crop Phenology Detection Studies. 5th Computer and Computing Technologies in Agriculture (CCTA), Oct 2011, Beijing, China. pp.135-150, 10.1007/978-3-642-27278-3_16. hal-01360974

HAL Id: hal-01360974

https://inria.hal.science/hal-01360974

Submitted on 6 Sep 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Evaluation of temporal resolution effect in remote sensing
based crop phenology detection studies

Hu Zhao<sup>1,2</sup>, Zhengwei Yang<sup>3</sup>, Liping Di<sup>2</sup>, Zhiyuan Pei<sup>1</sup>

¹ Agriculture Resource Monitoring Station, Chinese Academy of Agriculture Engineering, 41 Maizidian Street, Chaoyang, Beijing 100125, China
zhaohu@whu.edu.cn, peizhiyuan@tom.com

² Center for Spatial Information Science and Systems, George Mason University, 6301 Ivy Lane, Suite 620, Greenbelt, MD 20770, USA
ldi@gmu.edu

³ National Agricultural Statistics Service, United States Department of Agriculture, 3251 Old Lee Highway, Room 305, Fairfax, VA 22030, USA
Zhengwei_Yang@nass.usda.gov

**Abstract.** Remote sensing based phenology detection method has been employed to study agriculture, forestry and other vegetations for its potential to reflect the variations in climate change. These studies usually utilized time series Normalized Difference Vegetation Index (NDVI) generated from various sensors through a Maximum Value Compositing (MVC) process, which minimized the contamination from cloud and simultaneously introduce degradation of temporal accuracy. In this study, we assess the impact of temporal resolution on crop phenology derivation researches by comparing three different Moderate Resolution Imaging Spectroradiometer (MODIS) datasets: daily surface reflectance, 8 day composited surface reflectance and 16 day composited NDVI. The surface reflectance data were first filtered by employing auxiliary data which contained quality and viewing geometry information, and then used to calculate NDVI with specific date. A least square method was taken to fit the survival data points to double logistic function. And finally, seven time-related metrics were obtained and matched with field observation crop phenology stages. These remote sensing derivate phenology dates were compared to National Agricultural Statistics Service (NASS) weekly crop progress reports to evaluate the capability of these datasets in temporal sensitive studies. The results illustrated that daily surface reflectance datasets were the most accurate source for time-sensitive studies. However, extra ancillary datum and appropriate denoising techniques should be applied to reconstruct the time series curve. Phenology matching process is a necessary step before detecting phenological information from remote sensing imagery for specific land cover type since same phenological stages of different crop types might have different counterparts on time series curve.

**Keywords.** Crop phenology; temporal resolution evaluation; least square; double logistic function fitting
# 1. Introduction

Remote sensing based phenology detection techniques offer great opportunity for many climate change related researches, such as land surface vegetation dynamics (Beck et al., 2006), agriculture management (Garrigues et al., 2008), forestry monitoring (Ahl et al., 2006) and global warming (Parmesan, 2007). Phenology is the study of the timing of cyclical events in nature such as springtime vegetation budburst or seasonal bird migrations (Badeck et al., 2004), while crop phenology refers to the developments, differentiation and initiation of organs of a crop (Hodges, 1991). Many researchers have studied phenological change using field-collected data, which is prohibitive for large scale monitoring due to the need for repeat data collection in order to monitor change from previous years (Thayn and Price, 2008). The situation has been totally changed after satellite sensor imagery, which record multi-scale (spatial or temporal) information of ground surface, was applied to monitor land vegetation dynamics.

A number of research works have been conducted to monitor 'generic or mixed' vegetation phenological events (Justice et al., 1985; Reed et al., 1994; Zhang et al., 2003; Fisher et al., 2006; White and Nemani, 2006; Stöckli et al., 2008), and many efforts have been made to develop a general procedure for remote sensing based crop phenology detection. Sakamoto et al. (Sakamoto et al., 2005) generated planted, heading and harvest dates of rice paddy in Japan using time series MODIS Enhanced Vegetation Index (EVI) datasets. Vina et al. (Vina et al., 2004) presented a visible atmospherically resistant index (VARI) and accumulated growing degree days (AGDD) to detect maize physiological transitions that undetectable from NDVI. Duchemin et al. (Duchemin et al., 2006) monitored wheat phenology through the analysis of relationship among evapotranspiration, crop coefficients, leaf area index (LAI) and NDVI. Islam and Bala (Islam and Bala, 2008) made use of NDVI and LAI to evaluate potato phenological characteristics by modeling its essential phenological metrics. Each of these pioneering works has its special contribution to remotely-sensed crop phenology estimation research.

However, the majority of above-mentioned studies utilized multi-temporal NDVI/EVI/LAI values processed with the MVC technique (Holben, 1986), which effectively removed the cloud-contaminated data points in original observations series. The maximum value of the composited period (usually range from 8 to 16 days) was assigned to this multiday period rather than the date when this value was captured (Thayn and Price, 2008). This processing step imbedded a temporal error of 8 (or 4) days in average and 15 (or 7) days maximum because a maximum value in 16 (or 8) days is selected and assigned to all of 16 (or 8) days. This could not be neglected in time-sensitive studies, especially for crop phenology detection. Therefore, it is essential to assess the effect of this composited process on crop phenology estimation accuracy and seek out an appropriate strategy for the remote sensing based crop phenology estimation.

The objectives of this paper are to: a) check the feasibility of employing daily surface reflectance datasets and its corresponding ancillary data to generate NDVI
curve and consequently calculate crop phenology dates; b) make a quantitative
comparison between different temporal resolutions (daily, 8 day and 16 day
composited) with fixed spatial scale (250m*250m / pixel) about their effects on
phenological stages detection studies; c) attempt to develop a framework for
matching ‘specific’ phenology stages with curve-generated Julian dates. Two crop
types (corn and soybean) have been tested on 50 sites distributing over US Corn
Belt in this experiment. For each test site, seven Julian dates related potential
phenology stages have been generated and matched with specific crop
phenological events.

## 2. Data Description

2.1. Study sites

In this study, the State of Iowa (range from 40° 23′ N to 43° 30′ N and 90° 8′ W to 96° 38′ W) in U.S. Corn Belt is selected for experiment because it is a predominant corn and soybean state and it is one of the most important corn and soybean production areas in United States. Iowa locates in the middle of American and its terrain is mainly dominated by plain. All these facts will help us to simplify the study and to keep focus on the method for phenological information extraction from the MODIS products. In this study, 50 sites (25 sites each for corn and soybean) within Iowa State, as shown in figure 1, were selected for method and data experiments.

**2.2. Datasets description and preprocessing**

Three MODIS products have been used in this experiment: Surface Reflectance Daily L2G Global 250m (MOD09GQ), Surface Reflectance 8-Day L3 Global 250m (MOD09Q1) and Vegetation Indices 16-Day L3 Global 250m (MOD13Q1). The MOD09GQ product provides an estimate of the surface spectral reflectance as it would be measured at ground level in the absence of atmospheric scattering or absorption. Low-level data can be corrected for atmospheric gases and aerosols, yielding a level-2 basis for several higher-order gridded level-2 (L2G) and level-3 products. Bands 1 (RED) and 2 (near infrared (NIR)) are provided at a 250-meter resolution in a daily gridded L2G product in the Sinusoidal projection. This product is meant to be used in conjunction with the Surface Reflectance Daily L2G Global 1km and 500m (MOD09GA) where important quality and viewing geometry information is stored. MOD09Q1 is a 8-day gridded level-3 product that stores an estimate of the surface spectral reflectance. Each MOD09Q1 pixel contains the best possible L2G observation during an 8-day period as selected on
the basis of high observation coverage, low view angle, the absence of clouds or
cloud shadow, and aerosol loading. MOD09GQ and MOD09Q1 are first
calculated NDVI through the following formula: NDVI = (NIR-RED) / (NIR+RED). MOD13Q1, unlike MOD09GQ and MOD09Q1, is designed to
provide consistent spatial and temporal comparisons of vegetation conditions.
Blue, red, and near-infrared reflectance are used to determine the MODIS daily
vegetation indices. The MODIS NDVI complements NOAA's Advanced Very
High Resolution Radiometer (AVHRR) NDVI products and provides continuity
for time series historical applications. Although MODIS also includes EVI that
minimizes canopy background variations and maintains sensitivity over dense
vegetation conditions, we only extract NDVI values for comparison purpose
(USGS, 2009). All these datasets are captured between March 31 and December
11 from 2006 to 2008. This period includes 256 yeardays and totally covers the
growing season for both corn and soybean in each year according to the statistical
information published by United State Department of Agriculture (USDA).

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Distribution of test sites in Iowa State, United States. Filled symbol pentagram (★) and triangle (▲) represented corn and soybean respectively. The background color gradient indicates the general condition of elevation.</figcaption></figure>

Many criterions, besides cloud status information recorded in quality flag, have
been presented previously to remove noisy data points affected by thick clouds,
gaseous and aerosol scattering as well as bidirectional reflectance distribution
function (BRDF) distortion. For example, Sakamoto et al. (Sakamoto et al., 2005)
considered the pixels whose blue band reflectance was greater than 10% or whose
sensor zenith angle was more than 32.25° should be removed as abnormal data.
Points with a random NDVI increase greater than 0.4 during 20 days are also
rejected as misleading data (Chen et al., 2004). We adopt these extra criterions to
reserve the pixels in time-series of a test site for above-mentioned three datasets:
a) the cloud and cloud shadow status should be ‘clear’; b) the blue band value
should be lower than 0.1; c) the sensor zenith angle is smaller than 30° and d)
difference between survival data point and its Inverse Distance Weighted (IDW) estimation should be small than 0.1(Lu et al., 2007). These standards could ensure that the survival sensor recorded reflectance values are not influenced by cloud or other noises.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Time-series datasets of MOD09GQ (open circle ‘o’), MOD09Q1 (dotted line) and MOD13Q1 (solid line) after preprocessing. Many MOD09GQ cloud-contaminated data points are removed and replaced with NAN so the time series looks like discrete.</figcaption></figure>

The National Agricultural Statistics Service (NASS), the USDA data collection arm, conducts hundreds of surveys every year and prepares weekly crop progress report containing crop planting and emergence percentage estimates. Crop progress data to be used for validation were acquired from the weekly 2006-2008 Crop Progress reports produced by the NASS for Iowa. (These reports are publicly available at http://www.nass.usda.gov/Statistics_by_State/Iowa/Publications/Crop_Progress_&_Condition/index.asp.) According to these reports, there are nine districts in Iowa State: northwest (NW), north center (NC), northeast (NE), west center (WC), center (C), east center (EC), southwest (SW), south center (SC) and southeast (SE). Each of these districts has its own validation information (planted, emerged, silking, etc.) for different kinds of crops. Each stage was described as a period as its percentage covers from 0 to 100. Generally, the 50 percent corresponding Julian date of each stage could be selected as final validation information (Wardlow et al., 2006). For corn, the report offered nine phenology stages: planted, emerged, tasseled, silking, milk, dough, dent, mature and harvest. For soybean, there are seven stages: planted, emerged, blooming, setting pods, leaf turning color, mature and harvest. All phenological stage definitions could be accessed at http://www.nass.usda.gov/Charts_and_Maps/Crop_Progress_&_Condition/Terms/index.asp. The Cropland Data Layer (CDL) program (http://www.nass.usda.gov/research/Cropland/SARS1a.htm), also offered by NASS, provides the information about what kind of crop has been planted on a specific field with coordinate information. USDA NASS online products are the
main source of our experiment validation information.

# 3 Methodologies

## 3.1. Curve fitting strategy for daily data

After preprocessing steps, the time series are usually not equal-spaced or continuous since some abnormal pixels have been removed (as shown in Figure 2). Meanwhile, even those survival points are still contain some undefined noise (soil background, water content, moisture, saturation, pre-crop vegetation, etc.). Curve fitting is proved to be one of the most effective strategies to eliminate those noises (Bradley et al., 2007).

Various techniques have been implemented to model the survival data points and to reconstruct the removed values, such as Best Index Slope Extraction (BISE) (Viovy et al., 1992), Fourier analysis (Roerink et al., 2000; Geerken et al., 2005), Savitzky-Golay filter and its improvement (Savitzky and Golay, 1964; Chen et al., 2004), Support Vector Machine (Sun et al., 2006), asymmetric Gaussian function (Jönsson and Eklundh, 2002), mean-value iteration (MVI) (Ma and Veroustraete, 2006), wavelet transformation and its improvements (Sakamoto et al., 2005; Lu et al., 2007) and logistic function (Zhang et al., 2003; Jönsson and Eklundh, 2004; Beck et al., 2006; Beck et al., 2007). Hird and McDermid (Hird and McDermid, 2009) made a model-based empirical comparison of some frequently-used NDVI time series denoising techniques and concluded that double logistic function fitting strategy demonstrated the superiority over other methods. The formula could be expressed as equation (1), where $wNDVI$ stands for minimum NDVI value (often arose in winter) of the final curve fitting time series. In this test, $wNDVI$ represents bare soil or pre-crop vegetation NDVI value. This value is often represented as geometric mean of NDVI values in non-growing season (winter or bare soil background)(Beck et al., 2006). According to the NASS crop progress reports; there is no crop or plant available in the corn or soybean fields before April from 2006-2008. This means that we could set the non-growing season as the date before the 120th day of the year. $mNDVI$ denotes the maximum NDVI and its difference with $wNDVI$ represents the amplitude of NDVI time series; $mS$ and $S$ are variables relating to the inflection point in vegetative growing season, while $mA$ and $A$ are referred to senescence season. Their relationship as shown in Figure 3 could be expressed (Beck et al., 2007).

$$NDVI_t = wNDVI + (mNDVI - wNDVI) \times \left( \frac{1}{1 + e^{-mS(t-S)}} + \frac{1}{1 + e^{-mA(t-A)}} - 1 \right) \quad (1)$$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Example of the double logistic function $NDVI_t$ (equation (1)) used to model the yearly NDVI time series. It is defined by six parameters: the minimum and maximum NDVI ($wNDVI$ and $mNDVI$), two inflection points ($S$ and $A$), and the rate of increase or decrease in NDVI at $S$ and $A$, respectively ($mS$ and $mA$) (Beck et al., 2007). Feature points G, P, Q and E stand for four characteristic curvature values, while M is the point when the NDVI reaches its maximum (or the central points of NDVI saturation period in case of maximum value continues during a period).</figcaption></figure>

Least square method is commonly applied in function fitting to select the ‘best fit’ parameters (Steinier et al., 1972; Sun et al., 2006; Markwardt, 2008) and often expressed as follows:

$$ \chi^2 = \sum_{i=1}^{n} \left[ w_i (NDVI_i - NDVI_{ti}) \right]^2 \quad (2) $$

Where $w_i$ is the weight of ith NDVI value. When $\chi^2$ reaches its minimum, the parameters are optimal. Beck et al. (Beck et al., 2006) presented a two-step schema to determine weight, which was supposed to get the upper envelop of the survival data. First, all weights were set equally as 1 to obtain the reference curve line. In the second step, the weights were reassigned different values according to relative relationship between original values and the reference curve line: if the original value is above (below) the reference curve line, it gets a higher (lower) weight. $NDVI_i$ is the ith originally calculated NDVI value, and $NDVI_{ti}$ is the corresponding estimated value. Since there are contaminated data in daily NDVI, we have made a little adjustment to this schema: set all contaminated points as zero and then applied Beck’s two-step weighting strategy. The weights of those contaminated data were set to zero throughout the whole weighting process. This ensures that those contaminated points will not affect the phenology estimation. To get more accurate phenology estimation, the composited products have been linearly interpolated to reconstruct contaminated daily data with ENVI Interactive Data Language (IDL) routine ‘interpol’. However, all interpolated values are weighted as zero during the weighting process to keep the ‘original’ survival
points' affect on resultant NDVI curve.

## 3.2. Phenology stage estimation based on NDVI curve feature points

As shown in Figure 3, a NDVI curve could be characterized with seven curve feature points: G, S, P, M, Q, A, E. The inflection points S and A could be obtained from the maximum and minimum of the first derivate of estimated NDVI curve. The point M represents the date of the maximum NDVI value. If there exist more than one date, M will be the central point of these maximum-value dates. Points S, P, M, Q, and A divided the estimated NDVI curve into four parts. In each of these four parts, the maximum curvature could be obtained and marked as G, P, Q and E according to the following formula:

$$ \kappa = \frac{|NDVI'_{i}|}{(1 + (NDVI'_{i})^2)^{3/2}} \quad (3) $$

Julian date $t$ is the independent variable. These seven temporal points will be used to match with actual phenological dates.

After calculating seven feature point dates for each test site, these dates should be matched with statistical stages to identify the relationship between curve-derived dates and specific statistical phenological event dates (treated as actual phenology stages). Before matching, these remote sensing derived dates could be grouped as set G, set S, set P, set M and set Q. We also marked the statistical validation phenology dates (planted, emerged, silking, milk, mature and harvest, etc.) as phenology stage A, stage B, stage C, stage D, etc. By calculating root-mean-square error (RMSE) between each set and stage, an error matrix was produced through formula (4), in which $ε_{i,j}$ was the RMSE between ith set $Set_{i,k}$ and jth phenology stage Stagej (includes m test sites, k was the serial number of test sites).

$$ \varepsilon_{i,j} = \sqrt{\frac{\sum_{k=1}^{m} (Set_{i,k} - Stage_{j})^2}{m}} \quad (4) $$

All RMSE matrices for nine districts were first computed with test sites that falls into each district, and the state-level RMSE matrix was then obtained by weighting every district matrix. The weight of each district-level matrix is the proportion of district test site numbers to the total test sites. Figure 4 illustrates the whole process of final matching pair generation and error analysis.

Here are three basic principles to determine the intermediate match pairs for each dataset: 1) smaller RMSE, the minimum $ε_{i,j}$ was selected to determine the first match pair, for example, Setp and Stageq; 2) time line consistency, the sets
before (or after) Setp will only be matched with stages before (or after) Stageq, which is proposed to keep the consistency over time line. The final matching pair is generated by considering daily, 8-day and 16-day intermediate matching results, which should abide by the 3) majority rule: first we pick out the matching pairs in which these three datasets have three same intermediate results. Then we choose the matching pairs that two of these three datasets have same matching pairs and the other dataset has different result which will be re-decided same as the former two counterparts. If one stage refers to three different datasets, the result will be determined by the minimum RMSE value. The second rule, time line consistency, is simultaneously applied during this process.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Flowchart of matching pairs determination and error analysis. It includes 5 steps, step 1: calculating district level RMSE matrix, step 2: integrating 9 district results into state level RMSE matrix; step 3: determining intermediate matching pairs for daily, 8-day and 16-day datasets; step 4: reaching a final matching pair for corn and soybean and step 5: conducting error analysis by using matching pair results in step 3 and 4.</figcaption></figure>

### 3.3. Evaluation methods

We used two complementary approaches to evaluate the performance of each data set: phenology matching accuracy and phenology estimation error. First, we define three intermediate matching results and the final matching pair as: Imp01, Imp08, Imp16 and Fmp. They record the matching pairs as phenology “milk” to metric “M”, phenology “mature” to metric “A”, etc. We interpreted phenology matching accuracy (εᵢⱼ) as ratio between number of final matching pairs (Fmp) and the number of intersection of intermediate matching pairs (Imp01, Imp08 or Imp16) and the final matching pairs (Fmp). The matching accuracy reflects the degree of matching between intermediate matching results and the final matching
pairs, εmi equal to 1 means the intermediate result is the same as the final result. The bigger the εmi is, the worse they coincide with each other. The formula could be expressed as:

$$ \epsilon_{mi} = \frac{N(Fmp)}{N(Fmp \cap Imp_i)} \quad (5) $$

where εmi stands for the phenology matching accuracy of datasets i, N is the function that extract the number of a set. Impi represents Imp01, Imp08 or Imp16 and symbol ∩ means the intersection of the two datasets.

The phenology estimation error (εi,j) refers to the summation of final matching pair estimation errors, which is determined by the ratio of RMSE and the mean of the matching pairs' RMSE. For phenology estimation error εi,j, it happened between different datasets and could be listed as the ratio of final matching RMSE value RMSEij to the mean of its counterparts RMSEij (of different datasets). Formula 6 defines phenology estimation error.

$$ \epsilon_{ei} = \sum_{j=0}^{m} \frac{RMSE_{ij}}{\frac{1}{3} \times \sum_{i=0,08,16} RMSE_{ij}} \quad (6) $$

The final error is the summation of matching error and estimation error. And it can be expressed as:

$$ \epsilon_i = a \cdot \sum \epsilon_{mi} + b \cdot \sum \epsilon_{ei} \quad (7) $$

Where a and b are the weight for phenology matching error and phenology estimation error. Here we consider these two kinds of error have the same contribution to the final total error, which sets both a and b as 1.0.

## 4. Results and Discussion

According to the basic rules described in section 3.2, we got the final matching pairs for these three different temporal resolution datasets in 2006-2008 (figure 5). For corn, the final matching pairs in 2007 and 2008 are the same: field observation phenology stage 'planted' corresponds to remote sensing derived first metric 'G'; emerged to 'S'; tasseled to 'P'; milked to 'M'; dough to 'Q'; mature to 'A' and harvest to 'E'. There is no metrics corresponded to phenology stage corn silking and dent in these two years. However, in 2006, the matching pairs seem slightly different: the stage emerged refers to metric 'G' and there is no stage matching with metrics 'S'; and for metric 'A', the correspondent phenological stage is dent, not dough. Other matching pairs are the same with the following two years. According to the statistical observation information obtained
from NASS, the maize dough phenology stage is recorded 5 (South center of IOWA, 2007) to 15 (North center of IOWA, 2007) days before dent stage. When comparing to the RMSE error between remote sensing derived phenological dates and the validation information (in table 1 and 2), the duration between dough and dent is somewhat negligible. This indicates that the matching results are not varying greatly across different years since phenology changes every year. Comparing with vegetation types, comprehensive human being participation (seed selection, field management, irrigation status, etc.) is the most important factor of crop land cover development, which is sensitive to both natural and manmade microenvironment.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Final matching pairs for corn (left) and soybean (right) in 2006-2008.</figcaption></figure>

<table><caption>**Table 1.** State level RMSE result for corn in 2007 (days)</caption><thead><tr><th></th><th>Planted</th><th>Emerged</th><th>tasseled</th><th>Silking</th><th>milk</th><th>Dough</th><th>Dent</th><th>Mature</th><th>Harvest</th></tr></thead><tbody><tr><td>01G</td><td>7.6</td><td>11.3</td><td>64.9</td><td>69.4</td><td>84.2</td><td>96.2</td><td>108.6</td><td>127.4</td><td>171.0</td></tr><tr><td>01S</td><td>26.4</td><td>16.6</td><td>39.8</td><td>44.3</td><td>59.1</td><td>71.1</td><td>83.5</td><td>102.4</td><td>144.4</td></tr><tr><td>01P</td><td>53.1</td><td>43.3</td><td>15.3</td><td>19.3</td><td>33.3</td><td>45.1</td><td>57.3</td><td>76.1</td><td>118.0</td></tr><tr><td>01M</td><td>78.0</td><td>68.2</td><td>15.4</td><td>12.3</td><td>11.4</td><td>21.0</td><td>32.9</td><td>51.4</td><td>93.2</td></tr><tr><td>01Q</td><td>99.3</td><td>89.4</td><td>34.5</td><td>30.1</td><td>16.3</td><td>9.9</td><td>15.5</td><td>31.4</td><td>72.4</td></tr><tr><td>01A</td><td>133.3</td><td>123.4</td><td>68.0</td><td>63.6</td><td>48.9</td><td>37.1</td><td>25.3</td><td>11.5</td><td>39.3</td></tr><tr><td>01E</td><td>127.9</td><td>125.7</td><td>94.6</td><td>90.1</td><td>75.4</td><td>63.6</td><td>51.6</td><td>33.7</td><td>22.8</td></tr><tr><td>08G</td><td>13.8</td><td>9.6</td><td>55.8</td><td>60.3</td><td>75.0</td><td>87.0</td><td>99.3</td><td>118.2</td><td>150.2</td></tr><tr><td>08S</td><td>35.3</td><td>25.6</td><td>31.8</td><td>36.2</td><td>50.9</td><td>62.8</td><td>75.1</td><td>94.0</td><td>135.9</td></tr><tr><td>08P</td><td>61.1</td><td>51.3</td><td>11.5</td><td>13.9</td><td>26.7</td><td>38.1</td><td>50.1</td><td>68.7</td><td>110.5</td></tr><tr><td>08M</td><td>85.2</td><td>75.4</td><td>21.4</td><td>17.4</td><td>9.8</td><td>15.7</td><td>26.5</td><td>44.6</td><td>86.3</td></tr><tr><td>08Q</td><td>103.7</td><td>93.7</td><td>38.7</td><td>34.3</td><td>20.2</td><td>10.6</td><td>10.6</td><td>26.5</td><td>67.7</td></tr><tr><td>08A</td><td>135.6</td><td>125.7</td><td>70.2</td><td>65.7</td><td>51.0</td><td>39.1</td><td>27.0</td><td>10.5</td><td>36.2</td></tr><tr><td>08E</td><td>140.8</td><td>134.5</td><td>95.1</td><td>90.6</td><td>75.9</td><td>64.0</td><td>51.9</td><td>33.7</td><td>19.7</td></tr><tr><td>16G</td><td>8.2</td><td>14.1</td><td>68.7</td><td>73.2</td><td>88.1</td><td>100.1</td><td>112.5</td><td>131.4</td><td>152.0</td></tr><tr><td>16S</td><td>22.0</td><td>12.8</td><td>44.4</td><td>48.9</td><td>63.7</td><td>75.8</td><td>88.1</td><td>107.0</td><td>149.0</td></tr><tr><td>16P</td><td>46.7</td><td>36.9</td><td>20.9</td><td>25.1</td><td>39.6</td><td>51.5</td><td>63.7</td><td>82.5</td><td>124.5</td></tr><tr><td>16M</td><td>76.4</td><td>66.6</td><td>15.8</td><td>13.3</td><td>15.4</td><td>24.5</td><td>54.0</td><td>95.6</td><td>35.8</td></tr><tr><td>16Q</td><td>97.5</td><td>87.7</td><td>34.0</td><td>30.1</td><td>17.7</td><td>11.5</td><td>17.5</td><td>34.4</td><td>75.3</td></tr><tr><td>16A</td><td>129.5</td><td>119.6</td><td>64.2</td><td>59.7</td><td>44.9</td><td>33.1</td><td>21.2</td><td>7.8</td><td>42.2</td></tr><tr><td>16E</td><td>149.4</td><td>144.1</td><td>88.6</td><td>84.1</td><td>69.3</td><td>57.3</td><td>45.2</td><td>27.0</td><td>21.2</td></tr></tbody></table>
<table><caption>**Table 2.** State level RMSE result for soybean in 2007 (days)</caption><thead><tr><th></th><th>Planted</th><th>Emerged</th><th>Blooming</th><th>Set pods</th><th>Coloring</th><th>Drop leaf</th><th>Harvest</th></tr></thead><tbody><tr><td>01G</td><td>14.5</td><td>21.4</td><td>60.4</td><td>76.9</td><td>119.2</td><td>128.5</td><td>147.5</td></tr><tr><td>01S</td><td>20.5</td><td>11.0</td><td>32.7</td><td>49.2</td><td>91.6</td><td>100.9</td><td>120.0</td></tr><tr><td>01P</td><td>52.2</td><td>42.1</td><td>13.3</td><td>20.9</td><td>60.7</td><td>69.8</td><td>88.7</td></tr><tr><td>01M</td><td>76.6</td><td>66.3</td><td>25.3</td><td>11.6</td><td>35.0</td><td>44.2</td><td>63.2</td></tr><tr><td>01Q</td><td>99.7</td><td>89.4</td><td>48.4</td><td>31.9</td><td>16.9</td><td>22.6</td><td>40.7</td></tr><tr><td>01A</td><td>132.6</td><td>122.4</td><td>81.4</td><td>64.8</td><td>24.4</td><td>19.3</td><td>18.9</td></tr><tr><td>01E</td><td>173.0</td><td>148.9</td><td>107.9</td><td>91.4</td><td>49.8</td><td>41.0</td><td>28.5</td></tr><tr><td>08G</td><td>14.5</td><td>24.3</td><td>65.1</td><td>81.8</td><td>124.4</td><td>133.6</td><td>152.7</td></tr><tr><td>08S</td><td>14.5</td><td>4.9</td><td>37.6</td><td>54.2</td><td>96.9</td><td>106.2</td><td>125.3</td></tr><tr><td>08P</td><td>44.8</td><td>34.7</td><td>14.0</td><td>26.3</td><td>67.6</td><td>76.8</td><td>95.8</td></tr><tr><td>08M</td><td>71.7</td><td>61.5</td><td>21.5</td><td>10.6</td><td>40.8</td><td>49.9</td><td>68.8</td></tr><tr><td>08Q</td><td>93.7</td><td>83.5</td><td>43.3</td><td>27.8</td><td>23.7</td><td>31.1</td><td>48.9</td></tr><tr><td>08A</td><td>128.1</td><td>117.8</td><td>76.9</td><td>60.5</td><td>22.3</td><td>19.3</td><td>22.3</td></tr><tr><td>08E</td><td>129.4</td><td>120.0</td><td>103.6</td><td>87.1</td><td>45.8</td><td>37.3</td><td>28.8</td></tr><tr><td>16G</td><td>18.7</td><td>28.6</td><td>69.4</td><td>86.1</td><td>128.7</td><td>137.9</td><td>157.0</td></tr><tr><td>16S</td><td>12.0</td><td>3.0</td><td>39.6</td><td>56.3</td><td>99.0</td><td>108.3</td><td>127.4</td></tr><tr><td>16P</td><td>53.1</td><td>43.2</td><td>20.9</td><td>27.2</td><td>61.3</td><td>70.3</td><td>89.0</td></tr><tr><td>16M</td><td>73.2</td><td>63.2</td><td>24.1</td><td>15.2</td><td>41.0</td><td>49.8</td><td>68.5</td></tr><tr><td>16Q</td><td>94.9</td><td>84.7</td><td>44.3</td><td>28.6</td><td>24.0</td><td>29.6</td><td>47.4</td></tr><tr><td>16A</td><td>126.4</td><td>116.2</td><td>75.2</td><td>58.7</td><td>20.1</td><td>17.2</td><td>21.4</td></tr><tr><td>16E</td><td>128.8</td><td>141.4</td><td>100.3</td><td>83.7</td><td>41.8</td><td>32.9</td><td>22.4</td></tr></tbody></table>

For soybean, year 2006 and 2008 have the same matching pairs: planted to 'G'; emerged to 'S'; blooming to 'P'; setting pods to 'M'; leaves turning color to 'Q' and harvest to 'A'. Remote sensing derived metric 'A' refers to dropping leaf rather than harvest in 2007. The inter-annual matching error often happened at beginning or end of growing season. We believe that, besides the different cultivate conditions mentioned above, the limitation of NDVI at low vegetation density might be responsible for this experiment matching results. Another matter of concern was that, even using the same definition, stage 'harvest' is matched with metric 'E' for corn and 'A' for soybean. This indicates that soybean was harvested when NDVI descended dramatically. This is consistent with the reality: soybean is harvestable after it dries completely (soybean has less than 15% moisture). At this stage, soybean plant turns completely brown. However, when corn reaches harvestable stage, its stem is usually still green. For corn and soybean, the 2007 state level RMSE results are listed in table 1 and 2 respectively, in which normal italic RMSE values (with or without underlines) are intermediate matching pairs for each datasets, while final matching results are highlighted by underlines. The RMSE range for corn is from 7.6 (between planted and 'G', daily surface reflectance) to 25.6 (between emerged and 'S', 8-day composited surface reflectance) days. The average difference between NASS observation and remote sensing results for these three datasets are 13.5, 14.4 and 13.9 days. For soybean, the RMSE ranges from 3.0 (between emerged and 'S', 16 day composited vegetation product) to 23.9 (between turning color and 'Q', 16 day composited vegetation product) days. And the average difference between statistical phenology stages and NDVI curve based estimates for these datasets are 14.4,
14.5 and 16.5 days. MODIS Daily surface reflectance improves the remote sensing estimation results in 1-2 days when compared with MODIS MVC products.

The total errors for three datasets are displayed in figure 6. Generally the daily surface reflectance product performs better than other two datasets when being applied to crop phenology estimation studies in 2006-2008. The 8-day composited surface reflectance dataset works better than MODIS vegetation production (MOD13Q1) overall except a little irregularity for corn in 2007. Moreover, there is no obvious conclusion could be reached for crop-specific difference between corn and soybean.

<figure><img src="image_8.png" /><figcaption>Fig. 6. RMSE result for corn and soybean calculated from different datasets.</figcaption></figure>

**5. Conclusion**

Remote sensing based phenology detection is highly related to the temporal resolution of datasets. We have evaluated the effect of different datum which has different temporal resolution for two crop types: corn and soybean. USDA NASS crop progress reports have been used as ground truth data to assess the performance of different datum. Several conclusions could be drawn from above analysis: 1) Daily surface reflectance dataset could be employed to detect phenology information and had slight advantage to the composited data products, however, extra ancillary data and appropriate denoising and reconstruct techniques should be applied to eliminate the contaminated points. Double logistic function fitting method would be highly recommended; 2) The MVC process brought additional temporal error to final time series trend curve by removing contaminated points. This had a negative effect on time-sensitive remote sensing applications, for example, crop phenology estimation; 3) phenology matching is a necessary step when detecting specific stage for specific
land cover type. Even some phenological stage have the same name, for example,
emerged or harvest, they might corresponded to different segments of the final
time series curve because of differences in cropping practice. Furthermore,
phenology is a multi-factor dependent index to describe plant growth; NDVI has
its own limitations to work as the only factor, other factors, like precipitation,
temperature and soil moisture, et al. should be considered as part of the
phenology model.

Reference

Ahl, D.E., Gower, S.T., Burrows, S.N., Shabanov, N.V., Myneni, R.B. and Knyazikhin, Y. 2006, Monitoring spring canopy phenology of a deciduous broadleaf forest using MODIS, *Remote Sensing of Environment*, *104*(1),pp. 88-95

Badeck, F.W., Bondeau, A., Bottcher, K., Doktor, D., Lucht, W., Schaber, J. and Sitch, S. 2004, Responses of spring phenology to climate change, *New Phytologist*, **162**(2),pp. 295-309

Beck, P.S.A., Atzberger, C., Hgoda, K.A., Johansen, B. and Skidmore, A.K. 2006, Improved monitoring of vegetation dynamics at very high latitudes: a new method using MODIS NDVI, *Remote Sensing of Environment*, **100**(3),pp. 321-334

Beck, P.S.A., Jönsson, P., Högda, K.-A., Karlsen, S.R., Eklundh, L. and Skidmore, A.K. 2007, A ground-validated NDVI dataset for monitoring vegetation dynamics and mapping phenology in Fennoscandia and the Kola peninsula, *International Journal of Remote Sensing*, 28(19),pp. 4311 - 4330

Bradley, B.A., Jacob, R.W., Hermance, J.F. and Mustard, J.F. 2007, A curve fitting procedure to derive inter-annual phenologies from time series of noisy satellite NDVI data, *Remote Sensing of Environment*, **106**(2),pp. 137-145

Chen, J., Josson, P., Tamura, M., Gu, Z., Matsushita, B. and 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**(3-4),pp. 332-344

Duchemin, B., Hadria, R., Erraki, S., Boulet, G., Maisongrande, P., Chehbouni, A., Escadafal, R., Ezzahar, J., Hoedjes, J.C.B. and Kharrou, M.H. 2006, Monitoring wheat phenology and irrigation in Central Morocco: On the use of relationships between evapotranspiration, crops coefficients, leaf area index and remotely-sensed vegetation indices, *Agricultural Water Management*, **79**(1),pp. 1-27

Fisher, J.I., Mustard, J.F. and Vadeboncoeur, M.A. 2006, Green leaf phenology at Landsat resolution: Scaling from the field to the satellite, *Remote Sensing of Environment*, **100**(2),pp. 265-279

Garrigues, S., Allard, D. and Baret, F. 2008, Modeling temporal changes in surface spatial heterogeneity over an agricultural site, *Remote Sensing of Environment*, **112**(2),pp. 588-602

Geerken, R., Zaitchik, B. and Evans, J.P. 2005, Classifying rangeland vegetation type and coverage from NDVI time series using Fourier Filtered Cycle Similarity, *International Journal of Remote Sensing*, 26(24),pp. 5535-5554

Hird, J.N. and McDermid, G.J. 2009, Noise reduction of NDVI time series: An empirical comparison of selected techniques *Remote Sensing of Environment*, **113**(1),pp. 248-258

Hodges, T. (Ed.), 1991. *Predicting crop phenology*, pp. CRC Press

Holben, B.N. 1986, Characteristics of maximum-value composite images from temporal AVHRR data, *International Journal of Remote Sensing*, 7(11),pp. 1417-1434

Islam, A.S. and Bala, S.K. 2008, Assessment of Potato Phenological Characteristics Using MODIS-Derived NDVI and LAI Information, *GIScience & Remote Sensing, 45*(4),pp. 454-470
Jönsson, P. andEklundh, L. 2002, Seasonality extraction by function fitting to time-series of satellite sensor data, *IEEE Transactions on Geoscience and Remote Sensing*, 40(8), pp. 1824-1832

Jönsson, P. andEklundh, L. 2004, TIMESAT-A program for analyzing time-series of satellite sensor data, *Computers and Geosciences*, 30(8), pp. 833-845

Justice, C.O., Townshend, J.R.G, Holben, B.N. andTucker, C.J. 1985, Analysis of the phenology of global vegetation using meteorological satellite data, *International Journal of Remote Sensing*, 6(8), pp. 1271-1318

Lu, X., Liu, R., Liu, J. andLiang, S. 2007, Removal of Noise by Wavelet Method to Generate High Quality Temporal Data of Terrestrial MODIS Products, *PHOTOGRAMMETRIC ENGINEERING AND REMOTE SENSING*, 73(10), pp. 1129-1139

Ma, M. andVeroustraete, F. 2006, Reconstructing pathfinder AVHRR land NDVI time-series data for the Northwest of China, *Advances in Space Research*, 37(4), pp. 835-840

Markwardt, C.B. 2008, Non-linear Least Squares Fitting in IDL with MPFIT, *Astronomical Data Analysis Software and Systems XVII ASP Conference Series, XXX*

Parmesan, C. 2007, Influences of species, latitudes and methodologies on estimates of phenological response to global warming, *Global Change Biology*, 13(9), pp. 1860-1872

Reed, B.C., Brown, J.F., VanderZee, D., Loveland, T.R., Merchant, J.W. andOhlen, D.O. 1994, Measuring phenological variability from satellite imagery, *Journal of Vegetation Science*, 7(3), pp. 703-714

Roerink, G.J., Menenti, M. andVerhoef, W. 2000, Reconstructing cloud free NDVI composites using Fourier analysis of time series, *International Journal of Remote Sensing*, 21(9), pp. 1911-1917

Sakamoto, T., Yokozawa, M., Toritani, H., Shibayama, M., Ishitsuka, N. andOhno, H. 2005, A crop phenology detection method using time-series MODIS data, *Remote Sensing of Environment*, 96(3-4), pp. 366-374

Savitzky, A. andGolay, M.J.E. 1964, Smoothing and differentiation of data by simplified least squares procedures, *Analytical chemistry*, 36(8), pp. 1627-1639

Steinier, J., Termonia, Y. andDeltour, J. 1972, Comments on Smoothing and differentiation of data by simplified least square procedure, *Analytical chemistry*, 44(11), pp. 1906-1909

Stöckli, R., Rutishauser, T., Dragoni, D., O'Keefe, J., Thornton, P.E., Jolly, M., Lu, L. andDenning, A.S. 2008, Remote sensing data assimilation for a prognostic phenology mode, *Journal of Geophysical Research*, 113(11)

Sun, J., Zhou, Y., Bai, Y. and Luo, J. 2006, Nonlinear Noise Reduction of Chaotic Time Series Based on Multi-dimensional Recurrent Least Squares Support Vector Machines, *LECTURE NOTES IN COMPUTER SCIENCE*, 423-2900-906

Thayn, J.B. andPrice, K.P. 2008, Julian dates and introduced temporal error in remote sensing vegetation phenology studies, *International Journal of Remote Sensing*, 29(20), pp. 6045-6049

USGS, 2009, Land Progresses Distribution Active Archive Center, Available online at: https://lpdaac.usgs.gov/lpdaac/products/modis_products_table, (accessed 23rd, July, 2009, 2009)

Vina, A., Gitelson, A.A., Rundquist, D.C., Keydan, G., Leavitt, B. andSchepers, J. 2004, Monitoring maize (Zea mays L.) phenology with remote sensing, *Agronomy Journal*, 96(4), pp. 1139-1147

Viovy, N., Arino, O. andBelward, A.S. 1992, The Best Index Slope Extraction (BISE): A method for reducing noise in NDVI time-series, *International Journal of Remote Sensing*, 13(8), pp. 1585-1590

Wardlow, B.D., Kastens, J.H. andEgbert, S.L. 2006, Using USDA crop progress data for the evaluation of greenup onset date calculated from MODIS 250-meter data, *PHOTOGRAMMETRIC ENGINEERING AND REMOTE SENSING*, 72(11), pp. 1225-1234

White, M.A. andNemani, R.R. 2006, Real-time monitoring and short-term forecasting of land
surface phenology, *Remote Sensing of Environment*, 104(1),pp. 43-49

Zhang, X., Friedl, M.A., Schaaf, C.B., Strahler, A.H., Hodges, J.C.F., Gao, F., Reed, B.C. and Huete, A. 2003, Monitoring vegetation phenology using MODIS, *Remote Sensing of Environment*, 84(3),pp. 471-475