Developing a Remote Sensing-Based Combined Drought Indicator Approach for Agricultural Drought Monitoring over Marathwada, India

The increasing drought severities and consequent devastating impacts on society over the Indian semi-arid regions demand better drought monitoring and early warning systems. Operational agricultural drought assessment methods in India mainly depend on a single input parameter such as precipitation and are based on a sparsely located in-situ measurements, which limits monitoring precision. The overarching objective of this study is to address this need through the development of an integrated agro-climatological drought monitoring approach, i.e., combined drought indicator for Marathwada (CDI_M), situated in the central part of Maharashtra, India. In this study, satellite and model-based input parameters (i.e., standardized precipitation index (SPI-3), land surface temperature (LST), soil moisture (SM), and normalized difference vegetation index (NDVI)) were analyzed at a monthly scale from 2001 to 2018. Two quantitative methods were tested to combine the input parameters for developing the CDI_M. These methods included an expert judgment-based weight of each parameter (Method-I) and principle component analysis (PCA)-based weighting approach (Method-II). Secondary data for major types of crop yields in Marathwada were utilized to assess the CDI_M results for the study period. CDI_M maps depict moderate to extreme drought cases in the historic drought years of 2002, 2009, and 2015–2016. This study found a significant increase in drought intensities (p ≤ 0.05) and drought frequency over the years 2001–2018, especially in the Latur, Jalna, and Parbhani districts. In comparison to Method-I (r ≥ 0.4), PCA-based (Method-II) CDI_M showed a higher correlation (r ≥ 0.60) with crop yields in both harvesting seasons (Kharif and Rabi). In particular, crop yields during the drier years showed a greater association (r > 6.5) with CDI_M over Marathwada. Hence, the present study illustrated the effectiveness of CDI_M to monitor agricultural drought in India and provide improved information to support agricultural drought management practices.


Introduction
Droughts are one of the most devastating natural hazards over the globe, affecting millions of individuals in multiple ways (e.g., food security, economic losses, and migration). Megacities such as In the United States, Wu and Wilhite [33] derived an operational drought monitoring technique for crops like soybean and corn. They considered the crop-specific drought index (CSDI) and SPI for their analysis. In recent years, the vegetation drought response index (VegDRI) was established and implemented over the continental US [32,34], which has also been developed in Korea [35]. Sepulcre-Canto et al. [36] developed an integrated agricultural drought Index for Europe by considering the SPI, SM, and photosynthetically active radiation (fAPAR) as input variables. An expert judgment-based integrated drought index was developed for Morocco [37]. During their research, SPI, LST, ET, and NDVI were utilized as contributing parameters to derive the new index. To support the humanitarian aid organizations, a weighted, dependent-combined drought indicator was also implemented over Ethiopia [38].
In Indian scenarios, Patel and Yadav [28] integrated the modified rainfall anomaly index (MRAI) and the VCI to understand crop loss and its association with droughts over the Bundelkhand region. Zhang et al. [39] developed a concurrent drought analysis method for India. This technique considered the SPI, SSI, SRI, and VCI to observe the changes in wheat production over various parts of India. On a country scale, Shah and Mishra [7] have executed the integrated drought monitoring approach and identified the highest severity drought years (1965, 1987, 2002, and 2015) over multiple river basins. Other researchers [40][41][42] have also reported the potential usefulness of applying an integrated drought index technique in India. Marathwada is among the severely drought-affected parts of India, where more than 1500 farmers committed suicide during the recent, severe drought years of 2015-16 [43]. Swain et al. [44]; Kulkarni et al. [16]; and Purandare [45] have studied frequent droughts in Marathwada using rainfall as an input variable for analyzing a specific drought year. Still, most of these studies did not combine multiple agro-climatological factors. Over the last few decades, the Marathwada semi-arid region of India is experiencing more frequent and extreme drought events, which have unprecedented water scarcity and agricultural crises [43]. The global drought monitoring studies have also projected increasing drought severities over India in the near future [46]. These studies have marked the semi-arid regions of India as one of the major worldwide hot spots to experience the more frequent and severe droughts. Therefore, a thorough understanding of drought scenarios over the Indian semi-arid regions has a local as well as global importance.
Most of the drought monitoring studies in India and elsewhere are still lingering with the single input parameters and worked on the country or subcontinent level. But at a regional scale, limited efforts have been taken; even though region-specific studies are crucial for developing policies and drought mitigation strategies. Many researchers have never implemented grid-based weight assignment methods to combine different input variables for a better assessment of droughts. Moreover, a single combining approach is widely used in previous studies, but limited endeavors have been made so far in testing two or more combining approaches at a time. Area-specific studies always have a crucial part in the policy management and designing the region-wise drought mitigation and adaptation plans, but most of the drought monitoring studies in India are still lagging. Therefore, the present study attempts to overcome the above-mentioned gaps in the current operational drought monitoring methods of India.
In this study, four primary drought indicator parameters (i.e., LST, SPI, NDVI, and SM) that are representing mainly the agricultural-field conditions (SM and NDVI), meteorological (SPI), and land atmospheric interaction (LST) were used to develop the combined drought index for Marathwada (CDI_M). Two approaches were also tested to combine these input parameters: (1) an expert judgment method (fixed weights, based on in-situ relationships and expert knowledge) and (2) a statistical (PCA) method. For the first time in India, the PCA technique is implemented through this study for the region-specific development of CDI. Previously, this approach was used by Hazaymeh and Hassan [47] and Bayissa et al. [48] over Jordan and Ethiopia, respectively. Implementing both the subjective (fixed weights) and objective (PCA) approaches and the comparative analysis of droughts results from these methods in this study stand as a unique contribution to the CDI-based drought assessment research. Hence, this newly developed CDI_M can be used for a comprehensive analysis of Remote Sens. 2020, 12, 2091 4 of 25 droughts and to better understand the relationship between major crop yields and drought events over Marathwada, India.

