Analysis of spatiotemporal variations of drought and its correlations with remote sensing-based indices via wavelet analysis and clustering methods

Drought is considered as a natural disaster, effects of which accumulate slowly over a long period. During the past decades, frequent and severe droughts caused considerable ecological, agricultural, and financial loss around the world (Van Loon 2015; Saghafian & Sanginabadi 2020). Drought is generally difficult to define and understand due to its invidious nature (Vicente-Serrano 2007; Dadson et al. 2019; Li et al. 2020). The type and severity of drought vary spatially, and the resulting impacts depend on the vulnerability of affected sectors (Zhou et al. 2018; Kisi et al. 2019). Monitoring the spatiotemporal variation of drought makes it possible to reduce vulnerability and establish adaptation strategies (Nischitha et al. 2014). There are different techniques for quantitative analysis of drought phenomenon, and among them drought indices have been widely applied for this aim. For instance, Zhou & Zhou (2021) investigated drought duration, severity, intensity, affected area, and centroids during 1960–2015 for the Pearl River basin based on the Standardized Evapotranspiration Deficit Index and three-dimensional clustering algorithm and then revealed how these drought characteristics have affected net primary productivity. Li et al. (2020) used the Standardized Precipitation Evaporation Index time series for classifying drought grades in the Guanzhong area based on the meteorological and remote sensing datasets. Among drought indices, the Standardized Precipitation Index (SPI) is an ideal tool in drought severity characterizing (Li et al. 2015). Nowadays, this index is known as one of the most common methods for drought monitoring (Kwak et al. 2016). The SPI index has statistical consistency and is capable of showing short- and long-term drought impacts on precipitation anomalies. Although the drought indices provide useful information in drought assessment, they are derived from gauge observations and have several limitations, such as low spatial resolution and discontinuity. The satellite-based products can provide alternative datasets with higher spatial resolution, which are temporally and spatially continuous (Yilmaz et al. 2008). Therefore, monitoring of drought conditions is possible in ungauged or inaccessible areas via the satellite datasets (Skakun et al. 2016). Although the satellite-based products represented have been widely used in hydrological simulations, raw satellite datasets have some disadvantages, such as the amount of cloudiness which affects the accuracy of the results. Also, the data quality of these products will inevitably be affected by observational bias, spatial scale, and retrieval method. Therefore, satellite-based datasets need to be examined and verified by comparison with ground-based data. Satellite products have been shown to play an important role in assessing drought-related vegetation conditions, especially in semi-arid and arid regions. Owing to the close relationship between available soil moisture and vegetation vigor, two indices, the Normalized Difference Vegetation Index (NDVI) and the Land Surface Temperature (LST) index, have been widely used to assess drought conditions. The NDVI is an indicator of green biomass, leaf area index, and pattern of production. According to Zhou et al. (2012), the NDVI principle is that the internal mesophyll structure of healthy green leaves strongly reflects near infrared radiation and leaf chlorophyll, and other pigments absorb a large proportion of the red visible red radiation. This becomes reverse in the case of unhealthy or water-stressed vegetation. Also, the LST is used to show the energy balance at the surface of earth, one of the important variables in the physics of land-surface processes on regional and global scales.

On the contrary, within the past few decades, the time-series pre-processing approaches are commonly applied for hydrometeorological data analyses, given the seasonal fluctuations and non-stationary features of signals (Roushangar et al. 2021a). The Discrete Wavelet Transform (DWT) is a common approach used by many researchers for time-series decomposition (Agarwal et al. 2016). The modified version of the DWT is termed as the Maximal Overlap Discrete Wavelet Transform (MODWT). In this latter approach, time series should be projected onto a set of non-orthogonal basis functions in order to generate the wavelet coefficients. Another technique used for signal decomposition is the Ensemble Empirical Mode Decomposition (EEMD) method. According to Rezaie-Balf et al. (2019), this method is an intuitive, empirical, and self-adaptive data processing tool. The EEMD method is especially appropriate for analysis of nonlinear and non-stationary time series, which comes out as a superiority of the EEMD. In addition, clustering methods are used to identify homogenous drought regions. Clustering can be defined as the classification of a large number of data into a smaller number of groups so that data belonging to the same cluster are as similar as possible, while they have heterogeneity with the data belonging to other groups (Barton et al. 2016). The clustering approaches regard the spatiotemporal domains of components; so, the obtained outcomes from them indicate the local characteristics and temporary variations. The K-means and self-organizing map (SOM) methods are among the most common clustering methods, which were adopted in the present study.

