Mapping of soil properties is an important operation as it plays an important role in the knowledge about soil properties and how it can be used sustainably. The study was carried out in a local government area in Bangladesh in order to map out some soil properties and assess their variability within the area. From the study area, a total of 92 soil samples (0–20 cm) were collected from different cropping patterns at an interval of 2.2 × 2.2 km2 on a regular grid design. A portable global positioning system (GPS) was used to collect coordinates of each sampling site. Then, soil properties, that is, pH, electrical conductivity (EC), soil organic carbon (SOC), total nitrogen (Total N), and soil available nutrients (P, K, and S) were measured in the laboratory. After the normalization of data, classical statistics were used to describe the soil properties, and geostatistical analysis was used to illustrate the spatial variability of the soil properties by using kriging interpolation techniques in a GIS environment. Results show that the spatial distribution and spatial dependency level of soil properties can be different even within the small or large scale. According to cross-validation results, for most soil properties, the kriging interpolation method provided the least interpolation error. The generated maps of soil properties that indicate soil nutrient status over the study region could be helpful for farmers and decision-makers to enhance site-specific nutrient management. Also, these prototype maps would be helpful for future nutrient and fertilizer applications management, including a site-specific condition to not only reduce the cost of input management but also prevent any environmental hazard. It also demonstrates that the effectiveness of geostatistics and GIS techniques provided a powerful tool for this study in terms of regionalized nutrient management.

1. Introduction

Bangladesh, one of the most densely settled countries in the world whose economy, is largely based on agriculture which contributes 13.31% of Gross Domestic Product (GDP) at current prices [1]. The land is the principal resource which employs around 40% of the total labor force and feeds about 164.6 million of its population [2, 3]. The total cultivable land is estimated to be 9.10 million hectares with an average cropping intensity of 179 percent per year. High population growth (1.37 per year) with low growth in agricultural productivity adversely affects the living standard in the country [3, 4]. The country has one of the lowest land-person ratios in the world, estimated at 0.088 ha per person [3]. The number of agricultural farm households is estimated at 1.66 million, which accounts for 46.61% of total households [4]. There is huge pressure on the land to produce more crops to ensure self-sufficiency in food. According to the agricultural statistics database of the Food and Agriculture Organization of the United Nations [5], total cropland extent at the global scale, computed as the sum of arable land and permanent crop area, is about 15.3 million km2. These statistics account for all cropland used at least once in five years but neglect areas with long fallow periods. The total harvested crop area reported in the same database is 11.8 million km2 yr−1, indicating a global average cropping intensity of 0.77 crop harvests per year. However, the extent of fallow land is larger than the difference between global cropland extent and global harvested crop areas because many areas are harvested more than once per year.

Soil properties vary in different spatial areas due to the combined effect of biological, physical, and chemical processes over time [6] and can vary within farmland or at the landscape scale [7, 8]. Different land use and management practices greatly impact soil properties [9], and knowledge of the variation in soil properties within farmland use is essential in determining production constraints related to soil nutrients. With the introduction of commercial fertilizers, the physical properties of the soil were thus seriously neglected by both the farmers and soil scientists. This is especially true for Bangladesh, where the focus of soil science has mainly been concerned with soil chemistry and soil fertility. However, better crop varieties and increased use of fertilizers will fail to increase yield in the long run. The contribution of these inputs can only be realized if soil physical properties are properly managed [10]. Land use without adequate planning leads to soil impoverishment and low crop yield, which results in a decline in the socioeconomic and technological level of the rural population [11, 12]. Therefore, to undertake soil planning for use and management purposes, it is important to evaluate how the chemical and physical properties of the soil are distributed in a determined area.