Study Area
The Marathwada meteorological subdivision is situated in the central and southeastern parts of Maharashtra state, India. This region falls under the tropical semi-arid zone, which is recently known for its farmers' distress and suicides in the severe drought years 2014-15 [43]. Administratively, the Marathwada region is known as the 'Aurangabad' division of Maharashtra state, which consists of eight districts, named as Osmanabad, Jalna, Parbhani, Nanded, Aurangabad, Beed, Hingoli, and Latur. The geographical latitudinal expansion of Marathwada is 17 • 35 N to 20 • 41 N, and the longitudinal extent is 74 • 40 E to 78 • 16 E (Figure 1). The entire administrative expanse of the Marathwada region is 64,590 sq. Km, which covers 21% of the overall spread of Maharashtra state. Physiologically, Marathwada is part of the Deccan plateau and is mostly covered with deep black cotton soil. This area consists of three main geographical regions based on its altitude variations. The northern and southern uplands of Marathwada are known as the Ajantha mountain ranges and Balaghat hills, whereas the central part contains the Godavari and Manjara river basins. Marathwada belongs to the rain shadow region of the Sahyadri mountain ranges, making this area vulnerable to rainfall variabilities.

Study Area
The Marathwada meteorological subdivision is situated in the central and southeastern parts of Maharashtra state, India. This region falls under the tropical semi-arid zone, which is recently known for its farmers' distress and suicides in the severe drought years 2014-15 [43]. Administratively, the Marathwada region is known as the 'Aurangabad' division of Maharashtra state, which consists of eight districts, named as Osmanabad, Jalna, Parbhani, Nanded, Aurangabad, Beed, Hingoli, and Latur. The geographical latitudinal expansion of Marathwada is 17 0 35′ N to 20 0 41′N, and the longitudinal extent is 74 0 40′E to 78 0 16′ E (Figure 1). The entire administrative expanse of the Marathwada region is 64,590 sq. Km, which covers 21% of the overall spread of Maharashtra state. Physiologically, Marathwada is part of the Deccan plateau and is mostly covered with deep black cotton soil. This area consists of three main geographical regions based on its altitude variations. The northern and southern uplands of Marathwada are known as the Ajantha mountain ranges and Balaghat hills, whereas the central part contains the Godavari and Manjara river basins. Marathwada belongs to the rain shadow region of the Sahyadri mountain ranges, making this area vulnerable to rainfall variabilities. The average normal rainfall in Marathwada is 826 mm, ranging from 700 mm to 850 mm among the districts. The majority of annual precipitation (80%) is reliant on the southwest monsoon rains, which occurs in June, July, August, and September (J, J, A, and S) [49]. Depending on these monsoon rains, rainfed agriculture has become the dominant agricultural system in this region. Therefore, a failure of the southwest monsoon rains directly or indirectly leads to devastating drought scenarios and agrigrain crises in Marathwada. In this area, the highest temperature (45 0 C) occurs in May, and the coldest temperature (11 0 C) is observed in December ( Figure 2) [49]. Kharif (sowing is at the starting of the monsoon) and Rabi (seeding is at the end of the rainy season) are two major cropping seasons in Marathwada. Cotton, jowar, bajra, wheat, maize, pulses, groundnut, and sugarcane are the main crops grown in this region. The total share of food grains is 52% to the gross cropped yield, followed by cotton, which contributes 38% in overall agricultural output [50]. Among food grains, jowar accounts for 25 to 40% of the total yield in each district of Marathwada. The variations in rainfall The average normal rainfall in Marathwada is 826 mm, ranging from 700 mm to 850 mm among the districts. The majority of annual precipitation (80%) is reliant on the southwest monsoon rains, which occurs in June, July, August, and September (J, J, A, and S) [49]. Depending on these monsoon rains, rainfed agriculture has become the dominant agricultural system in this region. Therefore, a failure of the southwest monsoon rains directly or indirectly leads to devastating drought scenarios and agrigrain crises in Marathwada. In this area, the highest temperature (45 • C) occurs in May, and the coldest temperature (11 • C) is observed in December ( Figure 2) [49]. Kharif (sowing is at the starting of the monsoon) and Rabi (seeding is at the end of the rainy season) are two major cropping seasons in Marathwada. Cotton, jowar, bajra, wheat, maize, pulses, groundnut, and sugarcane are the main crops grown in this region. The total share of food grains is 52% to the gross cropped yield, followed by cotton, which contributes 38% in overall agricultural output [50]. Among food grains, jowar accounts for 25 to 40% of the total yield in each district of Marathwada. The variations in rainfall Remote Sens. 2020, 12, 2091 5 of 25 and temperature follow the nonuniform relief features in this region. Approximately 70% expanse of the overall geographic area is in agricultural production, whereas only 3.55% Marathwada region is under natural vegetation [50].
Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 26 and temperature follow the nonuniform relief features in this region. Approximately 70% expanse of the overall geographic area is in agricultural production, whereas only 3.55% Marathwada region is under natural vegetation [50].

Data
The CDI_M developed in this study is based on four different remotely sensed and model-based input parameters. Monthly time-series data of two agricultural drought indicators, i.e., SM and NDVI, and two climatological parameters, i.e., SPI and LST, were used. For accurate spatial assessment of agricultural droughts, only the spatial extent of croplands was considered and extracted from the MODIS (MCD12Q1) Land Use Land Cover (LULC) datasets for Marathwada. This LULC datasets shows the human intervene land operation, and vegetations and man made constructions over the study region. CDI_M was compared with crop yield information from the study area to assess the agricultural drought monitoring capabilities of this method. The selection of specific input datasets was based on long-term data availability, reliability, high spatial resolution, accessibility, and use in the worldwide operational agro-climatic monitoring systems. Details of these datasets are presented below, and the detailed description of datasets and their sources is given in Table 1.

SPI (CHIRPS Data)
Rainfall deficit is an essential component of drought. Variations in rainfall and the south-north deviation of the intertropical conversion zone (ITCZ) alter the precipitation conditions and drought in particular [51]. Therefore, the commonly used rainfall-based SPI [52] was chosen as one of the meteorological components for this study. By normalizing the rainfall, the SPI was calculated on various time intervals, like 1, 2, and 3 months for the 2001 to 2018 study period. The SPI was computed using IMD-derived SPI ranges from higher negative values, representing extremely dry (-2.00 or less) conditions to the positive values of extremely wet (+2.00 or more) conditions. To maintain the seasonal relationship between rainfall and vegetation response, the SPI-3 (3-month SPI) was used for this study. The SPI values were calculated by using long-term (1981-2018) precipitation data. Climate Hazards Group Infrared Precipitation with Station data (CHIRPS) [53] generated by the U.S Geological Survey were incorporated to derive the time-series of SPI. CHIRPS data were used for this study, due to these main reasons: 1) data availability on a longer temporal scale (1981-present); 2) higher spatial resolution (0.05°) when compared to other available rainfall datasets; and 3) its

