RICE CROP MONITORING AND YIELD ASSESSMENT WITH MODIS 250 m GRIDDED VEGETATION PRODUCT : A case study in Sa Kaeo Province , Thailand

Billions of people in the world depend on rice as a staple food and as an income-generating crop. Asia is the leader in rice cultivation and it is necessary to maintain an up-to-date rice-related database to ensure food security as well as economic development. This study investigates general applicability of high temporal resolution Moderate Resolution Imaging Spectroradiometer (MODIS) 250m gridded vegetation product for monitoring rice crop growth, mapping rice crop acreage and analyzing crop yield, at the province-level. The MODIS 250m Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) time series data, field data and crop calendar information were utilized in this research in Sa Kaeo Province, Thailand. The following methodology was used: (1) data pre-processing and rice plant growth analysis using Vegetation Indices (VI) (2) extraction of rice acreage and start-of-season dates from VI time series data (3) accuracy assessment, and (4) yield analysis with MODIS VI. The results show a direct relationship between rice plant height and MODIS VI. The crop calendar information and the smoothed NDVI time series with Whittaker Smoother gave high rice acreage estimation (with 86% area accuracy and 75% classification accuracy). Point level yield analysis showed that the MODIS EVI is highly correlated with rice yield and yield prediction using maximum EVI in the rice cycle predicted yield with an average prediction error 4.2%. This study shows the immense potential of MODIS gridded vegetation product for keeping an up-to-date Geographic Information System of rice cultivation.


INTRODUCTION 1.1 Background
Rice is the second most harvested staple food in the world and the leading staple food in Asia.Rice would be contributes heavily to the alleviation of hunger and poverty since millions of small farmers grow millions of hectares of rice in the Asian region.In addition, there are millions of landless workers who generate some income by working on these farms.Moreover, 60% of the global population and 90% of the world rice production is derived from the Asian continent (Claessens, 2013).
Rice monitoring and mapping are thus very important for optimizing food security, environmental sustainability, water security, greenhouse gas emission reduction and general economic prosperity.Most countries in the Asian region use limited survey methods to collect rice paddy data from community level to national level.These data sources do not meet the needs of science and policy researchers, who need geo-spatial databases of rice agriculture with updated spatial and temporal resolution (Xiao et al., 2006).
Remote Sensing is becoming an essential tool for monitoring and mapping rice growing over large areas, at repeated time intervals (Son et al., 2012).According to a review article of Remote Sensing of rice crop areas, Remote Sensing combined with Geographical Information System (GIS) (Kuenzer and Knauer, 2013), can provide reliable information for the following purposes related to rice farming: • Mapping and monitoring the extent of rice growing ecosystems * Corresponding author.
• Monitoring and assessment of rice growth and health status • Assessment of crop pattern and crop system efficiency • Estimation of crop growth related parameters • Precision agriculture • Improvement, and extension of, crop growth and yield models High temporal and medium spatial resolution optical Remote Sensing Vegetation Indices (VI) known as Normalized Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) and Land Surface Water Index (LSWI) are utilized to map and monitor rice agriculture on a small scale, and analyzed in a global, continental, and subcontinental context (Sakamoto et al., 2005); (Xiao et al., 2006).Moderate Resolution Imaging Spectroradiometer (MODIS) sensor data was the major data source for those studies.