A tool often used to analyze how soil properties are spatially distributed in an area is geostatistics. It is effective for understanding the magnitude and structure of the spatial variability of the physical and chemical properties [13]. The study of spatial variability of soil chemical and physical properties is important for agriculture because it aims to minimize the effects of variability on crop yield, optimizing the agricultural production systems [14, 15]. Meanwhile, the soil is a nonrenewable resource [16], and the perception of soil nutrient valuation becomes extremely significant for better agrarian productivity and the economic development of each nation. Correspondingly, inappropriate fertilizer rate application [17] and improper land management [18] contribute to the low nutrients contents in these soils [1921]. This had resulted in tenaciously muted crop yields. For improving crop yields and increasing the income of smallholder farmers, we need to increase the productivity of these soils. Also, the spatial distribution of the major soil nutrients needs to be mapped, and based on the outcomes of the distribution map, applicable fertilizer requirements could be recommended for localized intervention. Subsequently, the mapped outcomes of the soil nutrient valuation could then be used for effective monitoring of changes that might occur between cropping schemes and seasons over time. Then, monitoring of the nutrient concentrations will enable shareholders to evaluate soil fertility enhancement or else in such localities. Consequently, the costly and tedious conventional methods required to acquire soil nutrient information will also be abridged when nutrient levels are mapped since those conventional techniques are no more affordable [22, 23]. Therefore, mapping of the nutrient concentrations will provide spatial soil nutrients information that can be used as a decision support tool. Thus, developing spatial distribution maps of soil nutrients levels is incredibly necessary for the “diversified cropping” regions of which Jessore district of Bangladesh is one [24] because it will help improve agricultural management tractrices and refine sustainable resource use together with providing a baseline against which future soil nutrients can be recommended at site-specific localities [25, 26]. However, the mapping of soil nutrients concentrations, especially nitrogen (N), phosphorus (P), potassium (K), and sulfur (S), would also accelerate appropriate monitoring and review of recommended agricultural tools at localities from time to time. This may also help in the valuation of the influence of a specific technology at a specific time (e.g., every 10 years) in a specific place depending on the evaluation of the soil quality [27]. Likewise, due to the promising development of precision agriculture knowledge [25, 28, 29], scientists and decision-makers in the field of soil science would be in a better position to implement site-specific technologies, if accurate levels of soil nutrients in particular locations in the territory are mapped. This tactic will enhance soil fertility management outcomes and increase the interest of small-holder farmers to invest more in the agricultural sector.

The lack of any previous study makes it hard to assess the nutrient profile of the land (specific area) and suggest proper management protocols to farmers. Many smallholder farmers reside in the Jessore district in Bangladesh and have poor access to agricultural inputs, including organic and inorganic fertilizers. Therefore, the objective of this study was to generate appropriate soil properties models in order to generate a spatial distribution map of major soil nutrient (N, P, K, and S) contents across the Jessore district of Bangladesh. The results of the study would, therefore, reveal the spatial variation and pattern of distribution of N, P, K, and S nutrient contents across the study area. Also, their evaluation would help to make the appropriate decision to the application of inorganic fertilizer into the soil for sustainable crop production.

2. Materials and Methods

2.1. Study Area

The research was conducted in the Jessore region (23.1681 N to 89.2042 E), located in the western part of the Ganges River floodplain which is predominantly highland and medium highland. Most areas have a complex relief of broad and narrow ridges and interridge depressions, separated by areas with smooth, broad ridges, and basins. There is an overall pattern of olive-brown silt loams and silty clay loams on the upper parts of floodplain ridges and dark grey mottled brown, mainly clay soils on ridge sites and in basins. In this region, major cropping patterns are Rice-Mustard-Rice, Rice-Mustard-Rice/Pulse, Rice-Vegetables, Wheat-B.Aus/Jute-Fallow, Wheat-B.Aus/Jute-T.Aman, Wheat-Pumpkin-Orchard, Lentil-Sesame-T.Aman and Sugarcane-Boro-T.Aman, and so on. Soil fertility and productivity change over time and this change is in a negative direction because of intensive cropping with modern varieties, improper and imbalanced use of fertilizer and manure [30]. Again, crops grown in different cropping patterns and environments responded differently to fertilizer nutrients. Mineral fertilizer inputs are the crucial factors to the overall nutrient balance in intensive cropping systems [31, 32].

2.2. Soil Sampling and Analysis

