Utilizing TVDI and NDWI to Classify Severity of Agricultural Drought in Chuping, Malaysia

Abstract

Agricultural drought is crucial in understanding the relationship to crop production functions which can be monitored using satellite remote sensors. The aim of this research is to combine temperature vegetation dryness index (TVDI) and normalized difference water index (NDWI) classifications for identifying drought areas in Chuping, Malaysia which has regularly recorded high temperatures. TVDI and NDWI are assessed using three images of the dry spell period in March for the years 2015, 2016 and 2017. NDWI value representing water content in vegetation decreases numerically to −0.39, −0.37 and −0.36 for the year 2015, 2016 and 2017. Normalized difference vegetation indices (NDVI) values representing vegetation health status in the given area for images of years 2015 to 2017 decreases significantly (p ≤ 0.05) from 0.50 to 0.35 respectively. Overall, TVDI in the Chuping area showed agricultural drought with an average value of 0.46. However, Kilang Gula Chuping area in Chuping showed a significant increase in dryness for all of the three years assessed with an average value of 0.70. When both TVDI and NDWI were assessed, significant clustering of spots in Chuping, Perlis for all the 3 years was identified where geographical local regressions of 0.84, 0.70 and 0.70 for the years 2015, 2016 and 2017 was determined. Furthermore, Moran’s I values revealed that the research area had a high I value of 0.63, 0.30 and 0.23 with respective Z scores of 17.80, 8.63 and 6.77 for the years 2015, 2016 and 2017, indicating that the cluster relationship is significant in the 95–99 percent confidence interval. Using both indices alone was sufficient to understand the drier spots of Chuping over 3 years. The findings of this research will be of interest to local agriculture authorities, like plantation and meteorology departments to understand drier areas in the state to evaluate water deficits severity and cloud seeding points during drought.

1. Introduction

Drought in agriculture occurs when there is insufficient water available for a specific crop to grow at a specific time. Agricultural drought typically happens after meteorological drought, an indication of dropped rainfall, but before a hydrological drought, in which river, lake and reservoir water level decline. Drought in agriculture is dependent not only on the amount of rainfall, but also on how efficiently the available water was utilized. Droughts have different effects on irrigated and non-irrigated agriculture land, which should be acknowledged [1]. The impacts of drought in irrigated areas are normally less severe since water supplies in reservoirs are restricted. As a result, even if there is no rain, these crops will receive the water needed. In non-irrigated areas, crops depend heavily on precipitation. If the precipitation rates are reduced, crops will suffer from water scarcity. In such conditions and circumstances, drought stress is the most common environmental factor limiting crop productivity. It is reported that the frequency of severe drought conditions is increasing in accordance with global climate change [2].

A good description of agricultural drought takes into account crop vulnerability at various stages of growth, from emergence to maturity [3]. There are mechanisms that plants experience during agriculture drought that reduce crop yield due to soil water deficit [4]. It is suggested that in nature, plants have evolved to endure drought stress with a range of morphological, physiological and biochemical adaptations [4]. Drought escape, drought avoidance, drought tolerance and drought resistance are the four mechanisms of drought resistance that help plants thrive in moisture-deficient conditions [5]. The reproducibility of drought stress treatments is large and difficult to carry out, significantly obstructs research on plant drought tolerance [6]. Traditional breeding efforts and the use of modern genetic methods to increase the drought tolerance of crop plants have been hindered by a slow speed in revealing drought tolerance mechanisms [7]. Oil palm crop is the dominant crop cultivated in Malaysia. The decline of oil palm yield is a direct effect of inhibition of the photosynthetic rate of oil palm. Agricultural drought should be identified early as to avoid yield loss and to reduce the impact of global warming as water deficit will become more widespread in the coming years, including in the field of oil palm plantations. Drought affected areas have climate irregularities such as ‘El Nino’ causing drought in southeast Asia [6]. Drought causes the stomatal conductance of oil palm leaves to decline rapidly because the stomata tend to close [8]. Oil palm yield is highly dependent on water availability during the differentiation sex of its inflorescences, which occurs approximately 28 months prior to bunch harvest [9]. Drought affects oil palm crop yield by reducing fruit bunch number through changes in the ratio of female to male inflorescences and inflorescence abortion rate [6]. These changes occur up to two years before fruit harvest. Thus, there may be differences between existing photosynthetic production and demand from developing bunches [10]. The current bunch numbers are highly dependent on drought effects in previous years [11].