MODIS Gridded Vegetation Product
The MODIS land science developed Global MODIS VI to provide consistent spatial and temporal comparisons of vegetation conditions.Visible red (620 670nm), near infrared (NIR) (841 876nm) and visible blue (459 479nm) surface reflectance of MODIS sensor data are used to make global NDVI and EVI products with 250m, 500m, 1000m, and 5600m spatial resolution, and sixteen-day, or one-month temporal resolution.MODIS sensor contains both a Terra and Aqua satellite platform, which were used to make vegetation products MOD13XX and MYD13XX.
(The O refers to Terra platform, the Y refers to Aqua platform, and XX refers to one letter and one number.For example, if XX is Q1, then a sixteen-day composite image with 250m resolution is shown.If XX is A3, then a one-month composite with 1000m resolution is presented).Each year, one satellite sensor provides 23 images of sixteen-day composites.Production of Terra and Aqua products can be improved using the following temporal phase shifts: Terra 16-day period starting day 001, with Aqua 16-day period starting day 009) (Didan et al., 2010).(Shao et al., 2011).
NDVI defined by equation ( 1) is a normalized difference measure, comparing the near infrared and visible red bands defend by equation 1 where ρN IR (841-876nm) and ρRed (620-670nm) are the bi-directional surface reflectance for their respective MODIS bands.The NDVI is closely correlated with paddy area Leaf Area Index (LAI) (Xiao et al., 2006).
The MODIS Land VI QA provides quality assurance data for each pixel of the product as a bit field value.Composite dayof-year layer offers information for each pixel, representing the Julian day during the sixteen-day compositing period when reflectance was recorded (Pringle et al., 2012).
The noticeable difference between MODIS NDVI and EVI is the dynamic range.In high bio mass areas, MODIS NDVI appears to be saturated, whereas EVI shows more sensitivity (Shao et al., 2011).

Research Problem and Objectives
National and regional level rice crop monitoring and mapping is not enough to address food security, water management, and environmental stability problems in detail.Provincial or district level crop monitoring and mapping will be an effective way to address the above problems.
A number of studies have involved high level image processing techniques to rice crop monitoring, using satellite data.However, these methodologies were not applied to develop an active upto-date rice crop monitoring system at the provincial-or districtlevel, which can be easily used by government officials as well as administrators to manage and plan rice cultivation.
The specific objectives of this study were to investigate general applicability of high temporal resolution MODIS 250m gridded vegetation product for rice crop growth monitoring, rice crop acreage mapping, and yield analysis at the province level.This study addresses the following three primary research questions concerning rice crop monitoring and mapping: I. How is the rice plant growth illustrated by the Remote Sensing VI? II.How can you apply Threshold (Beurs and Henebry, 2010) method to determine rice crop information (area, phenology) using MODIS data?III.How do MODIS NDVI and EVI time series data represent a rice crop and its yield?
To address these questions, a Remote Sensing image analysis was performed combination with statistical analysis using MODIS 250m gridded vegetation product data from May to December 2013 across the Sa Kaeo Province in Thailand.

Step Processing
The whole study can be divided into the following four parts: (1) data pre-processing, (2) rice plant growth analysis using VI, (3) extracting rice crop parameters (e.g.acreage, start-of-season mapping from VI time series and accuracy assessment), and (4) rice yield analysis with MODIS VI.
2.3.1 Data pre-processing: Vegetation index layers (NDVI and EVI), VI QA layer, spectral bands (red, NIR, and MIR), and composite day-of-year layer were extracted from each MODIS HDF product as Geotiff images and then projected to WGS84 geographic projection.Eight days of NDVI and EVI images were stacked as a time series and detailed beyond Sa Kaeo province were masked out.The main element of pre-processing was the smoothing of the time series image stack.The VI time series was smoothed by using modified Whittaker Smoother with negatively biased noise of VI QA values and composite day-of-year values (Atzberger and Eilers, 2011); (Mattiuzzi and Lobo, 2013).Figure 2 shows selected time series before and after smoothing.

Extracting Rice Cultivation Parameters:
There are several spatio-temporal methods for extracting rice acreage and phenology from VI time series (Kuenzer and Knauer, 2013).This research work utilizes the simple threshold method (Beurs and Henebry, 2010).First, VI time series graphs were plotted for each land cover type (forest, water, urban, and paddy) to identify patterns and determine threshold values.After analyzing each land cover time series graph, an unique pattern was observed in  If the time difference was more than 80 days, those pixels were considered as representing a double rice paddy.There are two reasons for this condition: (1) a typical rice cycle is about 120 days, and thus cannot be much less than that, and (2) cassava and maze are cultivated with rice as a double crop pattern in some areas of the Sa Kaeo province.ALOS forest cover data was used to mask out the forest area in the study area.
In general, Sa Kaeo province rice cultivation cycle varies from 120 days 140 days.Therefore the rice planting date was occurred 60-70 days before the heading date.The start-of-season for each rice cultivation pixel was calculated by subtracting 70 days from heading date, and start dates were grouped monthly.
Accuracy of rice acreage was measured by comparing provincelevel and district-level rice cultivation area statistics and calculating the confusion matrix and kappa coefficient using ground truth from the field data.

