Validation of synthetic daily Landsat NDVI time series data generated by the improved spatial and temporal data fusion approach

Accepted Manuscript

Validation of Synthetic Daily Landsat NDVI Time Series Data Generated by the Improved Spatial and Temporal Data Fusion Approach Mingquan Wu , Wenjiang Huang , Zheng Niu , Changyao Wang , Wang Li , Bo Yu PII: DOI: Reference:

S1566-2535(17)30367-6 10.1016/j.inffus.2017.06.005 INFFUS 882

To appear in:

Information Fusion

Received date: Revised date: Accepted date:

14 September 2016 10 April 2017 11 June 2017

Please cite this article as: Mingquan Wu , Wenjiang Huang , Zheng Niu , Changyao Wang , Wang Li , Bo Yu , Validation of Synthetic Daily Landsat NDVI Time Series Data Generated by the Improved Spatial and Temporal Data Fusion Approach, Information Fusion (2017), doi: 10.1016/j.inffus.2017.06.005

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights 

Daily 30m data were generated by combining MODIS and Landsat data.



A temporal validation method for synthetic daily 30m data was proposed and tested. Vegetation phenology was monitored using synthetic daily Landsat NDVI

CR IP T



data. 

Phenology was monitored at 30 m spatial resolution while the MODIS product

A method for applying ISTDFA when land cover type changes was proposed.

AC

CE

PT

ED

M



AN US

is 500m.

1

ACCEPTED MANUSCRIPT

Validation of Synthetic Daily Landsat NDVI Time Series Data Generated by the Improved Spatial and Temporal Data Fusion Approach

Bo Yub

a

CR IP T

Mingquan Wua, c, *, Wenjiang Huangb, Zheng Niua, Changyao Wanga, Wang Lia and

The State Key Laboratory of Remote Sensing Science, Institute of Remote Sensing

Chaoyang, Beijing 100101, China b

AN US

and Digital Earth, Chinese Academy of Sciences, P.O. Box 9718, Datun Road,

Key laboratory of Digital Earth Science, Institute of Remote Sensing and Digital

M

Earth, Chinese Academy of Sciences, P.O.Box 9718, Datun Road, Chaoyang, Beijing

c

ED

100101, China

USDA Agricultural Research Service, Aerial Application Technology Research

PT

Unit, 3103 F & B Road, College Station, TX 77845, USA *Corresponding author:

CE

Mingquan Wu

AC

Tel: +86-10-6480-6258; Fax: +86-10-6480-6258 E-Mail: [email protected]

2

ACCEPTED MANUSCRIPT

Abstract: Correlation analysis has been widely used to validate the accuracy of synthetic medium-resolution data generated by spatial and temporal fusion methods. However, as the temporal resolution of Landsat data is 16 days, this method only can be used for the validation of multi-temporal Landsat data. This means that the fusion

CR IP T

accuracy of most images in a synthetic daily time series have not been validated. Furthermore, the fusion accuracy of each image in a synthetic Landsat time series is different, because there is a negative correlation with the time interval length between

AN US

the fusion date and the base image date. Therefore, there is a need for temporal validation of synthetic daily Landsat time series data. We propose a suitable validation method in this paper. The improved spatial and temporal data fusion approach

M

(ISTDFA) was applied to generate synthetic daily Landsat Normalized Difference Vegetation Index (NDVI) time series, which were then validated for both spatial and

ED

temporal dimensions using actual MODIS NDVI time series. For temporal validation,

PT

the correlation coefficients (R) between the actual and synthetic 500 m NDVI time series were calculated by pixel-by-pixel to generate imagery with an R value for each

CE

pixel. For spatial validation, R between MODIS NDVI imagery and synthetic Landsat

AC

NDVI imagery was calculated day by day to generate an R time series. This method was tested and validated in two locations (Bole and Luntai) in Xinjiang Province, China. The results show that, for temporal validation, the R values of 86.08% pixels in Bole and 94.71% pixels in Luntai are higher than 0.9, and in spatial validation, R values are higher than 0.8 on most days. Synthetic daily Landsat NDVI data was used to monitor the phenology of vegetation at a spatial resolution of 30 m successfully, 3

ACCEPTED MANUSCRIPT

while the MODIS product is limited to 500 m. Keywords: Temporal validation; spatial and temporal data fusion; remote sensing; MODIS; Landsat

CR IP T

1. Introduction Phenology, the study of seasonal natural phenomena, provides some of the most direct indicators of global climate change [1]. Long-term phenological monitoring helps monitor climate change and the responses of plant to it [2]. Satellites can obtain

AN US

periodical observations of the earth’s surface, making satellite-based remote sensing one of the main methods used for phenology monitoring. Systems such as the Advanced Very High Resolution Radiometer (AVHRR) [3], Systeme Pour

M

l’Observation de la Terre (SPOT), Vegetation (VGT) [4], and Moderate Resolution

ED

Imaging Spectroradiometer (MODIS) [5, 6] provide coarse-resolution data suitable for this purpose [7, 8]. These satellite systems can image the entire earth every 1 - 2

PT

days, creating high temporal resolution data that is highly suitable for phenology

CE

monitoring. The AVHRR system, for example, has provided daily global observations

AC

since July 1981 and constitutes an excellent long-term dataset for both global and local monitoring of phenology [9, 10]. However, the spatial resolution of these data is limited to 250 m (i.e. each pixel represents a square area of land measuring 250 x 250 m) and, hence, they cannot infer land cover types at smaller scales. Thus, they are unsuitable for high spatial resolution phenology monitoring. Medium spatial resolution sensors, such as the Thematic Mapper (TM), the Enhanced Thematic 4

ACCEPTED MANUSCRIPT

Mapper Plus (ETM+), and the Operational Land Imager (OLI) sensors on Landsat satellites provide other types of suitable data [11, 12]. Their spatial resolution is much finer (30 m), allowing them to be used for medium spatial resolution phenology monitoring [13-15]. However, their imaging frequency is much lower, at 16 days. Due

CR IP T

to cloud interference, the probability of obtaining clear Landsat images can be as low as 4% [16]. Thus, there is a shortage of valid Landsat data for the monitoring of phenology [17]. Overall, there is a shortage of high resolution satellite data for the

AN US

monitoring of phenology in both spatial and temporal dimensions.

A solution to this problem is to combine high temporal resolution data with high spatial resolution data to simulate high resolution data for both dimensions. Recently,

M

several image fusion methods were proposed to generate this kind of data. Gao et al. [18] introduced the spatial and temporal adaptive reflectance fusion model (STARFM)

ED

for blending MODIS and Landsat imagery. Several studies have applied STARFM for

PT

urban environment variable extraction, vegetated dry land ecosystem monitoring, public health studies, and generation of daily land surface temperatures [19-23].

CE

