SPATIAL PATTERN OF TEMPORAL TREND OF CROP PHENOLOGY MATRICES OVER INDIA USING TIMESERIES GIMMS NDVI DATA ( 1982-2006 )

NOAA-AVHRR bi-monthly NDVI data of 8×8 km for the period of 1982-2006 were used to analyze the trend of crop phenology matrices over Indian region. Time series principal component analysis of NDVI was performed to produce six calibration zones for fitting equations of temporal NDVI profile. Savitzky-Golay filter with different seasonality parameters, adaptation strengths and window sizes for different calibration zones were use to smoothen the NDVI profile. Three crop phenology matrices i.e. start of the growing season (SGS), Seasonal NDVI amplitude (AMP), Seasonally Integrated NDVI (SiNDVI) were extracted using TIMESAT software. Direction and magnitude of trends of these crop phenology matrices were analyzed at pixel level using Mann-Kendall test. Further the trends was assessed at meteorological subdivisional level using “Field significance test”. Significant advancement of SGS was observed over Punjab, Haryana, Marathwada, Vidarbha and Madhya Maharashtra where as delay was found over Rayalaseema, Coastal Andhra Pradesh, Bihar, Gangetic West Bengal and Sub-Himalayan West Bengal. North, West and central India covering Punjab, Haryana, West & East Uttar Pradesh, West & East Rajasthan, West & East Madhya Pradesh, Sourastra & Kutch, Rayalaseema, Marathwada, Vidarbha, Bihar and Sub-Himalayan West Bengal showed significant greening trend of kharif season. Most of the southern and eastern part of India covering Tamilnadu, South Interior Karnataka, Coastal Andhra Pradesh, Madhya Maharashtra, Gujarat region, Chhattisgarh, Jharkhand and Gangetic West Bengal showed significant browning trend during kharif season.


INTRODUCTION
Crop phenology is a study of recurring pattern of crop growth and development.Changes of crop phenology are recognized as being possibly the most important early indicator of impact of climate change on ecosystem.Variation of crop phenology has broad impacts on terrestrial ecosystems and human societies by altering global carbon, water and nitrogen cycles, crop production, duration of pollination season and distribution of diseases etc (Penuelas and Fiella, 2001;Soudani et al., 2008).Therefore, phenology has emerged recently as an important focus for ecological and climatic research (Menzel et al., 2001;Cleland et al., 2007).In situ observations, bioclimatic models and remote sensing constitute the three possible ways for monitoring the variations of vegetation phenological events (Schaber and Badech, 2003;Fisher et al., 2007).Remote sensing has the advantages of being the only way of sampling at low-cost with good temporal repeatability over large and inaccessible regions.Phenology from satellite observation is aggregate information at coarse spatial resolution that relates to the timing and rates of greening (growth) and browning (senescence), timing of maximum photosynthetic activities and duration of active growth phase at seasonal and inter annual time scales.Crop phenology is usually addressed through temporal monitoring of the normalized difference vegetation index (NDVI) time series, available since the launch of the earth observation satellites.The NDVI is calculated as the normalized difference between near-infrared and red bands, and quantifies the photosynthetic capacity of plant canopies, thus indicating the amount of vegetation present in one place.Time series NDVI data sets of NOAA/AVHRR, TERRA/MODIS and SPOT/VGT have been used in numerous studies to estimate the shift of vegetation phenological parameters and its relationships with different climatic and anthropogenic factors (Menzel et al., 2001;Zhou et al., 2001;Tao et al., 2006;Stockli and Vidale, 2004;Fontana et al., 2008;van Leeuwen et al., 2010;Melezdez et al., 2010;Kariyeva et al., 2011).Most of these studies were carried out over North America, Europe, China.Central Asia, Sahel and Sudan region.The present study is an attempt to analyze the direction and magnitude of trend of few crop phenology matrices over Indian region.

METHODOLOGY 2.1 Study area
The study is carried out during kharif season over 25 meteorological subdivisions of India Meteorological Department (IMD) spread across 16 states of India (fig.1).The selected meteorological subdivisions have large contiguous areas of agricultural land and crop phenology could easily be indentified using coarse resolution time series satellite data.There are wide diversity in crop type and its cultivation practices across these meteorological subdivisions.During kharif season, rice is the dominating crop throughout India with some areas of cotton, sugarcane, maize and vegetable cultivation.