Rice yield analysis with Remote
Sensing VI: Chlorophyll content is related to the production of yield of the rice plant and it is reflected by Remote Sensing VI.The VI data was therefore used to determine if there is a correlation between MODIS VI and rice yield, and how the correlation changes with NDVI and EVI.Yield assessment was done using field data gathered from Site 2 and Site 4. The correlation coefficient and liner regression analysis of yield was evaluated against maximum VI of the season, cumulative mean VI of the season, and VI of the flowering date (for both NDVI and EVI).Regression analysis results from each VI were compared to determine which VI was associated with the highest rice yield.The yield of the Site 3 was predicted from the most highly correlated VI model and compared to the actual yield.Figure 6 illustrates rice plant height time series in each field site.

RESULTS & DISCUSSION
According to the time series graph, rice plants grew more than 1.2m height and plants grew with approximately uniform speed.

Plant height and vegetation indices
We can observe how VI changes with plant height from the graphs of Figure 7.When plant height reached between 0.7m and 1.0m,    Table 3: Accuracy assessment of EVI-based rice acreage mapping 138795ha, and total acreage determined by EVI was 152944ha.These estimates deviate from the published statistics by +3529ha (NDVI) and +17678ha (EVI), producing relative error of 2.6% (NDVI), and 13.1% (EVI).At the district-level, rice acreage accuracy was calculated by comparing the VI-determined acreage with the statistics reported by Watthana Nakhon District Agriculture Office.Rice cultivation area within Watthana Nakhon District is 33174ha.However, the district rice acreage determined by NDVI was only 29491ha, and the district acreage determined by EVI was only 31916ha, producing relative error of 11.1% (NDVI) and 3.8% (EVI).These statistics suggests rice acreage can be esti-mated by NDVI and EVI with in 14% maximum deviation.Confusion matrix analysis shows that rice acreage estimation using NDVI and EVI has an accuracy of more than 75%.

Yeild Analysis
Rice yield and NDVI/EVI correlation analysis was done in three different ways.First, a correlation and linear regression analysis was performed using maximum VI of the time series and the yield.For the yield with maximum NDVI and maximum EVI, the correlation coefficient were 0.45 and 0.98, respectively.Analysis shows that the maximum EVI has the highest correlation with yield.The R-square values from linear regression analysis are 0.2 (NDVI) and 0.95 (EVI).
Second, the correlation between yield and flowering date VI was investigated.The average flowering date was founded by a survey done on farmers and it was on the 14th October, 2013.The R 2 values for linear regression are 0.72 (NDVI) and 0.92 (EVI).The results of all three yield analyses show that yield has a higher correlation with EVI than NDVI.The R 2 values of regression analysis with yield are 0.95 (for maximum EVI), 0.93 (for flowering dates EVI), and 0.88 (for cumulative mean EVI).
The highest correlated EVI was utilized to predict yield for Site 3, which was then compared with actual yield.The most precise prediction is given by the maximum EVI of the time series, with an average error of 4.21%.The flowering data EVI also give better results (5.62% error) than cumulative mean EVI (10.87%).
Figure 13 compares the actual yield with EVI-predicted yield.
Figure 13: Actual yield and predicted yield from EVI by different method in Site 3)