However, it has displayed insufficient accuracy in complex areas. Thus, Zhu et al. [24]

AC

enhanced STARFM (ESTARFM) for use in complex heterogeneous regions. The STARFM system has also been improved by other methods [25-27]. Wu et al. [28-30] proposed a high resolution spatial and temporal data fusion approach (STDFA) based on unmixing theory. It predicts Landsat data using the time-series mean reflectance of each land cover class, disaggregated from MODIS time-series using unmixing methods. This method was applied to the estimation of a high spatial and temporal 5

ACCEPTED MANUSCRIPT

resolution leaf area index [31] and land surface temperature data [32]. It also was used to combine 30 m resolution Huanjing satellite charge-coupled device (HJ CCD) data and 16 m Gaofen satellite no. 1 wide field-of-view camera (GF-1 WFV) data with MODIS data [16]. It was optimized by the improved spatial and temporal data fusion

CR IP T

approach (ISTDFA) to avoid the weaknesses of sensor differences and the spatial variability of STDFA [30].

Most data fusion methods (e.g. STDFA, STARFM and ESTARFM) have

AN US

validated their simulated images against actual Landsat images and quantified the differences with statistical parameters such as the correlation coefficients (R), variance, bias, and RMSE [18, 24, 29, 30]. Several studies have also used this method

M

to validate multi-temporal synthetic Landsat data time series [31-34]. However, due to the 16-day temporal resolution of actual Landsat data, it can only be used for the

ED

validation of limited synthetic data on those days with Landsat acquisitions. When

PT

synthetic daily Landsat time series were generated using high spatial and temporal data fusion methods, their accuracy was not validated due to the low temporal

CE

resolution of actual Landsat data. Furthermore, there was a negative correlation

AC

between Landsat time intervals and overall fusion accuracy [34, 35]. This means that the fusion accuracy of each image in a synthetic Landsat time series was different. The large time interval between the generation of fusion and base images resulted in large errors. So, when a synthetic daily Landsat time series is generated, there is a need for temporal validation. To address this problem, the objectives of the present study are: (1) to propose a method for the temporal validation of synthetic daily 6

ACCEPTED MANUSCRIPT

Landsat data generated by ISTDFA, and (2) to test the suitability of this data for monitoring phenology.

2. Study area

CR IP T

Two counties located in Xinjiang Province, China, were selected as the study areas (Fig. 1). The first county, Bole, is located at 44° 36′ – 44° 58′ N, 82° 05′ – 82° 40′ E,

in northwestern Xinjiang. It has a continental arid semi-desert and desert

AN US

climate. The main landforms are plains, which are largely comprised of agricultural land and desert. Main agricultural crops are cotton, winter wheat and corn. The second county, Luntai, Bayinguoleng Mongolian Autonomous Prefecture, is located at 41° 39′ – 41° 56′ N, 84° 00′ – 84° 21′ E, in the southern region of Tianshan

M

and the northern Tarim Basin. It has a warm temperate continental arid climate.

ED

Again, the main landforms are plains that support agriculture, cities and deserts. Crops include cotton, corn and winter wheat.

PT

The two countries have different spatial heterogeneity and are very

CE

representative of Xinjiang province in both land cover type and agriculture. The main

AC

difference between the counties is in the size of their crop plots, which are large in Bole and small in Luntai. Therefore, they provide suitable study sites for comparing the performance of ISTDFA in different land cover heterogeneity scenarios. This information may also aid crop monitoring in Xinjiang Province.

7

ACCEPTED MANUSCRIPT

3. Data and Pre-processing Six Landsat-5 TM data sets, seven Landsat-8 OLI data sets, and 228 500 m MODIS surface reflectance data sets (MOD09GA) acquired during clear sky conditions and of good quality were used in this study (Table 1). Six Landsat-5 TM

CR IP T

data sets and 101 500 m MODIS surface reflectance data sets acquired in 2011 were used for the Bole study area. Seven Landsat-8 OLI data sets and 117 500 m MODIS surface reflectance data sets acquired in 2013 were used for Luntai. Only the red and

AN US

near infrared (NIR) bands of Landsat and MODIS data were used for calculation of the Normalized Difference Vegetation Index (NDVI).

M

3.1 Landsat data and pre-processing

Six Landsat-5 TM data sets from the Bole area and seven Landsat-8 OLI data

ED

sets from the Luntai area were surface-reflectance products and were corrected using

PT

the 6S algorithm [36, 37]. The path/row of Landsat and horizontal/vertical tile of

CE

MODIS data are given in Table 1. All data were provided by the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center. All

AC

Landsat data were georeferenced using a second-order polynomial warping approach based on the selection of 56 ground control points (GCPs) using a 1:10,000 topographic map and the nearest neighbor resampling method with positional error within 0.56 Landsat pixels. The NDVI data, which are the most widely-used vegetation index for phenology, were then calculated from the Landsat time series. Finally, those Landsat NDVI data were used as the base image to generate synthetic 8

ACCEPTED MANUSCRIPT

daily Landsat NDVI data by combining Landsat NDVI and MODIS NDVI data with the ISTDFA method.

3.2 MODIS data and pre-processing

CR IP T

One hundred and one MOD09GA datasets from Bole and 117 MOD09GA datasets from Luntai were used. The MOD09GA dataset is a level 2G MODIS product that had been atmospherically corrected. No additional atmospheric

AN US

correction was needed. The MODIS images were reprojected from the native sinusoidal projection to a UTM-WGS84 reference system and resized to the selected study area using MODIS Reprojection Tool (MRT) software. The NDVI values were then calculated from the MODIS time-series data. The MODIS NDVI data were

M

georeferenced by a second-order polynomial warping approach based on the selection

ED

of 43 GCPs on 500 m Landsat images, with a nearest neighbor resampling method. The position error was within 0.69 500 m Landsat pixels. The 500 m Landsat images

CE

method.

PT

were resized from georeferenced Landsat images using the pixel aggregate resampling

AC

4. Method

4.1 Generating synthetic daily Landsat data using ISTDFA 4.1.1 Introduction of STDFA and ISTDFA

The STDFA model generates synthetic daily Landsat data based on the assumption that the temporal variation of each finer resolution pixel within class c is 9

ACCEPTED MANUSCRIPT

the same as the mean temporal variation of class c, which can be described as follows [29]:

r ( x, y, tk )  r ( x, y, t0 )  r (c, tk )  r (c, t0 ) ,

(1)

CR IP T

where r ( x, y, t0 ) and r ( x, y, t k ) are the reflectance of finer resolution pixels at location (x,y) at times t0 and tk; r (c, tk ) and r (c, t0 ) are the mean reflectance of land cover c at time t0 and tk. Values for r ( x, y, t0 ) can be obtained from a Landsat image at time t0. Values for r (c, tk ) , and r (c, t0 ) can be calculated from MODIS time-series data