Satellite data pre processing
Global Inventory Modeling and Mapping Studies (GIMMS) AVHRR 8 km 15-days composite NDVI for the period of 1982 to 2006 (Tucker et al. 2005) were analyzed to monitor changes in crop phenology over Indian region.The dataset is corrected against residual sensor degradation, inter-sensor calibration differences, distortions caused by persistent cloud cover globally, solar zenith angle and viewing angle effects due to satellite drift and volcanic aerosols.A 15 days maximum value composite has been used to reduce cloud cover while maintaining temporal resolution for accurate phenological estimates.
The software TIMESAT was used to generate smooth time series of NDVI as well as to estimate the vegetation phenology for the study area.TIMESAT software has options to fit a smooth continuous curve to time series of NDVI data using Savitzky-Golay (SG), Asymmetrical Gaussion (AG), or Double logistic (DL) functions, and an adaptive upper envelope to account for negatively biased noise such as cloud (Jonsson & Eklundh. 2002, 2004).Indian agriculture differs greatly in terms of its time of sowing, cropping pattern, cropping intensity etc.Thus a single equation to fit NDVI temporal curve may not hold good across the area.To resolve this issue, Indian region was divided into different calibration zones, based on NDVI time series principal component analysis (TS-PCA), to select smoothening techniques and its associated parameters settings (Heumann et al., 2007).Potential kharif agricultural mask is prepared from multi-year land use land cover map using AWiFS data.To cope up with the difference in spatial resolution between GIMMS NDVI (8 km) and AWiFS (60m), majority rule approach was adopted to aggregate the potential kharif agricultural area mask to coarse GIMMS NDVI.The mask is applied to subset kharif agricultural area over India and processed further to extract information of crop phenology.NDVI images of the first fortnight of each month of the year 1984, 1995 and 2006 (probable normal year from each decade) were selected to perform TS-PCA.First three TS-PCA components explained 96.8% of variance.Six zones for calibration were identified iteratively by isodata clustering technique using these first three TS-PCA components.These six calibration zones have six distinct temporal responses of NDVI (figure 1).Savitzsky-Golay filter (Chen et al., 2004) was used with different seasonality parameters, envelope adaptation strengths and window sizes to smoothen of temporal NDVI profile over the calibration zones (Table 1).Smoothened NDVI and original NDVI of total time series was compared with least square techniques and majority of the area showed R 2 > 0.9.All the pixels were found to be significant (p = 0.05) under F-test.

Computation of crop phenology matrices
There are several phenology matrices that can be extracted from time series NDVI data using TIMESAT software (Fig 2).The present study analyzed three important phenology matrices i.e. start of the growing season (SGS), seasonal NDVI amplitude (AMP), seasonal greenness or seasonally integrated NDVI (SiNDVI) during kharif season.In general, there are two approaches for determining the SGS from the seasonal variation in NDVI.The first method uses a NDVI threshold to identify the beginning of photosynthetic activities in monsoon /kharif period.
The second one identifies the period of large increase in NDVI as the beginning of growing season.A limitation associated with the threshold approach is that various land use types requires the use of different thresholds.
Since most of the agricultural area of India varied greatly in the soil & crop type, time of showing etc, determining the optimal threshold value for whole India would be difficult.Thus, the second approach was used in the present study.The SGS were estimated as the time point at which the value of the fitted function has first exceeds 20% of the distance between minimum and miximum at the rising limb.During calibration, it was found that estimation of SGS using values less than 20% had increased error due to noise effects.AMP is defined as the difference between the base and maximum NDVI values.SiNDVI is defined as the area under the vegetation curve during the growing season and above the mean of the base value.

