Principal component and cluster analyses based characterization of maize fields in southern central Rift Valley of Ethiopia

2.1 The study area

The study areas are located between Lakes Shala and Abaya in southern part of central Rift Valley of Ethiopia. It encompasses area coverage of 1,021,332 ha residing in over 30 districts located in Oromia region (West Arsi), SNNPR (Halaba, Hadiya, Kembata, Wolayta, and Gedeo) and Sidama region (Figure 1).

Figure 1

Study area and weather stations.

2.2 Data source

A continuous time series of 43 available weather stations located within and in close proximity to the study area were used as source to collect daily rainfall data from January 1991 to December 2020 (Figure 1). Geo-referenced daily minimum and maximum temperatures, rainfall, relative humidity, sunshine hours, and altitude of 30 years were acquired from the National Meteorology and Climate Change Institute (the then National Meteorology Agency), Hawassa Branch Office. Rainfall amount for the months of each growing season and year were summed to form a series of accumulated rainfall data for each season October, November, December, and January (ONDJ), FMAM, and JJAS. Data standardization was applied to the entire database by ranking or counting events based on departures from the local climatology in units of standard deviation using the equation (1) [37]

where X is value, μ _reff is mean, and σ _reff is standard deviation. Subsequently, L-moments, dispersion, and central tendency were computed for all physiographic and weather variables.

2.3 Independence of samples, homogeneity of variance, and equality of means

Normality test for rainfall across locations and years was carried out by a procedure set by Mishra et al. [38]. Bartlett’s test was applied to rainfall data to verify the homogeneity of variance between the months of the year (12 variances – temporal variability) at each site and among the sites (43 variances – spatial variability) at each of the 12 months of the year following a procedure developed by Steel et al. [39].

2.4 PCA

PCA, a dimensionality reduction tool in multivariate analyses, was used for examining multidimensional data to reveal patterns between objects that would not be apparent in a univariate analysis. The PCA with Varimax rotation was aimed at finding the relative influence of each variable explaining the variance of the system in each separate cluster in ONDJ, FMAM, and JJAS seasons separately [40]. For analysis of the data over the southern central Rift Valley, let A be a(t × i) anomaly matrix of rainfall data over a series of t = 30 years and i = 43 stations. A system of t × i numbers (elements) arranged in a rectangular arrangement along r rows and s columns and bounded by the brackets [] or () was called an r by s matrix as in equation (2).

The elements of the matrix A represented departure of rainfall values from their mean values. PCA was performed to examine the coupling between predictors and predicand. The determinant |A − λI| when expanded would give a polynomial, which is called as characteristic polynomial of matrix A. If, A = [ a i j ] was a square matrix, then the characteristic equation of A was |A − λI| = 0. The roots of the characteristic equation were called eigen values of A. Since, A = [ a i j ] was a square matrix of order “n,” then there existed a non-zero vector X, where is given by equation (3):

such that AX = λX, then the vector X was called an eigen vector of A corresponding to the eigen value λ. λ (Lamda) was an eigen value of A and x is its corresponding eigen vector. The covariance or correlation matrix of the dataset was used to derive the coefficients a _jj , which are the eigen vectors. Alternatively, x becomes the eigen vector pertaining to the eigen value λ, and vice-versa. If A _it is the actual rainfall data at stations i (i = 1, 2…s) in the year t (t = 1, 2…r), the mean value is given in equation (4).

The centered data are given in equation (5) as follow:

The covariance matrix is shown below in equation:

The eigen values {λ _j} and the eigen vectors {ϕ _ij } of this symmetric matrix were extracted. PCs were a linear combination of the original variables used as predictor variable in multiple linear regressions [40] and composed in general as in equation (7):

This transforms the original time series d _jt into the new time series PC _jt , orthogonal variables. First few PC _jt ’s were generally sufficient to explain the variation in the original data. The eigen vectors {ϕ _ij } represented spatial patterns, while (PC) _jt could be used to understand temporal variability in the data. Determination of the number of components to be used in PCA was performed using eigen values measuring more than 1 [41,42].

2.5 Cluster analysis

For cluster analysis, a partitioning method was adopted using two algorithms i.e., first hierarchical method was used to divide a set of N data vectors into K non-overlapping clusters (K ≤ N) [43] using identified variables from PCA. Second, each data vector was associated with a centroid by similarity. The Ward’s agglomerative hierarchical algorithm [44,45] was used as dissimilarity measure [46] (equation (8)):