Data
The CDI_M developed in this study is based on four different remotely sensed and model-based input parameters. Monthly time-series data of two agricultural drought indicators, i.e., SM and NDVI, and two climatological parameters, i.e., SPI and LST, were used. For accurate spatial assessment of agricultural droughts, only the spatial extent of croplands was considered and extracted from the MODIS (MCD12Q1) Land Use Land Cover (LULC) datasets for Marathwada. This LULC datasets shows the human intervene land operation, and vegetations and man made constructions over the study region. CDI_M was compared with crop yield information from the study area to assess the agricultural drought monitoring capabilities of this method. The selection of specific input datasets was based on long-term data availability, reliability, high spatial resolution, accessibility, and use in the worldwide operational agro-climatic monitoring systems. Details of these datasets are presented below, and the detailed description of datasets and their sources is given in Table 1.

SPI (CHIRPS Data)
Rainfall deficit is an essential component of drought. Variations in rainfall and the south-north deviation of the intertropical conversion zone (ITCZ) alter the precipitation conditions and drought in particular [51]. Therefore, the commonly used rainfall-based SPI [52] was chosen as one of the meteorological components for this study. By normalizing the rainfall, the SPI was calculated on various time intervals, like 1, 2, and 3 months for the 2001 to 2018 study period. The SPI was computed using IMD-derived SPI ranges from higher negative values, representing extremely dry (−2.00 or less) conditions to the positive values of extremely wet (+2.00 or more) conditions. To maintain the seasonal relationship between rainfall and vegetation response, the SPI-3 (3-month SPI) was used for this study. The SPI values were calculated by using long-term (1981-2018) precipitation data. Climate Hazards Group Infrared Precipitation with Station data (CHIRPS) [53] generated by the U.S Geological Survey were incorporated to derive the time-series of SPI. CHIRPS data were used for this study, due to these main reasons: (1) data availability on a longer temporal scale (1981-present); (2) higher spatial resolution (0.05 • ) when compared to other available rainfall datasets; and (3) its significant correlation (r = 0.7) with the IMD-derived station-based data (as CHIRPS blends the station data, only independent stations were selected for correlation analysis).

LST (NOAH Data)
Land surface temperature is a key indicator of land surface processes, which reflects the energy balance at the surface of the Earth [54]. Being a fundamental linking parameter between SM and evapotranspiration, increasing LST precedes the onset of droughts. Hence, we have considered LST as an input for the CDI_M. The monthly LST data for the period 2001 to 2018 were retrieved from the Global Land Data Assimilation System (GLDAS), NOAH land surface model of NASA [55]. This LST data (spatial resolution-0.25 • ), is a part of 36 land surface fields of monthly data product version 2, which was forced using a combination of observations and modeled data [55]. NDVI measures the vigor and greenness of vegetation. NDVI deals with vegetation health, based on the differences between its reflection in the red and near-infrared spectrum. Unhealthy vegetation highly reflects in the red spectrum and fewer in near-infrared. Stressed and unhealthy vegetation results in low NDVI values, which is useful to identify and measure the agricultural drought intensities [2,56]. Therefore, to understand the agrarian stress and its relation with overall agricultural drought conditions, this study utilized NDVI data from 2001 to 2018. Monthly NDVI anomalies were collected from Moderate Resolution Imaging Spectroradiometer (MODIS) with a spatial resolution of 0.01 • (1 km). The benefits of the data obtained from MODIS are better radiometric calibration and its high temporal resolution, which provides continuous satellite data for drought monitoring.

SM (NOAH Data)
Soil moisture has the potential to indicate agricultural drought and level of water storage. The root zone soil moisture is one of the critical parameters in the overall vegetation growth. It is also responsible for water availability during the transpiration process, which directly indicates the agricultural drought condition [57]. Hence, we have used SM as a significant input parameter to identify agrarian droughts. Monthly SM data, for the period 2001 to 2018, were derived from the GLDAS, NOAH land surface model of NASA [54] for this study. Like LST, this SM data were also one of the parts of 36 land surface fields data product, version-2, and available at 0.25 • spatial resolution. Concerning land use type, root zone soil moisture at middle depth (10 to 40 cm) was used as an input for the SM. The primary objective of the current research was to monitor the spatio-temporal variations and changes during the drought conditions for agrarian activities. Hence, for the accurate assessment of the CDI_M results, the historical time-series of data for each input variable was masked to retain only pixel locations associated with agricultural activities within the study area. Yearly LULC data (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) were acquired from MODIS, with a spatial resolution of 1 km. These LULC datasets (MCD12Q1) were generated by a supervised classification method, which has divided Marathwada into 17 different classes (forest, water bodies, barren, urban and built-up land, etc.) [58]. Most of the study area is covered with croplands (95.93%), followed by urban and built-up lands (1.72%), water bodies (0.83%), etc. (Figure 3) ( Table 2). For masking and data extraction, only image pixels associated with the croplands (class 12) and croplands/vegetation (class 14) were retained for CDI_M development and assessment, while pixel locations associated with the other nonagricultural classes were eliminated. The LULC-masked input datasets were then used for the CDI_M estimation. agricultural drought condition [57]. Hence, we have used SM as a significant input parameter to identify agrarian droughts. Monthly SM data, for the period 2001 to 2018, were derived from the GLDAS, NOAH land surface model of NASA [54] for this study. Like LST, this SM data were also one of the parts of 36 land surface fields data product, version-2, and available at 0.25 o spatial resolution. Concerning land use type, root zone soil moisture at middle depth (10 to 40 cm) was used as an input for the SM.

LULC (MODIS Data)
The primary objective of the current research was to monitor the spatio-temporal variations and changes during the drought conditions for agrarian activities. Hence, for the accurate assessment of the CDI_M results, the historical time-series of data for each input variable was masked to retain only pixel locations associated with agricultural activities within the study area. Yearly LULC data (2001-2018) were acquired from MODIS, with a spatial resolution of 1km. These LULC datasets (MCD12Q1) were generated by a supervised classification method, which has divided Marathwada into 17 different classes (forest, water bodies, barren, urban and built-up land, etc.) [58]. Most of the study area is covered with croplands (95.93%), followed by urban and built-up lands (1.72%), water bodies (0.83%), etc. (Figure 3) ( Table 2). For masking and data extraction, only image pixels associated with the croplands (class 12) and croplands/vegetation (class 14) were retained for CDI_M development and assessment, while pixel locations associated with the other nonagricultural classes were eliminated. The LULC-masked input datasets were then used for the CDI_M estimation.    In this study, we have considered season-based, highly cultivated, and maximum yield producing crops from Marathwada that include bajra, jowar, maize, cotton, wheat, and rice over the year 2001 to 2018. These crops were grouped and analyzed as per their growing seasons, i.e., Rabi (November to April) and Kharif (June to October). The station-based crop area and crop production data were incorporated from the national level agricultural authority and open government data platform [59]. These datasets were further arranged as per seasonal yield production, using the following formula, Yield = Tonnes/Hectare (1)