Spatial trends of crop phenology matrices
Anomalies of the crop phenology matrices from its long term mean were calculated and trend of it were analyzed using nonparametric statistical approach.The Mann-Kendall test was employed to detect presence of trends of phenology matrices, while Sen's method was used to compute the magnitude of the change.The advantages of this method over parametric statistics are that the missing values in time series data are allowed and the data do not require to conform any particular distribution.The Mann-Kendal test searches for a trend in a time series without specifying whether the trend is linear or nonlinear.Mann-Kendall method firstly calculates S statistic, which indicates the sum of the difference between the data points indicated in following equation.
Where, x i is observed value at time j, x k is observed value at time k, j is time after time k and n is the length of the data set.The sign of the value is defined as follows: The value of S indicates the direction of trend.A negative (positive) value indicate falling (rising) trend.Mann-Kendall have documented that when n 8 , the test statistics S is approximately normally distributed with mean and variance as follows: Where n indicates the number of data points, q indicates the number of tied groups (a tied group is a set of sample data having the same value) and t indicates the number p of data points in the p th group.
International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XXXVIII-8/W20, 2011 ISPRS Bhopal 2011 Workshop, 8 November 2011, Bhopal, India The statistical significance of S is checked using a test statistic (or Z score).The standardized Mann-Kendall statistics Z follows the standard normal distribution with zero mean and unit variance.
The Mann-Kendall method test the null hypothesis, H 0 , is that data indicate no distinct trend against the alternative hypothesis H i is that data indicate a trend.Two tailed test is used to decide whether or not the null hypothesis should be rejected in favour of hypothesis test.The null hypothesis, H 0 is tested by the Z test statistic value.A decreasing trend is identified if Z is negative and an increasing trend is identified if the Z is positive.In the present study, H is rejected at p = 0.1 if the absolute value of Z is greater than Z tab , where Z tab is taken from the standard normal cumulative distribution tables.
Linear regression method is generally used to estimate the slope of possible linear trend in time series but for noisy (outliers or missing data) short time series data it may grossly miscalculate the magnitude of trend.Sen's slope estimator is robust nonparametric trend operator and is highly immuned to gross data error (Sen., 1968).The break down bound of Sen's estimator is approximately 29%.Thus the estimated trends have to have persisted for more than 29% of the length of the time series.In the present study, magnitude of the trend over time is estimated according to Sen et al., 1968.It is calculated by determining the possible slope between all possible data pairs and then finding the median value as under; Where, Q i is slope between data points, x i is data measure-ment at time j, x k is data measurement at time k and j > k.

Field significance test
For each grid cell, and separately for SGS, AMP and SiNDVI anomalies data, Mann-Kendall test (at p= 0.1) for the null hypothesis of no trend were performed.While analyzing phenology matrices during kharif season over a large region like India to identify trends, it is also important to ascertain the computed trend across that region using the field significant test (Livezey and Chen, 1983).Field significant test is performed over each meteorological subdivision (field or global).First the trend is computed at each individual grid level and the number of grid having trends (increasing/decreasing) are used as an input to the field significant test for checking the hypothesis that the respective meteorological subdivision has any trend (increasing/decreasing).
In terms of the probability theory, a collection of N independent significance tests (trend analysis of N grid points) is perfectly analogous to N tosses of a loaded coin.Instead of Heads or Tails, with odds based on the coin load, the outcome for a 0.1 significance level test is passed (probability equal to 0.1) or test failed (probability equal to 0.9).This can be modeled with a binomial distribution as under; ‫ݎܲ‬ ሺ‫ܯ‬ሻ ൌ ൫ ே ெ ൯ ሺͲǤͳሻ ெ ሺͳ െ ͲǤͳሻ ேିெ (7) Where M = 0, 1, 2, .....The field significant test is considered to be passed if the number (M) of individual tests passing (out of N tests), with a significance level 0.1, exceeds a threshold number (M 0 ).As an example, a meteorological subdivision having 40 grids is tested for trends using Mann-Kendal test at p=0.1.As per binomial distribution (equation 7), the probabilities are exactly 0.01478, 0.06569, 0.014233, 0.20032, 0.20588, 0.16471 and 0.10676 for 0, 1, 2, 3, 4, 5 and 6 passed tests out of 40 grids.Thus, the cumulative probabilities are 0.20627 for five or more passed test and 0.099517 for six or more.Therefore, if field significant test of a meteorological subdivision (having 40 grids) is performed, minimum six grids (threshold number, M 0 ) must pass (significant at p = 0.1) to guarantee at least 90% (in this case 0.90048%) significance.Following this procedure, threshold number (M 0 ) of grids having significant trend were computed for each meteorological subdivision at 90% level of significance.The meteorological subdivisions were declared to have significant trend (increasing/decreasing) if total grids with significant trend (increasing/decreasing) is equal or more than the threshold number (M 0 ).