AN US

using the ordinary least squares technique according to linear mixing theory. According to this theory, the response of each coarse spatial resolution pixel is assumed to be a linear combination of the responses of each land cover class

k

 f ( x, y, c)  r (c, t )   ( x, y, t ) c

ED

R( x, y, t ) 

M

contributing to the mixture, expressed as [38]:

(2)

c 0

k

PT

Constrained:  f c ( x , y ,c )  1, and 1  f c ( x , y , c )  0 for all,

CE

c 0

where f c ( x, y, c) is the fractional cover of class c in coarse pixel (x,y) at time t,

AC

usually assumed not to change over time, and  ( x, y, t ) is the residual error term. The fractional cover values were extracted from a high spatial resolution land cover map. In the STDFA model, the high spatial resolution land cover is mapped with the IsoData method, using two days of finer resolution data, one acquired at the beginning of the period and the other acquired at the end of the period. The ISTDFA model was proposed to avoid the weaknesses of the STDFA 10

ACCEPTED MANUSCRIPT

method, which include sensor differences and spatial variability [30]. The sensor difference was adjusted by a linear-regression method and the influence of spatial variability was removed by a weighted linear mixed model. By considering the influences of sensor difference and spatial variability, Equation (1) can be modified as

CR IP T

follows:

r( x, y ,t k )  a  ( s( i , j ,t k )  r( c, s ,t k )  s( i , j ,t 0 )  r( c, s ,t 0 ))  r( x, y ,t 0 ) , s( i , j ,t k )  r( c ,i ,t k ) / r( c , s ,t k )

(3)

(4)

AN US

where r( c , s ,tk ) is the mean surface reflectance value of class c in subset S on time tk; s( i , j ,tk ) and s( i , j ,t 0 ) are the spatial-variability adjustment factors of MODIS pixels j for target MODIS pixel i of land cover class c at times tk and t0; a is the slope

M

of the linear regression equation between the mean value disaggregated from MODIS

ED

and the actual mean value calculated from a finer Landsat image at time t0. In the ISTDFA method, the mean value was disaggregated from MODIS using a weighted

PT

Equation (2), which is as follows:

k

AC

CE

w( i , j ,t )  R( x , y ,t )  w( i , j ,t )   f c ( x , y ,c )  r( c ,t )  w( i , j ,t )  ( x , y ,t ) (5) c 0

k

Constrained:  f c ( x , y , c )  1 , and 1  f c ( x , y , c )  0 for all, c 0

where w(i,j,t) is the weight of MODIS pixel j included in subset S for disaggregation of the mean reflectance of target MODIS pixel i. The value w(i,j,t) was determined by the difference between the actual MODIS surface reflectance and the MODIS surface reflectance predicted by STDFA, and the distance between pixel j and pixel i. 11

ACCEPTED MANUSCRIPT

4.1.2 Application of ISTDFA

Several studies have demonstrated that STDFA and ISTDFA can be used to fuse Landsat NDVI data with MODIS NDVI data [16, 29-31]. Thus, ISTDFA was applied to fuse the Landsat NDVI and MODIS NDVI time series. The inputs of the ISTDFA

CR IP T

model included two datasets: a medium resolution time series dataset and a coarse resolution time series dataset. The medium resolution time series dataset must contain two medium resolution data. The first medium resolution data was used as the base

AN US

image. The second medium resolution data was used to calculate spectral similarity. The coarse resolution time series dataset must contain at least two periods of data. The first set of coarse resolution images must be acquired on the same date as those

M

for the base image, while the second set must be acquired on the date for which we are predicting. The high spatial resolution land cover data, mapped using two sets of

ED

medium resolution data, was used to extract fractional cover data for disaggregation

PT

of the time-series coarse resolution data.

CE

4.2 Temporal validation of synthetic daily Landsat NDVI time-series

AC

To evaluate the accuracy of synthetic daily Landsat NDVI time series data, a

validation method for both spatial and temporal dimensions is proposed. High temporal resolution MODIS NDVI time series data were used to validate the synthetic Landsat NDVI time series. There were three reasons to choose MODIS time series data for validation. First, there were no finer spatial resolution daily images that could be used for the validation of the synthetic Landsat time series. Second, there was a 12

ACCEPTED MANUSCRIPT

lack of daily field-measured NDVI data in the study area which could also have been used for the validation. Third, coarse-resolution data acquired by sensors like AVHRR and MODIS represented the only high temporal resolution data and had been successfully used in the validation of synthetic medium resolution time series in

CR IP T

previous studies [34, 39]. The proposed validation method for both spatial and temporal dimensions included two steps. First, we resampled the 30 m synthetic medium resolution data to

AN US

500 m resolution using the pixel aggregate resampling method, which is often used for generating coarse resolution images from finer resolution data according to linear mixing theory. This is because the spatial resolution of the synthetic medium

M

resolution data was 30 m, while the spatial resolution of MODIS data was 500 m. Thus, these two kinds of data could not be compared directly. Second, the actual and

ED

synthetic 500 m NDVI time series data were compared using the correlation analysis

PT

method for both spatial and temporal dimensions. For the validation of temporal dimensions, the R between the actual and synthetic 500 m NDVI time series were

CE

calculated by pixel-by-pixel and we attained an image with an R value for each pixel.

AC

If a pixel had a high R, there was a high correlation between the actual and synthetic 500 m NDVI time series in that pixel. If there were more pixels with high R values in the R image, the higher the accuracy of the synthetic NDVI time series. A similar validated method has been reported by Tran et al. [34, 39]. For the validation of spatial dimensions, R between an actual MODIS NDVI data and a synthetic 500 m NDVI data was calculated day by day and we attained an R time series. If there was a 13

ACCEPTED MANUSCRIPT

high R between an actual MODIS NDVI data and a synthetic 500 m NDVI data in one day, there was high accuracy in that day. The more high R values in an R time series, the higher the accuracy of the synthetic NDVI time series. Similar validation methods have been reported by Lu et al. [39].

CR IP T

4.3 Phenology extraction using synthetic daily Landsat time series 4.3.1 Phenology extraction

AN US

Generally, land surface phenology (LSP) [40, 41] recorded by satellite time series data is defined according to the seasonal growing cycles of vegetation. A typical seasonal cycle includes a “green-up” phase, a maturity phase, a senescence

M

phase and a dormant phase [42]. These four phases can be characterized by their starting dates. The onset date of greenness increase (OGI), the onset date of maximum

ED

greenness, the onset date of greenness decrease (OGD), and the onset date of

PT

minimum greenness are called greenup onset, maturity onset, senescence onset, and dormancy onset, respectively.

CE

The hybrid piecewise logistic model (HPLM) algorithm is one of the most

AC