On the other hand, for the paddy crop, water scarcity hugely impacts the harvest output and leads to lower yield and earnings. Agriculture tends to be highly vulnerable to short-term weather changes, resulting in fluctuating crop yields. Rice yield decreased due to the dryness of a region due to inadequate rainfall and gradual rise in temperature [3]. Drought stress causes early senescence, which causes a variety of changes in rice traits, including tillering, leaf expansion and midday photosynthesis [12]. Drought has a greater negative impact on rice during the reproductive process, such as the blooming period, filling stage and maturity, as plant growth increased [13].

Sugarcane is an exceptional crop with the ability to accumulate large amounts of sugar. In Malaysia, areas such as in Chuping, Perlis is a huge contributor to the sugarcane production for the country. Sugarcane is a high-water-demanding crop and water scarcity, which is a major abiotic stress affecting sugarcane productivity, affects its development [14]. Sugarcane plant morphological and physiological responses differ depending on genotype, stress period (rapid or gradual), stress intensity (severe or mild) and type of affected tissue when experiencing water stress [15]. A significant reduction in the growth and net assimilation rate was observed in sugarcane under heat stress [16]. Under heat stress, Ebrahim et al., 1998, identified a reduction in inter-nodal duration and biomass accumulation, as well as early leaf senescence in sugarcane [17].

In relation to that, the latest trend of satellite technology sensors can refine many agricultural and climate change variations through spectral radiance measures. Some of these manipulation steps take into account vegetation indices, which are responsive to both the rate and the amount of plant growth. Moisture stress causes changes in vegetation and certain indices are responsive to those changes as seen in Table 1. One of the most used indices is normalized difference vegetation indices (NDVI) [18]. Researchers have reported the use of NDVI for vegetation monitoring, assessing the crop cover, drought monitoring and agricultural drought assessment [19]. On the other hand, water stress can be detected using remote sensing sensors in visible, near-infrared (NIR), shortwave infrared (SWIR) [20] and thermal infrared (TIR) bands [21]. Table 1 shows vegetation condition index (VCI), temperature condition index (TCI), vegetation health index (VHI) and standardized vegetation index (SVI) used to understand the vegetation and water deficits that relate to each other to characterize the vegetation health. Using the mentioned region of the spectrum, normalized difference water indices (NDWI) are sensitive to the water content of the vegetation canopy [22]. On the other hand, the NDWI is often used to view drought assessments with a focus on vegetation moisture content [23,24]. Chlorophyll and water content of plants are regularly used as significant indicators of plant stress in this context. Chlorophyll content decreases in stressed vegetation, resulting in an overall difference in light absorption by leaf pigments [25]. Drought research has taken a detour since the implementation of NDWI in terms of approaches to drought detection and control. In the relevant context, the temperature vegetation dryness index (TVDI) is an alternative drought index that uses parameters dependent on the sensitivity of the near-infrared and visible spectrum to identify vegetation and stress conditions associated with water shortages. Combining both surface temperature and normalized difference vegetation indices (NDVI) forms a virtual triangle where TVDI was successfully used as an indicator of drought change application [26]. The TVDI considers the relationship and changes between the NDVI and Land Surface Temperature (LST) in depth as seen in Figure 1. Moreover, when using the added temperature information to monitor soil moisture, temperature is a time-sensitive parameter that serves as a predictor of water stress and can be influenced by vegetation coverage using NDVI [27]. With the ability to composite information from visible, near-infrared and thermal infrared bands of light, TVDI evaluates soil moisture by combining vegetation indexes and surface temperature.

agronomy-11-01243-t001_Table 1 Table 1

Influential satellite-based indices useful in agricultural drought assessment.