Methodology
The input parameters used for this study were available at different spatial resolutions (0.25 • and 0.01 • ). Therefore, the first step of preprocessing was applied to spatially resample all the data sets to a 1-km resolution on a monthly temporal scale. This downscaling of the input data was performed by adopting the inverse distance weighted (IDW) technique [60]. This IDW method is easy to define, simple to understand, commonly used [61,62], and widely accepted in the scientific community and hence, implemented for this study. The general assumption behind IDW is points in close proximity are highly related than the distinct points. In the next preprocessing step, all the grid values were standardized by 'Z-score' statistics and further extracted using MODIS-LULC data sets.
Z-score values were computed using the following formula, where, X = particular value of a parameter from the group µ = long-term mean δ = long-term standard deviation This study deals with two different methods to combine the input parameters for the CDI_M analysis. Method I is a fixed percentage of weights were subjectively assigned to each input variable, whereas, method II determines the spatio-temporal weights of a parameter based on the PCA technique. Details of these two methods are described below.

Expert Judgment-Based Percentage Weights (Method-I)
The four input parameters used in this study have their own individual contributions to explain the collective drought conditions that may vary in importance. Therefore, in this first method of combining input variables for the CDI_M, weights were assigned using historical records of each parameter and its correlation with the actual in-field chronological drought conditions reported by the IMD. These weights and their contribution were also discussed and approved by the domain drought experts. In this analysis, the rainfall-based variable (SPI) was recognized as a key element to represent the ground reality of droughts; hence, it is assigned with the highest weight (40%). The remaining parameters of LST, NDVI, and SM showed relatively equal ability to represent in-field drought scenarios, and therefore were assigned equal weightings of 20% ( Figure 4). The monthly time series of CDI_M, based on this method was computed by the following formula: where Wt SPI , Wt LST , Wt NDVI , and Wt SM stand for the percentage weight values for individual parameters. z and i represent a particular year and the specific month (January-December).

PCA-Based Weights (Method-II)
In hydrologic and atmospheric studies, principle component analysis (PCA) is commonly used to detect a dominant pattern in the datasets [63,64]. Based on this capability, the PCA method was incorporated into this study as an alternative means to quantitatively determine input variable weightings. In the first step, based on all the input parameters and their temporal values, the correlation coefficient matrix (p*p, where p indicates the number of input variables) was constructed for each grid cell. In general, eigenvectors help to describe the relationship between data and the principal components. Therefore, as a second step, eigenvectors were computed using a matrix, developed in stage one and further employed for the orthogonal principal components (PCs). It is observed that the first PC (PC1) reflects the maximum variability of input data [65]; hence, we have considered PC1 for the development of the CDI_M.
In the next step, the percentage contribution (square of eigenvectors) of each variable was estimated using eigenvector values. These percentage contributions were generated and used to avoid the peaks and unwanted extremes in the CDI_M values.
Based on the above-mentioned steps, 48 spatially weighted maps for each month (12) and all the input variables (4) were created. These PCA-based percentage weights were further used to evaluate the CDI_M by using the following formula: where Wt SPI , Wt LST , Wt NDVI , and Wt SM stand for the PCA-based weight values for individual parameters, i.e., SPI, LST, NDVI, and SM, respectively. The z and i represent a particular year and the specific month (January-December) of the standardized input data.
After the computation of the CDI_M using both methods, all results were normalized by considering its standard deviation. The normalization was carried out using the formula mentioned below, where the stdCDI_M (z, i) is a combined drought anomaly for the z year and ith month. The δ i represents the standard deviation value for the month i, over all the years. After the evaluation of the CDI_M, all the drought events were classified based on IMD derived drought ranges and categories (Table 3).

Evaluation of the CDI_M Based on Crop Yield
Drought variability and its spatio-temporal distribution are highly responsible for the fluctuations in crop productivity and overall agricultural conditions. Hence, part of this study deals with the drought conditions and variations detected by the CDI_M and its relationship to the yields of principal crops cultivated in Marathwada. The primary step of preprocessing involves detrending the crop yield datasets to account for historical crop yield improvements due to improved plant genetics and farming technologies over time. Yearly crop production is not only a result of climatic variations but also depends on other nonclimatic factors that are accounted through data detrending. The detrending technique was used to remove the effect of variables other than agro-climatological parameters in the overall increasing trends in crop yield. This analysis is based on an annual yield data for the five significant crops and two major crop growing seasons (Kharif and Rabi) of Marathwada. The results of the correlation were produced based on the crop cycle (sowing and harvesting) and its spatio-temporal visualization through maps. The detrending technique was used to remove the effect of variables other than agro-climatological parameters in the overall increasing trends in crop yield. This analysis is based on an annual yield data for the five significant crops and two major crop growing seasons (Kharif and Rabi) of Marathwada. The results of the correlation were produced based on the crop cycle (sowing and harvesting) and its spatio-temporal visualization through maps.

Results and Discussion
Monthly maps of the CDI_M in Figure 5 show the spatio-temporal variations of drought frequencies over Marathwada. In the past 18 years (2001-2018), both of the CDI_M methods detected the years 2002, 2009, 2015, and 2016 as severe drought years, with conditions ranging from mild to extreme drought intensities across the cropping season. For the study area as a whole, the highest drought frequencies (number of extreme drought years) were observed during the monsoon (Jun, Jul, Aug, Sep) and postmonsoon (Oct and Nov) months. For the agriculture sector, in particular, these observed months are crucial for the overall cropping cycle of the Kharif season, whereas increasing drought frequencies over these months may also directly impact the following Rabi season's sowing activities in Marathwada. Between the years 2001-2018, the eastern parts of Jalna and central Parbhani experienced maximum drought frequencies (3/4 extreme drought years) in June and October, whereas the Osmanabad district had the highest number of drought frequencies (4 extreme drought years) in August ( Figure 5). For most of the remaining months ( Figure A1), only one to three severe to extreme drought events occurred over various parts of the Marathwada during the past 18 years. In general, Parbhani, Beed, and parts of Hingoli witnessed a higher number of drought frequencies during the monsoon and postmonsoon season, which can also be marked as continuous drought-affected expanses within Marathwada ( Figure 5). Over any region, three to four occurrences of severe to extreme drought events in 18 years are of great concern in the face of future uncertainties.