widely used models for extracting phenology information. It has been successfully used to determine vegetation phenology using data from AVHRR [43, 44], MODIS [10, 42, 45, 46] and Landsat [17]. The HPLM method is described as follows: NDVI (t ) 

c1 1  e a b t 1

NDVI (t ) 

1

c 2  dt 1  ea

2  b2 t

 NDVI b

(6)

 NDVI b ,

(7)

14

ACCEPTED MANUSCRIPT

where t is the day of the year (DOY), a is related to vegetation growth time, b is associated with the rate of leaf development, c is the amplitude of NDVI variation, d is a vegetation stress factor, and NDVIb is the background NDVI value. Equation (6) is the HPLM method used for favorable growth conditions. Equation (7) is the HPLM

CR IP T

method used for vegetation under stress. Several studies have quantified that crop phenology can be mapped with an HPLM method, using synthetic 30 m NDVI time series generated by a spatial and temporal data fusion approach (e.g. STARFM and

AN US

STDFA) [16, 47].

Landsat and MODIS data resulting from clouds and snow cover were not used. To remove data noise, the synthetic daily 30 m Landsat NDVI data were smoothed

M

using a Savizky-Golay filter [10, 42, 45, 46]. Then, vegetation phenology was extracted in three steps using the smoothed synthetic daily 30 m Landsat NDVI data.

ED

First, the vegetation green-up and senescence phases were separated using the first

PT

derivative value of the NDVI time series curve. When the first derivative of the NDVI time series curve shifts from positive to negative, this point is assumed to be the peak

CE

of the vegetation growing cycle. The days before the peak date are considered to be

AC

the green-up phase. The days after the peak date are considered to be the greenness decrease phase. The HPLM algorithm is applied for these two stages. Because cloud-masked and snow-covered data are not used, the minimum values in these two stages are used as background NDVI values. The data from prior to the NDVIb date are set as the NDVIb values in the green-up phase. The data from after the NDVIb date are set as the NDVIb values in the greenness decrease phase. Finally, the OGI and OGD of 15

ACCEPTED MANUSCRIPT

each vegetation pixel are calculated. The date with the maximum first derivative curve in each green-up phase is set to OGI. The date with the minimum first derivative curve in each senescence phase is set to OGD.

CR IP T

4.3.2 Phenology extraction results validation

The accuracy of extracted phenology data can be evaluated with two kinds of data: via extraction from MODIS data and from field measurements of vegetation.

MODIS is a sensor onboard the Terra and Aqua satellites, both of which can

AN US

image the earth once per day. Due to its high temporal resolution, MODIS data is some of the most widely used for vegetation phenology monitoring [8, 15, 42, 48]. The MODIS Global Vegetation Phenology product (MCD12Q2) was provided by the

M

Land Processes Distributed Active Archive Center (LP DAAC) managed by the

ED

NASA Earth Science Data and Information System (ESDIS) project. Several studies have suggested that phenology information extracted from MODIS data which is

PT

validated using the field measurements of vegetation is of a high accuracy [32,

CE

35-37]. Thus, phenology information extracted from pure MODIS data can be used to evaluate the accuracy of phenology information extracted from synthetic daily

AC

Landsat data. This was conducted by the following steps. First, phenology information for the main land cover types (Bole: corn, cotton and trees; Luntai: winter wheat-summer corn, cotton and trees) in the two study areas was extracted separately from the synthetic Landsat NDVI time series and MODIS NDVI time series data according to the method described in Section 4.3.1. Then, 30 pure pixels of each land

16

ACCEPTED MANUSCRIPT

cover type in both synthetic Landsat and MODIS NDVI images, in the same locations, were selected by field survey. If the length and width of a plot were greater than 500 m, it was a pure pixel in MODIS images and several pure pixels in Landsat images. The coordinates of the four corners of this plot were measured using a GPS.

CR IP T

These pixels were used to calculate mean phenology curves for each land cover type in the synthetic Landsat and MODIS NDVI time series data. The reason for only using pure pixels in the validation was because there was a deviation in coarser

AN US

phenology and finer phenology in mixed pixels due to the influence of scaling effects [49]. Finally, R between mean phenology curves in the synthetic Landsat and mean phenology curves in MODIS data were calculated using correlation analysis for

M

comparing the phenology curves of each land cover type, which were plotted from synthetic Landsat and MODIS NDVI time series data. Higher R values indicate

ED

greater similarity between two curves and, therefore, higher accuracy of the

PT

phenology information extracted from the synthetic daily Landsat time series. Furthermore, the OGI and OGD of each land cover type were also calculated. For a

CE

given land cover type, if their OGI and OGD values extracted from the two kinds of

AC

data were very close (less than 3 weeks), then the phenology information extracted using the two kinds of data was similar [47]. In order to evaluate the correlation between the phenology parameters (OGI and OGD) of MODIS and Landsat, R between these two data was calculated using correlation analysis. The higher the R value was, the greater was the similarity between MODIS phenology and Landsat phenology. A similar method was used to evaluate the remote sensing phenology with 17

ACCEPTED MANUSCRIPT

ground observed phenology [47]. Phenology measurements made by flux tower are another widely-used method for validating remote sensing data [48, 50]. However flux tower measurements are only valid at the tower’s location. As there was no flux tower data available for the

4.4 Landsat time gap intervals assessment

CR IP T

two study areas, it could not be used for validation in this study.

AN US

It has been reported that the accuracy of STARFM will be reduced when there is a gap between the dates of the base and synthetic images [35]. This kind of gap also influences the accuracy of the ISTDFA method. The larger an interval, the more likely it is to include land cover change. In addition, environmental conditions (e.g.

M

rain) and management measures (e.g. irrigation, fertilization) are more likely to be

ED

different. Thus, the larger the interval between the base and synthetic Landsat images, the lower was the accuracy of the ISTDFA method. This influence was evaluated

PT

through three steps for the Bole study site. First, intervals were calculated between the

CE

dates of the synthetic Landsat NDVI images (June 22 - July 10, 2011), and the

AC

acquisition date of the base image (July 11, 2011). Then, synthetic NDVI time series data and 500 m MODIS time series data were calculated using the same pure vegetation pixels. Next, relative errors between the two datasets were calculated. Finally, correlation analysis determined the relationship between relative error and gap interval.

18

ACCEPTED MANUSCRIPT

5. Results 5.1 Validation of synthetic daily Landsat time series data By combining MODIS NDVI and Landsat NDVI data using ISTDFA, two

CR IP T

synthetic daily Landsat NDVI data sets were generated. There were 101 synthetic Landsat NDVI data sets for Bole and 117 for Luntai. Actual daily MODIS NDVI time series were used to validate these synthetic NDVI time series for both spatial and

AN US

temporal dimensions using the method described in Section 4.2.

5.1.1 Validation for temporal dimensions using actual MODIS NDVI time series

The R between the actual and synthetic 500 m NDVI time series was calculated