Agricultural DroughtIndices	Formula	Advantages	Limitation
NDVI	$NDVI=\frac{{\mathsf{\rho }}_{nir}-{\mathsf{\rho }}_{red}}{{\mathsf{\rho }}_{nir}+{\mathsf{\rho }}_{red}}$[28]	To interpret vegetation details, a widely used index is used to understand variations and changes in green leaves from plants, as well as canopy spectral characteristics [18].	In complicated stand structures and closed canopy forests, data saturation related NDVI is reported as a major constraint or cause of uncertainty in biomass assessment [28,29].
VCI	$VC{I}_i=100\times \frac{NDV{I}_i-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}$[30]noting that NDVIi represents the continuous mean NDVI of a certain period i e.g., a month, NDVImax represents the maximum NDVI and NDVImin represents the minimum NDVI.	Used for assessment of vegetation in drought situations affecting agriculture [31]. Mainly used when covering a large spatial and continuous data is provided [32].	Using sole VCI was found to have disagreement and showed to be insufficient for drought analysis accurately [33,34].
TCI	$TCI=100\times \frac{B{T}_{max}-B{T}_j}{B{T}_{max}-B{T}_{min}}$[30]where the smoothed weekly brightness temperature, multi-year maximum and multi-year minimum, respectively, at each grid cell are BTj, BTmax and BTmin.	Developed to assist the VCI in assessing vegetation stress relative to temperature and evaluating stress caused by extreme wetness and extracted from the thermal band that is translated to brightness temperature [35].	TCI to be coupled with VCI as a tool to monitor both drought and excessive wetness for better evaluation, however, both indices have no direct relationship with precipitation [36].
VHI	$VHI=\propto \times VCI+\left(1-\propto \right)\times TCI$[37] $\mathrm{in} \mathrm{reference} \mathrm{to} \propto$ is the coefficient defining the contribution of the two indices at 0.5. A VHI value of less than 40 indicates the presence of vegetation stress, while a value of greater than 60 indicates that the vegetation is in good condition.	Combination of VCI and TCI data representing the relationship between temperature and vegetation growth as a proxy index characterizing vegetation health [38].	The optimum weights usually used 0.5, not known typically estimated by assuming equal weights of 0.5 to VCI and TCI [39]. Assuming as 0.5 because of the unknown conditions of soil moisture, temperature and precipitation [40].
SVI	$Z_{ijk}=\frac{V{I}_{ijk}-{\mu }_{ij}}{{\sigma }_{ij}}$VI whereby both either NDVI or EVI can be utilized [41]. Adding on $Z_{ijk}$ is the z-value for pixel i during week j for year k.	It is based on the NDVI and improves the start, number, severity and duration of vegetation stress. The SVI has the capacity to work in tandem with standard drought indices, as well as other weather and supplementary data, to help drought-response decisions [42].	It is vital to note that during SVI data analysis, each pixel is only compared to the data contained in that single pixel, which is the mean over time [43].
NDWI	$NDWI=\frac{{\mathsf{\rho }}_{nir}-{\mathsf{\rho }}_{swir}}{{\mathsf{\rho }}_{nir}+{\mathsf{\rho }}_{swir}}$[22]	The amount of water in the vegetation canopy has an impact on the NDWI, where satisfactory performance to estimate canopy water content with a lower error despite saturation at high level [23]. It is often used for drought assessment interpretation with an emphasis on moisture content in vegetation [24].	Poor performance due to spectral similarities between built-up regions where specific types of land cover the region, such as sand and exposed bedrock [44,45].
TVDI	$TVDI=\frac{T_s-T_{S_{min}}}{a+bNDVI-T_{S_{min}}}$ [46]	It is also important to underline that estimation of soil moisture using the TVDI gives root-zone moisture which is available to the plants [47]. The use of TVDI is sufficient for deployment over vast areas [48].	Scale or resolution of data is reported to be one aspect that can enhance uncertainty in soil moisture extraction [49,50], additionally TVDI depends on time scale and location [51].

The goal of this research is to use NDWI and TVDI to determine the agricultural drought. These indicators will subsequently be used to categorize locations within the research region that have been severely impacted by the drought. Landsat 8 imagery were chosen to support the TVDI index as they combined optical-thermal infrared remote sensing data given with cloud-free satellite imagery of the study site. This study used Malaysian Meteorological Data and Global Agriculture Monitoring—JASMIN, GCOM to choose three different Landsat 8 satellite imagery dated 18 March 2015, 20 March 2016 and 23 March 2017 based on the drier seasons. This research is limited to the mentioned data and the availability of resources as a case study.

2. Materials and Methods

Chuping, Perlis is in the northwestern part of peninsular Malaysia, and it is mostly covered in oil palm, paddy, rubber trees and sugarcane. As shown in Figure 2, the region that borders southern Thailand is estimated to be about 102 km². The mentioned area is between 6°27′30″ N and 6°35′00″ N latitude and 100°15′00″ E and 100°22′30″ E longitude.

This research paper used meteorology data from the Malaysian Meteorology Department and Global Agriculture Monitoring satellite to recognize the dry season occurrence and capture the rainfall at ground level. In relation to that, remote sensing data available from National Aeronautics and Space Administration (NASA) where three images of Landsat 8 Operational Land Imager (OLI) data, dated 18 March 2015, 20 March 2016 and 23 March 2017 were collected to derive NDWI, NDVI, LST and TVDI. Statistical test was evaluated before the spatial correlation was explored for all the indices derived. The spatial correlation using Moran’s I was used to further evaluate the mentioned indices to classify significant spots of the study area in relevance to drought.