A total of 92 soil samples were collected from the surface (0–20 cm) at an approximate interval of 2.2 × 2.2 km2 on regular grid design (Figure 1) with the help of a handheld global positioning system (GPS) over the entire Jessore Sadar (total area 435.22 km2) of Bangladesh. Regular sampling is proved to be more accurate in predicting spatial distribution than random sampling [33]. For each soil sample, two to three replicates were taken in spring 2017 under different cropping patterns (e.g., wheat, rice, corn, pulse, sesame, maize, lentil, groundnut, watermelon, vegetable, ladies finger, chili, pumpkin, orchard, and fallow) within a distance of 100 m and the samples were combined and homogenized by hand mixing. Nineteen (19) more points were sampled for cross-validation (Figure 1) purposes. The collected soil samples were air-dried by spreading on a separate sheet of paper and were transported to the laboratory. After drying, the larger aggregates were broken gently by crushing them in a wooden hammer. A portion of the crushed soils was passed through a 2.0 mm sieve to separate the coarse (>2 mm) and fine (<2 mm) fractions. The sieved soils were then preserved in a plastic container and labeled properly. These were later used for various chemical analyses. Afterward, soil pH was determined electrochemically with the help of a glass electrode pH meter keeping up the soil-to-water ratio 1 : 2.5 as recommended by Jackson [34]. The EC (electrical conductivity) of the soil was measured at a soil-to-water ratio of 1 : 5 with the help of an EC meter [35]. Total nitrogen of the soils was determined by the colorimetric method [36] following sulfuric acid (H2SO4) digestion as recommended by Jackson [37]. Available phosphorus was extracted from the soil with the help of 0.5 M NaHCO3(Olsen’s method) at pH 8.5 and the molybdophosphoric blue color method of examination was used for determination [38]. The available potassium was extracted from the soil with 1 N NH4OAc at pH 7.0 [34] and analyzed by a flame analyzer at 589 nm for determination [37]. Organic carbon of samples was determined by Walkley and Black’s wet oxidation method as outlined by Jackson [34].

2.3. Geostatistical Analysis

Descriptive statistical analysis was carried out first to determine the mean, maximum, minimum, standard deviation, and coefficients of variation (CV) of the variables of the data. In geostatistical analysis, the semivariogram was calculated for each soil variable as follows [39, 40]:where z(xi) is the value of the variable z at the sampled location xi, h is the distance lag in meters, and N(h) is the number of pairs of sample points separated by h. For irregular sampling, it is rare for the distance between the sample pairs to be exactly equal to h. Therefore, h is often represented by a distance interval. For the distance lag h, the semivariance is γ(h).

An experimental semivariogram plot was obtained by computing semivariances at different distance lags. In general, there are three important parameters useful for characterizing the spatial dependence of soil variables in the semivariogram plot: nugget (C0), partial sill (C), and range. Partial sill reflects the amount of spatial structural variance. Range expresses the distance at which the semivariogram stabilizes around a limiting value, and the nugget is defined as the variability at a scale smaller than the sampling interval and/or sampling and analytical error. The experimental semivariogram was then fitted with a suitable theoretical model: spherical, exponential, and Gaussian [4143]. The models provide information about the spatial structure as well as the input parameters for kriging interpolation. Kriging is considered as the optimal spatial interpolation for making the best linear unbiased estimates of regionalized variables at unknown locations. The spatial prediction of the value of a soil variable z at an unknown point x0is calculated as a weighted average followed [39, 40]:where is the value to be estimated at the location x0, z(xi) is the known value at the sampling site xi,, and λiis the weight. There are n sites within the search neighborhood around x0used for estimation, and the magnitude of n will depend on the size of the moving search window and on user definition. Kriging differs from other methods (such as inverse distance-weighted), in that the weight function λiis no longer arbitrary, being calculated from the parameters of the best-fitted variogram model under the conditions of unbiasedness and minimized estimation variance for interpolation. In this study, ordinary kriging was deployed for the interpolation of soil variables on a grid with a spatial resolution of 10 m.

2.4. Statistical Analysis

The descriptive statistics (mean, median, standard deviation and standard error of mean, skewness, and kurtosis) were calculated for all soil parameters using SPSS v.16 classical statistics. Before the construction of variograms, the data were tested for normality (Kolmogorov–Smirnov test). The pairwise Pearson’s correlation between the soil chemical properties was calculated in RStudio v1.1.463. Geostatistical analysis, including the development of sample variograms and kriging, was performed with the help of a geostatistical software, ArcMap 10.3. The degree of spatial dependency of each variable was decided with geostatistical methods using semivariogram analysis and kriging [44].

2.5. Cross-Validation of the Predicted Model

In cross-validation, the values estimated by ordinary kriging were compared to the values observed in the nineteen (19) sampling points in the study area. Subsequently, different indices, mean absolute error (MAE), root mean squared error (RMSE), and coefficient of determination (R2) between observed and predicted soil variables of the cross-validation were used to assess the performance of the kriging interpolation.

3. Results and Discussion

3.1. Descriptive Statistics of Major Soil Properties in the Study Area