M

with pixel-by-pixel. An image with an R value for each pixel was generated. Figure 2

ED

shows the R images for Bole and Luntai. Figure 2 illustrates that the R values of 86.08% pixels in Bole and 94.71% pixels in Luntai were higher than 0.9, while only 4.95%

PT

pixels in Bole and 2.08% pixels in Luntai had R values lower than 0.7. So the

CE

synthetic Landsat NDVI time series were highly correlated with actual MODIS NDVI time series. R values higher than 0.95 were found in arable land. Low R values were

AC

found in areas where land cover types were changeable, for example swamps and areas around lakes and rivers which swell and shrink with rainfall. Figure 3a shows the NDVI time series of one randomly selected pixel in vegetation and non-vegetation. Figure 3b shows the scatter plot between MODIS NDVI and synthetic Landsat NDVI time series. It can be observed in Figure 3 that the synthetic Landsat NDVI time series were very similar to the MODIS NDVI time series. This was mainly caused because 19

ACCEPTED MANUSCRIPT

ISTDFA predicted synthetic Landsat NDVI data from the base Landsat NDVI imagery using the time-series mean NDVI of each land cover class. These were disaggregated from MODIS time-series using unmixing methods. Thus, the synthetic Landsat NDVI data could monitor the phenology exceptionally well and, hence, was

CR IP T

suitable for vegetation phenology monitoring.

5.1.2 Validation for spatial dimensions using MODIS time series data

Figure 4 shows the correlation between actual MODIS, and synthetic MODIS

AN US

NDVI time series data generated from synthetic Landsat data for Bole and Luntai. There was a high similarity, with R > 0.6 on most days. Note that lower R values were found at the beginning and the end of the year. Similar results were also found in the

M

validation results of synthetic Landsat time series generated by STARFM [39]. This

ED

was caused by the longer gap intervals between the fusion date and the base image date, and the influence of land cover changes which were analyzed in Sections 6.1 and

PT

6.2. Figure 4 illustrates the R between synthetic and actual Landsat NDVI data. It was

CE

found that R values between synthetic and actual Landsat NDVI data were higher than between synthetic and actual MODIS NDVI data, except on April 13, 2013 in Luntai.

AC

The mean R values on the days between synthetic and actual Landsat NDVI data were 0.064 and 0.078 higher than between synthetic and actual MODIS NDVI data in Bole and Luntai , respectively (except data on April 13, 2013). Thus, the synthetic Landsat NDVI data was more similar with actual Landsat NDVI than it was with actual MODIS NDVI. This was for two reasons: 1) ISTDFA predicted synthetic Landsat data

20

ACCEPTED MANUSCRIPT

from the base Landsat imagery. Thus it was more similar to the actual Landsat NDVI imagery in spatial dimensions; 2) The sensor difference and spatial resolution difference reduced the R value between synthetic Landsat NDVI imagery and actual

5.2 Phenology extraction and validation

CR IP T

MODIS NDVI imagery.

Figure 5 shows the averaged phenology curve of the main crops in Bole and

AN US

Luntai modeled with the HPLM algorithm using synthetic daily 30 m Landsat NDVI data. It shows that the phenology of different land cover types can be simulated using the HPLM algorithm. In addition, different land cover types have different phenology curves. Different crop patterns can also be distinguished. If only one crop was planted

M

(cotton in Bole and Luntai, or corn in Bole) there was one peak (Figs. 5a and 5b). If

ED

two crops were planted (winter wheat and summer corn in Luntai), there were two peaks (Fig. 5b). Because the spatial resolution of the synthetic data was 30 m, land

PT

cover types could be classified at higher spatial resolution than with MODIS data. Our

CE

previous studies showed that some land cover types (e.g. trees along the road) that

AC

cannot be classified with MODIS data can be classified with Landsat data [16]. If daily Landsat time series are generated, the phenology of land cover types can also be extracted at 30 m spatial resolution, while the resolution of the MODIS phenology product is limited to 500 m. For example, the phenology of trees along the road can be determined from synthetic Landsat NDVI data (Figs. 5a and 5b), which is hard to do with MODIS data. 21

ACCEPTED MANUSCRIPT

Table 2 and Figure 6 show the OGI and OGD values obtained using simulated phenology curves. Different land cover types have different OGIs and OGDs. Figure 7a shows the averaged phenology curves of field-measured plots of the main crops in Bole and Luntai, modeled with the HPLM algorithm using a 500 m

CR IP T

pure MODIS NDVI time series. A comparison of Figs. 5 and 7a indicates that the phenology curves modeled using synthetic daily 30 m NDVI data are very similar to those obtained using pure 500 m MODIS NDVI data. Figs. 7b and 7c show the scatter

AN US

plots for the phenology curves obtained using pure 500 m MODIS NDVI time series data and synthetic daily 30 m NDVI data. There was a high correlation between them, with R2 > 0.949. Table 3 shows the OGI and OGD values calculated using the MODIS

M

phenology curve. A comparison of Tables 2 and 3 indicates that the OGIs and OGDs of the land cover types calculated from MODIS and synthetic Landsat phenology

ED

curves are very similar. The maximum interval between those two data sets is 12 days.

PT

Figure 8 shows the result of correlation analysis between the MODIS phenology parameters and the Landsat phenology parameters. There was high correlation

CE

between these two data with an R2 value of 0.9954. The length of plots for corn in

AC

Bole and winter wheat/corn in Luntai were smaller than 100 m. Due to the low spatial resolution of MODIS NDVI data, there were no pure MODIS pixels for corn in Luntai and winter wheat/corn in Bole. Therefore, the phenology of those two land cover types were unable to be validated.

22

ACCEPTED MANUSCRIPT

6. Discussion 6.1 Influence of land-cover change Most of the spatial and temporal fusion methods (e.g. STARFM and STDFA)

CR IP T

combine MODIS and Landsat data based on the assumption that the land cover types do not change over time. As a result, these methods were not specifically designed for mapping land cover change and showed higher errors in areas where land cover types

AN US

changed. Hilker et al. found that STARFM failed to map the newly occurring forest disturbance events especially when it occurred in the sub-pixel range of MODIS data [19]. Weng et al. found that prediction errors of the Spatio-temporal Adaptive Data Fusion Algorithm for Temperature mapping (SADFAT), which had been improved

M

from STARFM, had twice the errors in changed land cover areas than it in the areas

ED

where the cover type had not changed [23]. Wu et al. found that land-cover changes in 20.25% of the pixels could produce a reduced correlation coefficient R of 0.295 of

PT

STDFA and ISTDFA in the blue band [23, 28, 30]. Furthermore, land cover type

CE

changes lead to changes in NDVI time series curves.

AC

To determine how to best use ISTDFA when land cover changes, an ISTDFA