2.1. Meteorology Data

In the study area, the early months of the years 2015, 2016 and 2017 showed very low precipitation rates. This is especially in the months of February and March where it is considered to be the driest period of the year based on the meteorological data of satellite Global Agriculture Monitoring, GCOM-JASMIN Sensor from Japanese Aerospace Agency, Japan as seen in Figure 3. Data from Malaysian Meteorological Department (MMD) was acquired for the mentioned time showed that the research site was under a dry spell as seen in Table 2.

2.2. Landsat Image Processing

NDVI is preferred for global vegetation monitoring by using the previously described satellite data because it compensates for changes in lighting conditions, surface slope, exposure and other external factors [19,40]. This index is computed using the near-infrared (NIR) and red channels which are obtained at 30 m resolution with a wavelength of 0.85–0.86 µm and 0.64–0.67 µm respectively [28]. The near-infrared band reflectance,${\mathsf{\rho }}_{nir}$ and the red band reflectance, ${\mathsf{\rho }}_{red}$ are in Equation (1). In this study, NDVI was measured using Landsat 8 Band 4, the red band reflectance, and Landsat 8 Band 5, the near-infrared band reflectance.

(1)$NDVI=\frac{{\mathsf{\rho }}_{nir}-{\mathsf{\rho }}_{red}}{{\mathsf{\rho }}_{nir}+{\mathsf{\rho }}_{red}}$

The NDWI is a measure of liquid water molecules in plant canopies that interact with incoming solar radiation and it was designed to estimate soil moisture and canopy water content [52]. As it includes a short-wave infrared band, the NDWI is frequently a function of local climate and soil qualities affecting water availability and is sensitive to changes in liquid water [53,54]. NDWI derived from Landsat uses two near-IR channels with 30 m spatial resolution, where NIR wavelength 0.85–0.86 µm and SWIR-1 wavelength is 1.57–1.65 µm as seen in Equation (2). The simplicity of NDWI [22] used the Landsat 8 Band 5 as near-infrared band reflectance and Landsat 8 Band 6 as the short wave near-infrared band reflectance.

(2)$NDWI=\frac{{\mathsf{\rho }}_{nir}-{\mathsf{\rho }}_{swir}}{{\mathsf{\rho }}_{nir}+{\mathsf{\rho }}_{swir}}$

In this study, TVDI was derived from NDVI and LST using Landsat 8 OLI images. TVDI was calculated using Equation (3), [46]

(3)$TVDI=\frac{T_s-T_{S_{min}}}{a+bNDVI-T_{S_{min}}}$

in which parameter a is estimates of the wet range while b is the dry range of surface moisture content from bare soil to fully vegetated surfaces, respectively; T_s is the radiative surface temperature observed in a given pixel in Kelvin unit; T_smin is the minimum surface temperature in Kelvin unit. The steps to derive the TVDI using Landsat 8 OLI imagery follows:

2.3. Steps in Derivation of LST

Step 1: The Landsat 8 OLI data was calibrated in the 16-bit unsigned integer format for the images used in this report. As a result, the Landsat 8 Band 10 of wavelength from 10.69–11.19 µm with the spatial resolution of 100 m was converted to TOA Spectral Radiance (L) using the radiometric rescaling coefficients given in the metadata register, such as band-specific multiplicative rescaling factor (ML) and band-specific additive rescaling factor (AL) using the Equation (4) referred from Landsat 8 Data Users Handbook, United States Geological Survey (USGS) [55].

(4)$L_{\mathsf{\lambda }}=M_L\times {\mathrm{Q}}_{cal}+A_L$

Step 2: As seen in Figure 4, using the spectral radiance, the thermal band 10 was converted to TOA brightness temperature (BT) using Equation (5). In the mentioned equation, K1 and K2 were taken from the metadata file of the respective images which are the thermal constants of Band 10. The equation will readily minus the 273.15 to convert from Kelvin to Celcius [55].

(5)$BT=\left(\frac{K2}{\ln \frac{K1 }{L}+1}\right)-273.15$

Step 3: Using Equation (2), NDVI is calculated for each respective image where the proportional vegetation (P_v) can be calculated using Equation (6). Hence, the area under the vegetation and bare soil is proportional to pure NDVI pixels [56].