Drought is a quickly increasing risk that can lead to hardly repairable or unrepairable damages to the nature as well as socio-economy, especially for arid and semi-arid regions. Analyzing the spatiotemporal pattern of drought can provide a better understanding of drought characteristics in these regions. Meteorological data mainly originate from a ground station network, and the quality and stability of them cannot be guaranteed due to the scarcity of stations, uneven spatial distributions, limited time scales, and vulnerability to environmental and human factors. Since meteorological data, such as precipitation and land air surface temperature, collected by surface observation stations often possess poor spatial resolution, especially in remote regions with difficult access and in some developing countries, satellite-based data and the remotely retrieved NDVI and LST data may provide a valuable source of information for monitoring drought. Therefore, in this study, a multi-spatiotemporal method based on both in situ observations and satellite datasets was proposed to analyze the drought correlations with the satellite sensor datasets in the northwest Iran. In this regard, first, the 3-month SPI series were calculated via (1) precipitation data obtained from 29 sites located in the selected area and (2) the CPC Merged Analysis of Precipitation (CMAP) satellite covering the period of 1990–2020. In the second step, the MODWT decomposition was applied for obtaining the temporal features of the signals, and further decomposition was performed via the EEMD. This step is generally termed as the pre-processing of the data. The time-series decomposition into various periodicity scales with the serial application of MODWT and EEMD results in more stationary subseries. Then, the energy values of each subseries were computed, and the mean energy values of each station were used for spatial categorization (i.e. clustering) of the area. Two clustering techniques, namely the K-means and SOM, were used for this aim.

In the next step, three satellite-derived NDVI and LST products, namely the Vegetation Condition Index (VCI), the Temperature Condition Index (TCI), and the Vegetation Healthy Index (VHI), were used to evaluate drought conditions. The NDVI and LST data were obtained from the MODIS satellite datasets MOD13Q1 and MOD11A2, respectively. The Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data series used for the present study covered the period of 2000–2020. From the remote sensing technologies, MODIS, Long-term Data Record, and Advanced Very High Resolution Radiometer (AVHRR) datasets have been extensively used to assess vegetation conditions and drought all over the world (El-Vilaly et al. 2018). Among these datasets, MODIS is more popular because of its heritage sensors in terms of spatial, spectral, and temporal resolutions. A major advantage of using the MODIS imagery is the availability of a suite of products ranging from raw images to highly processed products such as vegetation indices. The advantage of MODIS data over National Oceanic and Atmospheric Administration (NOAA)–AVHRR data is that the MODIS data provide continuity for historical applications of time series. The MODIS NDVI data are computed from atmospherically corrected bi-directional surface reflectance that have been masked for water, clouds, heavy aerosols, and cloud shadow. Therefore, in this study, the MODIS satellite datasets were used for calculating the NDVI and LST time series. On the other hand, the appropriate correlation of precipitation data of the CMAP satellite with the ground-based data has been proved in a previous study of the authors (Roushangar et al. 2021b). Finally, for further evaluating the relationships of SPIs with the NDVI and LST products, any possible links between the SPIs of the central stations of clusters and satellite-based indices were investigated using the wavelet coherence (WTC) method. To the best of the authors’ knowledge, there are no known studies dealing with the relationships between the SPIs and satellite-based indices in the selected area. The novelty of this study is that for the first time, we are able to determine the statistical relationships and possible links among the applied drought indices in the northwest Iran. This study can improve the insight into the relationships between the SPI and the NDVI, LST, VCI, and TCI, which helps decision-makers manage droughts and provide guidance for drought-related studies across other regions.

In this study, the spatiotemporal variations of drought were assessed for the northwest Iran using both in situ observation and remote sensing-based datasets. As shown in Figure 1, the procedures of the multi-spatiotemporal model could be summarized as follows:

• Obtaining the precipitation data from the in situ stations and the CMAP satellite during the period of 1990–2020.

• Calculating the SPIs time series.

• De-noising the SPIs time series.

• Decomposition of the de-noised series into several subseries via the MODWT and further decomposition of subseries with the EEMD (the pre-processing).

• Calculating the energy values of each subseries.

• Zoning the selected region based on the energy amounts of the subseries and the mean SPI values and investigating the relationship between them.

• Spatially identifying the homogenous drought regions via the K-means and SOM clustering techniques.

• Validating the outcome of spatial clustering methods using the following three indices: the Davies–Bouldin Index (DBi), the Dunn Index (DI), and the Silhouette coefficient (SC).

• Selecting the central station of each cluster.

• Obtaining the NDVI and LST time series from the MODIS satellite covering the period of 2000–2020.

• Calculating the VCI, TCI, and VHI indices from the NDVI and LST time series.

• Determining the possible links between the central stations of clusters (ground-based data) and satellite-based drought indices. The applied methodology is explained in more detail in the following sections.

Owing to the high variability of annual precipitation in the northwest part of Iran, drought is more frequent in this area; therefore, this area is selected for assessing the spatiotemporal variation of drought. The precipitation datasets of 29 sites obtained from the meteorology portal of Iran (http://www.weather.ir) were used for this aim. Beside the in situ observation data, the CMAP remote sensing product was employed to compute the SPI time series for the selected stations. According to Xie & Arkin (1996), the CMAP product is a collection of precipitation dataset, which is generated using remote sensing-based precipitation estimations and gauge data analysis. Also, the NDVI and LST series were obtained from the MODIS on the National Aeronautics and Space Administration (NASA) Earth Observation Satellite (EOS) Terra (https://modis.gsfc.nasa.gov/). Figure 2 shows the geographical locations of selected sites.

The empirical mode decomposition (EMD) is another efficient method for pre-processing of signals. In the EMD, the signal is broken down into several inherent mode functions (IMFs), which makes it possible to process nonlinear time series. The capability to specify the instantaneous frequency of the signal is one of the benefits of EMD. The EEMD is the modified version of the EMD. According to Wu & Huang (2009), the EEMD is used for solving the EMD mode mixing issue. Both the EEMD and the MODWT facilitate incredible vision in frequency and time domain of data to capture the seasonal and nonlinear properties.

Clustering is the task of dividing a set of patterns or data points into a number of groups, such that data points in the same groups are more similar to other data points in the same group than those in other groups. In simple words, the aim is to segregate groups with similar traits and assign them into clusters. Among different clustering techniques, the K-means generally applied for data-mining cluster analysis and vector quantization (Pham et al. 2005). Using the K-means, n observations are classified into k clusters and in this process, each entity (observation) is assigned to the cluster with the nearest mean (Pham et al. 2005). Based on Likas et al. (2003), the following steps should be performed in data clustering through the K-means method:

• Selecting an initial number of clusters and defining the centers;

• Assigning each observation to a cluster with the closest center considering its Euclidean distance;

• Reordering the positions of the centers after assigning all observation to one cluster;

• Repeating steps 2 and 3 until the memberships of cluster do not vary.

Given the ability of K-means in clustering large datasets as well as its implementation in a short time, it is one of the most popular clustering methods. Additionally, the SOM is another commonly used approach in exploring and extracting the inter-relationships of higher dimensional datasets. This method is useful in data clustering and anticipating in a widespread range of fields (Cereghino & Park 2009). The SOM has high ability in extracting implicit patterns from high-dimensional dataset and grouping the captured patterns into a low-dimensional output layer (Hsu & Li 2010). In this process, similar inputs remain close together in the output neurons, while the structure of data is preserved. Since the neurons of the output layer are usually ordered in two-dimensional grids, the produced topology can be visualized to provide an insight into the considered system (Kalteh et al. 2008; Iwashita et al. 2018).

The outcomes of spatial clustering were evaluated based on the following three validation metrics: DI, DBi, and SC. The DI goal is to distinguish a category of clusters that are well-set, with a small variance among components of the cluster, and well-detached, where the averages of the various clusters are adequately far apart when compared to the within-cluster variance. A higher DI value shows a better clustering outcome as it shows a well-compacted cluster (Agarwal et al. 2016). Computational cost increases when the number of clusters and dimensionality increases, which is a disadvantage for the DI. The DBi is a widely applied internal evaluation criterion, which is applied to distinguish the number of optimal clusters that are well-detached and well-set based on the content and specification of dataset. A lower DBi value represents better clustering results. On the other hand, the DBi has a disadvantage in that best information detection cannot be implied by a good reported DBi value (Kasturi et al. 2003). Also, the values of SC indicate the degree of similarity for samples within the cluster. The range of SC varies between −1 and 1, and higher values indicate that the samples are well-matched to their own clusters (Hsu & Li 2010).

For assessing the spatiotemporal variations of the drought time series in the considered region, first, the MODWT and EEMD methods were applied to decompose the de-noised SPIs into several subseries. Since different internal and external factors may affect the quality of data during the data collection process, the wavelet transform was applied to de-noise the original SPIs datasets. To decompose a time series via the MODWT, at first, an appropriate wavelet function (a function with similarities to the original signal as much as possible) should be selected. Therefore, different wavelet functions (i.e. Daubechies, Symlets, Haar, Coiflets, Meyer, Bior, and Rbio) were tested. On the other hand, the optimum decomposition level (L) was determined via a trial-and-error process. For this purpose, the best wavelet filters for L = 4 were determined and then, by changing the value of L parameter, its optimized amount was sought. The root-mean-square error criterion was used to select the best case. From Figure 3, it is seen that the db4 resulted in the best outcome among the applied wavelet filters. Also, the optimized decomposition level was found to be L = 4. Khan et al. (2020) used wavelet-based hybrid models for meteorological drought modeling and stated that among different mother wavelets, the debauches (db) type led to better results. This wavelet function has full scaling and translational orthonormality properties and also signifies that the wavelets have non-zero basis functions over a finite interval.

After signal decomposition, five subseries (i.e. four wavelet coefficients (W_i) and one scaling coefficient (V₄)) were extracted. In Figure 4, the decomposed subseries of the SPIs obtained from the in situ observations and CMAP satellite data are shown for the S18 station. As seen in this figure, the lower wavelet levels expressed higher frequency components, while the higher levels showed components with lower frequencies (including the trend of signal). Each wavelet coefficients represent a specific time period. For example, for monthly datasets, the W₁, W₂, W₃, and W₄ show 2¹ = 2, 2² = 4, 2³ = 8, and 2⁴ = 16-month periods, respectively. In the second step, the first two subseries of wavelet coefficients (i.e. W₁ and W₂) were further decomposed via the EEMD, in order to have more stationary subseries. Using the EEMD, the signal is decomposed into several IMFs and one residual subseries. The sum of these subseries will be the same original signal. Owing to the increase in the number of inputs after signal decomposition, for decreasing the number of inputs, the energy values of subseries were computed and applied as inputs of the clustering techniques.

In the following step, the selected region was zoned via the obtained mean energy values to investigate the relationship between energy values and the SPIs time series. Figure 5 shows the results of the area zoning using the mean energy values and the mean drought index for both in situ observations and the CMAP satellite datasets. It should be noted that in the zoning process, the inverse distance weight (IDW) was used for spatial interpolation. The interpolation methods can be classified into the following two major groups: deterministic and geostatistical. Deterministic interpolation techniques create surfaces from measured points, based on either the extent of similarity or the degree of smoothing. A deterministic interpolation can either force the resulting surface to pass through the data values or not. On the other hand, geostatistical (or stochastic) methods capitalize on the spatial correlation between neighboring observations to predict attributed values at unsampled locations. The IDW is a determinist method and estimates the values of an attribute at unsampled points using a linear combination of values at sampled points weighted by an inverse function of the distance from the point of interest to the sampled points. The assumption is that sampled points closer to the unsampled points are more similar to it than those located further away. This method has been widely used in many different fields, such as hydrology and earth science (Jeong et al. 2020). According to Figure 5, for both datasets, the southeast part suffered from severe drought, while drought was slighter in the southwest and northwest parts of the study area. In the southeast part of the selected region, the SPI amounts were higher than the other parts, and in the southwest part the SPI values reduced and reached the lowest value. The lowest energy values were obtained for the central parts, where there were more arid regions. Compact counter lines were observed on the northwest parts of the selected area, meaning that energy amounts varied rapidly on these parts which were located on rainy regions of the area. Spatial changes of energy became smoother for the eastern zones which were semi-arid regions. From these results, it is clear that there is an inverse relationship between the drought index and energy values, i.e. the smaller the energy values, the higher the intensity of the drought index is.

To cluster 29 stations in the selected area, the mean energy values of each station were used as inputs of the K-means clustering technique. According to Figure 6, three SC, DBI, and DI evaluation criteria were applied for determining the optimal number of clusters (NC). The DI is defined as the ratio of the minimal intra-cluster distance to the maximal inter-cluster distance. This index is an internal evaluation scheme for evaluating clustering algorithms, where the result is based on the clustered data itself. A large value of the DI is preferred as it represents a better separation. The DBi determines an optimal cluster number, such that clusters are well separated and compact. Since the objective is to obtain clusters with minimum intra-cluster distances, a small DBi value is important. Also, if most members have a high SC value, then the formation of the clustering is suitable. On the contrary, if many members have a low or negative SC value, accordingly the clustering formation may have too many or too few clusters. The number of clusters to adequately represent the dynamic properties of drought data was determined to be between 2 and 8. Following the application of the K-means clustering technique with a training program which ran through 1,000 trials, the results obtained for various numbers of energy-based clusters were compared and the optimal number of clusters (Figure 6) was determined through the comparison of three validity indices. In the clustering process via the K-means approach, it could be seen that for the in situ dataset, a cluster number of 4 (–) yielded values of 0.89, 1.36, and 1.82 for the SC, DBI, and DI, respectively, and for the CMAP data, the NC = 4 yielded values of 0.77, 1.36, and 1.76 for the SC, DBI, and DI, respectively, indicating a better performance in determining the homogenous regions compared with other NC values. Therefore, the NC = 4 was selected as being the optimum value for the clustering of stations. Also, the SOM method was applied for clustering the 29 stations into a visible two-dimensional topology map. In this regard, different map sizes were tested (i.e. 2 × 2 to 10 × 10). The results indicated that in the case of in situ dataset, the NC = 5 led to better outcomes in which values for the SC, DBI, and DI were 0.78, 1.44, and 1.78, respectively. While, in the case of satellite dataset, the NC = 6 performed better in determining homogenous regions due to having 0.68, 1.45, and 1.57 values for the SC, DBI, and DI indices, respectively. However, based on the results, a failure has occurred in the SOM outcome with NC = 5. In this case, there were three clusters with only one station as membership of the cluster, and two clusters with >10 memberships, which is an undesired result. Considering both performance criteria and classification of stations outcomes, it was proved that the K-means performed better than the SOM; thus, the K-means method was selected for further analysis.

Decomposed subseries in the pre-processing step show different monthly scales; however, some of these components have no strong correlation with the original SPIs. Therefore, mean energy values of the subseries were computed and considered as inputs of the K-means for performing spatial categorizing. The dynamic properties of time series based on their energy values can be used to improve the efficiency of the clustering approaches. The geographical location of stations based on clustering through energy values is shown in Figure 7(a). In Figure 7(b), the members of each cluster and their center stations are presented based on the SC index. From Figure 7, it can be seen that some of the stations located in the same cluster are spread across the study area, indicating that the geographic contiguity is not the sole basis of clustering given the differences in energy values calculated for each station.

In this section, the WTC analysis was applied for assessing the relationships between in situ observations and satellite-based drought indices in the study area. The WTC was used to find out whether there was a significant oscillatory period between the selected drought indices. Since there were 29 stations in the selected region, the use of wavelet dependencies on all stations and all variables was not possible; so, the central stations of each cluster were considered as representative of all stations. After clustering the area, a number of stations are assigned to each cluster as memberships of that cluster, and one of them is the central station. The central station should have the lowest distance from the cluster's center and also have the highest SC value compared to other stations. These members are marked in yellow in Figure 7(b) for all clusters. Finally, four stations, namely S12, S13, S10, and S21, were selected as representatives of each cluster. The results of WTC analysis between the drought indices are shown in Figure 10. A relative phasing of two time series is provided on each of the panel in this figure via arrows for indicating the cause and effect relations between applied indices. Cold parts outside the significant areas (bluish regions) showed frequencies and time with no dependence in the time series. The causality relation between two time series can be assessed using the phase patterns. The direction of the importance of arrows can be described as the following:

• X(t) (SPI) and Z(t) (second time series, i.e. VCI, TCI, or VHI) are in-phase (positively related) when arrows point to the right.

• X(t) and Z(t) will be anti-phase (negatively related) when the direction of the arrows is to the left.

• Arrows pointing in other directions show leads or lags between two time series. For instance, arrows with the straight down direction denote that Z(t) leads X(t), while arrows pointing up signify that X(t) leads Z(t).

For having a correct sight about the WTC graph, only significant sections of the graph should be taken in account and the regions beyond the cone of influence should not be explicated (Araghi et al. 2017). From the results, the WTC obtained for the NDVI–LST revealed some important in-phase patterns and the LST led the NDVI at 8–16 months band. Also, anti-phase patterns could be seen at 16–32 months band, and the NDVI led the LST in this band. At 8–16 months band, both in-phase and anti-phase relationships could be observed. The WTC of the NDVI–TCI showed anti-phase patterns at 8–16 and 16–32 months bands, and the LST led the NDVI. Some small areas with anti-phase relationships were detected at 4–8 months bands, in which the NDVI led the LST. The WTC of the SPI–VCI showed that the arrows were in-phase in most cases and there were positive relationships between the SPI and the VCI. The SPI–TCI wavelet analysis showed that in the 8–16 months band, the SPI led the TCI in all stations. For the 4–8 months band, in the S10 and S13 stations, there were in-phase relationships of the SPI and the TCI; while, for the S12 and S21 stations, an anti-phase relationship was found. Also, at the 32–64 months band, in-phase patterns could be seen for S10, S12, and S13, in which the SPI led the TCI. For the SPI–VHI, some areas with in-phase relationships were observed at the 4–8, 8–16, and 32–64 months bands. At the 8–16 months band, the VHI led the SPI and at the 32–64 months band, the SPI led the VHI. In general, it can be stated that the applied satellite-based drought indices have significant effects on the SPI drought index in the study area.

Drought as one of the severe and frequently occurring natural hazards has negative impacts on human life. In this study, the spatiotemporal pattern of drought indices based on in situ observations and satellite datasets was investigated for the selected area. For meeting the objectives of this study, 29 in situ gauges as well as CMAP and MODIS satellite products were used. For having a correct vision of time-series decomposition via the MODWT, different types of wavelet filters were tested, and an optimal level for time-series decomposition was selected. The K-means and SOM clustering techniques were applied to regionalize the selected area. It was shown that the method based on the combined MODWT–EEMD–Energy–K-means technique is a robust approach for drought regionalization. The results showed that the stations located in a specific cluster were spread in different areas and the proximity of the stations to each other was not the sole basis of clustering. The satellite-derived indices including the NDVI, LST, VCI, TCI, and VHI time series were applied in order to investigate drought conditions and their inter-correlations. The results showed a negative relationship between the LST and the NDVI and between the TCI and the VCI and a positive relationship between the VHI and SPI time series. From the VHI time series point of view, the moderate and mild droughts occurred in the region during the selected time period. In the final step, the relationships between SPIs of the central stations of clusters and the NDVI and LST products were assessed via the WTC method. Based on the results obtained for the WTC analysis, the satellite-based indices have significant statistical correlations with the ground-based drought index within a certain period. The proposed method provides valuable references for drought assessment and management, specifically at the regional level. However, future work in the area can be done using a larger dataset that might yield better results. For the purpose of having more detailed and thorough investigations, it is also suggested to take other variables into consideration; for instance, soil moisture, other remotely sensed parameters, and large-scale processes. At some regions, the parameters that are correlated might be different based on the general precipitation cycle. Certain weather phenomenon such as the El Nino may also need to be taken into consideration.

All relevant data are included in the paper or its Supplementary Information.