RESULTS AND DISCUSSIONS
Grid wise analysis of trends of kharif SGS, AMP and SiNDVI over different meteorological subdivisions are shown in Table 2.The numbers of grids having trend (positive/negative) or no trend are presented as N in

Table 1
. Variable settings of parameters of Savitzky-Golay filter of TIMESAT over different calibration zones.(CZ:Calibration Zones, SP: Seasonality Parameters, AS: Adaptation Strength)

Table 2
International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XXXVIII-8/W20, 2011 ISPRS Bhopal 2011 Workshop, 8 November 2011, Bhopal, India India.Southern and western part showed negative trend of SiNDVI.Statistically significant positive trend of SiNDVI (0.02 to 0.05 year -1 ) was found over Punjab, Haryana, West & East Uttar Pradesh, West & East Rajasthan, West Madhya Pradesh, Sourastra & kutch, Rayalaseema, Marathwada, Vidarbha, Bihar and Sub-Himalayan West Bengal.Statistically significant negative trend of SiNDVI (-0.02 to -0.04 year -1 ) was found over Tamilnadu, South Interior Karnataka, Coastal Andhra Pradesh, Madhya Maharashtra, Chhattisgarh, Gujarat and Gangetic West Bengal.East Madhya Pradesh passed field significant test both for positive and negative categories.Northern part of East Madhya Pradesh showed significant positive trend where as southern part shown significant negative trend of SiNDVI.Spatial pattern of median trend of kharif season NDVI amplitude or crop vigour (fig.3c)showed positive trend over north, east, central regions of India, where as southern part shown negative trend.Statistically significant positive trend of AMP (0.002 to 0.007 year -1 ) is observed over Punjab, Haryana, West & East Rajasthan, East Uttar Pradesh, Bihar, Sub-Himalayan West Bengal, West Madhya Pradesh, Sourastra & Kutch, North Interior Karnataka, Rayalaseema.On the other hand, significant negative trend of AMP (-0.002 to -0.003 year -1 ) is found over Tamilnadu, South Interior Karnataka, Coastal Andhra Pradesh, Marathwada, Madhya Maharastra, Orissa, Gujarat and Jharkhand.Field significant test was found to be passed both for negative and positive trend in case of Telangana, Vidarbha, Chhattisgarh, Gangetic West Bengal, East Madhya Pradesh, West Uttar Pradesh.Parts of these meteorological subdivisions showed both significant positive and negative trend of kharif seasonal NDVI amplitude.Figure 1.Study area with different calibration zones.Meteorological subdivisions under the present study: Punjab(1), Haryana Chandigarh Delhi (2), West Uttar Pradesh (3), East Uttar Pradesh (4), West Rajasthan (5), East Rajasthan (6), West Madhya Pradesh (7), East Madhya Pradesh (8), Gujarat region Daman DNH (9), Sourastra Kutch Diu (10), Chhattisgarh (11), Jharkhand (12), Bihar (13), Sub-Himalayana West Bengal (14), with increase in AMP was observed over Punjab, Haryana, West & East Uttar Pradesh, West & East Rajasthan, West & East Madhya Pradesh, Bihar, Sub-Himalayan West Bengal, Sourashtra & Kutch and Rayalaseema.Whereas, increasing SiNDVI along with decreasing of AMP was found over Marathawada, Vidarbha implying increase of the length of the growing.Significant browning trends during kharif season were observed in most of the south and eastern part of India covering Tamilnadu, South Interior Karnataka, Coastal Andhra Pradesh, Madhya Maharashtra, Gujarat region, Chhattisgarh, Jharkhand and Gangetic West Bengal.Shifts in the crop phenology matrices are a composite effect of climate and anthropogenic interventions such as introduction of new crop /crop varieties or management practices or irrigation facilities etc.Further research work is required to analyze the functional relationship between different causative factors (anthropogenic and climatic) and shifts in crop phenology matrices.