Monitoring Changes in Croplands Due to Water Stress in the Krishna River Basin Using Temporal Satellite Imagery

Remote sensing-based assessments of large river basins such as the Krishna, which supplies water to many states in India, are useful for operationally monitoring agriculture, especially basins that are affected by abiotic stress. Moderate-Resolution Imaging Spectroradiometer (MODIS) time series products can be used to understand cropland changes at the basin level due to abiotic stresses, especially water scarcity. Spectral matching techniques were used to identify land use/land cover (LULC) areas for two crop years: 2013–2014, which was a normal year, and 2015–2016, which was a water stress year. Water stress-affected crop areas were categorized into three classes—severe, moderate and mild—based on the normalized difference vegetation index (NDVI) and intensity of damage assessed through field sampling. Furthermore, ground survey data were used to assess the accuracy of MODIS-derived classification individual products. Water inflows into and outflows from the Krishna river basin during the study period were used as direct indicators of water scarcity/availability in the Krishna Basin. Furthermore, ground survey data were used to assess the accuracy of MODIS-derived LULC classification of individual year products. Rainfall data from the tropical rainfall monitoring mission (TRMM) was used to support the water stress analysis. The nine LULC classes derived using the MODIS temporal imagery provided overall accuracies of 82% for the cropping year 2013–2014 and 85% for the year 2015–2016. Kappa values are 0.78 for 2013–2014 and 0.82 for 2015–2016. MODIS-derived cropland areas were compared with national statistics for the cropping year 2013–2014 with a R2 value of 0.87. Results show that both rainfed and irrigated areas in 2015–2016 saw significant changes that will have significant impacts on food security. It has been also observed that the farmers in the basin tend to use lower inputs and labour per ha during drought years. Among all, access to water is the major driver determining the crop choice and extent of input-use in the basin.


Introduction
Spatial information on the distribution of irrigated areas is limited to district level crop statistics that are published by the state or national governments in different parts of the world.Although data that has been collected by both irrigation and agriculture departments are often available, there are often differences in the extent of irrigated area reported between the sources [1,2].Water availability in major command areas (areas which can be irrigated under a reservoir and are fit for cultivation) vary frequently due to rainfall fluctuations which leads to stress during critical crop growth stages.Major command areas of the Krishna Basin often suffer from inequitable distribution of water due to overuse in head reaches, which is partly caused by farmers' preference for water-intensive crops such as rice and sugarcane [3,4].Decisions on the allocation of water in the river basin and command area require timely availability of location specific information.Monitoring and understanding LULC is important for planning to overcome problems, such as loss of agriculture land, decreasing environmental quality, water resources management etc. [5].Land use, as defined here, is how people use the landscape whether for development, conservation or mixed uses (e.g., agriculture, built-up lands).Land use data is important to analyze environmental processes and problems for sustainable development.Land cover is how much area is covered by naturally occurring features such as water bodies and natural vegetation (forest, grass lands, shrublands).Continuous monitoring of changes in land use and water use in the command areas can be cumbersome, labor intensive, and costly.Satellite images available in the public domain at coarse and medium spatial resolution (such as AVHRR, MODIS, Landsat and ASTER) are proving to be a viable option to continuously monitor land use [6][7][8].
Accurate information on the basin-wise extent of cropland is critical for food security assessment, water allocation decisions and yield estimation.This information will also help decision makers monitor the dynamic changes in landscapes, such as agricultural lands, fallow croplands and land cover such as forests, water bodies and wetlands.Ex-ante assessments of the impacts of different land use changes will facilitate socially, economically and ecologically sustainable land use planning.Moreover, the departments of agriculture and revenue within countries will need such spatial information at the village level to send advisories to farmers on timely access to inputs and adopt best agricultural management practices.
Croplands in South Asia are frequently stressed due to abiotic stress such as recurrent droughts.The year 2015 was declared the hottest year on record by the World Meteorological Organization (WMO).Very high temperatures over both land and ocean in 2015 were accompanied by many extreme weather events such as heatwaves, flooding and severe drought [9].Several studies have been conducted on land use mapping and changes using different satellite imagery.However, the main purpose of analyzing changes in agricultural land use is to monitor shifts in cropping patterns and cropland changes [10][11][12][13][14][15][16].These analyses rely heavily on agricultural statistics (e.g., area extent).Besides being reported at the district and state or provincial level, there are often discrepancies in the statistics reported by agricultural and irrigation agencies.Variations in land use at such a large scale are insufficient to fully understand their effect on the river basin.On the other hand, remote sensing with satellite imagery can give detailed maps of land use and identify where cropping pattern changes have occurred significantly in response to variations in rainfall [17,18] and other parameters.Remote sensing utilizing satellite imagery has been used to quantify water use and productivity in irrigation systems [18], but is less frequently used to identify how irrigated area changes in command areas in response to variations in rainfall and water supply.Time series data of the normalized difference vegetation index (NDVI) have been used for mapping land use changes [19] and irrigated areas [20,21].Time series data have also been used for detecting changes in irrigated areas in major river basins [4,22].
Given this background, the present study analyzed land use changes in the Krishna Basin including source wise irrigated areas.Water deficit for rainfed crops is also understood based on low rainfall which is reflected in the NDVI imagery during 2013-2014 and 2015-2016.The area under land use for each year and the changes between years were estimated.Both land use and land cover change estimates were compared against ground survey data and secondary sources such as published statistics on rice area.The study focused on the areas where significant changes had occurred in the cropping pattern.