Results and Discussion
Monthly maps of the CDI_M in Figure 5 show the spatio-temporal variations of drought frequencies over Marathwada. In the past 18 years (2001-2018), both of the CDI_M methods detected the years 2002, 2009, 2015, and 2016 as severe drought years, with conditions ranging from mild to extreme drought intensities across the cropping season. For the study area as a whole, the highest drought frequencies (number of extreme drought years) were observed during the monsoon (Jun, Jul, Aug, Sep) and postmonsoon (Oct and Nov) months. For the agriculture sector, in particular, these observed months are crucial for the overall cropping cycle of the Kharif season, whereas increasing drought frequencies over these months may also directly impact the following Rabi season's sowing activities in Marathwada. Between the years 2001-2018, the eastern parts of Jalna and central Parbhani experienced maximum drought frequencies (3/4 extreme drought years) in June and October, whereas the Osmanabad district had the highest number of drought frequencies (4 extreme drought years) in August ( Figure 5). For most of the remaining months ( Figure A1), only one to three severe to extreme drought events occurred over various parts of the Marathwada during the past 18 years. In general, Parbhani, Beed, and parts of Hingoli witnessed a higher number of drought frequencies during the monsoon and postmonsoon season, which can also be marked as continuous drought-affected expanses within Marathwada ( Figure 5). Over any region, three to four occurrences of severe to extreme drought events in 18 years are of great concern in the face of future uncertainties. Therefore, a timely assessment of these dry spells is a basic need for the socio-economic and environmental well-being of Marathwada. Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 26  Figure 6 represents the spatial association between the CDI_M results from both methods. The greater correlation (0.71 < r < 1.00) in Method-I and Method-II was found in March (average = 0.92) and May (average = 0.88) over most of the study area. The lowest values (0.24 < r < 0.55) of relationship were found in the Nanded, Aurangabad, and Parbhani districts in January (average = 0.71), and April (average = 0.78). In this study, Method-I is based on fixed assigned weights for calculating the CDI_M, where the highest weight (40%) was given to the rainfall parameter. For most of the months ( Figure  A2), the correlation between both of the methods was low when Method-II (PCA) has assigned lower weights to rainfall in calculating the CDI_M (e.g., in January, Method-I allotted a 40% weightage to SPI, whereas the PCA method considered only a 9% contribution from SPI). Months with comparatively similar weight in both of the methods for all input parameters (or higher weight to rainfall in the PCA-based method) illustrates a greater correlation in CDI_M results between Method-I and Method-II ( Figure A2).   Figure 6 represents the spatial association between the CDI_M results from both methods. The greater correlation (0.71 < r < 1.00) in Method-I and Method-II was found in March (average = 0.92) and May (average = 0.88) over most of the study area. The lowest values (0.24 < r < 0.55) of relationship were found in the Nanded, Aurangabad, and Parbhani districts in January (average = 0.71), and April (average = 0.78). In this study, Method-I is based on fixed assigned weights for calculating the CDI_M, where the highest weight (40%) was given to the rainfall parameter. For most of the months (Figure A2), the correlation between both of the methods was low when Method-II (PCA) has assigned lower weights to rainfall in calculating the CDI_M (e.g., in January, Method-I allotted a 40% weightage to SPI, whereas the PCA method considered only a 9% contribution from SPI). Months with comparatively similar weight in both of the methods for all input parameters (or higher weight to rainfall in the PCA-based method) illustrates a greater correlation in CDI_M results between Method-I and Method-II ( Figure A2).

CDI_M-Based Spatio-temporal Drought Analysis-Method-I
Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 26 Therefore, a timely assessment of these dry spells is a basic need for the socio-economic and environmental well-being of Marathwada.  Figure 6 represents the spatial association between the CDI_M results from both methods. The greater correlation (0.71 < r < 1.00) in Method-I and Method-II was found in March (average = 0.92) and May (average = 0.88) over most of the study area. The lowest values (0.24 < r < 0.55) of relationship were found in the Nanded, Aurangabad, and Parbhani districts in January (average = 0.71), and April (average = 0.78). In this study, Method-I is based on fixed assigned weights for calculating the CDI_M, where the highest weight (40%) was given to the rainfall parameter. For most of the months ( Figure  A2), the correlation between both of the methods was low when Method-II (PCA) has assigned lower weights to rainfall in calculating the CDI_M (e.g., in January, Method-I allotted a 40% weightage to SPI, whereas the PCA method considered only a 9% contribution from SPI). Months with comparatively similar weight in both of the methods for all input parameters (or higher weight to rainfall in the PCA-based method) illustrates a greater correlation in CDI_M results between Method-I and Method-II ( Figure A2).