(6)$P_v=\left(\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}\right){ }^2$

Step 4: Land surface emissivity (ε) is calculated by using the average emissivity of an element surface of the earth estimated from proportional NDVI values using Equation (7) [57].

(7)$\epsilon =0.004\times P_v+0.986$

Step 5: Lastly, to calculate the land surface temperature (LST) in Celsius using the TOA brightness temperature (BT), wavelength of emitted radiance (λ) and land surface emissivity (ε) in Equation (8) were used where is ρ is 1.438 × 10⁻² mK [58].

(8)$LST=\frac{BT}{1+\left( \frac{\mathsf{\lambda }\times BT}{\mathsf{\rho }}\right)} x Ln \left(\epsilon \right)$

2.4. Statistical Analysis and Derivation of Global Moran’s Index

In analyzing the derived indices using Landsat 8 mentioned earlier, NDVI, NDWI, TVDI and LST images were evaluated for basic raster statistics. Raster statistic was calculated for mean values. The mean values of each sub-area based on Figure 2 was calculated in ArcMap software using zonal statistic to determine the values in the respective years of 2015, 2016 and 2017. Using the calculated mean values, ANOVA statistical testing was carried out to understand the significant values from each of the indices derived to emphasize the difference in the sub-areas in the respective years.

In order to further evaluate the spatial difference in the study area, Moran’s I Index was quantified. Moran’s I approach is one of the earliest methods for assessing spatial autocorrelation [59] and most commonly used approach [60]. Moran’s I is defined as a measure of the correlation among neighboring observations in a pattern [61]. Moran’s I approach uses correlation that can determine the spatial correlation between feature locations and values [62]. With the given Equation (9), the approach will evaluate whether features and associated variables had a clustered, dispersed or random distribution pattern [63]. Moran’s I procedure is statistically represented as follows [59]:

(9)$I=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{W}_{ij} \left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{s^2\times {\sum }_{i=1}^n{\sum }_{j=1}^n{W}_{ij} }$

where n is the number of observations,$x_i$ is the value at a point i, $x_j$ is the value at point i’s neighbor j, $\overline{x}$ is the observed mean of all observation,$s^2$ is the variance of x value and $W_{ij}$ is the spatial weight of location i and j. For the measurement of the spatial autocorrelation a binary weight matrix,$W_{ij}$ is vital to show the connectivity among the observed value. Adding on,$z$ or Z score is calculated in Equation (10) where the E (I) is the expected value shown in Equation (11) and the $\sqrt{Var \left(I\right)}$ represent the variance [64]. In z-distribution with mean zero and variance one, p-values were obtained by comparison with the standard score or Z score. Z score is an excellent technique to determine where an observation falls into the distribution as a whole.

(10)$z=\frac{I-E\left(I\right)}{\sqrt{Var \left(I\right)}}$

(11)$E\left(I\right)=\frac{-1}{n-1}$

In general, a positive value for I indicates that a feature is surrounded by other features that have similar high or low attribute values; this feature is regarded to be a member of a cluster. A negative value for I indicates that an attribute contains adjacent features with separate values; this attribute is an outlier [65]. In either situation, the p-value for the feature must be small enough for the cluster or outlier to be considered statistically significant. The z-scores and p values indicate the statistical significance of the index values obtained. In our study area, the spatial autocorrelation was explored to understand the correlation between NDWI and TVDI where spatial significant clustering spots were generated using Equations (9) and (10) with input of 300 points of localized regression, where the stronger relationship is represented by p values and Z scores. Moran’s I equation was evaluated to produce Z score and p values embedded in the Figure 2 boundaries in ArcMap software environment for 2015, 2016 and 2017 years to identify classifying clusters in relation to agriculture drought.

3. Results and Discussions

3.1. NDWI Derivation

NDWI was used in this analysis to assess the reflection of water content in the soil and on plant surfaces. Rather than focusing exclusively on a spectral band whose reflection strength is largely determined by the phase of chlorophyll in leaves, adding short wave near-infrared (SWIR) emphasized light absorption by water [66,67]. The NDWI unit has no dimensions and ranges between −1 and +1 depending on the amount of surface water present. A higher NDWI value indicates a high-water content in the plant. Low NDWI values indicate low vegetation water content and thus the NDWI rate decreases during periods of water stress [68]. At this point, it is frequently stated that the NDWI can be used to determine the degree of wetness or dryness. As precipitation increases, the NDWI value increases [69]. In general, the NDWI for the Chuping district did not indicate substantial improvements over the reference years as shown in Figure 5. The average NDWI value representing the water content of vegetation decreases numerically, −0.39, −0.37 and −0.36 for the given Chuping area images from 2015, 2016 and 2017. However, when compared to other areas in the respective years, the values for the Kilang Gula Chuping region were significantly different. NDWI values of −0.46, −0.44 and −0.45 were observed in the respective years, as shown in Figure 6a. On the other hand, the Kubang Perun registers the lowest value in 2017 at −0.28. NDWI alone indicated a drier area in Kilang Gula Chuping.