usage experiment was conducted. Winter wheat was harvested in Bole in July 2011. This land cover change can be observed in a Landsat image acquired on July 11, 2011. In addition, the change in the NDVI time series can be monitored in the synthetic Landsat time series generated by the ISTDFA method, using a Landsat image acquired on July 11, 2011 as the base image (Figure 9). To determine whether 23

ACCEPTED MANUSCRIPT

this change can be monitored by ISTDFA, the ISTDFA method was used to generate synthetic Landsat images from May 24 - July 27, 2011, without using the Landsat image acquired on July 11, 2011. Figure 9 shows that the synthetic NDVI time series generated by the ISTDFA method without using the image acquired on July 11, 2011

CR IP T

cannot monitor the change effectively. STDFA and ISTDFA were based on the assumption that the change difference during the fusion period is homogeneous for each component [29, 30]. Zhang et al. asserted that the change rate component map

AN US

from the Landsat difference image should give a better performance than a real land cover map extracted from the base Landsat image [51]. A similar result in the unmixing-based, spatio-temporal reflectance fusion model, has been confirmed by

M

Huang and Zhang [26]. Thus, when land cover types change, a new base image must be selected to divide the fusion period into two smaller fusion periods in which the

ED

land cover types or change rate component map do not change.

PT

6.2 Influence of gap intervals

CE

Figure 10 shows the relationship between relative error and gap interval. The

AC

relationship was linear and positive - the longer the gap interval, the greater the relative error. This correlation exists because, in the period from June 22 - July 11, 2011, the winter wheat crop was harvested. Because harvest dates differed across plots, there were several changes in land cover type over the period. For dates earlier than July 11, 2011, the land cover difference is larger.

24

ACCEPTED MANUSCRIPT

6.3 Advantage and limitations Compared with MODIS data, the synthetic daily Landsat NDVI data have the advantage of being able to observe phenology at finer scales. However, there are also limitations. First, actual MODIS NDVI time series were used in both the fusion and

CR IP T

validation. In future work, it needs further validation using independent sample like field measured NDVI times series. Second, in the validation method using the MODIS time series, synthetic MODIS time series were generated by the pixel

AN US

aggregate resampling method. This is usually used to downscale high resolution images according to linear mixing theory which assumes that the reflectance of mixed pixels is a linear combination of the reflectance of each contributing pixel. This holds

M

true when the land surface is bare or covered with snow. However, a lot of studies

ED

have found that non-linear mixing is common in vegetation [52, 53]. Thus, a non-linear mixing validation method is an important direction for future studies.

PT

Furthermore, the spatial resolution of synthetic MODIS was set to 500m. However,

CE

the spatial resolution of actual MODIS was not always 500m due to the footprint and spatial response function of the MODIS pixel [51, 54]. Finally, as this paper aimed to

AC

validate synthetic daily Landsat time series data, only one year of data was used. However, interannual variation is an important aspect of the phenology of vegetation. Therefore, it will be useful for future studies to validate the use of several years of synthetic daily Landsat data for monitoring interannual variation in phenology.

25

ACCEPTED MANUSCRIPT

7. Conclusion In this study, a method for the validation of synthetic daily Landsat NDVI data was proposed and tested for two study areas (Bole and Luntai) located in Xinjiang, China. By comparing the results with MODIS NDVI time series data, we found the

CR IP T

following:

1) The ISTDFA method can be used to generate synthetic daily Landsat NDVI data. Validation using actual MODIS NDVI time series for both spatial and

AN US

temporal dimensions demonstrated that the synthetic daily Landsat NDVI data were highly correlated with actual MODIS NDVI time series data in temporal dimensions, while more similar with actual Landsat NDVI imagery in spatial

M

dimensions. In temporal dimensions, R values of 86.08% pixels in Bole and

ED

94.71% pixels in Luntai were higher than 0.9. In spatial dimensions, most R values between actual and synthetic Landsat NDVI data were higher than 0.8.

PT

2) Synthetic daily Landsat NDVI data can be used to monitor temporal variation in

CE

vegetation NDVI and phenology. Phenology at higher spatial resolutions can be monitored using this medium-resolution synthetic time series. By comparing it

AC

with the 500 m MODIS phenology product, 30 m-resolution phenology information can be extracted using a synthetic Landsat NDVI time series.

3) When land cover types change, the accuracy of ISTDFA is lower. In this situation, a new base image must be selected to divide the fusion period into two smaller periods in which land cover type does not change. 4) The accuracy of the ISTDFA method is influenced by data gap intervals (time 26

ACCEPTED MANUSCRIPT

between the dates of the base image and the generated image). The longer the interval, the lower the fusion accuracy.

CR IP T

Acknowledgments This work was supported by the Youth Innovation Promotion Association CAS (2017089), the National Natural Science Foundation of China (41301390), the National Science and Technology Major Project (2014AA06A511), and the Major Basic

Research

Development

Program

of

China

AN US

State

(2013CB733405,

2010CB950603). The funders had no role in choosing the study design; in collection, analysis, and interpretation of data; in writing the report; and in the decision to submit

AC

CE

PT

ED

M

the article for publication.

27

ACCEPTED MANUSCRIPT

References: [1] M.L. Parry, Climate change 2007-impacts, adaptation and vulnerability: Working group II contribution to the fourth assessment report of the IPCC, Cambridge University Press, 2007. [2] X. Zhang, B. Tan, Y. Yu, Interannual variations and trends in global land surface phenology derived from enhanced vegetation index during 1982–2010, International journal of biometeorology, 58 (2014) 547-564. [3] A.P. Cracknell, Advanced very high resolution radiometer AVHRR, CRC Press, 1997. [4] P. Maisongrande, B. Duchemin, G. Dedieu, VEGETATION/SPOT: an operational mission for the Earth

CR IP T

monitoring; presentation of new standard products, International Journal of Remote Sensing, 25 (2004) 9-14.

[5] V.V. Salomonson, W. Barnes, P.W. Maymon, H.E. Montgomery, H. Ostrow, MODIS: Advanced facility instrument for studies of the Earth as a system, IEEE Transactions on Geoscience and Remote Sensing, 27 (1989) 145-153.

[6] J. Verbesselt, R. Hyndman, A. Zeileis, D. Culvenor, Phenological change detection while accounting

AN US

for abrupt and gradual trends in satellite image time series, Remote Sensing of Environment, 114 (2010) 2970-2980.

[7] C. Wu, A. Gonsamo, C.M. Gough, J.M. Chen, S. Xu, Modeling growing season phenology in North American forests using seasonal mean vegetation indices from MODIS, Remote Sensing of Environment, 147 (2014) 79-88.

[8] R. Cao, J. Chen, M. Shen, Y. Tang, An improved logistic method for detecting spring vegetation phenology in grasslands from MODIS EVI time-series data, Agricultural and Forest Meteorology, 200