CONCLUSION
This study illustrates the significant potential of using MODIS data for rice crop monitoring.MODIS 250m gridded vegetation NDVI and EVI data from both sensors (Aqua and Terra) increase the frequency of observation and improve rice monitoring.The use of the Whittaker smoothing technique, which includes a combination of quality assessment data and R program library, is recommended (Atzberger and Eilers, 2011).
The results obtained by this research confirm that there is an observable relationship with rice plant growth and VI changes in rice paddy land cover.MODIS NDVI and EVI rice paddy acreage mapping (using the Threshold method and crop calendar data) gives a minimum of 86% area accuracy and 75% classification accuracy.Rice crop cultivation starting month extracted via MODIS NDVI data matched the field data.MODIS NDVI data is better than EVI for rice acreage mapping and start-of-season mapping.
Pixel level yield analysis shows yield and VI have a noteworthy relationship.EVI is better for yield analyzing and Maximum VI throughout the crop season which indicates the highest correlation with yield.Yield was predicted for site 03 with 4.2% error using maximum EVI of the season.
It was proved that time series MODIS 250m vegetation index product data has immense potential for monitoring rice crops.
This study shows the usefulness of MODIS gridded vegetation product for the following techniques: (1) mapping rice cultivation area at the province-level, (2) generating an up-to-date database of rice acreage and start-of-season dates, and (3) predicting yield of rice crop before harvesting.Applying this novel methodology for implementing up-to-date GIS rice database is recommended for other areas where rice cultivation is prominent.

Figure 1 :
Figure 1: Map of the Sa Kaeo Province and Watthana Nakhon District plant growth monitoring with RS VI: Rice plant height data from field observation was plotted against extracted corresponding pixel values in VI (NDVI, EVI) data.Rice plant growth and VI pattern was observed.

Figure 2 :
Figure 2: VI time series before and after applying smoothing

Figure 3 :
Figure 3: NDVI Time series graphs of different land cover

Figure 5 :
Figure 5: Geo-tag photograph time series of rice paddy cultivation

Figure 6 :
Figure 6: Rice Plant Time Series

Figure 7 :
Figure 7: Graph of Rice plant height vs VI (NDVI and EVI)

Figure 10 :
Figure 10: Single rice crop area start-of-season in Sa Kaeo Province, using NDVI

Figure 12 :
Figure 12: Correlation coefficient changing graph of the Cumulative VI mean and yield)

Table 1 :
Filed visit and field data collection was performed at four paddy plots in Watthana Nakhon District over the course of one complete rice cycle, during which time the MODIS satellites were observing the crops.A summary of the field sites is presented in Table1.Plant height, water level in the field and canals and yield are collected during the field visit.Geotag photographs were taken in every visit to record the change of paddy growth.Province-level agriculture statistics were obtained from the Office of Agriculture Economics and district-level agriculture statistics were provided by the Watthana Nakhon District Agriculture Office.Summary of the field data collection sites 2.2 Data 2.2.1 Remote Sensing / GIS Data: MODIS gridded VI product with 250m resolution (MXD13Q1) h27v07 tile was used as the main Remote Sensing data.The VI product data from both Aqua and Terra Satellites (31 images from each platform) was taken over an eight month period (from 2013-05-01 to 2013-12-27).Advance Land Observation Satellite (ALOS) global forest/non forest data was used to identify the forest area in the study.2.2.2 Field Statistics:
Accuracy assessment of NDVI-based rice acreage mappingProvince-level rice acreage was calculated by comparing the VIdetermined acreage with the statistics provided by the Office of Agriculture Economics.Total acreage determined by NDVI was

Table 4 :
Table 4 summarizes the area-based accuracy.Summary of the area based accuracy assessment

Table 5 :
Table 5 summarizes regression analysis results of the above two methods.Regression analysis results of the highest correlation between cumulative mean VI and yield are described in Table6.R 2 values from linear regression analysis are 0.71 (NDVI) and 0.87 (EVI).
Regression for yield with maximum VI and flowering date VIA third yield analysis compared yield and cumulative mean of VI.Cumulative VI mean was calculated for the eight days after the start-of-season.The correlation coefficient of each cumulative VI mean and yield was also calculated.The cumulative VI with highest correlation coefficient was used in a regression analysis of the yield.Figure12illustrates how the correlation coefficient changes with cumulative NDVI/EVI and yield.Both NDVI and EVI exhibits a similar pattern of correlation coefficient; both VI are negatively correlated with yield at the start of the season, and then become positively correlated, with high correlation coefficient values.