3.2. NDVI Derivation

NDVI indices in the remote sensing field are used to show the relation between spectral variability and changes in the vegetation at growth rate. As the index is determined using a normalization process, NDVI values fall between 0 and 1, indicating a sensitive reactivity to green vegetation even in places with little vegetation cover. This index commonly applies to studies on regional and global vegetation evaluations and was not only related to the structure of the canopy but to the photosynthesis of the canopy as well [18]. In this study, the Chuping area has shown mean values of 0.50, 0.30 and 0.35 in the years 2015, 2016 and 2017 respectively in Figure 6b. For the year 2015, the NDVI values from the year are significantly higher (p ≤ 0.05) than the years 2016 and 2017. Generally, a declining trend of NDVI was observed in the area. The NDVI mentioned values were calculated using Landsat remote sensing images with a resolution of 30 m by 30 m which is moderate-resolution data that can represent the properties of vegetation coverage in a specific area. The choice of fine resolution will aid in not only increasing the accuracy of field features but also in improving the efficiency of crop and hydrological models [70], resulting in more conclusive outputs and planners in pertaining vegetation cover in this study.

In recent literature, there are discussions that NDVI values derived at non-irrigation and irrigation areas are different. Recently, NDVI derived from optical remote sensing data enabled the separation of irrigated and non-irrigated land usage in India and Turkey’s semi-arid and desert regions [71,72]. In relation to that, there is a study in Thailand in the tropical areas in which NDVI have reflected the actual conditions of crop growth and irrigation in the area [73]. It is recorded that NDVI provided accurate mapping of the irrigated lands at a regional scale to understand when there is a decrease of rainfall observed [74]. The changes in rainfall patterns in a warmer climate will have a strong influence on the hydrologic cycle of all agricultural lands. To estimate the likely effect of these NDVI values on these areas, accurate knowledge of the geographic distribution of irrigated areas at the national scale is required [75]. In order to evaluate this effect, it is suggested that the area needs to add irrigated agriculture as a separate class [75] with ground data to comprehend the effects of non-irrigated and irrigated crops. In this study, the values of NDVI taken as a single class were then limited to represent the overall vegetation covered in this research area.

3.3. TVDI Derivation

When using the TVDI approach, LST and NDVI coupled will provide information on vegetation and surface temperature. The wavelengths of the thermal and visible spectrum have been successful in monitoring vegetation conditions when water stress occurs [76]. Interpretations LST-NDVI relationship is based on vegetation crops, climate, at various scales [77,78,79]. Since the vegetation stays green after the onset of the stress, the NDVI alone as a water stress indicator is less reliable [80]. The LST, on the other hand, rises rapidly in response to water stress. [46]. The association between LST and NDVI is categorized by a cloud of dispersion space for an assumed dry area, with the highest values of LST corresponding to the lowest values of NDVI [81]. LST in the Chuping area generally shows a significant difference (p ≤ 0.05) of a declining trend from year image 2015 to year image 2017 with an average temperature of 32.70 °C, 28.50 °C and 26.99 °C as seen in Figure 7c. It was found out that the land use type has a close impression on the LST [82], e.g., oil palm trees have a mean of 21.65 °C when the rubber crop had a higher LST of 21.62 °C [83]. In contrast, as urban areas showed higher temperature records in peninsular Malaysia where an average of 23.17 °C to 26.45 °C was observed [84]. Adding on, Ulu Klang which is on the fringes of the city of Kuala Lumpur, recorded a higher LST of 26–30 °C from 1994 to 2019 [69]. However, in Chuping the urban area is limited and less developed when compared to other mentioned studies [48,55,56,57]. In this research, the entire Chuping recorded an average high LST of 29.39 °C from 2015–2017 in the driest month where more agriculture cultivation is present. Temperature vegetation dryness index indicates that all the 3 years of images of the Kilang Gula Chuping area are significantly different (p ≤ 0.05) as seen in Figure 7. As observed in Figure 7, it is evidently noted that Kilang Gula Chuping has higher values of TVDI indication the dryness level. It is noted that 0.0–0.2 shows most wet value, 0.2–0.4 representing generally wet, 0.4–0.6 showing normal or balance value of wet to dry, 0.6–0.8 as more dry where 0.8–1.0 shows extremely dry; based on the NDVI and LST values seen in triangular relationship proposed [46]. The index also shows that the area faces agricultural drought or dryness with values of 0.73, 0.66 and 0.71 as seen in Figure 6d, where else the entire Chuping is averaged at 0.48, 0.45 and 0.47 in the respective years mentioned.