M

(2015) 9-20.

[9] S. Kandasamy, R. Fernandes, An approach for evaluating the impact of gaps and measurement

ED

errors on satellite land surface phenology algorithms: Application to 20year NOAA AVHRR data over Canada, Remote Sensing of Environment, 164 (2015) 114-129. [10] X. Zhang, Reconstruction of a complete global time series of daily vegetation index trajectory

PT

from long-term AVHRR data, Remote Sensing of Environment, 156 (2015) 457-472. [11] J. Luo, X. Li, R. Ma, F. Li, H. Duan, W. Hu, B. Qin, W. Huang, Applying remote sensing techniques to monitoring seasonal and interannual changes of aquatic vegetation in Taihu Lake, China, Ecological

CE

Indicators, 60 (2016) 503-513.

[12] Y. Mo, B. Momen, M.S. Kearney, Quantifying moderate resolution remote sensing phenology of Louisiana coastal marshes, Ecological Modelling, 312 (2015) 191-199.

AC

[13] Z. He, J. Du, W. Zhao, J. Yang, L. Chen, X. Zhu, X. Chang, H. Liu, Assessing temperature sensitivity of subalpine shrub phenology in semi-arid mountain regions of China, Agricultural and Forest Meteorology, 213 (2015) 42-52. [14] J.I. Fisher, J.F. Mustard, M.A. Vadeboncoeur, Green leaf phenology at Landsat resolution: Scaling from the field to the satellite, Remote sensing of environment, 100 (2006) 265-279. [15] J.I. Fisher, J.F. Mustard, Cross-scalar satellite phenology from ground, Landsat, and MODIS data, Remote Sensing of Environment, 109 (2007) 261-273. [16] M. Wu, X. Zhang, W. Huang, Z. Niu, C. Wang, W. Li, P. Hao, Reconstruction of Daily 30 m Data from HJ CCD, GF-1 WFV, Landsat, and MODIS Data for Crop Monitoring, Remote Sensing, 7 (2015) 16293-16314. [17] V. Kovalskyy, D.P. Roy, X.Y. Zhang, J. Ju, The suitability of multi-temporal web-enabled Landsat data 28

ACCEPTED MANUSCRIPT

NDVI for phenological monitoring–a comparison with flux tower and MODIS NDVI, Remote Sensing Letters, 3 (2012) 325-334. [18] F. Gao, J. Masek, M. Schwaller, F. Hall, On the blending of the Landsat and MODIS surface reflectance: Predicting daily Landsat surface reflectance, IEEE Transactions on Geoscience and Remote sensing, 44 (2006) 2207-2218. [19] T. Hilker, M.A. Wulder, N.C. Coops, J. Linke, G. McDermid, J.G. Masek, F. Gao, J.C. White, A new data fusion model for high spatial-and temporal-resolution mapping of forest disturbance based on Landsat and MODIS, Remote Sensing of Environment, 113 (2009) 1613-1627. [20] H. Liu, Q. Weng, Enhancing temporal resolution of satellite imagery for public health studies: A

CR IP T

case study of West Nile Virus outbreak in Los Angeles in 2007, Remote Sensing of Environment, 117 (2012) 57-71.

[21] M. Schmidt, R. Lucas, P. Bunting, J. Verbesselt, J. Armston, Multi-resolution time series imagery for forest disturbance and regrowth monitoring in Queensland, Australia, Remote Sensing of Environment, 158 (2015) 156-168.

[22] J. Walker, K. De Beurs, R. Wynne, F. Gao, Evaluation of Landsat and MODIS data fusion products for analysis of dryland forest phenology, Remote Sensing of Environment, 117 (2012) 381-393.

AN US

[23] Q. Weng, P. Fu, F. Gao, Generating daily land surface temperature at Landsat resolution by fusing Landsat and MODIS data, Remote Sensing of Environment, 145 (2014) 55-67.

[24] X. Zhu, J. Chen, F. Gao, X. Chen, J.G. Masek, An enhanced spatial and temporal adaptive reflectance fusion model for complex heterogeneous regions, Remote Sensing of Environment, 114 (2010) 2610-2623.

[25] D. Xie, J. Zhang, X. Zhu, Y. Pan, H. Liu, Z. Yuan, Y. Yun, An Improved STARFM with Help of an

M

Unmixing-Based Method to Generate High Spatial and Temporal Resolution Remote Sensing Data in Complex Heterogeneous Regions, Sensors, 16 (2016) 207. [26] B. Huang, H. Zhang, Spatio-temporal reflectance fusion via unmixing: accounting for both

ED

phenological and land-cover changes, International Journal of Remote Sensing, 35 (2014) 6213-6233. [27] X. Zhu, E.H. Helmer, F. Gao, D. Liu, J. Chen, M.A. Lefsky, A flexible spatiotemporal method for 165-177.

PT

fusing satellite images with different resolutions, Remote Sensing of Environment, 172 (2016) [28] M. Wu, W. Huang, Z. Niu, C. Wang, Generating daily synthetic Landsat imagery by combining

CE

Landsat and MODIS data, Sensors, 15 (2015) 24002-24025. [29] M. Wu, Z. Niu, C. Wang, C. Wu, L. Wang, Use of MODIS and Landsat time series data to generate high-resolution temporal synthetic Landsat data using a spatial and temporal reflectance fusion model,

AC

J. Appl. Remote Sens., 6 (2012) 063507-063501-063507-063513. [30] M. Wu, C. Wu, W. Huang, Z. Niu, C. Wang, W. Li, P. Hao, An improved high spatial and temporal data fusion approach for combining Landsat and MODIS data to generate daily synthetic Landsat imagery, Information Fusion, 31 (2016) 14-25. [31] M. Wu, C. Wu, W. Huang, Z. Niu, C. Wang, High-resolution Leaf Area Index estimation from synthetic Landsat data generated by a spatial and temporal data fusion model, Computers and Electronics in Agriculture, 115 (2015) 1-11. [32] M. Wu, H. Li, W. Huang, Z. Niu, C. Wang, Generating daily high spatial land surface temperatures by combining ASTER and MODIS land surface temperature products for environmental process monitoring, Environmental Science: Processes & Impacts, 17 (2015) 1396-1404. [33] I.V. Emelyanova, T.R. McVicar, T.G. Van Niel, L.T. Li, A.I.J.M. van Dijk, Assessing the accuracy of 29

ACCEPTED MANUSCRIPT