Study Area
The Krishna River Basin is India's fourth largest river basin and covers 258,948 km 2 of southern India, covering the states of Karnataka (113,291 km 2 ), Andhra Pradesh (31,638 km 2 ), Telangana (51,845 km 2 ) and Maharashtra (69,028 km 2 ) [3,23].The river Krishna originates in the Western Ghat mountains, flows east across the Deccan Plateau, and discharges into the Bay of Bengal.The Krishna has three main tributaries that drain from the northwest, west, and southwest (Figure 1).The climate is generally semi-arid, with some dry, sub-humid areas in the eastern delta and humid areas in the Western Ghats.Annual precipitation averages 780 mm and decreases gradually from 850-1000 mm in the Krishna Delta to 300-400 mm in the northwest, and then increases to >1000 mm in the Western Ghats (Figure 1).In the extreme western parts of the basin, the Western Ghats have high annual precipitation (1500-2500 mm).Most of the rainfall occurs during the monsoon from June to October.Crops are cultivated in three seasons: kharif during the monsoon (June to mid-December), rabi in the post-monsoon dry season (mid-December to March), and in the summer season (April and May).In the irrigated areas there is double cropping of rice and other grains, single cropping of sugarcane, chile, cotton, and fodder grass.Areas with less irrigation grow corn, sorghum and sunflower.Rainfed crops include grains (sorghum, millet), pulses (red and green gram, chickpea), and oilseeds (sunflower, groundnut).

Study Area
The Krishna River Basin is India's fourth largest river basin and covers 258,948 km 2 of southern India, covering the states of Karnataka (113,291 km 2 ), Andhra Pradesh (31,638 km 2 ), Telangana (51,845 km 2 )and Maharashtra (69,028 km 2 ) [3,23].The river Krishna originates in the Western Ghat mountains, flows east across the Deccan Plateau, and discharges into the Bay of Bengal.The Krishna has three main tributaries that drain from the northwest, west, and southwest (Figure 1).The climate is generally semi-arid, with some dry, sub-humid areas in the eastern delta and humid areas in the Western Ghats.Annual precipitation averages 780 mm and decreases gradually from 850-1000 mm in the Krishna Delta to 300-400 mm in the northwest, and then increases to >1000 mm in the Western Ghats (Figure 1).In the extreme western parts of the basin, the Western Ghats have high annual precipitation (1500-2500 mm).Most of the rainfall occurs during the monsoon from June to October.Crops are cultivated in three seasons: kharif during the monsoon (June to mid-December), rabi in the post-monsoon dry season (mid-December to March), and in the summer season (April and May).In the irrigated areas there is double cropping of rice and other grains, single cropping of sugarcane, chile, cotton, and fodder grass.Areas with less irrigation grow corn, sorghum and sunflower.Rainfed crops include grains (sorghum, millet), pulses (red and green gram, chickpea), and oilseeds (sunflower, groundnut).