3.4. Correlation of TVDI and NDWI

In order to relate TVDI and NDWI and also to identify the stronger correlated regions in the mentioned region, local regression and spatial autocorrelation using Moran’s I was performed. Using the mentioned images, the relationship of the local regression of TVDI and NDWI was observed at R-squared at 0.8399, 0.6969 and 0.6974 for the respective years of 2015, 2016 and 2017. The correlation of TVDI and NDWI was observed to be very strong at the dryer regions in Chuping as seen in Figure 7., especially at the Kilang Gula Chuping for all the mentioned years. Similarly, NDWI showed dryer values in the respective region as seen in Figure 5. In learning the correlation to both the NDWI and TVDI, as seen in Table 3, the values of local regressions shown significant relation where p ≤ 0.05 and the clustered relationship was identified based on the Z score. Table 3 shows that the Z score is highest in the year 2015 with a score of 17.806, and lowest in 2017 with a score of 6.769, while in 2016, it had a score of 8.625 when the p value was 0.0001 for all the years. Spatial clustering can be determined when the p value is minimal and the absolute value of the Z score is large enough or higher than 2.58 to exceed the acceptable confidence threshold [63]. When the value of Z score is lower than the 1.65 with Moran’s I value of negative, the spatial relationship will be either random or dispersed [64]. On the mentioned dry season, areas in Kilang Gula Chuping and Felda Chuping show a significant hot spot of higher correlation where the agricultural drought is clustered using Z score from the Moran’s I spatial correlation (refer Figure 8). In Figure 8, the classification of the clustering relationship is shown based on the Z-score that relates with the p value to allow the understanding of the severity of drought in the area. The Z score of more than 1.65 values from p value interval shows that Kilang Gula Chuping region has the driest seasons in the mentioned years. The values that fall in the −1.65 to the 1.65 Z score are not significant when compared against surrounding the adjacent areas. In the Kubang Perun area, seen with lower Z score starting from −1.65, the representation of the values is seen to be disperse where distribution TVDI and NDWI is showing dispersed or random representation of pattern.

In Figure 8, the interval confidence of 99%, 90% and 95% was shown to emphasis that the values from all years, Kilang Gula Chuping and Felda Chuping are experiencing low NDWI values interpreting as low water content of soil and vegetation. Furthermore, a high TVDI value was recorded in reference to high temperature and dryness in the mentioned areas in regards to the proportional vegetation cover. When looking at the Kubang Perun, Guar Nangka and part of Sungai Buluh areas, represented in blue from Figure 8, the dispersed pattern was due to the sufficient water content in the areas. Furthermore, the area had moderately dry or normal TVDI values to show that the mentioned areas are moderately dry but sufficiently wet. This resulted in a spatial relation of a dispersed pattern.

In summary, using both indices of TVDI and NDWI, a better understanding of indices representing the surface water content utilizing the NDWI and the NDVI-LST relationship taking the vegetation cover by TVDI in Chuping was explored. The localized relationship for the mentioned area improved agricultural drought assessment provided a significant classification of very dry areas of the Chuping identifies using the Moran’s I spatial relationship. It is observed that Kilang Gula Chuping and part of Felda Chuping have experienced the most impact in all the dry seasons in the years 2015, 2016 and 2017.

4. Conclusions

This study has intended to identify the extent of agricultural drought over the Chuping, Perlis area using satellite-borne remote sensing data based on TVDI and NDWI. Both TVDI and NDWI were used in a spatio-temporal assessment for the month of March, which has regularly recorded high temperatures for the years 2015, 2016 and 2017. The Kilang Gula Chuping region, on the other hand, saw a substantial increase in TVDI values and was dry for all three years surveyed, with an average value of 0.70. When both the TVDI and the NDWI were examined, substantial clustering of spots was discovered in Chuping, Perlis for all three years, with local regressions of 0.8399, 0.6979 and 0.6974 for the years 2015, 2016 and 2017.

