DERIVING CROP CALENDAR USING NDVI TIME-SERIES

Agricultural intensification is defined in terms as cropping intensity, which is the numbers of crops (single, double and triple) per year in a unit cropland area. Information about crop calendar (i.e. number of crops in a parcel of land and their planting & harvesting dates and date of peak vegetative stage) is essential for proper management of agriculture. Remote sensing sensors provide a regular, consistent and reliable measurement of vegetation response at various growth stages of crop. Therefore it is ideally suited for monitoring purpose. The spectral response of vegetation, as measured by the Normalized Difference Vegetation Index (NDVI) and its profiles, can provide a new dimension for describing vegetation growth cycle. The analysis based on values of NDVI at regular time interval provides useful information about various crop growth stages and performance of crop in a season. However, the NDVI data series has considerable amount of local fluctuation in time domain and needs to be smoothed so that dominant seasonal behavior is enhanced. Based on temporal analysis of smoothed NDVI series, it is possible to extract number of crop cycles per year and their crop calendar. In the present study, a methodology is developed to extract key elements of crop growth cycle (i.e. number of crops per year and their planting – peak harvesting dates). This is illustrated by analysing MODIS-NDVI data series of one agricultural year (from June 2012 to May 2013) over Gujarat. Such an analysis is very useful for analysing dynamics of kharif and rabi crops.


INTRODUCTION
Agriculture is one of the most important sectors for India.It is necessary for our country to arrange enough food for the people of our country.Proper planning for this sector requires relevant and reliable information in timely manner.Planning commission of India recognizes the crucial role of remote sensing (RS) technology in generating quality crop statistics (Gupta and Rajan, 2009).One of the first requirements for drawing a development plan is to have an inventory of the resource.India is one of the few select countries which have a reliable system of generating forecast of crop production using RS and other collateral sets of data.Crop inventory includes knowing area under a crop as well as its place, time and quantity.Since crop is a very sensitive biological system which is affected by many biotic and abiotic factors, it is essential to have dynamic assessment crop condition so that crop growth cycle can be simulated and reliable forecasts of crop production can be generated and/or modulated.
The chronological sequence of the occurrence of different physiological stages of a crop in its growth cycle forms the crop calendar.Crop calendar is one way of summarizing the information of a crop cycle.Information about crop calendar is essential for proper management of agriculture.Information on a crop calendar is a prerequisite in all crop estimation surveys using remote sensing data.
Remote sensing is a valuable tool for generating crop statistics that saves resources like time and money and provides real time information to farmers, policy makers, etc. (www.fao.org/agriculture/seed/cropcalendar).The use of remote sensing data * Corresponding author for estimating crop area has reached an operational level (Navalgund et al., 1991;Parihar and Dadhwal, 2002;Parihar and Oza, 2006).But the dynamics of agriculture requires monitoring the area repetitively in a year in order to map the crop regions.This calls for a time series analysis of satellite images.This is made possible by the advent of moderate resolution satellite imagery with high revisit such as that of MODIS (Moderate-Resolution Image Spectroradiometer) data.
Vegetation indices (VI) computed from satellite images gives an indication of the presence of vegetation and its health.Several studies on remote sensing applications have proved that VI can be used effectively in crop monitoring as well as in characterizing the vegetation with phenology.Time series profiles of VI derived from satellite data are potential tools to interpret the dynamics and phenological development of vegetation in different areas.
The present study demonstrates the use of temporal NDVI data in extracting key phenological information calendar (i.e.planting and harvesting dates and date of peak vegetative stage) by interpreting temporal variations in spatial domain to generate crop calendar.

2.
STUDY AREA AND DATA USED

Study area
The study area was taken as whole Gujarat state of India (Figure 1).Gujarat is the main producer of tobacco, cotton, and groundnuts crops in India.Other major food crops produced in Gujarat are rice, wheat, mustard, jowar, bajra, maize, tur, and gram.Gujarat has an agricultural economy; the total crop area amounts to more than one-half of the total land area (www.agri.gujarat.gov.in).