where d _ij is the Euclidean distance, n is the number of vectors, X _ik is the input image vector and X _jk is the comparison image vector. Ward’s algorithm [44] creates groups, minimizing the dissimilarity, or minimizing the total sums of squares within groups, also known as sum of square deviations (SSD). At each stage of the procedure, groups were formed in such a manner that the resultant solution bears the lowest SSD within groups. At these stages, the joints of all possible pairs of groups were considered, and the two resulting smaller SSD were grouped together gathering alike individuals [46].

2.6 Data analysis

The presence of significant difference among clusters was tested using identified variables from PCA with one-way fixed effect ANOVA F test (equation (9)) and comparisons of group means among clusters were carried out by factor analysis provided in SPSS [47].

where SSB is sum of squares between treatments, SSW is the sum of squares within treatments, (k − 1) and (N − k) are the degrees of freedom between and within treatments, respectively [48]. Homogeneity of variance among clusters was tested using Levene test statistic, and robust tests of equality of means were checked for those with homogeneous variance using Welch and Brown–Forsythe statistic [49]. Post hoc comparison of means among clusters was done using Bonferroni method, which is a family contrast comparison method, used to compare with m number of post hoc comparisons to assure m = k(k − 1)/2 possible comparisons where k = the number of groups being considered [50–52]. Finally LSD test was done using equation (10):

where P _Bonf is the P value from Bonferroni test, P _LSD is the P value from LSD, and m is the group mean of a given cluster.

2.7 Cluster mapping

The clusters were mapped in ArcGIS environment version 10.3 using Inverse Distance Weighting method for better visualization of the results [53]. The clusters assume value of an attribute z at any un-sampled point as a distance-weighted average of sampled points lying within a defined neighborhood. The map of clusters and their altitude, rainfall, soils, texture, and geology was drawn using a 1:125,000 scale.

2.8 Principal component regression

Considering pre-processed or identified variables through PCA as main inputs determining the spatial and temporal variations in the study area, without loss of generality, the regression equation can be given as in equation (11) [54]:

where k is the number of independent variables, ßj ( j = 1, 2, …., k) is the partial regression coefficient, xi is the ith independent variable, y is the dependent variable, and ἐi is the error term corresponding to yi.

3.1 Correlation between predictor variables and distribution of predicand

Strong, positive, and significant association was observed between mean rainfall with latitude (R = 0.556, P < 0.001), STDEV (R = 0.754, P < 0.0001), CV(R = 0.529, P < 0.001), and maximum rainfall (R = 0.937, P < 0.0001) in ONDJ season. However, the association of mean rainfall with mean temperature (R = –0.388, P < 0.05) and minimum temperature (R = –0.432, P < 0.005) was negative and strong in ONDJ period (Table 1). Thus, areas that exhibit higher temperature experience lower mean rainfall in ONDJ period in the study area. Conversely, areas existing in higher latitude experience higher rainfall in ONDJ season. Consequently those variables uncorrelated to mean rainfall in ONDJ period are longitude, altitude, T _max, and RF_min. These uncorrelated variables were used in subsequent PCA analysis for ONDJ period. This result is in line with those of Machiwal and his colleagues who elaborated rainfall correlation with geographical factors and statistical parameters [43].

In FMAM season, the mean rainfall was positively and strongly associated with STDEV (R = 0.603, P < 0.0001), CV (R = 0.651, P < 0.0001), and maximum rainfall (R = 0.876, P < 0.0001); however, the association of mean rainfall was negative with latitude (R = –0.413*, P < 0.05) (Table 1 and Figure 2). The increment in mean rainfall with decrement in latitudes was due to the increased effects of ENSO in low tropical areas, which agrees with reports of Bates et al. [55]. The association of mean annual rainfall was negative with latitude (R = –0.580, P = 0.01), CV (R = –0.206, P = 0.05), T _min (R = –0.352, P < 0.05), and T _mean (R = –0.291, P < 0.05). Thus, small increment in latitude results in decrement in the mean annual and FMAM rainfall. Conversely, small increment in latitude results in increment of the mean rainfall in ONDJ season. JJAS season marks mainly rainy period for all places in the study area due to which there was no statistical differece in the amount of rain across clusters (Table 1 and Figure 2). Thus, the predictor variables uncorrelated with mean monthly rainfall were T _max, altitude, and longitude.