blending Landsat–MODIS surface reflectances in two landscapes with contrasting spatial and temporal dynamics: A framework for algorithm selection, Remote Sensing of Environment, 133 (2013) 193-209. [34] T.V. Tran, K.M. de Beurs, J.P. Julian, Monitoring forest disturbances in Southeast Oklahoma using Landsat and MODIS images, Int. J. Appl. Earth Obs. Geoinf., 44 (2016) 42-52. [35] D. Fu, L. Zhang, H. Chen, J. Wang, X. Sun, T. Wu, Assessing the Effect of Temporal Interval Length on the Blending of Landsat-MODIS Surface Reflectance for Different Land Cover Types in Southwestern Continental United States, ISPRS International Journal of Geo-Information, 4 (2015) 2542-2560. [36] J.G. Masek, E.F. Vermote, N.E. Saleous, R. Wolfe, A Landsat surface reflectance dataset for North America, 1990-2000, IEEE Geoscience & Remote Sensing Letters, 3 (2006) 68-72.

CR IP T

[37] E. Vermote, C. Justice, M. Claverie, B. Franch, Preliminary analysis of the performance of the

Landsat 8/OLI land surface reflectance product, Remote Sensing of Environment, 185 (2016) 46-56.

[38] J. Settle, N. Drake, Linear mixing and the estimation of ground cover proportions, International Journal of Remote Sensing, 14 (1993) 1159-1177.

[39] M. Lu, J. Chen, H. Tang, Y. Rao, P. Yang, W. Wu, Land cover change detection by integrating

object-based data blending model of Landsat and MODIS, Remote Sensing of Environment, 184 (2016) 374-386.

AN US

[40] K.M. De Beurs, G.M. Henebry, Land surface phenology, climatic variation, and institutional change: analyzing agricultural land cover change in Kazakhstan, Remote Sensing of Environment, 89 (2004) 497-509.

[41] M. Friedl, G. Henebry, B. Reed, A. Huete, Land surface phenology: a community white paper requested by NASA, in, 2006.

[42] X. Zhang, M.A. Friedl, C.B. Schaaf, A.H. Strahler, J.C. Hodges, F. Gao, B.C. Reed, A. Huete,

M

Monitoring vegetation phenology using MODIS, Remote sensing of environment, 84 (2003) 471-475. [43] Y. Julien, J. Sobrino, Global land surface phenology trends from GIMMS database, International Journal of Remote Sensing, 30 (2009) 3495-3513.

ED

[44] X. Zhang, D. Tarpley, J.T. Sullivan, Diverse responses of vegetation phenology to a warming climate, Geophysical Research Letters, 34 (2007).

[45] L. Liang, M.D. Schwartz, S. Fei, Validating satellite phenology through intensive ground

PT

observation and landscape scaling in a mixed seasonal forest, Remote Sensing of Environment, 115 (2011) 143-157.

CE

[46] X. Zhang, M.A. Friedl, C.B. Schaaf, Global vegetation phenology from Moderate Resolution Imaging Spectroradiometer (MODIS): Evaluation of global patterns and comparison with in situ measurements, Journal of Geophysical Research: Biogeosciences, 111 (2006).

AC

[47] F. Gao, M.C. Anderson, X. Zhang, Z. Yang, J.G. Alfieri, W.P. Kustas, R. Mueller, D.M. Johnson, J.H. Prueger, Toward mapping crop progress at field scales through fusion of Landsat and MODIS imagery, Remote Sensing of Environment, 188 (2017) 9-25. [48] Y. Liu, C. Wu, D. Peng, S. Xu, A. Gonsamo, R.S. Jassal, M.A. Arain, L. Lu, B. Fang, J.M. Chen, Improved modeling of land surface phenology using MODIS land surface reflectance and temperature at evergreen needleleaf forests of central North America, Remote Sensing of Environment, 176 (2016) 152-162. [49] X. Zhang, J. Wang, F. Gao, Y. Liu, C. Schaaf, M. Friedl, Y. Yu, S. Jayavelu, J. Gray, L. Liu, D. Yan, G.M. Henebry, Exploration of scaling effects on coarse resolution land surface phenology, Remote Sensing of Environment, 190 (2017) 318-330. [50] V. Rodriguez‐Galiano, J. Dash, P.M. Atkinson, Intercomparison of satellite sensor land surface 30

ACCEPTED MANUSCRIPT

phenology and ground phenology in Europe, Geophysical Research Letters, 42 (2015) 2253-2260. [51] H.K. Zhang, B. Huang, M. Zhang, K. Cao, L. Yu, A generalization of spatial and temporal fusion methods for remotely sensed surface parameters, International Journal of Remote Sensing, 36 (2015) 4411-4445. [52] Z. Malenovský, H.M. Bartholomeus, F.W. Acerbi-Junior, J.T. Schopfer, T.H. Painter, G.F. Epema, A.K. Bregt, Scaling dimensions in spectroscopy of soil and vegetation, Int. J. Appl. Earth Obs. Geoinf., 9 (2007) 137-164. [53] T.W. Ray, B.C. Murray, Nonlinear spectral mixing in desert vegetation, Remote Sensing of Environment, 55 (1996) 59-64.

CR IP T

[54] M.L. Campagnolo, Q. Sun, Y. Liu, C. Schaaf, Z. Wang, M.O. Román, Estimating the effective spatial resolution of the operational BRDF, albedo, and nadir reflectance products from MODIS and VIIRS,

AC

CE

PT

ED

M

AN US

Remote Sensing of Environment, 175 (2016) 52-64.

31

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Figure captions

32

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

33

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

34

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

35

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

36

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

37

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

38

ACCEPTED MANUSCRIPT

Table captions Table 1 Images used in this research. Landsat data sets MODIS data sets

Bole

Sensor

Acquisition datePath/rowNo.Horizontal and Vertical Tile

Landsat-5 TM

Landsat-8 OLI

146/29

04/22/2011

146/29

05/24/2011

146/29

07/11/2011

146/29

07/27/2011

146/29

09/13/2011

146/29

04/13/2013

144/31

05/15/2013

144/31

05/31/2013

144/31

08/19/2013

144/31 117

09/04/2013

144/31

10/06/2013

144/31

10/22/2013

144/31

101

h23v04

h24v04

M

AN US

Luntai

04/06/2011

CR IP T

Study Area

ED

Table 2 The OGI and OGD of each land use type extracted using synthetic daily 30 m NDVI data

Luntai

Bole

Planting patterns

Land Use

OGI (DOY)

OGD (DOY)

OGI (DOY)

OGD (DOY)

CE

PT

from Luntai and Bole

Tree

127

310

109

284

Cotton

177

295

171

276

171

235

AC

One crop

Wheat - Summer corn

Summer corn Wheat

129

181

Summer corn

204

253

39

ACCEPTED MANUSCRIPT

Table 3 The OGI and OGD for each land use type extracted from MODIS NDVI data for Luntai and Bole Land Use

Bole

OGI (DOY)

OGD (DOY)

115 175

317 292

OGI (DOY)

OGD (DOY)

163 172

277 228

AC

CE

PT

ED

M

AN US

CR IP T

Tree Cotton Summer corn

Luntai

40