MODIS data sets
The MODIS Terra 8-day composite Surface Reflectance Product (MOD09Q1) is used.Each 8-day composite includes estimates of ground spectral reflectance of the two spectral bands at 250-m spatial resolution.These data sets provided for this product include reflectance values for Bands-1(Red) and 2(NIR), and a quality rating.MOD09Q1 provides Bands 1 and 2 at 250-meter resolution in an 8-day gridded level-3 product in the Sinusoidal projection.Each MOD09Q1 pixel contains the best possible L2G observation during an 8-day period as selected on the basis of high observation coverage, low view angle, the absence of clouds or cloud shadow, and aerosol loading.In this study, we downloaded MOD09Q1-H24V06 tile index data for May 2012 to June 2013 (fifty-five data sets of 8day composites) from the USGS EROS Data Centre (www.edc.usgs.gov).
Figure 1.FCC of MODIS (Tile H24V06) covering the study area (Gujarat state) 2.3 NDVI Normalised Difference Vegetation Index (NDVI) utilises the absorptive and reflective characteristics of vegetation in the red and near infrared portions of the electromagnetic spectrum.Therefore, changes in the NDVI time-series indicate changes in vegetation conditions proportional to the absorption of photosynthetically active radiation (Sellers, 1985).Therefore, NDVI was used in the present study.For each 8-day composite, we calculated NDVI using surface reflectance values from the red (620-670 nm), NIR (841-875 nm) bands using following equation: (1) where, NIR and Red represent reflectance in NIR and Red region of electromagnetic spectrum, respectively.The NDVI is found to be an effective measure for crop growth assessment and monitoring.Satellite data are radiometrically normalized to compensate for changes in radiation caused by differences in illumination, atmospheric conditions, sun zenith angle, sensor and satellite characteristics, time of data acquisition, etc. (Jonna et al., 1994).

METHODOLOGY
The major steps of methodology are (i) NDVI image generation from geo-referenced multi-date MODIS data; (ii) overlaying of land use land cover (LULC) mask; (iii) two stage of NDVI profile and (iv) generation of crop calendar.The methodology adopted in this study is shown in figure 2.

NDVI generation
Even though the available NDVI data sets were corrected for gaseous and aerosol scattering, the thick clouds still remained as noise (Los et al., 1994).To remove data affected by thick clouds information on quality control flags in the MODIS file was used (Xiao et al., 2005).Crops exhibit distinctive seasonal patterns.The NDVI time series data can be used to dynamically reflect crop growth and track crop phenological metrics change, because of the positive correlation between NDVI and LAI (Leaf Area Index) (Jakubauskas et al., 2001).

LULC Mask
For further refinement LULC mask for Gujarat state provide by NRDB database was used.These LULC was made by IRS AWiFS data source and 1:250,000 scale agricultural layer were created for making agriculture mask for India.Using this agriculture mask non-agriculture pixels were omitted for further analysis.

Data filtering
The high level of noise present in MODIS data often makes it difficult to determine the number of annual seasons.For this reason, smoothed NDVI data is used instead of raw NDVI series.The noise may be due to instrumentation behaviour, changes in sensor angle, atmospheric conditions (clouds and haze) and ground conditions.
The filtering was achieved in two-steps.In the first step, the subsequent NDVI value is accepted only if it is within a certain range.In the present case a threshold of 20 percent was used (Viovy et al., 1992).If the subsequent value is beyond the threshold value then it is replaced by average of previous and next NDVI value.The idea is that sudden rises and falls in NDVI are not compatible with the gradual process of growth, but may be due to departure in environmental conditions or viewing geometry.This eliminates high frequency 'noise' related changes and allows genuine changes in NDVI to be represented.This is followed by 13-point Gaussian filtering.
Essentially it is a weighted low-pass filter in time domain.Thus, the two-step filtering of NDVI series adopted eliminates unrealistic changes and outputs a smooth NDVI series.Smoothed series is used in subsequent analysis.

Single/Double crop detection
A peak finder algorithm designed to reject small depressions as seasons reasonably finds the number of cropping seasons for a particular pixel.
From the smoothed data, all the local peaks and troughs were identified.Troughs just before and after this time were identified.This trough -peak -trough duration was considered as a valid crop cycle if and only if (i) duration is greater than 60 days and (ii) NDVI value is greater than 0.3, inherent assumption is that a crop would be of atleast 60 days and if it is a cropped area, than it should be greater than 0.3 NDVI value.This will eliminate features as water, desert, built up and barren land, if not eliminated by applying the LULC mask.If the spectral emergence is detected after Nov 01, then it was considered as rabi crop otherwise the crop was categorized as Kharif crop.
A program was written in C-language to perform this task.