Figure 2

Rainfall distribution in annual, ONDJ, FMAM, and JJAS seasons in the study area.

3.2 Eigen values of different components and their significance level

For ONDJ, three PCs explained 74.7% of magnitude of variance, whereas four variables explained 79.3% of variations for FMAM. For JJAS, three PCs explained 78.2% of variance in annual rainfall, while four variables explained 80.5% of variance (Table 2) using a cutoff point of eigen value 1. This was because the variables with higher eigen values were more important than those with smaller magnitude in explaining the variances and thus contain most of the information of the selected factors. These will be used subsequently as forecasting factors. The 25.3, 20.7, and 21.8% variability unexplained by the model in ONDJ, FMAM, and JJAS seasons meant that there were other factors apart from geographic, climatic, and statistical factors that could be used to elucidate the variation. Thus, the original 12 variables were reduced dimensionally to 3–4 new variables and were used subsequently to explain maximum amount of variance.

The maximum degree of variability in the dataset was expressed with the first four components in FMAM period. Thus, PC1, PC2, PC3, and PC4 explained 30.3, 27.8, 12.8, and 8.4% of variance, respectively (Table 2). FMAM is important because land preparation, planting, urea application, and weeding are carried for maize crop in the majority of the study areas. In all the seasons, PC1 was mean rainfall component with strong association with altitude and longitude, but strong negative association with T _min, T _max, and T _mean (Table 3) implying the dynamics of the ITCZ over the study area [56,57]. This was in turn due to a high elevation leading to low temperatures (an increase in 1,000 m in altitude leads to a decrease of 6.5°C in temperature). This result is in agreement with reports of Belay in the Beles basin of Ethiopia, who reported inverse relationship among mean temperatures and elevations [58]. PC2 was termed as temperature component and showed strong positive association with latitude but strong to moderate negative association with altitude. The variables that correlate most with PC1 are RF_mean (0.401), RF_max (0.478), and variance (0.400). The first PC is positively correlated with these three variables. Therefore, increasing the value of RF_mean, RF_max, and variance subsequently increases the value of PC1 (Table 3). The variables positively associated with PC2 were T _min (0.456), T _max (0.433), and T _mean (0.511). This summarizes the direction of eigen vectors.

3.3 Rotated factor loadings and communalities

After Varimax rotation, the PC results showed that T _mean (0.960), T _max (0.895), and T _min (0.784) had large positive loadings on factor 1 (Table 4). Hence, their increment is vital for homogenous clustering. Whereas altitude (–0.829) has large negative loading on factor 1 and thus acted as a base for heterogeneity. RF_mean (0.927) and RF_max (0.942) have large positive loadings on factor 2 (Table 4; Figures 3–6). Bilate, Bulbula, Aje, Alaba Kulito, Humbo, and Ropi had high loadings with PC2, whereas Dilla, Humbo, Seraro, and Shone locations had high loadings with PC1 (Figure 3). The lowest loadings with PC1 were Woteraresa and wondo genet whereas the lowest loadings with PC2 were Woteraresa, Tefer Kela, and Yirgalem (Figure 3). The mean rainfall in FMAM was also meager in places like Aje, Ropi, Seraro, and Abaya. Conversely Wolaita Sodo, Yirgalem, Tefer Kela, Hawassa, Wondoget, Abaro, Boditi, Shamana, Durame, and Maykote received quit good rain in FMAM season (Figure 3).

Figure 3

Score plot of observations for PC1 and PC2.

Figure 4

Bi-plot of variables used in PCA analysis.

Figure 5

Agglomeration schedule coefficients.

Figure 6

Dendrogram using ward linkage.

When the contribution of each variable to the linear combinations was considered, longitude and altitude had strong negative influence on PC2 (Figure 4). But T _mean had strong positive influence on PC2. Similarly, RF_min exerted weak influence on PC1 (Figure 4).

3.4 Cluster analysis