Satellite Data
MODIS data for the Krishna River Basin was available from the data archive at the NASA-USGS website (http://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.005).MODIS 2013-2014 to 2015-2016 every 16 days (Table 1) Terra sensor data were used for present study.It has two specific bands (band 1 (red) and band 2 (near infrared)) are processed for land applications as a MODIS vegetation index product (MOD13Q1).MOD13Q1 is computed from MODIS level 5 bands 1-2 (centered at 648 nm, and 858 nm).

Satellite Data
MODIS data for the Krishna River Basin was available from the data archive at the NASA-USGS website (http://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.005).MODIS 2013-2014 to 2015-2016 every 16 days (Table 1) Terra sensor data were used for present study.It has two specific bands (band 1 (red) and band 2 (near infrared)) are processed for land applications as a MODIS vegetation index product (MOD13Q1).MOD13Q1 is computed from MODIS level 5 bands 1-2 (centered at 648 nm, and 858 nm).MODIS imagery was used to map the spatial extent of LULC for the years 2013-2014 and 2015-2016.The process began with rescaling 16-day NDVI images, and these were later stacked into a single file data composite for each cropping year [24][25][26][27].MODIS 16-day composites were converted to NDVI monthly maximum value composite (MVC) (NDVI MVC) using Equation (1), where MVCi is the monthly MVC of the ith month and i1 and i2 are every 16-day data in a month: In the present study, monthly NDVI MVC were used for classification and a NDVI 16-day data set was used for identifying and labelling LULC classes including irrigated areas.The main purpose of using MVC was to avoid noise (cloud) in some of the areas [28].

Ground Survey Data
Ground survey data was collected during 13-26 September 2013 for 227 sample points, 21-30 September 2015 for 243 sample points and 19-21 October 2015 for 47 sample points covering about 7500 km of road travel in the Krishna River Basin (Figure 2).In 2013, out of 227 points, 73 were collected for class identification and labeling and the remaining 154 points were used to assess accuracy.Similarly, in 2015, out of 243 points, 81 were used for class identification and labeling and the remaining 162 points were used to assess accuracy.Further, an additional 47 points were used to assess accuracy for the 2015-2016 classification.Based on pre-classified output, Google Earth imagery and GPS tracking attached to the image processing software captured ground survey information while moving on the road.Detailed information was collected for class identification and labeling point locations.Point specific information was collected from 250 m × 250 m plots and consisted of GPS locations, land use categories, land cover percentages, cropping pattern during different seasons (through farmer interviews), crop types, and watering method (irrigated, rainfed).Samples were obtained within large contiguous areas of a particular LULC.Landsat 8 products were used as additional ground survey information in class identification.A stratified-systematic sample design was adopted based on road network or footpath access.Wherever possible, a systematic location of sites was done every 5 or 10 km along the road network by vehicle or on foot [27,29,30], which is detailed in a description of the ground survey methodological approach.

National Statistics & Secondary Data
Statistics on cultivated area, production, and price were obtained at the sub-national level (districts) from India-stat [31].
Daily rainfall data for the study area were downloaded from the tropical rainfall measuring mission (TRMM) [32] and processed as monthly rainfall to compare with 11 years (from 2005 to 2016) of monthly mean rainfall data, spatial resolution is 0.25 degrees × 0.25 degrees (covers ~27 km).The TRMM estimates are near real-time grids from which average monthly rainfall data for each month in the rainy season was extracted for the selected 11 years, including for the 2002 drought year.Due to the lack of official statistics for drought-affected areas [33], the rainfall data was not used to identify drought areas but as an accumulation of evidence of water stress occurrence.

Mapping Land Use/Land Cover
An overview of the methods (Figure 3) and details are described below.The process began with mapping LULC areas using MODIS 16-day time series data with spectral matching techniques and field-plot information.

National Statistics & Secondary Data
Statistics on cultivated area, production, and price were obtained at the sub-national level (districts) from India-stat [31].
Daily rainfall data for the study area were downloaded from the tropical rainfall measuring mission (TRMM) [32] and processed as monthly rainfall to compare with 11 years (from 2005 to 2016) of monthly mean rainfall data, spatial resolution is 0.25 degrees × 0.25 degrees (covers ~27 km).The TRMM estimates are near real-time grids from which average monthly rainfall data for each month in the rainy season was extracted for the selected 11 years, including for the 2002 drought year.Due to the lack of official statistics for drought-affected areas [33], the rainfall data was not used to identify drought areas but as an accumulation of evidence of water stress occurrence.

Mapping Land Use/Land Cover
An overview of the methods (Figure 3) and details are described below.The process began with mapping LULC areas using MODIS 16-day time series data with spectral matching techniques and field-plot information.
A time series of MODIS 16-day composite vegetation index images at 250 m resolution was obtained for 1 June 2013 to 31 May 2014 and 1 June 2015 to 31 May 2016 (MOD13Q1 data product).The 16-day composite images in the MOD13Q1 dataset are available in the public domain and are pre-calibrated (http://modis-sr.1tdri.org/html).The large scene size and daily overpass rate of MODIS makes it attractive for mapping large crop areas, and NDVI derived from MODIS has high fidelity with biophysical parameters [34].The 16-day NDVI images were stacked into a 23-band file for each crop year (two images per month).The monthly maximum value composites were created using 16-day NDVI MODIS data to minimize cloud effects [35].The 16-day composite images in the MOD13Q1 dataset are available in the public domain and are pre-calibrated (http://modis-sr.1tdri.org/html).The large scene size and daily overpass rate of MODIS makes it attractive for mapping large crop areas, and NDVI derived from MODIS has high fidelity with biophysical parameters [34].The 16-day NDVI images were stacked into a 23-band file for each crop year (two images per month).The monthly maximum value composites were created using 16day NDVI MODIS data to minimize cloud effects [35].
Unsupervised classification as described by Cihlar et al. 1998 was used to generate initial classes.The unsupervised ISOCLASS cluster algorithm (ISODATA in ERDAS Imagine 2014 TM ) run on the NDVI MVC generated an initial 40 classes, with a maximum of 40 iterations and convergence threshold of 0.99.Though ground survey data was available at the time of image classification, unsupervised classification was used in order to capture the full range of NDVI over a large area.Use of unsupervised techniques is recommended for large areas that cover a wide and unknown range of vegetation types, and where landscape heterogeneity complicates identification of homogeneous training sites [20,30,36].Identification of training sites is particularly problematic for small, heterogeneous irrigated areas.
LULC classes were identified based on NDVI temporal signatures along with ground survey data.We observed crop growth stages and cropping pattern from temporal signatures, such as (a) onset of cropping season (e.g., monsoon and winter); (b) duration of cropping season such as monsoon and winter; (c) magnitude of crops during different seasons and years (e.g., water stress and normal years); and (d) end of cropping season [26,30].
The process of labeling class identification was done based on spectral matching techniques (SMTs) [8,37].Initially, 40 classes from the unsupervised classification were grouped based on spectral similarity or closeness of class signatures.Each group of classes was matched with ideal spectra signatures and ground survey data and assigned class names [8,28].Classes with similar NDVI time series and land cover were merged into a single class.Classes showing significant mixing,   [20,30,36].Identification of training sites is particularly problematic for small, heterogeneous irrigated areas.
LULC classes were identified based on NDVI temporal signatures along with ground survey data.We observed crop growth stages and cropping pattern from temporal signatures, such as (a) onset of cropping season (e.g., monsoon and winter); (b) duration of cropping season such as monsoon and winter; (c) magnitude of crops during different seasons and years (e.g., water stress and normal years); and (d) end of cropping season [26,30].
The process of labeling class identification was done based on spectral matching techniques (SMTs) [8,37].Initially, 40 classes from the unsupervised classification were grouped based on spectral similarity or closeness of class signatures.Each group of classes was matched with ideal spectra signatures and ground survey data and assigned class names [8,28].Classes with similar NDVI time series and land cover were merged into a single class.Classes showing significant mixing, e.g., continuous irrigated areas and forest, were masked and reclassified using the same ISOCLASS algorithm.Some continuous irrigated areas mixed with forests in the Western Ghats were separated using a 90 m digital elevation model (DEM) from the shuttle radar topography mission (SRTM) and an elevation threshold of 630 m, Landsat imagery and ground survey data through spatial modeling techniques such as overlay matrix, recode and proximity analysis [3,38].This resulted in nine classes of LULC.While class aggregation could have been performed statistically using a Euclidean or other distance measure, we employed a user-intensive method that incorporates both ground survey data and high resolution imagery in order to avoid lumping classes that might be spectrally similar but have distinct land cover.The NDVI of some classes differed in only one or two months, which would have caused the classes to be merged if an automated similarity index were used.
Land 2017, 6, 72 7 of 18 Classification was done at 250 m spatial resolution.In the present study area, average land holding size is less than a pixel and there are different LULC classes with 250 m × 250 m pixel (6.25 ha).Full pixel areas are not an accurate representation of actual cropland areas.Cropland fraction was calculated using the methodology described in [3,26,37,39].Sub-pixel areas are important when a particular pixel was named as cropland but also contained different other LULC classes (e.g., grasses, trees, shrubs, etc.).
Ground survey points were used to assess the accuracy of the classification results, based on a standard procedure [40][41][42], to generate an error matrix and accuracy measures for each LULC map.Error matrices and Equation (2) 'Cohen's kappa coefficient (κ)' are commonly used for accuracy assessment.For example, these are useful when building models that predict discrete classes or when classifying imagery.κ can be used as a measure of agreement between model predictions and reality [43] or to determine if the values contained in an error matrix represent a result significantly better than random [42].κ is computed as: where, N is the total number of sites in the matrix, r is the number of rows in the matrix, x ii is the number in row i and column i, x +i is the total for row i, and x i+ is the total for column i [42].

Water Stress Mapping and Categorization
Water stress was measured based on NDVI signatures and water stress areas were classified into three categories based on intensity/crop condition observed from field survey data.Severe water stress areas were categorized as areas with no crop throughout the cropping season due to water shortage.Moderate and mild water stress areas were identified where a fraction of the plot was used to grow rice or partial damage had occurred and the rest was left fallow.This information was collected from 205 fields out of 584 locations.
NDVI signatures were extracted from the NDVI stack for the years 2013 and 2014 for the ground survey locations (conducted in 2015).Out of the 243 locations considered, 81 points were selected based on field plot information and spectral similarity with water stress.Farmers' responses were used along with the NDVI signatures to identify ideal spectral signatures (temporal signatures) for water stress and non-water stress areas.During the ground survey, we observed that the year 2015-2016 was one of the severe water stress year and large agricultural areas including irrigated command areas were left fallow.

Results
This section discusses LULC areas including major croplands, changes in irrigated area, water stress areas and accuracy assessment based on ground survey data and a comparison of irrigated areas from the present study and national statistics.

Spatial Distribution of LULC
The 16-day time series NDVI stack and the spectral signatures (temporal signatures) generated from it for sampled locations and class signatures obtained after classification were used to identify LULC.Nine LULC classes were mapped for the normal year (2013-2014; Figure 4 and Table 2) and water-deficit years (2015-2016; Figure 5 and Table 2).Classes identified were based on spectral matching techniques along with ground survey data include GPS referenced digital images and field observations.The cropland areas in the Krishna River Basin for 2013-2014 are shown in Table 2. Rainfed croplands were 4.8 Mha of full pixel area, which is dominated by rainfed croplands mixed with other land cover areas.In Class 1, the actual rainfed agriculture area was 4.1 Mha, and remaining 0.7 Mha was mixed with other LULC areas.The total irrigated area was 9.5 Mha which includes groundwater, tank irrigation, and major canal irrigation during 2013-2014.Groundwater irrigated areas alone covered 5.1 Mha, mainly located across the Krishna River Basin and major command areas.
matching techniques along with ground survey data include GPS referenced digital images and field observations.The cropland areas in the Krishna River Basin for 2013-2014 are shown in Table 2. Rainfed croplands were 4.8 Mha of full pixel area, which is dominated by rainfed croplands mixed with other land cover areas.In Class 1, the actual rainfed agriculture area was 4.1 Mha, and remaining 0.7 Mha was mixed with other LULC areas.The total irrigated area was 9.5 Mha which includes groundwater, tank irrigation, and major canal irrigation during 2013-2014.Groundwater irrigated areas alone covered 5.1 Mha, mainly located across the Krishna River Basin and major command areas.Note: GW = groundwater; SW = surface water; DC = double crop.The spatial extent of LULC in the Krishna River Basin for 2015-2016 is shown in Figure 5 and areas shown in Table 2. Altogether nine classes were delineated, of which four were croplands and five were other LULC.Rainfed croplands covered 2.8 Mha while total irrigated area was 7.6 Mha, which included groundwater irrigated, tank irrigated and major canal irrigated areas during 2015-2016.Groundwater irrigated area alone covered 4.0 Mha which was spatially distributed across the Krishna River Basin including major command areas.
The results show a significant decrease in total irrigated area (including surface and groundwater areas) in the Krishna Basin between 2013-2014 (9.The spatial extent of LULC in the Krishna River Basin for 2015-2016 is shown in Figure 5 and areas shown in Table 2. Altogether nine classes were delineated, of which four were croplands and five were other LULC.Rainfed croplands covered 2.8 Mha while total irrigated area was 7.6 Mha, which included groundwater irrigated, tank irrigated and major canal irrigated areas during 2015-2016.Groundwater irrigated area alone covered 4.0 Mha which was spatially distributed across the Krishna River Basin including major command areas.
The results show a significant decrease in total irrigated area (including surface and groundwater areas) in the Krishna Basin between 2013-2014 (9.due to lower water storage in the reservoir during 2015-2016.Figure 6 shows the spatial extent of crop land changes from crop lands converted to fallows due water stress.

Water Scarcity and Cropland Changes
Agricultural land uses decreased from 14.7 Mh in a normal year (2013-2014) to 11.5 Mha during the water stress year (2015-2016) (Figure 7 and Table 3).These changes are mainly due to reduced cropping intensity and conversion of crop lands to fallow due to water scarcity.Rainfed croplands in the Krishna Basin were reduced to half of the normal year to 2.8 Mha, and there was significant reduction in the irrigated croplands.Overall, 7.6 Mha of croplands were affected by water scarcity, including rainfed and irrigated croplands, and 7.1 Mha were not affected (Table 3).

Water Scarcity and Cropland Changes
Agricultural land uses decreased from 14.7 Mh in a normal year (2013-2014) to 11.5 Mha during the water stress year (2015-2016) (Figure 7 and Table 3).These changes are mainly due to reduced cropping intensity and conversion of crop lands to fallow due to water scarcity.Rainfed croplands in the Krishna Basin were reduced to half of the normal year to 2.8 Mha, and there was significant reduction in the irrigated croplands.Overall, 7.6 Mha of croplands were affected by water scarcity, including rainfed and irrigated croplands, and 7.1 Mha were not affected (Table 3).
Land 2017, 6, 72 10 of 17 due to lower water storage in the reservoir during 2015-2016.Figure 6 shows the spatial extent of crop land changes from crop lands converted to fallows due water stress.

Water Scarcity and Cropland Changes
Agricultural land uses decreased from 14.7 Mh in a normal year (2013-2014) to 11.5 Mha during the water stress year (2015-2016) (Figure 7 and Table 3).These changes are mainly due to reduced cropping intensity and conversion of crop lands to fallow due to water scarcity.Rainfed croplands in the Krishna Basin were reduced to half of the normal year to 2.8 Mha, and there was significant reduction in the irrigated croplands.Overall, 7.6 Mha of croplands were affected by water scarcity, including rainfed and irrigated croplands, and 7.1 Mha were not affected (Table 3).There were some visible changes in cropping patterns in 2015-2016 (Table 4).Long-duration irrigated cropland under sugarcane was converted to crops such as rice and maize (Figure 7, class 2).The area under double crop (monsoon rice + winter rice) also fell, with pulses included in the system.Intensively irrigated areas, which include groundwater irrigated areas, in 2013-2014 were under fallows in 2015-2016.Water scarcity in minor irrigation systems forced farmers to drill bore wells to irrigate.Note: GW = groundwater; SW = surface water; and DC = double crop.
A sudden reduction in food production could lead to unforeseen consequences in the market, leading to food insecurity.These results are especially relevant in a changing climate scenario and underline the need for appropriate adaptation strategies like growing water efficient crops during a water deficit year.Water scarcity is a frequent reality faced by farmers in the semi-arid tropics, even in the basin command areas.This leads to a chain of events affecting crop production, cropping pattern changes and eventually reduced in incomes to farmers.To maintain food security sustainably, large river basins like the Krishna should adopt new adaptation strategies such as the use of climate smart varieties, best bet management practices along with institutional mechanisms to remunerate the crops during such times.

Accuracy Assessment and Comparison with Other Published Data Sets
A quantitative accuracy assessment was done through error matrix [44] to examine whether a known LULC was identified as the same LULC or not.The ground survey data was based on an extensive field campaign conducted throughout the Krishna Basin during the kharif season for the crop years 2013-2014 and 2014-2015.Accuracy was performed on classified three products (LULC maps of 2013-2014 and 2015-2016).The 363 points that were not used in the classification were used to assess classification accuracy.Accuracy assessment was performed with independent datasets.
Tables 5-7 show the error matrix of each product.In LULC 2013-2014 of the first class (01.Rainfed croplands), ten ground survey reference points of the 14 matched with rainfed croplands; three points matched with 02.Irrigated-tanks/GW-croplands and one point matched with other LULC class.
The same assessment was carried out for all the classes shown in Table 5.For all the six classes, 125 points out of 152 matched with the same class of reference data.Accuracy for the final six classes of 2013-2014 was 82.4% with a Kappa value of 0.78 (Table 5).Table 6 shows the accuracy assessment for the six LULC classes for the year 2015-2016.Classification accuracy was performed with 162 sample points, of which 138 matched correctly with the present classification.The overall accuracy was 85% with a kappa value was 0.82 (Table 6).
The final MODIS-derived cropland areas statistics of the Krishna Basin districts were compared with district-wise statistics obtained from the Directorate of Economics and Statistics (DES), Andhra Pradesh and Telangana; Directorate of Economics and Statistics, Karnataka; and the of Agriculture, Maharashtra.This data was fractionalized based on the district-wise area covered in the Krishna Basin for a comparative study with the MODIS derived data.It revealed that most of the district level statistics from DES matched with the MODIS derived statistics.Difference in the statistical data and MODIS data varied between −35% and +35% (Figure 8) and R 2 value was 0.8686.Figure 9 illustrates the decrease in rainfall percentage from 17 years of normal rainfall to 2015 rainfall for monsoon and post-monsoon seasons (from June to November).Also, basin level inflows/outflows are a good indicator of water availability during crop production.Figure 10 clearly indicates that the inflows into and outflows from the sub-basins were at a very low level during 2015-2016 due to low rainfall [45,46].The changes in land use due to low rainfall and inflows are explained in Section 3.2.
Land 2017, 6, 72 12 of 17 assessment for the six LULC classes for the year 2015-2016.Classification accuracy was performed with 162 sample points, of which 138 matched correctly with the present classification.The overall accuracy was 85% with a kappa value was 0.82 (Table 6).
The final MODIS-derived cropland areas statistics of the Krishna Basin districts were compared with district-wise statistics obtained from the Directorate of Economics and Statistics (DES), Andhra Pradesh and Telangana; Directorate of Economics and Statistics, Karnataka; and the Department of Agriculture, Maharashtra.This data was fractionalized based on the district-wise area covered in the Krishna Basin for a comparative study with the MODIS derived data.It revealed that most of the district level statistics from DES matched with the MODIS derived statistics.Difference in the statistical data and MODIS data varied between −35% and +35% (Figure 8) and R 2 value was 0.8686.Figure 9 illustrates the decrease in rainfall percentage from 17 years of normal rainfall to 2015 rainfall for monsoon and post-monsoon seasons (from June to November).Also, basin level inflows/outflows are a good indicator of water availability during crop production.Figure 10 clearly indicates that the inflows into and outflows from the sub-basins were at a very low level during 2015-2016 due to low rainfall [45,46].The changes in land use due to low rainfall and inflows are explained in Section 3.2.

Discussion
This study identified the changes in irrigated areas including other cropland areas in a large river basin due to water scarcity.First, a baseline irrigated area map of the study area was produced for 2013-2014 with an estimation of cropland area under different classes.Accuracy was determined by correlating the MODIS-derived land use/land cover areas with field-plot information and subnational statistics obtained from the Ministry of Agriculture.MODIS imagery plays an important role in this type of study, where time series (composites of every 16 days) imagery not only helps in identifying a land use type by crop, based on its specific growing season, but also monitors the dynamics of such land use over time and space.The compositing of daily data to 16-day imagery has the advantage of minimal cloud contamination since the best of 16 days imagery is retained for analysis (MVC).The spectral matching technique used in this study is a very useful method to identify a specific land use where the spectral profile follows the phenology (Thenkabail et al., 2007).

Discussion
This study identified the changes in irrigated areas including other cropland areas in a large river basin due to water scarcity.First, a baseline irrigated area map of the study area was produced for 2013-2014 with an estimation of cropland area under different classes.Accuracy was determined by correlating the MODIS-derived land use/land cover areas with field-plot information and subnational statistics obtained from the Ministry of Agriculture.MODIS imagery plays an important role in this type of study, where time series (composites of every 16 days) imagery not only helps in identifying a land use type by crop, based on its specific growing season, but also monitors the dynamics of such land use over time and space.The compositing of daily data to 16-day imagery has the advantage of minimal cloud contamination since the best of 16 days imagery is retained for analysis (MVC).The spectral matching technique used in this study is a very useful method to identify a specific land use where the spectral profile follows the phenology (Thenkabail et al., 2007).

Discussion
This study identified the changes in irrigated areas including other cropland areas in a large river basin due to water scarcity.First, a baseline irrigated area map of the study area was produced for 2013-2014 with an estimation of cropland area under different classes.Accuracy was determined by correlating the MODIS-derived land use/land cover areas with field-plot information and sub-national statistics obtained from the Ministry of Agriculture.MODIS imagery plays an important role in this type of study, where time series (composites of every 16 days) imagery not only helps in identifying a land use type by crop, based on its specific growing season, but also monitors the dynamics of such land use over time and space.The compositing of daily data to 16-day imagery has the advantage of minimal cloud contamination since the best of 16 days imagery is retained for analysis (MVC).
The spectral matching technique used in this study is a very useful method to identify a specific land use where the spectral profile follows the phenology (Thenkabail et al., 2007).
This study applied spectral matching techniques introduced in the seminal paper by Thenkabail et al. (2007) for mapping irrigated areas using AVHRR time-series imagery along with intensive field-plot information.The purpose, focus, and advance in this paper is to map changes in cropping patterns in a region with two crop seasons per year.Mapping that spatial and temporal complexity in land use changes due to water stress is a significant new information.The methodology starts with unsupervised classification of NDVI MVCs.In the class identification and labeling process, we used 16-day as well as monthly MVCs.The main advantage of this approach during classification is the usage of monthly cloud-free or near-cloud-free images with consideration of maximum value composites.In addition, MODIS 16-day temporal resolution data easily identified cropping systems across different cropping pattern and changes (irrigated, rainfed, etc.).Spectral matching techniques were particularly successful in differentiating cropping patterns such as mild-water stress, moderate water-stress, severe water-stress areas.

Conclusions
The study highlighted the significant changes in agricultural land use that took place as a result of water stress during 2015-2016.Rainfed and irrigated areas were mapped with classification accuracy between 77 and 85% using MODIS 250 m time series images and spectral matching techniques (SMTs).The MODIS-based irrigated cropland statistics for the districts were highly correlated (R 2 value of 0.86) with the figures reported by the Directorate of Economics and Statistics.Though the current study focused on a large river basin, these methods and approaches are also applicable for large areas, such as countries and continents.
Estimates of the basin-level water inflow and outflow are useful in understanding the water balance in the basin.It also provides a realistic picture of the effects of water scarcity during a low rainfall year on crop production, and in turn on the food security of the population in the basin.The outcome of this research highlights the value of using MODIS time series 250 m data and advanced methods such as spectral matching to study changes in the agricultural cropland in large river basins.It also contributes significantly to the knowledge base of earth observation groups involved in monitoring irrigated areas.

Figure 1 .
Figure 1.Map of the study area with a location map.

Figure 1 .
Figure 1.Map of the study area with a location map.

Figure 2 .
Figure 2. Ground survey locations in the Krishna River Basin and representative photographs taken during the survey.

Figure 2 .
Figure 2. Ground survey locations in the Krishna River Basin and representative photographs taken during the survey.

Figure 3 .
Figure 3. Flow chart describing the analysis process.A time series of MODIS 16-day composite vegetation index images at 250 m resolution was obtained for 1 June 2013 to 31 May 2014 and 1 June 2015 to 31 May 2016 (MOD13Q1 data product).The 16-day composite images in the MOD13Q1 dataset are available in the public domain and are pre-calibrated (http://modis-sr.1tdri.org/html).The large scene size and daily overpass rate of MODIS makes it attractive for mapping large crop areas, and NDVI derived from MODIS has high fidelity with biophysical parameters[34].The 16-day NDVI images were stacked into a 23-band file for each crop year (two images per month).The monthly maximum value composites were created using 16day NDVI MODIS data to minimize cloud effects[35].Unsupervised classification as described byCihlar et al. 1998  was used to generate initial classes.The unsupervised ISOCLASS cluster algorithm (ISODATA in ERDAS Imagine 2014 TM ) run on the NDVI MVC generated an initial 40 classes, with a maximum of 40 iterations and convergence threshold of 0.99.Though ground survey data was available at the time of image classification, unsupervised classification was used in order to capture the full range of NDVI over a large area.Use of unsupervised techniques is recommended for large areas that cover a wide and unknown range of vegetation types, and where landscape heterogeneity complicates identification of homogeneous training sites[20,30,36].Identification of training sites is particularly problematic for small, heterogeneous irrigated areas.LULC classes were identified based on NDVI temporal signatures along with ground survey data.We observed crop growth stages and cropping pattern from temporal signatures, such as (a) onset of cropping season (e.g., monsoon and winter); (b) duration of cropping season such as monsoon and winter; (c) magnitude of crops during different seasons and years (e.g., water stress and normal years); and (d) end of cropping season[26,30].The process of labeling class identification was done based on spectral matching techniques (SMTs)[8,37].Initially, 40 classes from the unsupervised classification were grouped based on spectral similarity or closeness of class signatures.Each group of classes was matched with ideal spectra signatures and ground survey data and assigned class names[8,28].Classes with similar NDVI time series and land cover were merged into a single class.Classes showing significant mixing,

Figure 3 .
Figure 3. Flow chart describing the analysis process.
Unsupervised classification as described by Cihlar et al. 1998 was used to generate initial classes.The unsupervised ISOCLASS cluster algorithm (ISODATA in ERDAS Imagine 2014 TM ) run on the NDVI MVC generated an initial 40 classes, with a maximum of 40 iterations and convergence threshold of 0.99.Though ground survey data was available at the time of image classification, unsupervised classification was used in order to capture the full range of NDVI over a large area.Use of unsupervised techniques is recommended for large areas that cover a wide and unknown range of vegetation types, and where landscape heterogeneity complicates identification of homogeneous training sites

Figure 6
shows the spatial extent of crop land changes from crop lands converted to fallows due water stress.Land 2017, 6, 72 10 of 17

Figure 8 .
Figure 8.A comparison of district-wise cropland area from the MODIS classification and national statistics.

Figure 8 .
Figure 8.A comparison of district-wise cropland area from the MODIS classification and national statistics.

Figure 9 .
Figure 9. Decrease in rainfall percentage from 2013 to 2015 using the tropical rainfall measuring mission (TRMM) 0.25 degree rainfall grid.

Figure 9 . 17 Figure 9 .
Figure 9. Decrease in rainfall percentage from 2013 to 2015 using the tropical rainfall measuring mission (TRMM) 0.25 degree rainfall grid.

Table 1 .
MODIS-250 m Terra vegetation indices 16-day L5 product used in this study.

Table 3 .
MODIS-derived water stress areas across the Krishna River Basin.The table shows full-pixel area (FPA), crop area fraction (CAF), and sub-pixel area (SPA) or actual area.SPA = FPA × CAF.

Table 4 .
Changes in irrigated agricultural cropland from 2013-2014 (a normal year) to 2015-2016 (a water-deficit year) as a response to water stress, where FPA = full-pixel area; CAF = crop area fraction; SPA = and sub-pixel area or actual area.SPA = FPA × CAF.

Table 5 .
Accuracy assessments done using error matrix for LULC classes of 2013-2014.

Table 6 .
Accuracy assessments done using error matrix for LULC classes of 2015-2016.

Table 7 .
Accuracy assessments done using error matrix for water stress classes of 2015-2016.