RESULTS & DISCUSSIONS
There were forty six, 8-day NDVI values considered in a year.Since moving filter is used, the series was extended by taking four additional data on either side.By doing so, the effects in head and tail region are reduced considerably.After smoothing, these additional values were deleted from further analysis.
A graph indicating necessity of smoothing and improvement achieved by adopting suggested two-step procedure is shown in figure 3. Figure 3 represents the behavior of a typical raw NDVI of a location.It can be seen that there are unrealistic fluctuation which are reduced after first stage of (20 % threshold).After 13-point Gaussian filter the artifacts in the series are eliminated and dominant cyclic patterns emerge.The value of the following equation 46 V 1 = ∑ (y t -y t-1 ) 2 (2) t=2 which measures roughness in the series.Rough series has higher value where as smooth series has lower values of this measure.For the representative pixel plotted in figure 3, V 1 for raw series was 0.86 whereas after first stage of filtering (20 % gradient), this value was reduced to 0.47.After second stage of filtering (13-point Gaussian filter), this further reduced to 0.054.This indicates success in eliminating noise from the series.The filtered is now amenable to further analysis and extraction of desired information.From table 4, it can be seen that out of total cropped pixels, about 57.9% area is under single crop.For single crop pixels, about 21.6% is having crop length of 200 days or more.The area under category 120-200 crop length days is about 18.6% and less than 120 days is about 17.7%.Double cropped is followed in about 39% area.The first cycle of the double crop region is kharif crop and most of this category is under less than 120 days.While second cycle is rabi Is Quality flag acceptable ?
Calculate NDVI Ignore from further processing Applying Gaussian Filter

NDVI Image Stack
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-8, 2014ISPRS Technical Commission VIII Symposium, 09 -12 December 2014, Hyderabad, India This contribution has been peer-reviewed.doi:10.5194/isprsarchives-XL-8-869-2014crop and its major cropped cover under 120-200 days and less than 120 days.While more than two cropped area is about 2.9% area only.Overall crop seasons having less than 120 crop days is about 46.9%, 120-200 crop days is about 35.5% and more than 200 days is about 17.6%.Knowing area under single, double and multiple cropping intensity defined as gross cropped area divided by net sown area, can be calculated as (1*Single Crop + 2*Double Crop + 3*Multiple Crop) / Total Cropped Area.In the present study cropping intensity for Gujarat in 2012-13 is 1.45.
In Gujarat for the year 2011-12, total area under kharif crop was 10.35 M ha and rabi crop area was 4.79 M ha (DOA-Gujarat, 2011).Assuming that farmers growing crop in rabi season would have taken a crop in kharif season, the area under kharif season can be considered as a first approximation (in fact lower bound) of net sown area.Net sown area by remote sensing based method adopted in this study was 10.06 M ha, which is 97% of ground area reported.The cropping intensity calculated from ground data was 1.46 where as from present study it came out to be 1.45.Thus this can be considered that there is a good agreement between DOA-Gujarat, 2011 data with present study output data.

CONCLUSIONS
A methodology was developed to extract some key elements (i.e.number of crops / year and their planting -peakharvesting dates) of crop growth cycle in a region.The methodology consists of applying two-step smoothing followed by interpreting temporal variations.The methodology is illustrated by taking an example of temporal NDVI data over Gujarat for 2012-13.For this region, about 57.9% area was single crop and 39.2% was double crop.Such an analysis is very useful for analysing dynamics of kharif and rabi crops.

Figure 5 .
Figure 5. Spatial distribution of crops for Gujarat in 2012-13 (a) Number of Crops (b) Crop days within Single Crop (c) Crop days within Double Crop

Table 4 .
Category-wise distribution of pixels for crop cycles versus crop length in daysThe International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, VolumeXL-8, 2014   ISPRS Technical Commission VIII Symposium, 09 -12 December 2014, Hyderabad, IndiaFrom last column of table 4 it be seen that Gujarat is having 17.6% crop area with duration more than 200 days.Typically castor, tur, sugarcane and cotton in Gujarat have crop duration have more than 200 days.It can also be seen that 35.5% crop area is having crop duration of 120-200 days.Remaining 46.9% crop area is having crop duration less than 120 days.Predominant region of double crop area is in parts of North and Central Gujarat.There are pockets of multiple cropped areas in and around Banaskantha district in north Gujarat.The regions of double and multiple crops are regions where cash crops (e.g.tobacco, potato and cumin) and/or coarse cereals such pearl millet is grown.Figure5(b) shows distribution on duration of a crop within single crop area.Major area having long duration crops (i.e. more than 200 cropping days) can be seen in parts of north, central and south Gujarat regions which are regions of sugarcane, cotton and castor cultivation (Crop Weather Calendars for Gujarat, 2002).Figure5(c) shows distribution of crop duration (short, medium and long) within double crop.Large proportion of area having 120-180 cropping days region is seen in parts of north and central Gujarat districts.Long kharif and short rabi season area shown in parts of south Gujarat region.Thus it can be seen that, qualitatively the output of present analysis agrees well with the known existing agricultural patterns of Gujarat.