It is noted that TVDI index which uses both the NDVI and LST components, can successfully be used for recognizing the spatio-temporal range of agricultural drought. In addition, it can also be utilized to clarify drought severity classes in the exploration areas through merged analysis of both vegetation coverage and temperature of vegetation.

Moran’s I values showed that the research area had a high value of I for 2015 with a value of 0.6252, while the I values for 2016 and 2017 were 0.2962 and 0.2299 respectively, with corresponding Z scores of 17.806, 8.625 and 6.769 respectively for the years 2015–2017. This showed that the clustered relationship is significant in the 95–99% range of confidence interval. The combined TVDI and NDWI, which both correlate using the Moran’s Index, will lead to tracking the occurrence of agricultural drought as a warning system mechanism, where major areas affected can be detected. It is noted that when the precipitation level is recorded to have long dry spells, the use of TVDI will show how severe the agricultural drought is. Using the two indices, TVDI and NDWI alone were sufficient to understand and identify the drier spots of Chuping over 3 years. However, cloud-free satellite data will be useful in future studies to provide a better view of the excess water stress conditions during the rainy seasons in the area and sufficient rainfall data to support the approach. It is suggested that additional years with more frequent observations will further explore the method advocated in this research.

Author Contributions

Conceptualization, V.S., A.R.M.S. and Y.P.L.; Data curation, V.S.; Formal analysis, V.S.; Investigation, V.S. and A.R.M.S.; Methodology, V.S., A.R.M.S. and A.W.; Resources, A.R.M.S., A.W., M.R.K., Y.P.L. and W.T.; Software, V.S.; Supervision, A.R.M.S. and A.W.; Validation, A.R.M.S., A.W., M.R.K. and Y.P.L.; Visualization, A.R.M.S. and A.W.; Writing—original draft, V.S.; Writing—review & editing, V.S., A.R.M.S., A.W., M.R.K., Y.P.L. and W.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Japanese Aerospace Agency (JAXA), using the Earth Observation-Research Announcement (EO-RA2) collaboration with PI number of ER2A2N180 to obtained Global Agriculture Monitoring, JASMIN sensor rainfall data.

Institutional Review Board Statement

Not Applicable.

Informed Consent Statement

Not Applicable.

Data Availability Statement

JASMIN rainfall data can be obtained from Japanese Aerospace Agency upon request. Restrictions apply to Malaysian Metorological rainfall data that was obtained through a special request. Landsat 8 satellite data can be obtained in the National Aeronautics and Space Administration (NASA) website.

Acknowledgments

The authors would like to thank the Malaysian Meteorological Department for the rainfall data in the study site. The authors would like to thank National Aeronautics and Space Administration (NASA) for the Landsat 8 images used in the study.

Conflicts of Interest

The authors declare no conflict of interest in the process of writing this paper locally and internationally.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. A conceptual Land Surface Temperature (LST) with NDVI triangle [46] X-axis shows the increasing normalized difference vegetation index from 0 to 1. Y-axis shows the increase in land surface temperature (Celsius).

Figure 2. Study area shown using latitute and longitute derived from World Geodetic System 1984 (WGS 1984) expressed in degree unit.

Figure 3. Precipitation rates in Perlis per month for the year 2015, 2016 and 2017 and overall average.

Figure 4. Derivation of LST from Landsat 8/OLI.

Figure 5. NDWI Map for 2015 (A), 2016 (B) and 2017 (C).

Figure 6. Average values over 3 years (2015–2017) for Chuping, Perlis for the driest months: (a) NDWI; (b) NDVI; (c) LST where 2016 marked with “**” and 2017 marked with “***” shows that over 3 years, the LST significantly different from each other in all the areas and subarea.; (d) TVDI.

Figure 7. TVDI Map of Chuping for 2015 (A), 2016(B) and 2017 (C).

Figure 8. Regression analysis of NDWI against TVDI for the year 2015 (A), 2016 (B) and 2017 (C) in the month of March where clustered significant spots are identified with 90%, 95% and 99% accuracy are shown using Z Score.

Table 2

Information of length of dry Spell in Chuping using MMD data.

Year	Maximum Length of Dry Spell per Year (Days)
2015	20
2016	43
2017	22

Table 3

Spatial autocorrelation results to identify the clustered spot.

Year	Moran’s Index	Z Score	p Value	Variance
2015	0.6252	17.806	0.0001	0.0012
2016	0.2962	8.625	0.0001	0.0012
2017	0.2299	6.769	0.0001	0.0011