Spatio-temporal reconstruction of air temperature maps and their application to estimate rice growing season heat accumulation using multi-temporal MODIS data

Spatio-temporal reconstruction of air temperature maps and their application to estimate rice growing season heat accumulation using multi-temporal MODIS data 1 2 3 first-author †E-mail:flowerpapa@hotmail.comhjf@zju.edu.cn †E-mail:flowerpapa@hotmail.comhjf@zju.edu.cn flowerpapa@hotmail.com 1 2 3 corresp ‡Corresponding Author ‡Corresponding Author †E-mail:flowerpapa@hotmail.comhjf@zju.edu.cn †E-mail:flowerpapa@hotmail.comhjf@zju.edu.cn flowerpapa@hotmail.com 1 2 3 1 2 3 4 5 Institute of Remote Sensing and Information Application, Zhejiang University, Hangzhou 310058, China Institute of Remote Sensing and Information Application, Zhejiang University, Hangzhou 310058, China Key Laboratory of Polluted Environment Remediation and Ecological Health, Ministry of Education, College of Natural Resources and Environmental Science, Zhejiang University, Hangzhou 310058, China Key Laboratory of Polluted Environment Remediation and Ecological Health, Ministry of Education, College of Natural Resources and Environmental Science, Zhejiang University, Hangzhou 310058, China Key Laboratory of Agricultural Remote Sensing and Information System of Zhejiang Province, Hangzhou 310058, China Key Laboratory of Agricultural Remote Sensing and Information System of Zhejiang Province, Hangzhou 310058, China Meteorological Information Technology Support Center of Jilin Province, Changchun 130062, China Meteorological Information Technology Support Center of Jilin Province, Changchun 130062, China Institute of Remote Sensing and Earth Sciences, Hangzhou Normal University, Hangzhou 311121, China Institute of Remote Sensing and Earth Sciences, Hangzhou Normal University, Hangzhou 311121, China The accumulation of thermal time usually represents the local heat resources to drive crop growth. Maps of temperature-based agro-meteorological indices are commonly generated by the spatial interpolation of data collected from meteorological stations with coarse geographic continuity. To solve the critical problems of estimating air temperature ( T a ) and filling in missing pixels due to cloudy and low-quality images in growing degree days (GDDs) calculation from remotely sensed data, a novel spatio-temporal algorithm for T a estimation from Terra and Aqua moderate resolution imaging spectroradiometer (MODIS) data was proposed. This is a preliminary study to calculate heat accumulation, expressed in accumulative growing degree days (AGDDs) above 10 °C, from reconstructed T a based on MODIS land surface temperature (LST) data. The verification results of maximum T a , minimum T a , GDD, and AGDD from MODIS-derived data to meteorological calculation were all satisfied with high correlations over 0.01 significant levels. Overall, MODIS-derived AGDD was slightly underestimated with almost 10% relative error. However, the feasibility of employing AGDD anomaly maps to characterize the 2001–2010 spatio-temporal variability of heat accumulation and estimating the 2011 heat accumulation distribution using only MODIS data was finally demonstrated in the current paper. Our study may supply a novel way to calculate AGDD in heat-related study concerning crop growth monitoring, agricultural climatic regionalization, and agro-meteorological disaster detection at the regional scale. MODIS land surface temperature Air temperature estimation Reconstruction Heat accumulation Rice growing season Growing degree day (GDD) Introduction 1 Heat units or thermal time, introduced by Reaumut in 1730 to explain the relationship between the growth duration and temperature, has been successfully used in a variety of agricultural sciences (McMaster and Wilhelm, 1997 1997 ). The accumulation of heat units in a region usually represents the local heat resources. As a crucial driving factor to attain a certain phenological stage, those temperature-based agro-meteorological indices expressed in growing degree days (GDDs, °C·d), photothermal units (PTU, °C·d), and heliothermal units (HTU, °C·d), are quite useful in agricultural management decisions such as crop growth monitoring, yield prediction, and the detection of phenological developments (Gordon and Bootsma, 1993 1993 ). Advances in temporal and spatial resolution and ease of availability have made remotely sensed data derived from the moderate resolution imaging spectroradiometer (MODIS) instrument onboard the Terra and Aqua satellites the most commonly used data sets in agricultural and environmental applications at the regional and global scale. Numerous studies using MODIS products have developed a series of algorithms to carry out crop identification and mapping, phenophase detection, yield estimation (Sun et al., 2009 2009 ), and meteorological disaster monitoring. Studies have shown that temperature-based agro-meteorological indices such as GDD, have vastly improved the precision of phenological event identification (McMaster and Wilhelm, 1997 1997 ; Vina et al., 2004 2004 ; Boschetti et al., 2009 2009 ) and models of yield estimation (Teal et al., 2006 2006 ; Sarma et al., 2008 2008 ). Traditionally, temperature based thermal indices are available from the point calculation of field observations of near surface air temperature ( T a , usually 1.5 m higher than the ground) from a limited number of meteorological stations. A few previous studies paid attention to the methodology for rasterizing accumulated heat units based on the Geographic Information System (GIS) interpolation (Liu et al., 2011 2011 ; Ye et al., 2011 2011 ) at the regional scale with the conventional methods such as kriging and inverse distance weighting (IDW) and modified IDW method improved by digital elevation model (DEM) (Li et al., 2010 2010 ). However, these indices are usually irregularly and sparsely distributed or unavailable from areas in complex terrain without meteorological stations; traditional geospatial interpolation approaches perform unreliably for some big and continuous scale models. Additionally, their results contain lots of uncertainty since the quality of interpolation is very much a function of the method employed (Hassan et al., 2007a 2007a ). Alternatively, the thermal images from remotely sensed satellites, usually give better spatial results (Florio et al., 2004 2004 ). MODIS provides an abundant series of land surface temperature (LST) products with different spatial and temporal resolutions from both Terra and Aqua platforms. The spatiotemporal patterns of MODIS LST-based GDD maps can preserve microhabitat-like effects better than interpolation maps from meteorological data (Neteler, 2010 2010 ). There are some previous studies that try to carry out such mapping of agro-meteorological indices by remotely sensed LST data. Hassan et al. ( 2007a 2007a ) listed several previous studies providing long-term expression of GDD on the relationship between remotely sensed normalized difference vegetation index (NDVI) time series and GDD. Instead of this GDD-VI approach, Hassan et al. ( 2007b 2007b ) focused on the methods using MODIS-based LST to calculate spatially continuous GDD maps at different spatial scales. Zorer et al. ( 2011 2011 ) calculated accumulative Winkler Index (WI) or GDDs maps derived from aggregated daily MODIS LST of Terra and Aqua satellites to assess their suitability for viticulture of the region. However, this work largely neglects the problem that the temperature in GDD calculations is designed by T a but not LST, which probably accounts for the over-estimation of remotely sensed index values relative to the true ground values. Mapping GDD from remotely sensed data, in actuality, is equal to mapping maximum and minimum T a . At present, MODIS LST has been used, with some difficulty, to estimate minimum and maximum near-surface T a with several categories of methods such as the statistic regression method, artificial neural network models, the temperature-vegetation index (TVX) method, and a method based on surface energy balance. Among those methods, the statistical method, often linking LST and other geographical parameters such as latitude, longitude, and continentality (distance from the sea), has been a classical and most accurate approach to model maximum and minimum T a at a local scale (Mostovoy et al., 2006 2006 ; Cristóbal et al., 2008 2008 ; Vancutsem et al., 2010 2010 ). Another problem in calculating GDD by optical remotely sensed data is the missing values due to clouds and other reasons. The reconstruction of T a maps is necessary to ensure the accuracy of long time accumulations. Some spatial interpolation algorithms such as volumetric spline interpolation (Neteler, 2010 2010 ; Zorer et al., 2011 2011 ) and temperature gradient-based interpolation (Ke et al., 2011 2011 ) have been used to fill in the missing pixels in LST maps. A more complex method of filling gaps can be performed on the basis of the surface energy balance equation which estimates the cloudy pixel values from neighboring clear day pixels (Lu et al., 2011 2011 ). Actually, clouds move over time, which makes merging multi-temporal satellite (e.g., Terra and Aqua platform) data with high temporal resolution an effective solution for the problem of missing pixels (Coops et al., 2007 2007 ; Crosson et al., 2012 2012 ). The objectives of the current study are as follows: (1) to provide a methodology taking advantage of the complementarity of the Terra and Aqua platform to spatio-temporally reconstruct maximum and minimum T a maps derived from MODIS LST data; (2) to evaluate the feasibility of employing accumulative growing degree day (AGDD) anomaly maps to characterize the 2001–2010 spatial and temporal variability of heat accumulation during the rice growing season in Jilin Province, China, and estimating the 2011 heat accumulation distribution using only MODIS data. Materials and methods 1 General description of study region 2 The study region, Jilin Province, is situated in the centre of Chinese northeast region between 40°52′–46°18′ Latitude N and 121°38′–131°19′ Long E (Fig. 1 1 ). It has a temperate continental monsoon climate which features four distinct seasons with a hot rainy season, and receives an annual precipitation ranging from 400 to 600 mm. The annual mean air temperature is between 2–6 °C and an annual average sunshine duration of 2 259–3 016 h. The landforms vary from western plains and central low hills to the eastern high mountains, with altitudes ranging from 4 m along the Tumen River to 2 691 m on Changbai Mountain. The district’s vertical climate zone supports many types of vegetation including meadow, crops, and needle\broad-leaved forest. Rice is the main food crop grown in Jilin Province, with a growing season of April to September. Therefore, the current study only assesses the local heat condition during this period of time. The distribution of a total of 50 meteorological stations is labeled with small white squares Available data and pre-processing 2 In this paper, the details of available data sources with remotely sensed, meteorological and other geographic data were listed in Table 1 1 . The four categories of MODIS data, namely 8-d LST/emissivity 1 000 m products (MOD/MYD11A2.005), 8-d surface reflectance 500 m product (MOD09A1.005), and yearly land cover 500 m product (MCD12Q1.005) were from the NASA Terra and Aqua satellites. The data were downloaded from NASA FTP site spanning the years of 2001–2011 (2001–2010 for model and 2011 for verification). For each year a total of 22 8-d composites of all products (from the 105th day of year (DOY 105 ) to DOY 273 ) were used. The MODIS reprojection tool software (MRT V4.0) was used to mosaic and reproject all the MODIS data from HDF format with the original sinusoidal (SIN) to a TIFF file with the Lambert conformal conic projection (SD 1 =25, SD 2 =47, CM=105, WGS84 ellipsoid). In order to operate calculations at the pixel scale, all MODIS images need to be resampled into 1000 m with the same spatial resolution of LST data by nearest neighbor (NN) resampling, which can better maintain original brightness values of MODIS pixels and have a quicker processing time during reprojection and spatial resampling (Li et al., 2009 2009 ; Neteler, 2010 2010 ). 2 Data type Original resolution Utility Source 4 MODIS products MOD11A2 8 d, 1 000 m 3 Max/min T a estimations and GDD/AGDD calculations 4 http://lpdaac.usgs.gov/lpdaac/products http://lpdaac.usgs.gov/lpdaac/products MYD11A2 8 d, 1 000 m MOD09A1 8 d, 500 m MCD12Q1 Yearly, 500 m Water mask 2 Meteorological temperature data Maximum Daily 2 Precision validation 2 Jilin Provincial Meteorological Bureau Minimum Daily 2 Auxiliary geographic data ASTER GDEM 30 m Fitting models http://asterweb.jpl.nasa.gov/gdem.asp http://asterweb.jpl.nasa.gov/gdem.asp Boundaries Map http://nfgis.nsdi.gov.cn/ http://nfgis.nsdi.gov.cn/ The values of the MOD11A2/MYD11A2 LST were calculated by averaging all the valid pixels under clear-sky during each 8-d period. This method may introduce uncertainty due to using clear sky samples to represent actual value from all sky average (Shen and Leptoukh, 2011 2011 ) However, we selected MODIS 8-d LST to calculate the multi-year heat accumulations during the whole period of the rice growing season mainly for two reasons: (1) Comparing the daily products, 8-d products can minimize the loss of data due to clouds to a large extent, which is essential to retain a more reliable source for the estimation of T a . (2) There is a high consistency between the two estimations of heat accumulation from meteorological daily T a and 8-d composited mean T a with the slopes of linear regression around 1.0 and correlation coefficients over 0.99. It is reasonably assumed that the temporal resolution change of original data from daily to 8-d composition brings little error to the final heat accumulation estimation. Therefore, we chose to characterize the multi-year variation of rice growing season heat accumulation by using less data-intensive 8-d LST products. Enhanced VI (EVI) time series were first calculated by blue, red, and near infrared bands from MOD09A1 based on the Eq. ( 1 1 ). Then an algorithm proposed by Sun et al. ( 2009 2009 ) was implemented to remove the cloudy pixels which were identified by the information in the quality assessment (QA) layer. The MODIS land cover type product contains multiple classification schemes, which were derived from observations spanning a year of Terra and Aqua data. In this study, we extracted all water pixels labeled by the plant functional type (PFT) scheme in MCD12Q1 of each year and then aggregated multi-year water pixels into one water mask. All of the MODIS-derived values were extracted from the MODIS pixel where the meteorological stations located. No aggregation technique was used in the current paper since there is a complex elevation gradient in the eastern study area. The meteorological data for the Jilin Province (daily minimum and maximum temperature) were provided by the Jilin Provincial Meteorological Bureau for the period 2001–2011. Corresponding to the DOY 105–273 of MOD/MYD11A2 products, we employed the simple average method for daily maximum and minimum air temperature data at 8-d intervals from 50 meteorological stations to compose corresponding 8-d mean maximum air temperature ( T a-max ) and 8-d mean minimum air temperature ( T a-min ) for each station over the rice growing season of all 10 years. Additionally, we use the same periods of T a-max and T a-min data to calculate 8-d mean meteorological derived GDD (Meteo-GDD) based on Eq. ( 2 2 ), 8-d accumulated Meteo-GDD (Meteo-AGDD) based on Eq. ( 3 3 ), and Meteo-AGDD anomaly based on Eq. ( 4 4 ). All of the above meteorological indices were used for model construction and accuracy verification of corresponding remotely sensed estimated data. It is worth noting that this study found that when DEM data was used, a high correlation with correlation coefficient R =0.923 was observed between 1 000 m pixel resolution DEM and actual elevation of 50 meteorological stations. Therefore, the minimal variation of elevation values after being resampled into 1 000 m from the original 30 m can be totally ignored during the following model fitting and spatial interpolation. Thermal indices and validation index 2 For this work, annual time series of EVI from DOY 105 to DOY 273 were calculated at each pixel from 8-d MOD09A1 data using the formula developed by Huete et al. ( 2002 2002 ): $${\text{EVI}} = G\frac{{{\rho _{{\text{NIR}}}} - {\rho _{{\text{red}}}}}}{{{\rho _{{\text{NIR}}}} + {C_1}{\rho _{{\text{red}}}} - {C_2}{\rho _{{\text{blue}}}} + L}}$$ , where ρ blue , ρ red , and ρ NIR are blue, red, and near infrared bands in MOD09A1, respectively. For MODIS sensor, L (=1) is the canopy background adjustment, C 1 (=6) and C 2 (=7.5) are aerosol resistance coefficients, and G (=2.5) is a gain factor. Instead of NDVI, we used EVI data as an indicator for analyzing the impacts of vegetation coverage on the relationship between near surface air temperature and LST since EVI exhibits sensitivity across high vegetation densities and atmospheric effects, and shows better performance in the downscaling study of GDD maps (Hassan et al., 2007a 2007a ). The reason why we selected an 8-d surface reflectance product to calculate EVI rather than directly employing the 16-d composited EVI product (MOD13Q1) was because the former product has the same temporal resolution of 1 000 m with 8-d LST product. There are numerous modified expressions on the basis of the original concept of GDD. In recent years, an increasing number of scientists have defined GDD as the daily thermal-time above the base temperature for plant growth (McMaster and Wilhelm, 1997 1997 ; de Beurs and Henebry, 2004 2004 ; Vina et al., 2004 2004 ; Sarma et al., 2008 2008 ), and have taken AGDDs as the index representing heat summation during the whole growing period (Yang et al., 1997 1997 ; Wang et al., 2001 2001 ; Vina et al., 2004 2004 ; de Beurs and Henebry, 2010 2010 ). In the current paper, 8-d mean GDD (°C·d) was calculated by using the canonical GDD equation from McMaster and Wilhelm ( 1997 1997 ): $$\begin{gathered}{\text{GD}}{{\text{D}}_t} = ({T_{\max {\text{ - }}t}} + {T_{\min {\text{ - }}t}})/2 - {T_{{\text{base}}}}, \\{T_{{\text{base}}}} \leqslant {T_{\max }} \leqslant {T_{{\text{upper}}}}, \\{T_{{\text{base}}}} \leqslant {T_{\min }} \leqslant {T_{{\text{upper}}}}, \\\end{gathered}$$ where T max- t is reconstructed complete T max (8-d mean maximum air temperature estimated from the MODIS daytime LST image) at DOY t , T min- t is reconstructed complete T min (8-d mean minimum air temperature estimated from the MODIS nighttime LST image) at the same period. T base is the minimum temperature threshold for crop growth, which varies among species. T upper is the maximum temperature threshold for crop growth. Usually 10 °C is the base temperature for thermophilous crops such as rice (Islam and Sikder, 2011 2011 ); therefore 10 °C is set to be the base temperature and 40 °C as the upper limit in the current study. When calculating GDD, the following adjustments were carried out: T min is set to equal to T base if T min is lower than T base , and set to T upper if T min is higher than T upper , while T max is set to equal to T base when lower than T base and equal to T upper when higher than T upper . To simplify the equation, the quantity [( T max + T min )/2] is sometimes used instead of daily mean air temperature (Morrison et al., 1989 1989 ; McMaster and Wilhelm, 1997 1997 ). However, there was no significant variation between the two interpretations of GDD equation during the period of crop growth, and the equation with T max and T min used in this paper is the most commonly used for application in crop study since it can remove impacts on the crop growth from extreme high and low air temperature (McMaster and Wilhelm, 1997 1997 ). We accumulated each 8-d interval GDDs over the growing season by simple summation when GDD exceeded the base of 10 °C by the equation from de Beurs and Henebry ( 2004 2004 ): $${\text{AGD}}{{\text{D}}_t} = {\text{AGD}}{{\text{D}}_t}_{ - {\text{8}}} + {\text{8}} \times {\text{GD}}{{\text{D}}_t}$$ , where GDD t is the 8-d AGDDs at DOY t and AGDD t is the GDDs accumulated from the beginning of the time period till DOY t+ 7 . GDD calculated by T max and T min images at DOY t only means the average GDD for 8 d during the period of DOY t and t +7 . Therefore, every pixel in each GDD image needed to be multiplied by 8 first when it was accumulated into AGDD, which can express the meaning that AGDD is the summation of daily GDD. In this paper, we accumulate GDDs by beginning with DOY 105 and ending with DOY 280 . When GDD t was greater than zero, AGDD t was added, otherwise it was ignored. Based on the formula, if we choose 10 °C for rice physiology lower limit temperature as the base temperature in GDD, the value of AGDD is equal to the value of effective accumulative temperature spending 10 °C (de Beurs and Henebry, 2004 2004 ). AGDD anomaly (AGDD a ) indicates the departure from one year AGDD to the multi-year average level. The equation of AGDD anomalies is as given below: $${\text{AGD}}{{\text{D}}_{{\text{a}}i}} = {\text{AGD}}{{\text{D}}_i} - \overline {{\text{AGDD}}}$$ , where AGDD a i is AGDD anomaly at the year of i , AGDD i is accumulated GDDs from DOY 105 to DOY 273 at the same year, $$\overline {{\text{AGDD}}}$$ is the average of 10 years AGDD from DOY 105 to DOY 273 . The consistency between all the meteorological data and corresponding remotely sensed estimated data is evaluated by the mean bias error (MBE), mean absolute error (MAE), and root mean squared error (RMSE). The MBE is given by: $${\text{MBE}} = \frac{{\sum\limits_{i = 1}^n {({x_i} - {y_i})} }}{n}$$ , where x i and y i are the measured and estimated values from meteorological station and MODIS data, respectively, and n is the number of data pairs. In the study, MBE can explain the average over or under-estimation of estimated data from the true value. The MAE and RMSE are given by: $${\text{MAE}} = \frac{{\sum\limits_{i = 1}^n {\left| {{x_i} - {y_i}} \right|} }}{n}$$ , $${\text{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{\left({{x_i} - {y_i}} \right)}^2}} }}{n}}$$ . As the name suggests, the MAE is an average of absolute errors and the RMSE is a quadratic scoring rule which measures the average magnitude of the errors. In statistics, both quantities are used to measure how close predictions or estimation are to the true observation, with RMSE giving a relatively high weight to large errors. The parameters of these two equations are the same as those for MBE. Usually, lower values are desirable. Data processing steps 2 The details of the algorithm in this study are shown in Fig. 2 2 . The five main steps were hereinafter described in sequence. The first step is the filtering of low-quality quality control (QC) values due to clouds or other processing failures of each of the LST images. The QC SDSs for LST and emissivity provide quality information on algorithm results for each pixel, which can be extracted by reading the bits in the 16-bit unsigned integer (Wan, 2007 2007 ). We determined the good quality (valid) LST pixels according to the following criteria: mandatory QA flags=pixel produced, good quality; data quality flag=good; emissivity flag=average emissivity error ≤0.01; LST error flag: 00=average LST error ≤1K. Only those good quality pixels were used for model fitting. After quality filtering and converting LST values from Kelvin to Celsius (LST C =LST K ×0.02−273.15), to link the LST Day /LST Night , EVI, and elevation with actual maximum/minimum T a from meteorological stations by multiple stepwise regression analysis, the optimal regression model with biggest adjusted R 2 value was applied to estimate fitted daytime and nighttime T a maps from Terra and Aqua platforms, respectively. The third step was the temporal merging of fitted T a maps at each 8-d period. For the year of 2003–2010, the 8-d synthetic maximum air temperatures ( T max ) were merged from two fitted daytime T a maps ( T max-Terra and T max-Aqua ) of the same 8-d period. For example, the T max map of DOY 105 is the merging result of T max-Terra of DOY 105 and T max-Aqua ) of DOY 105 . The algorithm of temporal merging is as follows: for a given pixel, if both input maps have the fitted values, the value of the synthetic T max map is the average of T max-Terra and T max-Aqua . If only one of two fitted values was present in the maps, that single value was used. Otherwise, if neither of the two maps had the fitted values, it was kept as an invalid value and was set aside for filling by spatial interpolation in the next step. Similar algorithms were carried out to produce the synthetic T min map from T min-Terra and T min-Aqua images of the same 8-d period. It is worth mentioning that temporal merging from two platforms was not implemented for the T max and T min maps of 2001 and 2002 since we could only acquire the complete available data of Terra satellite before 2003; therefore, for 2001–2002, T max and T min were directly substituted by daytime fitted T max-Terra and nighttime fitted T min-Terra , respectively. The residual invalid pixels presenting in the synthetic T max and T min images were spatially interpolated with an algorithm inspired by Neteler ( 2010 2010 ) which is based on regression analysis between surface temperature and elevation in each sliding window of the original image. If there were more than 10% valid T a pixels per sliding window and existing a regression relationship built by those 10% valid T a with corresponding DEM, the remaining invalid pixels will be filled by this regression relationship with DEM, otherwise it skipped to the next sliding window. This process was implemented until all the invalid pixels of whole image were filled. Twenty-four pixels were taken as the best sliding window radius in the current paper after the repeated test. More algorithm details can be seen in the literature (Neteler, 2010 2010 ; Ke et al., 2011 2011 ). Finally, GDD and AGDD and AGDD anomaly were calculated by using synchronous completely reconstructed T max and T min for every 8-d period in the rice growing season of each year, and then verified against corresponding ground meteorological stations’ data. Results 1 Maximum and minimum T a estimation and model validation 2 The basic mechanism for estimating air temperatures with remotely sensed LST relies on the intense heat exchange between land surface and atmosphere. In this study, air temperature was modeled as a function of EVI, elevation, and MODIS LST. Given the seasonality effect of the relationship between T a and LST, models were established for four different remotely sensed LST data across three seasons during the rice growing season. Thus, a total of 12 regression models were created for each of LST Day-Terra , LST Night-Terra , LST Day-Aqua , and LST Night-Aqua for three different seasons, corresponding to data from 50 meteorological stations which were listed in Table 2 2 . On the whole, compared with Aqua, RMSE of maximum air temperature ( T max ) modeled by Terra daytime LST was lower, the same as minimum air temperature ( T min ) modeled by Terra nighttime LST. Models were assessed based on coefficient of determinations (R2) and RMSE of fitting and test data, respectively Significance at p < 0.0001 Referring to the season category reported by Crosson et al. ( 2012 2012 ), period of DOY of 105–145 (March to May) is considered ‘spring’, DOY of 153–233 (June to August) is ‘summer’, and DOY of 241–280 (September to November) is ‘fall’. Models based on Terra LST are estimated by all the data from 10 years (2001–2010) for each season. Models based on Aqua LST are estimated by all the data from eight years (2003–2010) for each season 2 Type 2 Season 2 Model 3 T a 2 n Adjust R 2 RMSE-fit RMSE-test Fit 60% Test 40% T max-Terra Spring T max =3.473+0.568LST+12.268EVI 0.573 ** ** 2.712 2.674 1 006 670 Summer T max =15.935+0.347LST+4.469EVI−0.002Elevation 0.367 ** ** 2.010 2.040 1 706 1 138 Fall T max =4.488+0.674LST+8.702EVI 0.709 ** ** 1.610 1.780 1 133 755 T min-Terra Spring T min =0.293+0.832LST+3.615EVI−0.003Elevation 0.777 ** ** 1.663 1.570 976 651 Summer T min =3.339+0.742LST+2.592EVI−0.002Elevation 0.723 ** ** 1.452 1.554 1 568 1 045 Fall T min =0.608+0.809LST+3.464EVI−0.003Elevation 0.792 ** ** 1.593 1.685 855 570 T max-Aqua Spring T max =4.508+0.495LST+10.588EVI 0.473 ** ** 2.920 3.017 612 409 Summer T max =19.642+0.237LST+2.666EVI−0.002Elevation 0.271 ** ** 2.106 2.162 1 054 703 Fall T max =7.452+0.519LST+8.240EVI−0.002Elevation 0.575 ** ** 1.885 1.958 728 485 T min-Aqua Spring T min =4.658+0.754LST+3.034EVI−0.003Elevation 0.762 ** ** 1.985 1.898 1 054 702 Summer T min =8.086+0.613LST+1.257EVI−0.004Elevation 0.667 ** ** 1.665 1.598 850 566 Fall T min =2.184+0.757LST+5.841EVI−0.004Elevation 0.743 ** ** 1.720 1.858 772 515 Some literature also pointed out that T a estimated from Terra MODIS LST data has a smaller deviation with ground observation than the estimation derived from Aqua MODIS LST data (Zhang et al., 2011 2011 ). For 40% of the test samples, seasonal statistics show the lowest RMSE value of 1.780 and 1.958 °C in Terra and Aqua fall models for T max estimation and 1.554 and 1.598 °C in summer models for T min estimation. Spring models showed the poorest performance in estimating both T max and T min from Terra and Aqua daytime/nighttime LSTs. Reconstruction of estimated T a maps and accuracy validation 2 Temporal merging of T max and T min from Terra and Aqua T a estimation 3 There is a large amount of temporal and spatial missing pixels accounted for by clouds and low quality in LST data (Crosson et al., 2012 2012 ). Here, we introduced a spatio-temporal algorithm to fill those gaps in estimated T a images with two main strategies: firstly, temporally merged images of two acquired times and secondly, spatially interpolated in the current merged image. The object of this two-steps reconstruction for estimated T a images over the entire study area was to deal with the significant problem of missing data in determining GDDs. The performance of the temporal reconstruction was illustrated in Fig. 3 3 . For the DOY 105 of 2003, the T max-Terra image shows the pattern of T max with large invalid areas (cloudy area and an area with a large LST retrieval error) in the southeast and smaller regions in the central area, while the T max-Aqua image displays the pattern of T max with a large invalid area in the northwest and smaller regions in the east. After merging T max-Terra (Fig. 3a 3a ) and T max-Aqua (Fig. 3b 3b ) into T max (Fig. 3c 3c ), there was a small area of concentrated invalid pixels allocated in the eastern region. We assure that the temporal merging of maps from Terra and Aqua platforms can introduce higher quality T a pixels reflecting true temperature values when the cloudy condition varies quickly between two satellite’s overpass times. After this temporal merging step, the remaining pixels present in the synthetic images were spatially interpolated with an algorithm based on the relationship between surface temperature and elevation in each sliding window on the original image (Neteler, 2010 2010 ), then a series of complete T max / T min images were reconstructed. The Fig. 3d 3d showed the eventual result after two-step reconstructions. Areas of water and missing data (invalid pixels) are in white We defined the percentage of valid air temperature pixels out of the complete number of pixels covering the study area as the coverage rate. Thus, a coverage improvement rate means an increased pecentage of valid pixels from either the Terra-based T a map or the Aqua-based T a map as compared to the temporal-merged T a map. The coverage improvement rate of temporally merged T max / T min images from Terra/Aqua T max / T min and eventual average coverage rate of merged T max / T min were displayed in Table 3 3 . Because the Aqua MODIS data of the rice growing season can only be completely acquired from 2003 onwards, we only compared the coverage improvement rate for the years of 2003–2010. Generally, data from satellite images taken in the afternoon are more frequently interrupted by cloud, therefore it is not surprising that Table 3 3 reveals a remarkably better improvement for Aqua platform than for Terra platform, with annual average coverage improvement rate ranges of 5.60%–7.22% for T max-Terra , 9.56%–14.53% for T max-Aqua , 4.84%–5.86% for T min-Terra , and 13.85%–17.72% for T min-Aqua images over eight years. Eventually, the average spatial coverage during the rice growing season for eight years could reach 84.07% for merged T max image and 85.54% for merged T min image. For seasonal analysis (Fig. 4 4 ), the improvement was the most significant for the summer season with the average spatial coverage increasing from 70.37% and 64.43% of T max-Terra and T max-Aqua respectively to 78.03% of merged T max image, 69.70% and 55.98% of T min-Terra and T min-Aqua to 76.50% of the merged T min image, 86.28% and 95.90% for spring merged T max and T min images, and 91.65% and 94.13% for fall T max and T min . After carrying out the temporal merging, it is easier to fill in the remaining invalid pixels by interpolation based on temperature gradient. CIR: coverage improvement rate; CR: coverage rate. T- T max , A- T max , and M- T max are T max-Terra , T max-Aqua, and temporally merged T max images, so T- T min , A- T min , and M- T min the corresponding T min images 2 Year 4 CIR (%) 2 CR (%) T- T max A- T max T- T min A- T min M- T max M- T min 2003 7.17 13.95 4.93 16.37 86.44 84.86 2004 5.81 9.56 4.84 13.85 84.26 85.65 2005 7.22 13.54 5.05 16.15 78.98 82.05 2006 6.48 11.77 5.12 15.19 83.10 84.86 2007 5.89 11.35 5.86 14.58 86.68 87.64 2008 6.51 14.53 5.18 16.12 85.42 90.23 2009 5.60 11.79 5.18 14.52 85.35 86.42 2010 7.09 12.55 5.04 17.72 82.30 82.59 8 Mean 6.47 12.38 5.15 15.56 84.07 85.54 T- T max , A- T max , and M- T max are T max-Terra , T max-Aqua, and temporally merged T max images, so T- T min , A- T min , and M- T min the corresponding T min images Precision validation of T max and T min estimation 3 We spatially interpolated the remaining pixels by adjacent good quality (clear-day) pixels. This method may introduce uncertainty since cloudy areas may have had lower solar irradiation and T a than clear areas. However, it is an inevitable defect of optical remote sensing that the accurate estimation of some pixels will be impossible due to cloud cover. In the current paper, we only focus on the supposition whether or not the reconstructed T a can be taken as a proxy for real T a to estimate the heat accumulation. To clarify how many errors were introduced in each different reconstruction steps, the precision validations of all pixels in each 8-d T max and T min image were calculated respectively for the following two categories: good pixels and invalid pixels. Good pixels are those T a estimation pixels from LST data after QC filter and temporal merging, which can be taken as the high quality T a pixels of clear sky. The invalid pixels are the remaining pixels from cloudy sky and low quality in the product and are filled by the spatial interpolation. It could be seen in Tables 4 4 and 5 5 that statistic regression and temporal merging brought a total underestimation for maximum T a with 10-year mean MBE in 1.48 °C and quite few underestimations for minimum T a with a mean MBE in 0.29 °C. While these underestimations further increased after spatial interpolation, the 10-year mean MBEs for invalid pixels of T max and T min were 1.82 and 1.97 °C, while RMSE increased from 2.87 to 3.55 °C and 1.75 to 2.63 °C, respectively. However, the reconstruction effects of T max and T min during the whole rice growth season were satisfied with the 2001–2010’s RMSE varying from 2.77 to 3.55 °C for T max and 1.84 to 2.19 °C for T min . 2 Year 3 Good pixel 3 Invalid pixel 2 All pixel MBE (°C) RMSE (°C) n MBE (°C) RMSE (°C) n MBE (°C) RMSE (°C) 2001 1.77 2.87 743 3.32 4.73 335 2.25 3.55 2002 1.56 3.01 675 1.91 3.80 403 1.69 3.33 2003 1.68 3.03 778 2.05 3.46 300 3.46 3.16 2004 1.57 2.68 761 1.85 3.19 317 1.65 2.84 2005 1.10 2.91 681 1.71 3.06 397 1.33 2.97 2006 1.66 3.11 766 1.78 3.58 312 1.70 3.25 2007 1.46 2.59 786 1.97 3.36 292 1.60 2.82 2008 1.40 2.90 804 1.82 3.36 274 1.51 3.03 2009 1.44 3.17 792 0.15 3.54 286 1.10 3.27 2010 1.18 2.39 727 1.59 3.43 351 1.31 2.77 9 Mean 1.48 2.87 1.82 3.55 1.76 3.10 2 Year 3 Good pixel 3 Invalid pixel 2 All pixel MBE (°C) RMSE (°C) n MBE (°C) RMSE (°C) n MBE (°C) RMSE (°C) 2001 0.48 1.65 725 2.16 2.71 353 1.03 2.06 2002 0.44 2.05 729 1.58 2.28 349 0.81 2.12 2003 0.30 1.77 820 1.66 2.34 258 0.62 1.92 2004 0.05 1.85 853 1.86 2.46 225 0.42 1.99 2005 0.61 1.75 783 2.49 3.08 295 1.12 2.19 2006 0.16 1.62 834 1.86 2.56 244 0.54 1.88 2007 0.21 1.72 863 1.99 2.60 215 0.56 1.93 2008 0.19 1.71 892 1.82 2.60 186 0.47 1.89 2009 0.25 1.63 849 1.77 2.45 229 0.57 1.84 2010 0.26 1.74 803 2.52 3.17 275 0.84 2.19 9 Mean 0.29 1.75 1.97 2.62 0.70 2.00 The comparisons between time series of meteorological and MODIS-derived maximum and minimum air temperature at each 8-d period are shown in Fig. 5 5 . To reflect the overall level of the study area, the values in the curves were the averages of T a values at each 8-d interval from 50 meteorological stations and the 50 MODIS pixels corresponding to the stations’ positions. Overall, the reconstructed estimation of minimum temperature fitting was better than the maximum temperature fitting. It can be seen from the Fig. 5b 5b that 10-year fitted T min from MODIS data basically had a consistent pattern with the T a-min measured by meteorological stations, except for a tiny discrepancy in the peak of every year. Obviously, T min estimations are quite good in spring and autumn and acceptable in the summer of each year; the average MAEs of all 10 years for spring, summer and fall season are 0.944, 1.231, and 0.917 °C, respectively. However, two almost overlapping curves showed the fact that T min was capable of actually reflecting T a-min . By contrast, difference curve from meteorological T a-max minus MODIS-derived T max had a larger fluctuation over the “0” value line (Fig. 5a 5a ). Ten-year average MAE between T max and T a-max ranges from 1.697 to 2.541 °C, with 2.592 °C for ‘spring days’, 1.648 °C for ‘summer’, and 2.446 °C for ‘fall’, respectively. The time series of T a-max and T max basically raised and declined synchronously during the whole 10 years, except when there is a sudden drop of actual weather change. The values in curve are the average of 50 meteorological stations at each 8-d intervals during the whole growing season which begins with DOY 105 and ends with DOY 273 . Difference curves (max_d and min_d) are the result of meteorological T a time series minus MODIS-derived T a time series. The little breaks in the curves represent the time gaps between every year’s growing seasons The complex effects of vegetation cover on the heat exchange between LST and near surface air temperature can probably account for the discrepancies between meteorological observation and MODIS estimation. Another potential explanation is the possible situation of one or many of the meteorological stations near residential areas, in which case they would be influenced by the urban heat island. Comparing meteorological and MODIS-derived GDDs and AGDDs 2 The differences between Meteo-GDD and MODIS-GDD at each 8-d period spanning the year of 2001–2010 are shown in the Fig. 6 6 to analyse the error sources from different 8-d periods. The gap was a comprehensive reflection of the difference between reconstructed T max / T min and ground T a observation. Fig. 6 6 displayed that the difference between meteorological and MODIS-derived 8-d GDD summation was relatively large in the fall seasons from 2001–2006 and low from 2007–2010 seasons. The “level off” of MODIS-derived from 8-d GDDs in summertime can be accounted for the “plateau” in estimation of T max from LST. Difference curve is the result of Meteo-GDD time series minus MODIS-derived GDD time series. The little breaks in the curves represent the time gaps between every year’s growing seasons A comparison between meteorological and MODIS-derived AGDD with the period of entire rice growing season was carried out to display the discrepancy at each 8-d period (Fig. 7 7 ). It is clearly that there was quite a good agreement during the earlier period of the growing season in spring and summer time. The MODIS-derived AGDDs become slightly lower than the meteorological-derived AGDDs since the fall with largest error at the end of the growing season with 10-year mean difference of 157.7 °C·d. The same results are found in the statistics in each year of the 2001–2010. The reason for this negative difference is under investigation. However, our result of the annual variation of AGDDs with the DOY is similar to the finding of Neteler ( 2010 2010 ), who also demonstrated that AGDDs from reconstructed MODIS LST were systematically slightly lower than those calculations from meteorological stations in the fall time. Except for the overall comparison, the correlation coefficient between meteorological and remotely sensed AGDDs estimation of 50 stations in every 8-d periods is also analyzed (Fig. 8 8 ). The result showed that the correlation coefficients of AGDDs in the spring stage were unstable among ten years with relatively low values. This is probably because the air temperature in spring time is apt to experience a recovering warm fluctuation of around 10 °C, which is difficult to capture with satellite data. The adjustment scheme for lower and upper temperatures in GDD calculation would account for these larger gaps in spring time. It has become constant with R values more than 0.7 since DOY 145 of all 8-d periods in ten years. From the above, our result gives a powerful demonstration of the capability to employ remotely sensed AGDD to assess the accumulative heat condition in later rice growing season. Using MODIS-derived indices to evaluate the rice growing season heat accumulation 3 Comparing 2001–2010 spatial and temporal variabilities of heat accumulation 4 AGDD at DOY 273 of 8-d period in this paper (equal daily GDD accumulation from DOY 105 –DOY 280 ) represents the accumulative temperature conditions during all stages of rice growth in Jilin Province. The detailed statistics of the heat accumulation for 10 years were displayed in Table 6 6 . The accumulated heat derived from MODIS data was basically consistent with the calculated value from T a from meteorological stations. Overall, 10-year correlation coefficient R value between two categories AGDDs ranged from 0.73 to 0.86, which were all above the 0.01 significant levels. This indicates that there was an underestimation of MODIS-derived AGDDs from meteorological AGDDs with 10-year average MBE of 156.75 °C·d, amounting to almost 10% of the relative error. This result can be explained by the underestimation of T max and T min in the later growing season, especially in the fall (Fig. 5 5 ). The RMSE of 10 years spans from 143.91 to 236.81 °C·d. However, the interannual variation trend of MODIS-derived AGDDs was almost the same with the heat accumulation calculated from meteorological data. We evaluated this interannual variation by AGDD a ; the departure from one-year AGDD to the average of 10 years (Fig. 9 9 ). MODIS-derived AGDD a and meteorological AGDD a both showed the lower heat accumulation of the years of 2002, 2003, 2005, 2006, and 2009 with negative values, while better heat conditions in 2001, 2007, 2008, and 2010 with positive values. Compared to AGDDs, the correlation between two AGDD a was less high (only seven years over the significant level) because the multi-year universal models adopted to estimate T max / T min were not ideal for all 10 years. We chose these universal models because they are better at comparing the interannual variation by removing the uncertainty caused by single models, and forecasting the accumulative heat resource in the later year without using meteorological data. Significance at p < 0.01 2 Year 5 AGDD (°C·d) 2 R MODIS- Meteo- MBE MAE RMSE 2001 1 392.83 1 606.07 212.6 215.79 236.81 0.83 ** ** 2002 1 316.88 1 475.8 157.46 163.6 185.69 0.83 ** ** 2003 1 324.22 1 486.85 161.91 163.21 183.38 0.86 ** ** 2004 1 378.04 1 495.24 116.66 125.51 148.17 0.83 ** ** 2005 1 342.71 1 473.63 129.74 140.34 166.43 0.74 ** ** 2006 1 347.43 1 444.71 95.8 117.59 143.91 0.73 ** ** 2007 1 430.78 1 619.27 187.92 187.92 212.79 0.84 ** ** 2008 1 392.51 1 560.83 166.38 168.85 193.76 0.75 ** ** 2009 1 332.6 1 474.02 140.7 145.74 170.28 0.80 ** ** 2010 1 402.99 1 601.72 198.28 199.18 222.14 0.77 ** ** 7 Mean 1 366.10 1 523.81 156.75 162.77 186.34 0.80 ** ** The dotted line stands for the “0” value line An illustration of AGDD anomaly maps at 1 km spatial resolution was produced for the 2001–2010 period in Jilin Province and is displayed in Fig. 10 10 . Based on this we could characterize the spatial and temporal variability in heat accumulation during the whole rice growing season (AGDD at DOY 273 ) at a large spatial scale year to year. The positive AGDD anomalies with warm colours (red, orange) indicated a greater heat accumulation in a specified year than the normal level based on the average of the recent 10 years. This is usually an indicator above average crop development. Vice versa, the negative AGDD anomalies represented in cool colours (light and dark blue) mean lower heat accumulation and probably indicate areas where below-average rice growth delayed harvest. A similar method of employing AGDD anomalies for corn growth monitoring was found by the North Dakota Agricultural Weather Network. Our results revealed a better overall heat condition for rice cropping in the years of 2001, 2007, and 2010 versus lower heat accumulation in the 2002, 2003, and 2009 years. This characterization was almost the same as the reports of annual analysis of Jilin Agro-meteorological conditions and their impacts on crop production which has been released by the Jilin Provincial Meteorological Bureau. The positive AGDD anomalies with warm colors (red, orange) indicate a greater heat accumulation of specified year than the normal level based on the average of the recent 10 years. Vice versa, the negative AGDD anomalies represented in cool color (light and dark blue) mean lower heat accumulation where there is probably a growth delay for rice reaching harvest Estimating GDDs and AGDDs by only using MODIS data 4 Only MODIS data for 2011 were chosen to evaluate the estimating ability of MODIS indices (GDD, AGDD) to reflect the true ground heat condition without T a data being taken to be modeled. A total of 44 maximum and minimum T a maps of each 8-d period of 2011 were calculated by employing the 12 T a -LST regression models built by 2001–2010 data sets (Table 2 2 ) and then spatio-temporal reconstruction. The average discrepancies between MODIS reconstructed T a , GDDs, AGDDs and corresponding meteorological values were compared in the forms of MBE, MAE, RMSE, and R (Table 7 7 ). The statistical results including all the data pairs from the 50 meteorological stations spanning the 22 8-d periods in rice growing season. Seasonal comparisons between MODIS-derived values and the meteorological estmation for T max , T min , 8-d summated GDDs and every 8-d intervals AGDDs of 2011 are listed in Table 7 7 . Almost all of the correlation of two categories of data for three seasons and four indices are over the 0.01 significant levels except for the T max in summer time ( R =0.32). The annual average MAEs of T min and T a-min in 2011 were 2.386 and 1.829 °C, respectively, slightly larger than those results of 2001–2010. Point-samples of MODIS-derived GDD and AGDD values were compared against the meteorological values; high R values were found with 0.91 for GDD and 0.99 for AGDD. 2 Season 4 T max (°C) 4 T min (°C) 4 GDD (°C·d) 4 AGDD (°C·d) MBE MAE RMSE R MBE MAE RMSE R MBE MAE RMSE R MBE MAE RMSE R Spring −0.05 2.65 3.37 0.67 0.55 1.55 2.11 0.88 0.82 11.67 15.36 0.73 −15.28 24.69 31.07 0.92 Summer 1.51 2.26 2.79 0.32 1.31 2.01 2.48 0.70 11.31 15.86 19.24 0.60 42.06 67.73 92.45 0.98 Fall 1.59 2.35 2.85 0.73 0.68 1.78 2.14 0.89 8.74 11.78 15.39 0.74 169.22 169.60 194.58 0.82 Year 1.10 2.39 2.96 0.86 0.96 1.83 2.31 0.94 7.86 13.79 17.33 0.91 55.32 79.15 114.11 0.99 The scatter plots of all data pairs in 50 meteoro-logical stations of 22 8-d period during the growing season (total 1 100 pairs) of 8-d GDDs and AGDD were displayed in Figs. 11a and 11b 11a and 11b , respectively. Compared with the correlation coefficient of 8-d GDDs with R 2 of 0.82, there was a much higher correlation between two AGDDs data sets with R 2 of 0.98. The discrepancies derived from remotely sensed data in real-time monitoring could be reduced to some extend in the long process of accumulation. In conclusion, methodology of estimating remotely sensed accumulative heat time from MODIS data proposed by this current paper is feasible to reflect the ground true heat condition. This would be an advanced replacement of raster maps of such agro-meteorological indices interpolated from meteorological data with GIS. Except for the quantitative analysis, the spatial pattern of AGDD of the rice growing season in 2011 was validated by the GIS interpolation of meteorological AGDD from 50 ground stations by employing the algorithm proposed by Li et al. ( 2010 2010 ) (Fig. 12 12 ). Nearly the same spatial distribution of heat accumulation of the rice growing season was shown in the two images. Compared with the traditional GIS interpolation method, remotely sensed accumulative heat estimation exhibited better performance in characterizing microclimates and spatial continuity. This was particularly evident in the northwestern and southeastern parts of the study region. (a) MODIS-derived; (b) GIS interpolation based on IDW method improved by DEM Discussion and conclusions 1 Estimation of air temperature and filling in missing values due to clouds and low quality pixels are the main critical issues in GDD mapping by satellite data. To solve the problems above, we first adopt the classical multi-regression model to estimate the maximum and minimum T a from LST images, and then implement a novel spatio-temporal algorithm that combines the temporal merging MODIS LST-based T a maps at four different acquired times from the Terra and Aqua platform with spatial interpolation based on the temperature gradient with elevation change. The reconstructed T a images were then used to map agro-meteorological indices expressed in GDDs and AGDDs. When compared with the traditional estimation of ground agro-meteorological indices from the limited number of meteorological stations, remotely sensed accumulative heat estimation exhibits superior performance in microclimates. In this paper, MODIS 8-d composite LST products (MOD11A2) were demonstrated to be a feasible approach to mapping the heat accumulation for the crop growing season at the regional scale. The statistical analysis reveals high correlations, over the 0.01 significant level of agro-meteorological indices (GDDs and AGDDs), between estimates from meteorological stations and MODIS LST data. However, the discrepancy between MODIS-derived AGDDs and ground-derived true values becomes a little large in the later stages of the growing season. GDD is designed to be calculated by daily maximum and minimum T a and AGDD is the summation of daily GDD. Therefore, 8-d accumulated GDD calculated with 8-d average GDD would lead to more errors in GDD and AGDD estimation comparing with the summation of eight actual daily GDDs. A crucial avenue for further research concerns the reconstruction of daily LST-based T a estimation for cloudy days, which to some extent can improve the accuracy of remotely sensed GDDs to satisfy the quantified needs of GDD-related remote sensing studies for agriculture and climatology. Aside from the limited success in applying AGDD and AGDD anomaly map to evaluating the spatio-temporal pattern of heat condition for a whole crop growing season, forecasting future MODIS-derived GDDs and AGDDs with the LST- T a models constructed by data from previous years was demonstrated to be feasible. Furthermore, it is suggested to employ our methodology to monitor crop growth and identify crop phenological stages by estimating GDDs reaching a certain threshold. Our study may supply a novel approach to heat-related study concerning crop growth monitoring, agricultural climatic regionalization, and agro-meteorological disaster detection at the regional scale. References M Boschetti , D Stroppiana , PA Brivio , S Bocchi . Multi-year monitoring of rice crop phenology through time series analysis of MODIS images . Int J Remote Sens , 2009 . 30 ( 18 ): 4643 - 4662 . DOI: 10.1080/01431160802632249 http://doi.org/10.1080/01431160802632249 . Multi-year monitoring of rice crop phenology through time series analysis of MODIS images N Coops , D Duro , M Wulder , T Han . Estimating afternoon MODIS land surface temperatures (LST) based on morning MODIS overpass, location and elevation information . Int J Remote Sens , 2007 . 28 ( 10 ): 2391 - 2396 . DOI: 10.1080/01431160701294653 http://doi.org/10.1080/01431160701294653 . Estimating afternoon MODIS land surface temperatures (LST) based on morning MODIS overpass, location and elevation information J Cristóbal , M Ninyerola , X Pons . Modeling air temperature through a combination of remote sensing and GIS data . J Geophys Res , 2008 . 113 ( D13 ): D13106 DOI: 10.1029/2007JD009318 http://doi.org/10.1029/2007JD009318 . Modeling air temperature through a combination of remote sensing and GIS data WL Crosson , MZ Al-Hamdan , SNJ Hemmings , GM Wade . A daily merged MODIS Aqua-Terra land surface temperature data set for the conterminous United States . Remote Sens Environ , 2012 . 119 315 - 324 . DOI: 10.1016/j.rse.2011.12.019 http://doi.org/10.1016/j.rse.2011.12.019 . A daily merged MODIS Aqua-Terra land surface temperature data set for the conterminous United States KM de Beurs , GM Henebry . Land surface phenology, climatic variation, and institutional change: analyzing agricultural land cover change in Kazakhstan . Remote Sens Environ , 2004 . 89 ( 4 ): 497 - 509 . DOI: 10.1016/j.rse.2003.11.006 http://doi.org/10.1016/j.rse.2003.11.006 . Land surface phenology, climatic variation, and institutional change: analyzing agricultural land cover change in Kazakhstan KM de Beurs , GM Henebry . Spatio-Temporal Statistical Methods for Modelling Land Surface Phenology . Phenological Research: Methods for Environmental and Climate Change Analysis , 2010 . : London Springer , 177 - 208 . DOI: 10.1007/978-90-481-3335-2_9 http://doi.org/10.1007/978-90-481-3335-2_9 . Spatio-Temporal Statistical Methods for Modelling Land Surface Phenology E Florio , S Lele , Y Chang , R Sterner , G Glass . Integrating AVHRR satellite data and NOAA ground observations to predict surface air temperature: a statistical approach . Int J Remote Sens , 2004 . 25 ( 15 ): 2979 - 2994 . DOI: 10.1080/01431160310001624593 http://doi.org/10.1080/01431160310001624593 . Integrating AVHRR satellite data and NOAA ground observations to predict surface air temperature: a statistical approach R Gordon , A Bootsma . Analyses of growing degree-days for agriculture in Atlantic Canada . Clim Res , 1993 . 3 169 - 176 . DOI: 10.3354/cr003169 http://doi.org/10.3354/cr003169 . Analyses of growing degree-days for agriculture in Atlantic Canada QK Hassan , C Bourque , FR Meng , W Richards . Spatial mapping of growing degree days: an application of MODIS-based surface temperatures and enhanced vegetation index . J Appl Remote Sens , 2007 . 1 ( 1 ): 013511 DOI: 10.1117/1.2740040 http://doi.org/10.1117/1.2740040 . Spatial mapping of growing degree days: an application of MODIS-based surface temperatures and enhanced vegetation index QK Hassan , CPA Bourque , FR Meng . Application of Landsat-7 ETM + and MODIS products in mapping seasonal accumulation of growing degree days at an enhanced resolution . J Appl Remote Sens , 2007 . 1 ( 1 ): 013539 DOI: 10.1117/1.2800284 http://doi.org/10.1117/1.2800284 . Application of Landsat-7 ETM+ and MODIS products in mapping seasonal accumulation of growing degree days at an enhanced resolution A Huete , K Didan , T Miura , EP Rodriguez , X Gao , LG Ferreira . Overview of the radiometric and biophysical performance of the MODIS vegetation indices . Remote Sens Environ , 2002 . 83 ( 1-2 ): 195 - 213 . DOI: 10.1016/S0034-4257(02)00096-2 http://doi.org/10.1016/S0034-4257(02)00096-2 . Overview of the radiometric and biophysical performance of the MODIS vegetation indices M Islam , S Sikder . Phenology and degree days of rice cultivars under organic culture . Bangladesh J Bot , 2011 . 40 ( 2 ): 149 - 153 . . Phenology and degree days of rice cultivars under organic culture L Ke , Z Wang , C Song , Z Lu . Reconstruction of MODIS LST time series and comparison with land surface temperature (T) among observation stations in the Northeast Qinghai-Tibet plateau . Prog Geogr , 2011 . 30 ( 7 ): 819 - 826 . (in Chinese) . Reconstruction of MODIS LST time series and comparison with land surface temperature (T) among observation stations in the Northeast Qinghai-Tibet plateau F Li , C Wang , J Zhao , J Zheng . The spatialization of multi-year average accumulated temperature in China . J Nat Resour , 2010 . 25 ( 5 ): 778 - 783 . (in Chinese) . The spatialization of multi-year average accumulated temperature in China J Li , X Wang , W Ma , H Zhang . Remote sensing evaluation of urban heat island and its spatial pattern of the Shanghai metropolitan area, China . Ecol Complex , 2009 . 6 ( 4 ): 413 - 420 . DOI: 10.1016/j.ecocom.2009.02.002 http://doi.org/10.1016/j.ecocom.2009.02.002 . Remote sensing evaluation of urban heat island and its spatial pattern of the Shanghai metropolitan area, China J Liu , Z Wang , M Chen , K Shi , T Liang , B Zhang , F Zhou , J Chang . Temporal and spatial variation characteristics of effective accumulated temperature in Northeast China during 1951-2007 . J Anhui Agric Sci , 2011 . 39 ( 25 ): 15655 - 15656 . 15680 (in Chinese) . Temporal and spatial variation characteristics of effective accumulated temperature in Northeast China during 1951-2007 L Lu , V Venus , A Skidmore , T Wang , G Luo . Estimating land-surface temperature under clouds using MSG/SEVIRI observations . Int J Appl Earth Obs Geoinform , 2011 . 13 ( 2 ): 265 - 276 . DOI: 10.1016/j.jag.2010.12.007 http://doi.org/10.1016/j.jag.2010.12.007 . Estimating land-surface temperature under clouds using MSG/SEVIRI observations GS McMaster , W Wilhelm . Growing degree-days: one equation, two interpretations . Agric For Meteorol , 1997 . 87 ( 4 ): 291 - 300 . DOI: 10.1016/S0168-1923(97)00027-0 http://doi.org/10.1016/S0168-1923(97)00027-0 . Growing degree-days: one equation, two interpretations M Morrison , P McVetty , C Shaykewich . The determination and verification of a baseline temperature for the growth of Westar summer rape . Can J Plant Sci , 1989 . 69 ( 2 ): 455 - 464 . DOI: 10.4141/cjps89-057 http://doi.org/10.4141/cjps89-057 . The determination and verification of a baseline temperature for the growth of Westar summer rape GV Mostovoy , RL King , KR Reddy , VG Kakani , MG Filippova . Statistical estimation of daily maximum and minimum air temperatures from MODIS LST data over the state of Mississippi . GISci Remote Sens , 2006 . 43 ( 1 ): 78 - 110 . DOI: 10.2747/1548-1603.43.1.78 http://doi.org/10.2747/1548-1603.43.1.78 . Statistical estimation of daily maximum and minimum air temperatures from MODIS LST data over the state of Mississippi M Neteler . Estimating daily land surface temperatures in mountainous environments by reconstructed MODIS LST data . Remote Sens , 2010 . 2 ( 1 ): 333 - 351 . DOI: 10.3390/rs1020333 http://doi.org/10.3390/rs1020333 . Estimating daily land surface temperatures in mountainous environments by reconstructed MODIS LST data A Sarma , TVL Kumar , K Koteswararao . Development of an agroclimatic model for the estimation of rice yield . J Ind Geophys Union , 2008 . 12 ( 2 ): 89 - 96 . . Development of an agroclimatic model for the estimation of rice yield S Shen , GG Leptoukh . Estimation of surface air temperature over central and eastern Eurasia from MODIS land surface temperature . Environ Res Lett , 2011 . 6 ( 4 ): 045206 DOI: 10.1088/1748-9326/6/4/045206 http://doi.org/10.1088/1748-9326/6/4/045206 . Estimation of surface air temperature over central and eastern Eurasia from MODIS land surface temperature H Sun , J Huang , AR Huete , D Peng , F Zhang . Mapping paddy rice with multi-date moderate-resolution imaging spectroradiometer (MODIS) data in China . J Zhejiang Univ-Sci A (Appl Phys & Eng) , 2009 . 10 ( 10 ): 1509 - 1522 . DOI: 10.1631/jzus.A0820536 http://doi.org/10.1631/jzus.A0820536 . Mapping paddy rice with multi-date moderate-resolution imaging spectroradiometer (MODIS) data in China R Teal , B Girma , K Freeman , K Arnall , DO Walsh , W Raun . In-season prediction of corn grain yield potential using normalized difference vegetation index . Agron J , 2006 . 98 ( 6 ): 1488 DOI: 10.2134/agronj2006.0103 http://doi.org/10.2134/agronj2006.0103 . In-season prediction of corn grain yield potential using normalized difference vegetation index C Vancutsem , P Ceccato , T Dinku , SJ Connor . Evaluation of MODIS land surface temperature data to estimate air temperature in different ecosystems over Africa . Remote Sens Environ , 2010 . 114 ( 2 ): 449 - 465 . DOI: 10.1016/j.rse.2009.10.002 http://doi.org/10.1016/j.rse.2009.10.002 . Evaluation of MODIS land surface temperature data to estimate air temperature in different ecosystems over Africa A Vina , AA Gitelson , DC Rundquist , GP Keydan , B Leavitt , J Schepers . Monitoring maize ( Zea mays L.) phenology with remote sensing . Agron J , 2004 . 96 ( 4 ): 1139 - 1147 . DOI: 10.2134/agronj2004.1139 http://doi.org/10.2134/agronj2004.1139 . Monitoring maize (Zea mays L.) phenology with remote sensing Z Wan . Collection-5 MODIS Land Surface Temperature Products Users’ Guide , 2007 . : Santa Barbara ICESS, University of California , . Collection-5 MODIS Land Surface Temperature Products Users’ Guide J Wang , K Price , P Rich . Spatial patterns of NDVI in response to precipitation and temperature in the central Great Plains . J Appl Remote Sens , 2001 . 22 ( 18 ): 3827 - 3844 . DOI: 10.1080/01431160010007033 http://doi.org/10.1080/01431160010007033 . Spatial patterns of NDVI in response to precipitation and temperature in the central Great Plains W Yang , L Yang , J Merchant . An assessment of AVHRR/NDVI-ecoclimatological relations in Nebraska, USA . J Appl Remote Sens , 1997 . 18 ( 10 ): 2161 - 2180 . DOI: 10.1080/014311697217819 http://doi.org/10.1080/014311697217819 . An assessment of AVHRR/NDVI-ecoclimatological relations in Nebraska, USA Q Ye , X Yang , Y Li , S Dai , J Xiao . Changes of China agricultural climate resources under the background of climate change. VIII. Change characteristics of heat resources during the growth period of double cropping rice in Jiangxi Province . Chin J Appl Ecol , 2011 . 22 ( 8 ): 2021 - 2030 . (in Chinese) . Changes of China agricultural climate resources under the background of climate change. VIII. Change characteristics of heat resources during the growth period of double cropping rice in Jiangxi Province W Zhang , Y Huang , Y Yu , W Sun . Empirical models for estimating daily maximum, minimum and mean air temperatures with MODIS land surface temperatures . Int J Remote Sens , 2011 . 32 ( 24 ): 9415 - 9440 . DOI: 10.1080/01431161.2011.560622 http://doi.org/10.1080/01431161.2011.560622 . Empirical models for estimating daily maximum, minimum and mean air temperatures with MODIS land surface temperatures R Zorer , D Rocchini , L Delucchi , 等 . Use of Multi-annual MODIS Land Surface Temperature Data for the Characterization of the Heat Requirements for Grapevine Varieties . 2011 . 2011 6th International Workshop on the Analysis of Multi-Temporal Remote Sensing Images, Multi-Temp 2011-Proceedings . Trento, Italy IEEE Computer Society , 225 - 228 . DOI: 10.1109/Multi-Temp.2011.6005089 http://doi.org/10.1109/Multi-Temp.2011.6005089 . Use of Multi-annual MODIS Land Surface Temperature Data for the Characterization of the Heat Requirements for Grapevine Varieties 10.1631/jzus.B1200169 Project supported by the National Key Technology R&D Program of China (No. 2012BAH29B02) and the PhD Programs Foundation of Ministry of Education of China (No. 200100101110035) 2012-06-14 2013-01-13 2012-10-15 2022-01-05 Copyright © Zhejiang University and Springer-Verlag Berlin Heidelberg 2013 2013 Journal of Zhejiang University-SCIENCE B (Biomedicine and Biotechnology) 02 14