CDI_M-Based Spatio-Temporal Drought Analysis-Method-I
This method of computing the CDI_M is based on predefined weights for all of the input parameters, where the meteorological parameter (SPI) has the highest weight (40%), and equal weights (20% each) were assigned to all the other input variables (LST, SM, and NDVI). The time series of CDI_M was computed to evaluate the spatio-temporal variations in droughts over Marathwada. This method highlighted the years of 2002, 2009, and 2015-2016 as severe drought years, whereas other years were marked with wet or low-intensity dry spells. In the exceedingly affected dry periods of 2002 (53%), 2009 (42%), and 2015 (82%), nearly 40 to 60% of the study area was under severe-to-extreme stress of droughts ( Table 4). The central and southern districts of Marathwada (Beed, Latur, and Osmanabad) depicted a consistently high magnitude of dryness over the severe drought years. The influence of long-term and extreme droughts often extend over adjacent years. Therefore, several months during 2003, 2008, 2014, and 2016 also experienced mild to moderate dry stress as an extension of extreme drought years (2002, 2009, and 2015). This method observes that, with SPI, the remaining three input parameters are also collectively responsible for the enhancement of drought severities. Figure 7 shows the spatial coverage of significant drought years over Marathwada, derived from individual SPI method (Figure 7a) of drought analysis and based on CDI_M (Figure 7b). The individual SPI technique was able to identify the drought conditions at some extent, but CDI_M gives the more detailed and precise nature of expanses of drought ( Figure 7). This method of computing the CDI_M is based on predefined weights for all of the input parameters, where the meteorological parameter (SPI) has the highest weight (40%), and equal weights (20% each) were assigned to all the other input variables (LST, SM, and NDVI). The time series of CDI_M was computed to evaluate the spatio-temporal variations in droughts over Marathwada. This method highlighted the years of 2002, 2009, and 2015-2016 as severe drought years, whereas other years were marked with wet or low-intensity dry spells. In the exceedingly affected dry periods of 2002 (53%), 2009 (42%), and 2015 (82%), nearly 40 to 60% of the study area was under severe-to-extreme stress of droughts ( Table 4). The central and southern districts of Marathwada (Beed, Latur, and Osmanabad) depicted a consistently high magnitude of dryness over the severe drought years. The influence of long-term and extreme droughts often extend over adjacent years. Therefore, several months during 2003, 2008, 2014, and 2016 also experienced mild to moderate dry stress as an extension of extreme drought years (2002, 2009, and 2015). This method observes that, with SPI, the remaining three input parameters are also collectively responsible for the enhancement of drought severities. Figure 7 shows the spatial coverage of significant drought years over Marathwada, derived from individual SPI method (Figure 7a) of drought analysis and based on CDI_M (Figure 7b). The individual SPI technique was able to identify the drought conditions at some extent, but CDI_M gives the more detailed and precise nature of expanses of drought (Figure 7).  In the CDI_M computation, Method-I assigned a 40% weight to the SPI, and in Marathwada, 80% of the annual rainfall comes from monsoon rains. Hence, compared to other months, this drought In the CDI_M computation, Method-I assigned a 40% weight to the SPI, and in Marathwada, 80% of the annual rainfall comes from monsoon rains. Hence, compared to other months, this drought Remote Sens. 2020, 12, 2091 13 of 25 assessment method is highly accurate for the monsoon season. The trends in CDI_M values over the monsoon months in Marathwada are graphically represented in Figure 8. This figure exhibits the increasing trend in CDI_M values during July, August, and September. The CDI_M values range from −2 to +2. Lower values of CDI_M (negative) indicate the drier nature of a surface. Therefore, upward trends in the CDI_M inferred the increasing drought intensities from 2001 to 2018 in Marathwada. The enhanced drought stress is more pronounced in August and July, which are the primary sowing months for the Kharif crops. Other than monsoon months, an increasing trend in CDI_M values was also noticed during February, April, and October within the study area.
Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 26 assessment method is highly accurate for the monsoon season. The trends in CDI_M values over the monsoon months in Marathwada are graphically represented in Figure 8. This figure exhibits the increasing trend in CDI_M values during July, August, and September. The CDI_M values range from −2 to +2. Lower values of CDI_M (negative) indicate the drier nature of a surface. Therefore, upward trends in the CDI_M inferred the increasing drought intensities from 2001 to 2018 in Marathwada. The enhanced drought stress is more pronounced in August and July, which are the primary sowing months for the Kharif crops. Other than monsoon months, an increasing trend in CDI_M values was also noticed during February, April, and October within the study area.

CDI_M-Based Spatio-temporal Drought Analysis-Method-II
In this method, the PCA technique is used to integrate all the input parameters for the development of the CDI_M. The PCA-based method takes care of monthly variations in input variables and their changing percentage contribution with respect to space and time. Hence, the PCAbased CDI_M method seems to be more effective for drought monitoring over Marathwada. In this method, the first PC value is considered to assign the specific weight for each input parameter (SPI, SM, NDVI, and LST). Table 5 indicates parameter wise average weights for every month over Marathwada, whereas Figure A3 represents the sample of the spatial distribution of contribution of each input parameter in June. It is noted that on an annual scale, and in both the crop harvesting seasons (Rabi and Kharif), soil moisture is the highest contributing (35%) parameter in the CDI_M followed by LST (25%), SPI (21%), and NDVI (18%). The spatial distribution maps of weights for the Rabi season (harvesting months) show the highest weight values for SM and NDVI, and the lowest values are for the SPI and LST, at the central part of Aurangabad, Beed, and Jalna district. During the Kharif season (harvesting months), in the extreme southern portion of Osmanabad district, SM is given the maximum weighting followed by eastern and central Aurangabad for SPI and LST.

