A phenology-based classification for crop in Great Mekong Subregion based on MODIS data

The difference and the complementarity on agriculture decides that there is a wide-ranging cooperation on agriculture between China and GMS. So the research on GMS’s agriculture has important significance. The objective of this research is to produce the phenology-based classification map for crop in GMS by investigating the spatial and temporal patterns of multiple cropping systems. Although the spatial resolution of MODIS-NDVI data is not very high, it is still useful for detecting very large-scale crop phenomenon. In this paper 16-day time-series MODIS data from 2001 to 2015 were used to exploer the imforamtion of crop phenology in GMS. After constructing a phenology-based approach, four types of cropping systems, namely single-cropped rain-fed crop, double-cropped irrigated crop, double-cropped rain-fed crop and triple-cropped crop, were classified in the region


Introduction
Agricultural land is one of the most important land uses in Earth especially in the Greater Mekong Subregion (GMS). The contribution of agriculture to gross domestic product (GDP) of the countries in the region has been declining over the past decade; however, the sector still plays a very important role in food security and providing employment for its population. Over the past several decades changes in land use and agricultural intensification have come at an environmental cost in GMS. Regular and timely monitoring of crop system at regional to local scale will result in efficient use of agricultural land resources, as well as provide valuable information in agricultural management to meet the needs for food while safeguarding the environment. Remote sensing offers a feasible tool for crop classification by providing valuable information regarding the extent, status of agricultural systems at various spatial and temporal scales. A lot of researches about regional vegetation phenology by using time series vegetation index (VI) data were . VI data acquired from high spatial resolution satellite images such as Landsat TM/ETM+, etc. do not provide sufficient information on seasonal change characteristics of cropland due to infrequent temporal coverage. So the generally used VI data for crop phenology research are available at coarse and medium spatial resolution produced by MODIS,AVHRR and SPOT-VGT data. The advantages of using coarse and medium spatial and high temporal resolution VI data are (i) increased ability to monitor the phenological change of crop plants, and (ii) the possibility of generating consistent large area crop cover maps. To date, many researchers have applied time-series analysis techniques to the NDVI time-series to obtain the information of seasonal vegetation changes characterized. In this paper the main objective According to the statistic by FAO, only the deltas and low-lying floodplains produced half of region's rice production although they constitute only 10% of the total land area.

Data
The satellite sensor and ancillary data used in this research were MODIS NDVI time-series. The NDVI data used in this paper are the Version 5 MODIS/Terra 250m resolution 16-day Normalized Difference Vegetation Index (MOD13Q1 V005) spanning from 2001 to 2005 and 2011 to 2015. They were acquired from Land Processes Distributed Active Archive Center (LP DAAC) at no cost. Statistic data on province-level were used in this paper to validate the result of the method.

Method
There are four main steps involved in classification scheme:1)Reconstructing of time-series MODIS NDVI data by applying Savitzky-Golay filter to remove deviations from the vegetation cycle due to residual cloud and atmospheric effects and so on;2) Masking non-crop area by calculating for each pixel how many times the NDVI values within each year were less than 0.4 or greater than 0.75 and then setting a threshold to separate the crop and non-crop area;3)Producing the cropping intensity map by checking how many peaks and heading date of NDVI profile on a pixel-by-pixel basis, which then be used to separate the crop area into single-cropped area, double-cropped area and triple cropped area and then separate the double-and triple-cropped area into rainy and dry growing season;4)Validation analysis.

Time series NDVI reconstructing.
Extraction of the phenology of the crop requires high quality time series data. Even though the constraint view angle -maximum value composite method has been applied on the 16-day composite NDVI time series, the data remain significant residual effects and noises levels. In order to reduce the impact of the noises and improve the accuracy of the extraction result it is important to reconstruct high quality NDVI data. The reconstruct processing in this paper consists of the following steps:  detecting of useless pixels (according to the Pixel Reliability flags in MDOIS VI products) and low reliable values (identified if they exceeded a change threshold in time series NDVI);  replacing bad data and low reliable points with boxcar averages of the specified width;  applying Savitzky-Golay(S-G) smoothing filtering to smooth NDVI time-series and making a comparison between the filtered NDVI time-series and the NDVI time-series acquired from step;  replacing the points in the original NDVI time-series that are below a user-defined threshold with the value of the filtered NDVI time-series;  next iteration is needed and again the NDVI time-series is searched for possible bad data until no new points are being found. 3) Cropland generally displays a distinctly higher annual standard deviation compared to forest due to high greenness during the growing season and extremely low greenness following harvest. According to these characteristics we firstly calculated the number of NDVI time-series whose value is less than 0.4 on pixel basis. If the number is less than 20 and the maximum value of NDVI is less than 0.5 within 23 periods(one year), that pixel is determined to non-vegetation. The forest area could be identified by calculating the number of NDVI time-series whose value is greater than 0.75 on pixel basis.If the number is greater than 16 within 23 periods(one year), that pixel is determined to forest. Then all other pixels were categorized into crop class.

Extracting cropping intensity.
NDVI time-series can provide a seasonal trajectory of vegetation phenology. For example, single cropping system exhibits a unimodal NDVI periodicity while double cropping system exhibits a bimodal NDVI periodicity. So by detecting the numbers and the date of the peaks of the smoothed NDVI profile on a pixel-by-pixel basis, cropping intensity and the date when maximum NDVI occurred can be acquired. The peak finding algorithm used to identify this parameters has been described in details in the paper(Lv Tingting, Liu Chuang,2010).  Figure3. Statistic chart of date when maximum NDVI occurred in the study area Figure 3. is the statistic chart about time when peaks appear in NDVI curves. The horizontal (category) axis, also known as the x axis, of a chart displays the cropping patterns. The legend '1' means single cropping system.The legend '2-1'means that two growing season were estimated while '-1'means the first growing season while '-2'means the second growing season. The y axis displays the temporal intervals of the 16-day NDVI data.We can see that in the study area the peaks on the single cropping line mostly appeared between July and September. The peaks for double cropping system in the rainy season mostly appeared from October to December while peaks for double cropping system in the dry season mostly appeared from from April to June.This phenomena can be attributed to that irrigated could have different planting time in the dry season while rain-fed crop's planting time is mainly determined by the time of onset of rainfall. The peaks for triple cropping system mostly concentrated in three periods: from March to May, July to October and from November to February the next year. The onset of the date of maximum NDVI of crop may vary inter-annually which corresponds to crop variety and agricultural activities scheduling. Based on analysis above a phenology-based approach was constructed to classify four types of cropping systems in the study area: single-cropped rain-fed crop, double-cropped irrigated crop, double-cropped rain-fed crop and triple-cropped crop. If a single cropping pattern had a peak between July to November, it was classified as single-cropped rain-fed rice . If one of peaks for double cropping system occurred between January to June, it was classified as double-cropped irrigated crop and the other double cropping system deemed as double-cropped rain-fed crop. It must be mentioned that triple cropping system general sparsely distributed in the study area which characterized by smaller crop field sizes. So the result of this cropping system derived from coarse-spatial-resolution MODIS data has higher uncertainty. Considering this the sub-class for triple cropping system did not detect.

Validation
Due to lack of ground reference data we used crop statistic data from 2001 to 2010 acquired on national level to validate the result. Figure 5 is a comparison of cropping intensity from MODIS-NDVI and statistical data on provincial-level. Because the boundary of study area does not match the administrative boundary of province in GMS, only statistic data of provinces which are total included in the rectangle area were used in this paper.From figure 6 we could seen that in most of provinces the estimated intensity by MODIS data demonstrates a good association with statistical data at countylevel.