EVALUATION OF DIFFERENT PHENOLOGICAL INFORMATION TO MAP CROP ROTATION IN COMPLEX IRRIGATED INDUS BASIN

Accurate information of crop rotation in large basin is essential for policy decisions on land, water and nutrient resources around the world. Crop area estimation using low spatial resolution remote sensing data is challenging in a large heterogeneous basin having more than one cropping seasons. This study aims to evaluate the accuracy of two phenological datasets individually and in combined form to map crop rotations in complex irrigated Indus basin without image segmentation. Phenology information derived from Normalized Difference Vegetation Index (NDVI) and Leaf Area Index (LAI) of Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, having 8-day temporal and 1000 m spatial resolution, was used in the analysis. An unsupervised (temporal space clustering) to supervised (area knowledge and phenology behavior) classification approach was adopted to identify 13 crop rotations. Estimated crop area was compared with reported area collected by field census. Results reveal that combined dataset (NDVI*LAI) performs better in mapping wheat-rice, wheat-cotton and wheat-fodder rotation by attaining root mean square error (RMSE) of 34.55, 16.84, 20.58 and mean absolute percentage error (MAPE) of 24.56%, 36.82%, 30.21% for wheat, rice and cotton crop respectively. For sugarcane crop mapping, LAI produce good results by achieving RMSE of 8.60 and MAPE of 34.58%, as compared to NDVI (10.08, 40.53%) and NDVI*LAI (10.83, 39.45%). The availability of major crop rotation statistics provides insight to develop better strategies for land, water and nutrient accounting frameworks to improve agriculture productivity.


INTRODUCTION
Crop rotation is the growth of different crops in successive seasons of a year (Panigrahy and Sharma, 1997).Information on spatial distribution of different crop rotations helps managers to perform different agriculture functions such as water/nutrient supply, crop production estimation, revenue generation in an effective way (Jones et al., 2017;Bégué et al., 2018).It optimizes the modeling of energy fluxes for food security and climate change studies in different agro-ecosystems (See et al., 2015;Leng and Huang, 2017).Accurate crop area is also the backbone of agriculture statistics.Pakistan being an agricultural country depends heavily on reliable and timely crop statistics for its food and economic policies.Complete enumeration and sample survey are the two traditional methods to generate crop area statistics (FAO, 1982).Complete survey is expensive and time consuming whereas sample survey is timely and cost effective.Nevertheless, sampling design and estimation procedure is very critical in sampling method (Cotter and Nealon, 1987;FAO, 1995).The survey area statistics than compiled at reasonable administrative (mostly district) unit.
Geographic location of different crop rotations and its variation within a district is not existed under traditional crop reporting system.This lack of spatial information may result in a less effective planning of valuable agriculture resources within a district.Advances in satellite remote sensing has proven to be an effective tool for crop mapping and made it possible to estimate crops area at coarser resolution of one kilometer to finer resolution within few meters (Zhou et al., 2013;Jiao et al., 2014;Qin et al., 2015;Belgiu and Csillik, 2018).For small study areas, high-resolution images are used to identify different crops based on supervised classification approach using field signatures (Blaschke, 2010;Yang et al, 2011;Huang et al., 2017).This approach is time dependent and additional training data is needed when applied to other periods (Zhong et al., 2016).Regional level crop mapping is done using low spatial but high temporal resolution satellite sensors data such as Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) (Gumma et al., 2014;Zhang et al., 2015;He et al., 2017).Use of phenology is a key component in crop mapping when using multi-temporal images (Knight et al., 2006;Geerken, 2009;Zhang et al., 2014).Remote sensing based phenology information has widely used for biomass monitoring (Xiao et al., 2004;Lewinska et al., 2018), farm management (Lobell et al., 2003;Suepa et al., 2016) and climate change analysis (Tateishi and Ebata, 2004;Wang et al., 2017).Vegetation indexes (VIs) are widely used parameters, such as normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) to develop phenological metrics for crop classification (Wardlow et al., 2007;Wardlow and Egbert, 2008;Shao et al., 2010;Yan et al., 2015).Yet, crop identification becomes difficult with low-resolution images in heterogeneous large regions and its accuracy remains questionable (Loveland et al., 2000;Chen et al., 2016).Subdivision of different feature clusters based on temporal space is also difficult for large regions having more than one cropping seasons.This can be minimized by image segmentation (Vintrou et al., 2012;He et al., 2017).However, this spatial division is time consuming and needs addition information about study area (Schultz et al., 2015).Mixing of non-vegetative features (e.g., bare soil and built-up) with cropland also hamper the crop area estimation if spatial segmentation is not adopted for unsupervised cluster identification.
Leaf area Index (LAI) is another dimensionless physiological parameter that characterize the structure of plant canopy (Bréda, 2008).It is the ratio of leaf area per unit ground surface area (Watson, 1947).There are many direct (e.g., allometry, harvesting and littering collection) and indirect methods (based on radiation transfer through plant canopy) to measure LAI in the field (Gower et al., 1999;Bréda, 2003).Plant canopy directly linked with spectral reflectance therefore, LAI retrieval from remote sensing data has stimulated many research studies in recent years (Zheng and Moskal, 2009;He et al., 2016).Positive correlation has been reported in many studies between NDVI and LAI (Wang et al., 2005;Maki and Homma, 2014).Most studies have related remotely sensed LAI with photosynthesis, transpiration, rainfall interception and energy fluxes in context of forest cover (Soegaard, 1999;Leuning et al., 2005;Chen et al., 2007).Wang et al., in 2017 concluded that phenology information extracted from LAI performs better than enhanced vegetation index (EVI) for croplands and evergreen forest regions.However, there is still a gap to evaluate LAI for mapping different crop types based on its temporal behavior.
This study aims to evaluate crop rotation mapping in large and heterogeneous irrigated area of Indus basin using different phenological information.NDVI, LAI and fusion of these two datasets are evaluated to map crop rotation in large and diverse region.Unsupervised classification is used to identify different clusters in the regions based on temporal space without spatial segmentation.Different clusters then classified based on curve behavior and spatial occurrence.Accuracy of estimated crop area is checked against area reported in statistical archives.