The results about the Kolmogorov–Smirnov test indicated normality for most of the variables (Table 1) which suggests that they were all normally distributed. Also, it indicates that the constraint values strongly varied among fields. Thus, the majority of soil properties measured in the study area were similar in terms of mean and median values, while mean values are usually somewhat higher than the median (Table 1), which indicates dominant measures of central tendency. Many other kinds of research also noted similar results, including Brejda et al. [45], Cambardella and Karlen [46], Cambardella et al. [47], Emadi et al. [48], and Young et al. [49]. The coefficient of variation for all the variables observed was very different ranging within 3.35–76.54% at 0–20 cm depth. The lowest coefficient of variation was observed in pH with a value of 3.35%, which could be as a result of the uniform conditions in the area such as little changes in slope and its direction leading to a uniformity of soil in the area [48, 5053], while electrical conductivity (76.54%) had the highest variation at a specific depth. The higher variability of soil properties in terms of coefficient of variation is SOC, electrical conductivity (EC), available phosphorus, and potassium (CV ≥ 35%). By contrast, the lower variability (CV ≤ 15%) was observed for soil pH, while moderate variability (CV = 15–34%) was found for total nitrogen and available sulfur according to the guidelines provided by Warrick [54] for the variability of soil properties. Table 1also indicates that the skewness and kurtosis indices of the soil variables deviated considerably from the standard values of 0 and 3, respectively, at 0–20 cm depth. These variations in chemical properties are mostly related to the different soil management practices carried out in the study area, parent material on which the soil is formed, the role of the depth of groundwater, and irrigation water quality [52, 53, 55].

3.2. Correlations among Selected Soil Properties

Pearson’s product-moment correlation coefficient was calculated for each property to characterize the relationship among selected soil chemical properties (pH, SOC, Total N, Av P, Av K, Av S, and EC) (Table 2). In general, moderately significant correlations were found between some variables in the study area as reported, for example, SOC with Av P (r = 0.501, ), SOC with pH (r = 0.342, ), Av K with pH (r = 0.419, ), Av K with Av P (r = 0.413, ), and Total N with EC (r = 0.358, ). Significant correlation can also be identified between pH and Av P (r = 0.290, ). It is clear that SOC was strongly negatively correlated with EC (r = −0.652, ) and moderately correlated with Total N (r = −0.209, ). Moreover, more significant and stronger correlation was observed between SOC and Av K (r = 0.558, ) than other soil properties. The soil EC was found to be responsible for variable distribution of different soil minerals, and it was found that Av P, Av K, and Av S were negatively correlated with EC with r = −0.174, −0.173 and -0.235 (), respectively. In general, a moderate association between different soil chemical properties was observed.

3.3. Geostatistical Analysis
3.3.1. Semivariogram and Spatial Dependency

The possible spatial structure of the different soil properties was identified by calculating the semivariograms and the best model that describes these spatial structures was identified. Model parameters for the best fit semivariogram models are presented in Table 3 and Figures (24). The model with the best fit was applied to each parameter, the accuracy of soil property values can be estimated through kriging (using modeled semivariogram parameters) at unsampled locations was tested using different error estimates. The parameters were noted including the nugget effect (Co), the sill (Co + C), and the range of influence for each soil. Also, the degree of autocorrelation between the sampling points was found to be related to spatial dependency (nugget: sill ratio) and expressed in percentages. The spatial dependent variables were classified as randomly spatially dependent if the ratio (nugget: sill) was >75% and moderately spatially dependent if the ratio is between 25 and 75%, while it is strongly spatially dependent if it is <25% [47, 5660]. A large range indicates that the observed values of a soil variable are influenced by other values for this variable over greater distances than soil variables, which have smaller ranges [61].

Various models are used for semivariogram analysis. The spherical model was used to estimate the hypothetical semivariogram parameter in ArcGIS 10.3. The nugget value denotes the random variation usually derived from the accuracy of measurement or variations of the properties that cannot be detected in the sample range [60]. The sill value is representing the upper limit of the fitted semivariogram model [42]. The nugget to sill ratio implies the spatial dependence of soil properties. The range of the semivariogram denotes the average distance through which the variable semivariance reaches its highest value. A small effective range indicates a distribution pattern composed of small patches [62].