CDI_M-Based Spatio-Temporal Drought Analysis-Method-II
In this method, the PCA technique is used to integrate all the input parameters for the development of the CDI_M. The PCA-based method takes care of monthly variations in input variables and their changing percentage contribution with respect to space and time. Hence, the PCA-based CDI_M method seems to be more effective for drought monitoring over Marathwada. In this method, the first PC value is considered to assign the specific weight for each input parameter (SPI, SM, NDVI, and LST). Table 5 indicates parameter wise average weights for every month over Marathwada, whereas Figure A3 represents the sample of the spatial distribution of contribution of each input parameter in June. It is noted that on an annual scale, and in both the crop harvesting seasons (Rabi and Kharif), soil moisture is the highest contributing (35%) parameter in the CDI_M followed by LST (25%), SPI (21%), and NDVI (18%). The spatial distribution maps of weights for the Rabi season (harvesting months) show the highest weight values for SM and NDVI, and the lowest values are for the SPI and LST, at the central part of Aurangabad, Beed, and Jalna district. During the Kharif season (harvesting months), in the extreme southern portion of Osmanabad district, SM is given the maximum weighting followed by eastern and central Aurangabad for SPI and LST.
The time series of CDI_M maps from this method have highlighted the monsoon and postmonsoon months of 2015 and the premonsoon season of 2016 as the most severe drought period over Marathwada during our study's period of record. These two years and in particular the months of July, August, October, and April were noted as one of the most life-threatening droughts of the 21st century [43]. In July 2015 and October 2016, nearly 92% of the Marathwada region experienced severe-to-extreme drought conditions, and the average CDI_M values for most of the districts were less than −1.92. This severity of the drought continued for almost eleven months, starting from 2015 July to mid-June Remote Sens. 2020, 12, 2091 14 of 25 2016 (Figure 9). During this period, the average spatial drought spread has varied from 30% to 95% with respect to each month. In general, the south, southeastern, and central region of the Marathwada, which covers a major part of the Osmanabad, Beed, Latur, and Parbhani districts, was frequently under rigorous dry spells (Figure 9). Higher elevation and less availability of surface water resources are the common characteristics for almost all of the regions mentioned above. These characteristics might have acted as additional parameters for the enhancement of drought severities and the reduction in agricultural productivity over the 2015-16 period. The time series of CDI_M maps from this method have highlighted the monsoon and postmonsoon months of 2015 and the premonsoon season of 2016 as the most severe drought period over Marathwada during our study's period of record. These two years and in particular the months of July, August, October, and April were noted as one of the most life-threatening droughts of the 21st century [43]. In July 2015 and October 2016, nearly 92% of the Marathwada region experienced severe-to-extreme drought conditions, and the average CDI_M values for most of the districts were less than −1.92. This severity of the drought continued for almost eleven months, starting from 2015 July to mid-June 2016 ( Figure 9). During this period, the average spatial drought spread has varied from 30% to 95% with respect to each month. In general, the south, southeastern, and central region of the Marathwada, which covers a major part of the Osmanabad, Beed, Latur, and Parbhani districts, was frequently under rigorous dry spells (Figure 9). Higher elevation and less availability of surface water resources are the common characteristics for almost all of the regions mentioned above. These characteristics might have acted as additional parameters for the enhancement of drought severities and the reduction in agricultural productivity over the 2015-16 period.  Like Method-I, this method also depicted 2002 and 2009 as drought years, but compared to the first method, the PCA-based technique showed higher variations in month wise spatial distribution of droughts over Marathwada. Figure 10 represents the monthly spatial spread of statistical significance (P values) of the CDI_M values with respect to years. In August, 66.4% of Marathwada noticed a significantly increasing trend (p ≤ 0.05, hence, accepted the alternative hypothesis, i.e., drought intensities are increasing with respect to time) in CDI_M intensities, followed by October (26.7% area, where p ≤ 0.05), February (23.4% area, where p ≤ 0.05), and September (13.3% area, where p ≤ 0.05) ( Figure 10). CDI_M values range from −2 to +2; therefore, the increasing trend means CDI_M values are shifting towards more negative values, overtime for the period 2001 to 2018. In all of the above-identified months, parts of Hingoli, Parbhani, and Beed districts have experienced a significant upwelling trend in CDI_M values, whereas Aurangabad and Osmanabad districts noticed it in August and October, respectively ( Figure 10). The results of spatio-temporal trends in the CDI_M emphasizes the necessity of a PCA-based drought assessment method for the region and month specific drought management strategies.
Like Method-I, this method also depicted 2002 and 2009 as drought years, but compared to the first method, the PCA-based technique showed higher variations in month wise spatial distribution of droughts over Marathwada. Figure 10 represents the monthly spatial spread of statistical significance (P values) of the CDI_M values with respect to years. In August, 66.4% of Marathwada noticed a significantly increasing trend (p ≤ 0.05, hence, accepted the alternative hypothesis, i.e., drought intensities are increasing with respect to time) in CDI_M intensities, followed by October (26.7% area, where p ≤ 0.05), February (23.4% area, where p ≤ 0.05), and September (13.3% area, where p ≤ 0.05) ( Figure 10). CDI_M values range from −2 to +2; therefore, the increasing trend means CDI_M values are shifting towards more negative values, overtime for the period 2001 to 2018. In all of the above-identified months, parts of Hingoli, Parbhani, and Beed districts have experienced a significant upwelling trend in CDI_M values, whereas Aurangabad and Osmanabad districts noticed it in August and October, respectively ( Figure 10). The results of spatio-temporal trends in the CDI_M emphasizes the necessity of a PCA-based drought assessment method for the region and month specific drought management strategies.

CDI_M and Its Relation with Crop Yield
Numerous studies have observed the significant negative effect of severe to extreme drought events on crop yields over the USA [66], Canada [66], China [67], Korea [68], Bangladesh [69], and other parts on the world. India has particularly shown the highest risk of reduction in maize crops under the drought years [66]. Intense dry spells in India always have a dire influence on the community and environment. Marathwada agriculture sector, in particular, bears 80% of the direct impacts of severe drought events. Therefore, in the present research, we have tried to comprehend the spatio-temporal association between major crop yields and cropping season drought incidences.

CDI_M and Its Relation with Crop Yield
Numerous studies have observed the significant negative effect of severe to extreme drought events on crop yields over the USA [66], Canada [66], China [67], Korea [68], Bangladesh [69], and other parts on the world. India has particularly shown the highest risk of reduction in maize crops under the drought years [66]. Intense dry spells in India always have a dire influence on the community and environment. Marathwada agriculture sector, in particular, bears 80% of the direct impacts of severe drought events. Therefore, in the present research, we have tried to comprehend the spatio-temporal association between major crop yields and cropping season drought incidences. To reduce the influence of other non-agro-climatological parameters on crop production, all the crop yield datasets were detrended, as shown in Figure 11. Further, these detrended crop yield values were compared with the monthly CDI_M intensities using the statistical correlation method. CDI_M values from both the CDI_M computation methods and for all the wet and dry years were compared to the detrended anomalies of crop yields; however, in the drought years, the outcomes from the PCA-based CDI_M method showed a greater correlation (r > 6.5). Hence, the relationship between the CDI_M from Method-II (PCA-based) and crop yield anomalies is discussed below. To reduce the influence of other non-agro-climatological parameters on crop production, all the crop yield datasets were detrended, as shown in Figure 11. Further, these detrended crop yield values were compared with the monthly CDI_M intensities using the statistical correlation method. CDI_M values from both the CDI_M computation methods and for all the wet and dry years were compared to the detrended anomalies of crop yields; however, in the drought years, the outcomes from the PCA-based CDI_M method showed a greater correlation (r > 6.5). Hence, the relationship between the CDI_M from Method-II (PCA-based) and crop yield anomalies is discussed below.    Figure 13 demonstrates a spatial distribution of correlation coefficient values among major crop yields and CDI_M drought intensities over the harvesting of the Rabi season. This figure has identified the eastern and southeastern districts of Marathwada with higher correlation values (r > 0.5), whereas parts of Beed, Aurangabad, and Osmanabad in the west had negative correlations. In the case of jowar and wheat crop, a significant relationship was marked over large portions of Hingoli (r = 0.82 over 78% area), Parbhani (r = 0.76 over 58%), and Nanded (r > 0.6 over 32% area) districts. For a detailed understanding of the temporal association between crop yield and CDI_M values, linear regression graphs were plotted, and regression equations were generated over randomly selected sites.
In a further step, two potential assumptions (i.e., CDI_M values may increase or decrease by 0.5 or 1) were made for the future CDI_M, and based on it, the future changes in crop yield were predicted [70]. Regression equations given in Table 6 suggest that in the future, if the values of CDI_M increases or decreases by 0.5, then crop yields for Wheat, Jowar, and Maize will rise or decline by 16.7%, 9.2%, 9.1%, respectively. If the CDI_M varies by + or -1, then the yield will alter by 33% (wheat), 18.4% (Jowar), and 18.2% (Maize) Tonnes/ Hectare.