Study area
Study region consists of 48 canal command areas that spread over four Provinces (Punjab, Sindh, Khyber Pakhtunkhwa (KPK) and Balochistan) of Pakistan.More than 90% irrigated area lies under Punjab (33 districts) and Sindh (23 districts) Province.Area lies between 24 • 05 ' to 34 • 31 ' N latitude and 67 • 17 ' to 74 • 44 ' E longitudes (Figure 1).The total area of the study region is 21.08 million hectares (mha).Two cropping seasons, summer (kharif) (May -Oct) and winter (rabi) (Nov -April) existed in the region.Wheat, rice, cotton and sugarcane are the cash crops contributing 2.2%, 0.7%, 1.4% and 0.7% respectively to the GDP (Gross Domestic Product) of the country (Usman, 2016).Wheat is the main winter crop with barley, oats and gram as minor crops.Rice and cotton are main summer crops with maize, sorghum, millet and sesamum as minor crops.There is no district level data available for minor crops therefore; minor crops referred as fodder crops in this study.Sugarcane is an annual crop that takes 10-12 months to get mature.Three dominate crop rotation, wheat-rice, wheat-cotton and wheat-fodder dominates the study region.

Crop area statistics
In Pakistan, final estimate of area sown under all major crops is collected through complete enumeration of all villages by field officers of revenue department twice a year (PBS).These acreage statistics compile at district level and publish by agriculture department in provincial and national statistical reports.Area statistics of crop calendar 2010-11 was used for this analysis.

MODIS data
Composite data product of Aqua and Terra satellite's MODIS sensor, MOD/MYD13A2 and MCD15A2, was used to get temporal information of NDVI and LAI respectively.Data has a spatial resolution of 1000 m with temporal resolution of 8 days.Total 46 images complete the crop calendar of 2010-11.LAI algorithm adopts radiative transfer approach by utilizing three, four or even seven MODIS Reflectance bands and some ancillary data on surface characteristics such as land cover type to quantify canopy structure (Knyazikhin, 1999).Whereas, NDVI is a normalized transform of the NIR to red reflectance ratio, designed to standardize vegetation index values between −1 and +1 (Didan et al., 2015).

Methodology
Two different phenological datasets NDVI, LAI and one derived dataset by multiplying NDVI with LAI (NDVI*LAI) were used to identify different crop rotation in study region based on temporal behavior of different pixels.Complete year phenology response for each dataset was stacked together for further analysis.An unsupervised to supervised approach was adopted for classification.Unsupervised classification was done using Iterative Self-Organizing Data Analysis Technique (ISODATA) clustering method to identify different temporal spaces.At first, each dataset was divided into two broad classes of vegetative and non-vegetative area.Vegetative area was further classified into 100 unsupervised classes.Temporal values of NDVI, LAI and (NDVI*LAI) for each class was extracted using zonal mean statistics and plotted against time.To make resulting curves smooth, a three period running average was taken so that beginning, peak and end of season in each curve could be identified easily.As a first step, different curves showing similar phenological behavior and occurred at same location were merge together.Crop calendar developed by Agriculture Information System (AIS) Pakistan was a useful information to design start and end of crop season in phenology matrix for different crop rotations in the region.Crop calendar can be found here: http://dwms.fao.org/~test/downs/docs/pak_crop_calendar.pdf.Seasonal duration and amplitude of phenology response were the key to find similar curves of different crop rotation.Supervised classification of different clusters were done using phenology information and expert knowledge of cropping rotation existed in area.After merging, clusters were reduced from 100 to 13, depicting different crop rotations.Flowchart of methodology is shown in Figure 2.After classification, area of major crops in each district were computed for three datasets.Estimated area of major crops per district was compared with area mentioned in statistical reports.Two statistical parameters, root mean square error (RMSE) and mean absolute percentage error (MAPE) were used to check the classification accuracy.Smaller the value of statistical parameter better is the estimation.The RMSE and MAPE can be calculated as: (2) Where   = crop area estimated from classification   = actual area in statistical reports  = total number of districts.