The agglomeration schedule coefficients end at 42 (right top) whereas the elbow of the graph breaks sharply at 38 (left). The difference between the top unit and elbow is four clusters (Figure 5). This result of agglomeration schedule coefficients was verified by employing dendrogram method that automatically generated four relatively homogeneous units in the study area (Figure 6). Locations in one cluster are more similar to one another than locations in another cluster within a given season with climatic, environmental, geological, and physiographic features. Thus, four clusters were formed namely cluster I (Shamana cluster), cluster II (Bilate cluster), cluster III (Hawassa cluster), and cluster IV (Dilla cluster) (Figures 5–7).

Figure 7

The location of the four clusters.

3.5 Map and attributes of clusters

Cluster II locations had largest area (388,193 ha) whereas cluster IV areas had smallest area (121,931 ha). The area coverage of clusters I (324,242 ha) and III (186,966 ha) were intermediate (Figure 7). The records of annual and JJAS rainfall were highest in cluster IV. The wettest cluster during ONDJ is also cluster IV with mean rainfall of 299.2 + 51.3 mm across the 4 months of ONDJ (Table 5). Cluster I areas include high altitude maize growing areas in Abaro, Duguna Fango, Damot Gale, Badawacho, Shala, Seraro, Haisawita, Shashamane, Telamokantise, Aleta Wendo, Woteraresa, and Wujegra. In these areas, tillage operation begins early in November and planting maize starts with some rainfall showers from mid-December to mid-January unlike the other clusters. The areas in cluster II include low-altitude areas in Humbo, Halaba, Seraro, Bedessa, Bilate Tena, Bilate, Alaba Kulito, Shala, and Bulbula. These areas obtained 758.7 ± 125.8 mm of rain on average every year. ONDJ period was critically dry in clusters II and III and hence, was unable to support agricultural operations including tillage practices. In the next FMAM season, all clusters receive fairly good rain, except cluster II (Table 5). Hence, planting maize is usually delayed up to the beginning of May in cluster II areas. The 3rd cluster also called Hawassa cluster includes mid-altitude areas of Badawacho, Hawassa, Seraro, Shashamane, Boricha, Balela, Halaba, Humbo, Wolaita Sodo, and Bedessa. The areas included in cluster IV were lowland to mid-altitude areas in Aleta Wendo, Dilla, Tefer Kela, Dara, and Chuko. Maize is usually grown as of March within the agro-forestry system in this area. Thus, a given administrative district had more than one clusters, and careful arrangements are required even within a district in planning extension service, input delivery, and agromet advisories. This finding is in agreement with elucidation of Befikadu et al. [59] who partitioned adjacent districts into Bilate lowlands, Wolaita Sodo midlands, and Boditi highland. The Bilate cluster represented maize production agro-ecology termed as semiarid lowlands of central Rift Valley or dry mid-altitudes by Abate et al. [1] and low moisture areas by Worku et al. [60]. The Hawassa cluster was one of the sub-moist mid altitude areas in central Rift Valley in which is also called mid-altitude sub-humid by Worku et al. [60]. The Shamana cluster stands for highland transitional moist areas with upper mid-altitudes category as reported by Abate et al. [1] or highland sub-humid areas as stated by Worku et al. [60]. Finally low to mid altitude areas adjacent to mountain chains of Sidama and Gedeo were grouped into Dilla cluster and represented as humid tropics by Shengu [61], moist lower mid-altitudes by Abate et al. [1], and as lowland sub-humid areas by Worku et al. [60].

The presence of diverse maize production and recommendation domains were also reported by Alemu et al. [62] who elaborated high altitude sub-humid, low altitude sub-humid, mid altitude sub-humid, and moisture stressed places in the country. Thus, the highland, mid-altitude, and lowland areas are grouped into four clusters having different agro-ecology, cropping season, and choice of varieties, and consequently require different input delivery and training arrangement.

3.6 Geology, parent materials, and soils

In Shamana cluster, soils are developed from volcanic sedimentary to lacustrine deposits and have relatively higher clay and dominantly belong to Vertic Luvisols, Eutric Cambisols, and Chromic Vertisols based on FAO soil classification [63] (Figure 8). However, in Bilate cluster, the most dominant soil is riverine or lacustrine alluvium derived from basalt, ignimbrite, lava, or ash with dark brown loamy coarse sands and sandy loams with often calcareous subsoil, and belong to Calcaric Fluvisol, Orthic Andosol, Ortic Phaeozems, and Chromic Luivisols. In Hawassa cluster, the dominant soils are Vitric Andosols, Eutric Cambisols, and Leptosols. The soils had lacustrine and pyroclastic deposits of sands and silts, interbedded with pumice (Figure 8). In Dilla cluster the soils are originated from basalt rocks, and had deep reddish to brown clayey to clay loam texture classified as Haplic Luvisols and Chromic Vertisols (Figure 8).