. Kharif Season
Monsoon rainfall is the primary resource for Kharif crops in Maharashtra. Hence, the fluctuations in the monsoon rains directly impact the dryness, soil moisture availability, and the overall crop productivity in the region. Figure 14 shows the spatial distribution of dry spells in the highest drought-affected years (2002,2009, over the sowing and harvesting months of Kharif. In 2002, eastern parts of Aurangabad and Nanded districts observed severe to extreme drought intensities during the starting of sowing and at the end of the harvesting season. 2015 was a crucial time over most of the cropping cycle, where more than 40% of the area of Marathwada continuously experienced extreme drought. In 2015, Kharif, Osmanabad, and Latur were the two most affected districts due to the life-threatening dry spells.

Kharif Season
Monsoon rainfall is the primary resource for Kharif crops in Maharashtra. Hence, the fluctuations in the monsoon rains directly impact the dryness, soil moisture availability, and the overall crop productivity in the region. Figure 14 shows the spatial distribution of dry spells in the highest drought-affected years (2002,2009, over the sowing and harvesting months of Kharif. In 2002, eastern parts of Aurangabad and Nanded districts observed severe to extreme drought intensities during the starting of sowing and at the end of the harvesting season. 2015 was a crucial time over most of the cropping cycle, where more than 40% of the area of Marathwada continuously experienced extreme drought. In 2015, Kharif, Osmanabad, and Latur were the two most affected districts due to the life-threatening dry spells. In the Rabi season, jowar and cotton crops had a stronger relationship (r > 0.5) with the CDI_M based drought intensities. At the district level, in particular, almost the entire area of Beed and some pockets of the Nanded showed a significant association (r > 0.7) in CDI_M and yield of jowar ( Figure  15). In the case of cotton, the entire Hingoli and northern horn of Osmanabad districts observed a significant correlation (r > 0.7) between crop yield and drought intensities ( Figure 15).
To get the general idea of future Kharif crop yield anomalies, based on assumed variations in the CDI_M values (increase or decrease by 0.5 or 1), the same regression equation technique was implemented as described for the Rabi crop. Regression equations given in Table 7 suggest that in the future, if the values of CDI_M increase or decrease by 0.5, then crop yields for Jowar and Cotton will rise or decline by 19.8%, and 11.3%, respectively. If the CDI_M varies by + or -1, then the yield will alter by 39.5 % (Jowar) and 22.6% (Cotton) Tonnes/ Hectare. In the Rabi season, jowar and cotton crops had a stronger relationship (r > 0.5) with the CDI_M based drought intensities. At the district level, in particular, almost the entire area of Beed and some pockets of the Nanded showed a significant association (r > 0.7) in CDI_M and yield of jowar ( Figure 15). In the case of cotton, the entire Hingoli and northern horn of Osmanabad districts observed a significant correlation (r > 0.7) between crop yield and drought intensities ( Figure 15).
To get the general idea of future Kharif crop yield anomalies, based on assumed variations in the CDI_M values (increase or decrease by 0.5 or 1), the same regression equation technique was implemented as described for the Rabi crop. Regression equations given in Table 7 suggest that in the future, if the values of CDI_M increase or decrease by 0.5, then crop yields for Jowar and Cotton will rise or decline by 19.8%, and 11.3%, respectively. If the CDI_M varies by + or −1, then the yield will alter by 39.5 % (Jowar) and 22.6% (Cotton) Tonnes/Hectare.

Limitations and Future Recommendations
In this study, a weighted combination of input parameters was used to compute the CDI_M. The weights for CDI_M were based on expert judgment (Method-I) and the principle component analysis (PCA) weighting approach (Method-II). However, several other methods like entropy [71,72], Copula [73] can also be implemented and compared in the development of combined drought indicators.
All the input parameters used for this study were initially available in different spatial resolutions. Some of the datasets were best available only with the coarse resolution (SM and LST-0.25 • ). Though these coarse resolution data still represented the spatial variability, this study recommends the use of high-resolution input data for future research.
The correlation coefficient maps of the crop yield with the CDI_M showed a weak correlation (r < 5.0) over some parts of the study area. Unavailability of grid-based or high-resolution crop yield data and shorter data length (18 years) can be responsible for the above-mentioned poor association between crop yield and CDI_M.
Current operational drought monitoring methods in the Marathwada, India, are based on a single input parameter, i.e., rainfall. Hence, the combined drought indicator developed in this study constitutes a new technique for agricultural drought monitoring over the study area. The drought assessment methods used in the framework of this study recommend its implementation over the various parts of the world by adopting the region-specific input variables and their relative weights.

Conclusions
During the past few decades, the global environment has been brutally affected by droughts [74]. In recent years, severe drought events and water crises have increased widely over the various parts of India, especially Marathwada. Hence, for the accurate assessment of droughts, this study has developed the combined drought index for Marathwada. The present study has used fixed weighted-based and PCA-dependent methods to develop a combined drought index for Marathwada. Both methods were capable of identifying the spatial extent of historical and current drought events over the study area. Two agricultural (NDVI and SM) and two meteorological (SPI and LST) variables were utilized in developing the CDI_M as an attempt to understand the agro-climatological drought severities over Marathwada in a better way. In the later stages of this study, both the CDI_M results were compared with the temporal trends in major crop yields in Marathwada. Over both crop growing seasons, the PCA-based CDI_M showed a higher correlation with most of the crop yields compared to the fixed weight-based CDI_M. The outcomes from this study conclude that, Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 26