Crop rotation phenology
Water availability to crops have significant effect on photosynthesis.Water stress conditions would reduce chlorophyll levels in plant leaves ultimately effecting photosynthesis and leaf development (Sanchez et al., 1983).Decrease in chlorophyll fluorescence will reduce NDVI and LAI values (Springer et al., 2017).Wheat, rice, cotton and sugarcane are cash crops for study area, which not only fulfill food and fabric needs but also a major source of economic benefit to grower.Therefore, these cash crops are cultivated in suitable growing conditions with ample supply of water and nutrient to get maximum yield.It is clear from Figure 3 that phenology response in all three datasets (NDVI, LAI and NDVI*LAI) for cash crops is higher than other minor crops of each season.
In all datasets, wheat crop has high phenology response than other minor crops of rabi season.Rice, Cotton and sugarcane curve response is also higher than maize and other kharif fodder crops.Cotton crop phenology curve is more elongated than rice crop because of longer growth period.Wheat-Rice crop rotation has the highest phenology curve response in all three datasets followed by Wheat-Cotton (in NDVI) and Fodder-Cotton (in LAI and NDVI*LAI).Wheat-Bare crop cycle has unusual high curve response in LAI dataset that can be explained by influence of adjacent non-vegetative (water/sand) pixels.Wheat-Bare rotation is mostly existed in area closer to rivers in flood plains.
During rabi season, reduced river discharge leaves behind high moisture and nutrient soil to be cultivated for only one season.In rabi season, for all major crop rotations, wheat crop has maximam NDVI, LAI and LAI*NDVI value range of (0.65-   It is also observed from Figure 3 that phenology response, in all datasets, of kharif crops is higher and more close to each other that rabi crops.This is because of high soil moisture, by either rainfall (monsoon events) or sufficient canal water supply due to increase in river flow caused by upstream snow melting.Base value before the start of crop growth is high in NDVI phenology matrix than in LAI and (NDVI*LAI).This high value is due to soil reflectance in NDVI dataset whereas, LAI base value is close to zero for non-vegetative area.This high base value leaves less value difference between start (tillering stage) and peak (full bloom) of phenology curve making curve separability difficult in NDVI profiles.On the other hand, curve separability is high in LAI and (NDVI*LAI) making crop rotation identification much easy.High base value in phenology matrix also makes difficult to separate desert shrubs from fodder crops.During curve analysis, it was also observed that crop clusters in Punjab region have overall high phenological behavior in all three datasets.It can be attributed to favorable crop growing conditions than other parts of study region.

Classification validation
Comparison of estimated wheat crop area with reported area reveals that derived dataset (NDVI*LAI) perform better with RMSE of 34.55 and MAPE of 24.56%.While, NDVI and LAI achieve RMSE of 38.35 and 35.37 and MAPE of 28.16% and 24.87% respectively (Figure 4a).Overall, NDVI and LAI has overestimated (↑) the wheat area by 9% and 2% respectively.Whereas, NDVI*LAI underestimate (↓) wheat area by 0.2%.
Average maximum error of 73%↑ and 43%↓ is estimated for Qamber Shadatkot and Jacobabad districts respectively.
For sugarcane crop cover, LAI dataset performs better by attaining RMSE and MAPE of 8.60 and 34.58%.Whereas, RMSE of 10.08 and 10.83 and MAPE of 40.53% and 39.45% was obtained by NDVI and NDVI*LAI datasets respectively (Figure 4d).Overall, each dataset exaggerated in estimating sugarcane area with least deviation of 1%↑ by LAI, followed by NDVI (2%↑) then NDVI*LAI (18%↑).Average maximum fluctuation of 68%↑ and 34%↓ is estimated for Sanghar and Mandi Bahauddin district respectively.

Spatial variation of crop rotation
Three provinces has wheat-rice rotation with, 86% in Punjab, 9% in Sindh and 4% in Balochistan.A total area of 1530 (000, ha) (11%) is under wheat-rice rotation of this irrigated basin.Nasirabad, Gujranwala and Larkana are major areas of wheatrice rotation in study region.Wheat-cotton rotation is only dominate in Punjab (78%) and Sindh (22%).Bahawalnagar and Sanghar are main districts of each province to have highest area cultivated under wheat-cotton rotation.A total area of 2503 (000, ha) (18% of study region) is under wheat-cotton rotation.Wheatfodder rotation distribution is 89% in Punjab, 9% in Sindh and 2% in KPK.Sargodha, Khairpur and D.I.Khan are the areas of dominant wheat-fodder rotation.A total area of 2723 (000, ha) (19%) is under wheat-fodder rotation.Sugarcane is cultivated on 871 (000, ha) (6%) land in which 68% is cultivated in Punjab, 31% in Sindh and 1% in KPK.Rahim yar khan, Badin and D.I.Khan has the highest area sown under sugarcane.Fraction of major crop rotations within each province is shown in Figure 5. Punjab has major crop rotation wheat-fodder (25%).Sindh and KPK has fodder-fodder crop rotation occupying 26% of cropland.
Figure 5. Fraction of different crop rotation in each province.
Pixel based distribution of different crop rotation in study region, derived from three phenological datasets, is shown in Figure 6.It is evident that mixing of bare soil and built-up with cropland is minimized by incorporating LAI with NDVI.This mixing was the main cause of overestimation different crop areas when using only NDVI dataset.LAI alone do not produce good result in areas with sparse cropland especially in desert or in delta regions.NDVI dataset also makes it difficult to distinguish fodder from shrubs due to high base value.This can be avoided by fusing NDVI and LAI together.

CONCLUSION
Individual analysis of each dataset reveals that LAI performs better in mapping land feature that contain high biomass (e.g., sugarcane, forest or orchards).Rice crop nursery is transplanted under standing water and NDVI dataset does contain high moisture information.Therefore, NDVI can be a good temporal phenological information to map single season rice crop.Nevertheless, analysis of crop rotation mapping concludes that fusion of NDVI and LAI perform better for mapping of wheatrice, wheat-cotton and wheat-fodder rotation.Derived data (NDVI*LAI) not only exclude non-vegetative land features but also facilitates in subdividing different crop cluster during unsupervised classification without image segmentation that cannot be achieved using NDVI and LAI individually.This split in clusters helps to classify different crop rotations more accurately.Combined data reduce the base value in phenology matrix and increase curve separability making crop identification easy and efficient.Phenology curve analysis reveals that cash crop rotations have high curve response indicating favorable crop growing conditions.These conditions are generally concentrated in Punjab province.Crop area analysis shows that Punjab has 71% irrigated cropland while Sindh, KPK and Balochistan has 25%, 3% and 2% of study region respectively.Cash crop rotation analysis indicates that wheat-cotton is the main cash crop rotation in Punjab and Sindh.KPK has wheatfodder as dominant cash crop rotation with maize as major fodder crop.Whereas, wheat-rice is major cash crop rotation in Balochistan.Sindh has highest fraction (18%) of sugarcane-cultivated area with respect to area sown under cash crops followed by KPK (17%) and Punjab (9%).
Major crop rotations of an area is a key information for planning of food security, resource management and rural development.
Water and nutrient demands vary for different crops.Therefore, crop rotation mapping helps policy analysts and managers to articulate better plans for effective utilization of limited resources.This data fusion offers a robust method to map crop rotation, without image segmentation, in a large and diverse irrigated basin having more than one cropping seasons.

Figure 1 .
Figure 1.Geographic location of study area Curves with two peaks indicate two cropping season.Curves with one peak represent area of single cropping season with other season as fallow.Multiple peak curves represent area with fodder cultivation that give multi cuts during cultivation period.Curves with one elongated peak represent sugarcane crop.Whereas, curves with almost constant temporal behavior above the base value represents forest or orchard land features.

Figure 3 .
Figure 3. Phenological behavior of different crop rotation using NDVI, LAI and NDVI*LAI.

Figure 3 .
Figure 3. Phenological behavior of different crop rotation using NDVI, LAI and NDVI*LAI.

Figure 4 .
Figure 4. Comparison of estimated crop area with reported area of (a) wheat, (b) rice, (c) cotton and (d) sugarcane.

Figure 6 .
Figure 6.Spatial distribution of different crop rotation in study region derived from NDVI, LAI and NDVI*LAI.