Figure 8

Geology, altitude, and textural class of soils in the four clusters.

Low bulk density and weak structure of soils in Bilate and Hawassa clusters render them vulnerable to erosion even on gentle slopes compared to those in Shamana and Dilla clusters. The altitude of the study area ranges from 1,174 m a.s.l. at Lake Abaya (Cluster II) to 3,160 m a.s.l. at Woteraresa (Cluster I) (Figure 8), which leads to dramatic variability in environmental conditions among clusters over relatively short distances.

3.7 Evaluation statistics among clusters for FMAM season

Each cluster was distinct and heterogeneous in climatic, physical, and statistical variables (Table 6). There were four clusters in FMAM season as dictated by altitude, RF_min, and T _min with 9, 14, 16, and 4 cluster members in clusters I, II, III, and IV, respectively. Cluster II had highest T _min (22.7°C), lowest RF_min (86.1 mm), and lowest altitude (1,701 m a.s.l.). Cluster I areas have lowest T _min (22.7°C) but highest altitude (2,180 m a.s.l.) (Table 6).

The ANOVA table showed that F-test was significantly (P < 0.001) different among clusters for T _min, RF_min, and Altitude (Table 7).

The mean separation based on Bonferroni technique showed that cluster I had significantly lower temperature by value of 3.76 and 2.05°C compared to clusters II and III, respectively. The difference in mean temperature between cluster I and cluster IV was not significantly different (P < 0.05). The significantly lower temperature in cluster I areas like Abaro, Bitena, Boditi, Haisawita, Mayokote, Shamana, Shone, Telamokantise, Woteraresa, and Wujegra could be due to higher rainfall received in the previous months that actually have cooling effect in the area. Cluster II areas with significantly higher T _mean compared to other clusters were Abaya, Alaba Kulito, Bilate, Bilate Tena, Bulbula, Chuko, Dilla, Felka, Hawassa, Humbo, Morocho, and Seraro. These were the areas extending from Lake Shala to Lake Abaya and lie in the central part of Rift Valley. The difference in minimum rainfall received in FMAM was significantly higher in cluster IV by 265, 267, and 246 mm in clusters I, II, and III, respectively. Similarly, significantly higher means of altitudes were measured in clusters I than cluster II (P < 0.001), III (P < 0.001), and IV (P < 0.05) by value of 479, 376, and 395 m, respectively (Table 8).

3.8 Seasonal rainfall prediction

As maize is planted in FMAM season and continues to grow in JJAS period, the two seasons are important for rain-fed maize production. Additionally, sea surface temperature and T _mean are not included in regression equation mainly because the former shows little variation across clusters [11] and the latter has well-established inverse relationship with altitude in Ethiopia [7,58], which necessitates reduction in input variables and associated redundancies. Hence, the percentage of the variance in mean rainfall is explained by local factors including altitude, latitude, and longitude. As depicted by regression of geographical and topographic parameters (predictors) on rainfall (predicand), coefficient of determination (R²) is greater than 0.50 for most of the seasons and clusters (Table 9). Latitude (ß₁) and longitude (ß₂) showed negative regression coefficients for most of the clusters across the seasons unlike altitude (ß₃) which showed positive regression coefficients. Hence, latitude showed significant and negative prediction on mean rainfall in third cluster in ONDJ and annual seasons, and also first cluster in annual season. Thus, areas in southern part of the study site (for instance, Dilla) would receive more rainfall compared to northern areas (for instance, Shala) (Figure 8). In FMAM season, cluster I areas lying in lower longitude and cluster IV areas located in lower altitude are going to receive higher rainfall (Table 9). Similarly, cluster II areas lying in relatively higher altitude are going to receive higher rainfall in JJAS season. Thus, there is a need to prioritize mitigation measures in cluster II areas lying in lower altitude and cluster I areas situated in higher longitude in the future. This agrees with findings of Teodoro et al. [64] who showed altitude followed by latitude as an important physiographic factor that influence the monthly rainfall behavior in dry, transitional, and rainy periods.