The nugget effect can be defined as an indicator of continuity at close distances. Soil properties with a lower nugget effect were generally defined by the spherical semivariogram model (Table 3). The soil parameters including pH, electrical conductivity (EC), soil organic carbon (SOC), and available potassium (Av K) follow a strong spatial distribution (except total nitrogen (Total N), available phosphorus (Av P), and available sulfur (Av S)) and clear patchy distribution all over the study area as the percent nugget value includes 7.33%, 22.92%, 1.47%, and 0.00%, respectively (Table 3 and Figures 24). Also, the value of the nugget effect for Av K and SOC was lower at a specific depth in the study area which suggests that the grid variance of variables is low. Thus, this result implies that the near and away samples have similar and different values, respectively. Therefore, the nugget effects that are small and close to zero indicate a spatial continuity between the neighboring points; this finding is similar to the result of Vieira and Gonzalez [63] and Jafarian and Kavian [64] which showed that semivariogram of Total N had very small nugget effect. The spatial variability of estimating soil properties varies within range and this model can calculate the unsampled within a neighboring distance of about 0.0383 degrees for pH, 0.0493 degrees for EC, 0.0688 degrees for SOC, 0.01734 degrees for Total N, 0.0402 degrees for Av P, 0.0414 degrees for Av K and 0.1593 degrees for Av S, respectively. Their spatial dependency might be controlled by both intrinsic variations of soil properties and extrinsic factors such as human-induced activities [65, 66].

3.3.2. Spatial Interpolation and Mapping of Soil Properties

Soil parameters maps obtained by ordinary kriging interpolation are presented in Figures 511. In addition, model parameters for the best fit semivariogram models are displayed in Table 3and Figures 24. The results showed that all the soil samples varied considerably all over the studied area. A high ratio (100%) of nuggets means that the slope of the semivariogram was close to zero, and the soil variable was considered nonspatially correlated (pure nugget). When the distribution of soil properties is moderately or strongly spatially correlated, the average extent of these patches is given by the range of the semivariogram. Again, if the ratio was <25% or the slope of the semivariogram was far from zero, that means a large part of the variance is introduced spatially, implying a strong spatial dependence of the variable. Their spatial dependency may be controlled by both intrinsic variations of soil properties and extrinsic factors such as human-induced activities [65, 66]. However, the unknown spatial dependence of the variable might exist at a lower scale even if a high nugget : sill ratio was obtained [67].

The map of different attributes (Figures 511) shows the spatial distribution map of soil properties generated using the best interpolation method (ordinary kriging interpolation). The maps of pH follow a strong spatial (Table 3) pattern and clear patchy distribution all over the studied area (Figure 5). The pH of the soil was lower in the southern part of the study area, while higher pH was found at the western site of the study area. Soils were usually characterized by having extensive areas with pH values between 7.35 and 7.82 (Figure 5), which indicates a mildly alkaline nature over the study area. Lake [68] reported that pH might be changed with current intensive agricultural practices and the relative absence of acidic cations which may raise pH above 7.

It was found that all over the study areas with a gradual decrease in electrical conductivity from northwest to southeast (Figure 6) site of the field. The electrical conductivity values of soil samples show a strong spatial dependency (Table 3) and an opposite trend like SOC. Soils were usually characterized by having extensive areas with concentrations from 0.39 to 1.38 dSm−1 (Figure 6) found almost all over the studied area. As the electrical conductivity of soils varies depending on the amount of moisture held by soil particles, thus its variability is distinct throughout the sampling area [69].

The map of SOC shows that almost two-thirds of the study area has a SOC content of less than 1%, which is the critical limit for organic carbon in most agricultural alkaline soils. The semivariograms and kriging interpolation map of SOC show a parallel (strong) spatial pattern trend as soil pH, while the opposite direction was observed for EC (Table 3). A lower percentage of SOC was found in the northwestern site of the study area, while a higher percentage of SOC was observed in the eastern and middle part of the study region (Figure 7). Also, soils were usually characterized by having extensive areas with SOC contents from 0.79 to 1.06% (Figure 7) in the study field. The major driving factors for low concentration of SOC over study area could be long-term cultivation practices [70], low water content, and high soil temperature [71]. The variability of SOC might be due to the land use pattern which might be the dominant factors of SOC in an area with the same parent material and the same climatic condition. Similarly, this variability may be based on landscape attributes including slope, aspect, and elevation [72].

Total N kriging interpolation map shows a moderate spatial (Table 3) dependency and clear patchy distribution pattern all over the studied area (Figure 8). The autocorrelation (semivariograms) analysis of the Total N displays a parallel (moderate) spatial scattering trend as Av S and Av P (Table 3). Also, soils were usually characterized by having extensive areas through Total N from 0.39% to 0.42% (Figure 8) noticed around the study area. Moreover, the higher contents of Total N were found only at the western site for this study. The moderate spatial dependence in Total N contents in diverse land-use patterns may be aided by the use or addition of NPK fertilizers, while grassland restoration is frequently encouraged by planting N fixing alfalfa [73].

According to the spatial distribution of Av P and Av K shows a clear patchy distribution all over the studied area (Figures 9 and 10). As expected from autocorrelation analysis (Table 3), the spatial distribution map of Av P shows moderate spatial continuity (small and high values of Av P occur near to each other). Also, soils were generally characterized by having extensive areas with available phosphorus (Av P) contents from 173.11 to 204.23 ppm (Figure 9). In contrast, the kriging map of Av K shows a strong spatial continuity of Av K with an increase along the northeast toward the center and southern part of the study area. Almost half the study region (west, northwest, and south) has an Av K less than 95 ppm (Figure 10), which is a critical limit for agricultural production [74]. In this study, the spatial variability of soil available potassium (Av K) may be influenced due to regional factors such as topography, climate, and soil matrix [39]. Also, we assume that an intensive cropping system and low application rate of K fertilizer could be the two major reasons causing severe potassium depletion in the observed area [7577]. The high and low amounts of Av P and Av K revealed that their spatial pattern might be also influenced by anthropogenic causes such as fertilization.

However, comparatively the lower content of available sulfur (Av S) was detected at the northern site for the present study area, while a relatively higher content of available sulfur (Av S) was detected at the southern site. The semivariogram of available sulfur (Av S) indicates that the moderate spatial correlation (Table 3) implies that a large part of the variance is introduced spatially. Also, soils were usually characterized by having extensive areas with available sulfur (Av S) contents between 95.45 and 121.19 ppm (Figure 11). In this study area, the lower contents of sulfur may be the consequence of predominantly washed-out sulfur in the form of sulfates, especially under leaching circumstances [69, 78].

3.4. Accuracy Assessment of the Predicted Model

The fit of the models of the semivariograms and kriging interpolation was based on the values of the coefficients of determination of the cross-validation near the coefficients of the straight line 1 : 1 (Figure 12), the lower values of MAE, and the higher values of R2 (Table 4), according to Robertson and Goovaerts [7981]. The spherical model has fitted to the variables pH, EC, SOC, Total N, Av K, and Av S, indicating that the behavior of these variables (soil properties) is less erratic on a small scale [82]. According to McBratney and Webster [83], these mathematical models best fit the soil properties. This also shows that semivariogram parameters obtained from fitting of the experimental semivariogram and kriging interpolation values were reasonable to describe the spatial variation of all the studied soil properties.

4. Conclusion

Assessing spatial variability and mapping of soil properties is an important prerequisite for precision agriculture because these maps will measure spatial variability and provide the basis for controlling it. It would also help in reducing the amount of inorganic (fertilizer) inputs being added to the soil in the form of supplements so as not to overburden the soil which can lead to pollution thereby degrading the land. Statistical analysis showed that the coefficient of variation for all the variables observed was very different ranging within 3.35–76.54% at 0–20 cm depth of soil. For geostatistical analysis of soil variables, the value of nugget: sill ratio ranges from 0% (Av K) to 64.10% (Total N), which indicates that internal (e.g., the soil-forming processes) factors were dominant over external (e.g., human activities) factors. However, the soil pH, EC, SOC, and Av K had a strong spatial dependency with a nugget : sill ratio of <25% since it was induced by structural factors. Meanwhile, other soil variables (Total N, Av S, and Av P) had moderately spatially dependency with nugget : sill ratio of 25–75% since these variables were mostly determined by both internal and external factors. The autocorrelation distances of all variables differed from 0.0383 degrees (pH) to 0.1743 degrees (Total N), which implies the sampling design is reasonable. Scattering maps, resulting from kriging interpolation, demonstrated that these studied areas were categorized by an opposite trend for Av S (higher in southeastern site) and EC (higher in northwestern site). By contrast, Total N and soil pH differed all over the studied area, whereas Av K and Av P had a high percentage at the center of the study area. The areas with low SOC contents (0.26–0.52%) and Av P (110.85–141.98 ppm) were mainly detected in the northwestern site of the study area. According to cross-validation results, kriging interpolation method provided the least interpolation error for most of the soil properties. It also demonstrates the effectiveness of GIS techniques in the interpolation of unsampled data. These results can be used to make recommendations of best agricultural management practices within the locality and also improve the livelihood of smallholder farmers.

Data Availability

The field and lab (raw) data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.