Delineation of Soil Texture Suitability Zones for Soybean Cultivation: A Case Study in Continental Croatia

Abstract

Soil texture is a vital criterion in most cropland suitability analyses, so an accurate method for the delineation of soil texture suitability zones is necessary. In this study, an automated method was developed and evaluated for the delineation of these zones for soybean cultivation. A total of 255 soil samples were collected in the Continental biogeoregion of Croatia. Three methods for interpolation of clay, silt and sand soil content were evaluated using the split-sample method in five independent random repetitions. An automated algorithm for soil texture classification based on the United States Department of Agriculture (USDA) in 12 classes was performed using Python script. Suitability classes for soybean cultivation per soil texture class were determined according to previous agronomic and soybean land suitability studies. Ordinary kriging produced the highest accuracy of tested interpolation methods for clay, silt and sand. Highly suitable soil texture classes for soybean cultivation, loam and clay loam, were detected in the northern part of the study area, covering 5.73% of the study area. The analysis of classification results per interpolation method indicated a necessity of the evaluation of interpolation methods as their performance depended on the normality and stationarity of input samples.

1. Introduction

Knowledge of naturally suitable locations for the cultivation of a particular crop type is a basis for effective natural resource management and sustainable agriculture [1]. The development of spatial models of cropland suitability analyses integrated with geographic information system (GIS) enabled decision-makers to develop a crop management system that increases yield and land productivity [2]. The selection of the relevant criteria is based on crop and site-specific requirements, typically compiled from previous studies or crop expert surveys in the study area. Land suitability analysis commonly integrates climate and soil criteria groups of specific land use [3]. Soil texture presents a vital criterion in the GIS-based soybean suitability analyses for the determination of locations with the highest potential yield increase [4,5]. Many soil physical and chemical properties are directly influenced by soil texture, affecting soil fertility and productivity [6,7]. Recent research proved the direct influence of soil texture on the important segments of the ecosystem, such as agricultural yield variability [8], soil microbial and structural recovery [9] and richness and diversity of soil bacterial communities [10]. Automation of classification algorithms allows more time-efficient and robust classification in comparison to the conventional manual approach [11]. Recently developed algorithms for soil classification cover a wide range of applications, most notably measurement of soil surface roughness [12], prediction of soil fertility [13] and estimation of degraded soils [14]. In addition, mapping of soil texture variability in both horizontal and vertical dimensions has been a subject of interest in previous studies [15,16]. The mapping of the GIS analysis results presents an effective tool in their distribution to the end-users, both for farmers on a large scale, as well as for land policy managers on a small (regional) scale. Such data representation does not require extensive knowledge of spatial data and therefore allows simple interpretation, especially when integrated into web GIS maps.

Interpolation overcomes the limited representation of the agricultural soil properties by discrete soil sampling. Instead, it incorporates these ground-truth data to predict soil properties and enable a continuous representation of the entire agricultural field [17]. Interpolation methods are generally divided into geostatistical and deterministic methods. According to [18], geostatistical interpolation methods model the uncertainty associated with a spatial estimation based on random function theory, while the deterministic methods use only the input sample values and mathematical functions for interpolation. Geostatistical methods also apply statistical functions to quantify the predictions’ uncertainty, in contrast to deterministic interpolation methods. The basic step in the implementation of the geostatistical interpolation methods is the variogram analysis for the fitting of the most suitable theoretical mathematical model to an empirical variogram [19]. The evaluation of the interpolation methods enables the selection of the optimal, most accurate interpolation method for a particular soil sample [20]. The interpolation accuracy of soil samples has a direct impact on the creation of land-use management plans and variable-rate applications in precision agriculture [21].

Soybean (Glycine max L.) is the main oilseed crop in the world, presenting a high-quality source of oil and protein [22,23]. The long-term projection by [24] predicts a 1.3% global annual increase of soybean yield, which does not meet the food demands in the world. The development of soybean suitability analyses is therefore necessary for the improvement of soybean yield, especially to land policy managers [25]. A study by [26] noted the importance of soybean introduction in the cropping system, allowing the effective exploitation of soil nutrients and the overall reduction of pesticides and herbicides in precision agriculture. Soil texture is proven to have a significant impact on soybean cultivation and yield [27,28]. Improvement of research regarding soybean yield enhancement and sustainability is essential to satisfy growing soybean global demands [29].

The purpose of this study was to delineate soil texture suitability zones for soybean cultivation using an automated soil texture classification algorithm in a GIS environment. The specific objectives of the study were to: (1) evaluate and select the most accurate method for the interpolation of clay, silt and sand from input samples, (2) automate the process of soil texture classification in a GIS environment and (3) provide a map of soil texture suitability zones for soybean cultivation at a regional scale.

2. Materials and Methods

2.1. Study Area

The study area is the Continental biogeoregion of Croatia, covering 30,863 km² and being one of the three national biogeoregions defined by the European Environment Agency [30] (Figure 1). The other two biogeoregions present in Croatia are alpine and Mediterranean biogeoregions. The study area is characteristically a moderately warm, rainy climate, as per Köppen’s classification. Soybean is gradually becoming one of the most produced crops in Croatia, with a 119% yield production increase between 2013 and 2016 as per Croatian Bureau of Statistics data. The terrain is dominantly lowland, at lower elevations. According to CORINE 2018 land cover data, agricultural areas represent the major land cover class with 54.0% of the study area. The second-largest land cover classes are forests and seminatural areas with 40.7% of the study area, followed by artificial surfaces (3.5%), water bodies (1.4%) and wetlands (0.5%).

2.2. Automated Classification for Delineation of Soil Texture Suitability Zones for Soybean Cultivation in a GIS Environment

The selected interpolation methods for evaluation in this research were ordinary kriging (OK), inverse distance weighted (IDW) and spline interpolation (SI). These methods were recommended for the interpolation of soil properties by [31]. The workflow of this research is shown in Figure 2. All spatial calculations in this study were conducted using open-source GIS software: SAGA GIS v7.4.0 for the samples pre-processing and spatial interpolation and QGIS v3.8.3 for interpolation results and soil texture map visualization. Python v3.7.4 was used for scripting as a base for the automation of soil texture classification. The selected coordinate reference system (CRS) was the Croatian Terrestrial Reference System (HTRS96/TM). The spatial resolution of interpolated rasters was set to 250 m.

Input soil sample data was downloaded using the Croatian Agency for the Environment and Nature soil data WFS service. Soil samples were collected in the field according to ISO 11272:2017 and were pre-treated according to ISO 11464:2006 [32]. A total of 255 samples were collected in the study area during the year 2016. The soil texture content was determined by sieving and the sedimentation method prescribed in ISO 11277:2009. Soil fractions (<2 mm) were reclassified to three classes regarding particle size: clay (<0.002 mm), silt (0.002–0.063 mm) and sand (>0.063 mm). Each soil sample collected in the field represented a set of clay, silt and sand data on the same location. These values were filtered by attributes from input samples into a separate vector point layer in a GIS environment. The split-sample method was applied for the creation of training and test datasets, randomly splitting input samples to training data (70%) and test data (30%). From the total of 255 collected soil samples in the study area, 179 training samples and 76 test samples were randomly split from the original set individually for each repetition. The split-sample method was performed in five independent repetitions, resulting in five training and five corresponding test sample sets. Descriptive statistics consisting of mean, standard deviation (SD), coefficient of variation (CV), skewness (SK) and kurtosis (KT) were calculated to quantify normality and stationarity of the training sets. These properties are regarded as necessary for the determination of interpolation parameters for each method [33]. Normality designated the distribution of sample values, quantifying its deviation to normal distribution. Stationarity represented the number of local variations in the sample values, where low local variations indicate higher stationarity. Input data normality and stationarity are regarded as necessary indicators for the selection of optimal interpolation parameters [34].

OK is the most commonly used geostatistical interpolation method and is thoroughly described in [19,35]. It is also considered as the best linear unbiased spatial predictor [36]. OK incorporates a variogram as a tool for the modelling of spatial dependence between the input sample values and the mutual distance of samples [37]. The variogram was modelled based on (1):

$$\mathsf{\gamma }\left(h\right)=\frac{1}{2N}{{\displaystyle \sum }}_{x=1}^N{[Z\left(x\right) – Z(x+h)]}^2,$$

(1)

where $\mathsf{\gamma }\left(h\right)$ is the semivariance, N is the number of sample pairs and $Z\left(x\right)$ and $Z(x+h)$ are the observation values at point x and at a point separated by h.

IDW is a deterministic interpolation method that predicts the values based on a weighted average of the input samples. The weights were assigned to each sample in an inversely proportional way to the distance from the predicted location [38]. The prediction principle of IDW is outlined in detail in [39].

SI is the deterministic interpolation method that performs the spatial interpolation by using the polynomial functions for describing the components of the interpolated area. Mutual fitting of those functions is performed, so that they join smoothly, consequentially resulting in the smooth representation of the interpolated area [17]. SI performs best in the case of gently varying surfaces, which is commonly associated with samples with high stationarity [31].

The three selected interpolation methods were applied to five randomly created training sets, resulting in 15 interpolation variations separately for clay, silt and sand soil content. The selection of interpolation parameters for all three interpolation methods was performed according to the evaluated normality and stationarity of the training sets. The existence and the range of the spatial variability of training samples and full samples were observed using Moran scatterplots, Moran’s I values and autocorrelograms [40]. All three tested interpolation methods were implemented in SAGA GIS software. OK was conducted with logarithmic transformation in the pre-processing due to lack of normal distribution of sample values for all soil parameters, as proposed by [41]. The starting parameters used for OK interpolation of training sets and final interpolation using full sample sets were a range of 170,000 m with 25 lags, with each lag covering a distance of 6800 m. Tested theoretical mathematical models for fitting to the empirical variogram were linear, square root, logarithmic, Gaussian and spherical. The selection of these models and corresponding ranges was performed based on the highest fitting coefficient of determination to the empirical variogram in ranges with existing spatial autocorrelation. The interpolation by the IDW method was performed by using the power parameter of 3 to prevent the emphasis of local variations in the interpolation results, due to moderate stationarity of training sets. SI was conducted by the thin plate spline process to ensure the smoothing of the interpolated surface to decrease the impact of higher local variabilities of the input samples [42]. The ranges of deterministic interpolation methods were selected based on observed spatial autocorrelation, amounting to maximum distances with continuous positive autocorrelation values.

Final clay, silt and sand soil content raster were interpolated with the full input sample data, using the most accurate interpolation method determined during accuracy assessment. Soil texture classification was performed in 12 classes according to the specifications of the United States Department of Agriculture (USDA) [43]. The classification process was automated using the algorithm realized in Python script (Figure 3, Appendix A). The automation in the GIS environment was developed to perform classification robustly and independently of the CRS, location extent, or spatial resolution of input rasters.

2.3. Accuracy Assessment of Interpolated Rasters

Accuracy assessment of the selected interpolation methods was conducted by the calculation of coefficients of determination (R²) and root-mean-square error (RMSE) for all interpolated rasters. These statistical values are considered as complementary for the evaluation of interpolation results for soil parameters [44]. The R² quantified fitting of interpolated values to the regression line, while RMSE was used to assess variation in data based on absolute fit. The R² values were calculated using a linear regression model. The higher R² and lower RMSE indicate a better model fitting. Formulas for the calculation of R² (2) and RMSE (3) were [45]:

$$R^2=1 – \frac{{{\displaystyle \sum }}_1^n{{(y}_i – {\widehat{y}}_i)}^2}{{{\displaystyle \sum }}_1^n{{(y}_i – {\overline{y}}_i)}^2},$$

(2)

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle {\sum }_1^n{(y_i – {\widehat{y}}_i)}^2}}$$

(3)

where $y_i$ is the measured (sampled) value, ${\widehat{y}}_i$ is the predicted value, ${\overline{y}}_i$ is the average of the measured value y, and n is the number of observations.

2.4. Determination of Soybean Suitability Zones According to Soil Texture Classes

Loamy soil is the most suitable for the majority of crops, generally containing more organic soil matter and retaining moisture and nutrients better than other soil textures [46]. Levels of the suitability of 12 soil texture classes for soybean cultivation were evaluated based on previous soybean-specific agronomic and land suitability studies. Soybean-specific agronomic studies that evaluated the effects of soil texture for various applications in soybean cultivation are shown in

Table 1. Soil texture suitability for various applications from previous soybean-specific studies.

Source	Soil Texture Classes in the Study	Application
[48]	loamy sand, sandy loam, silt loam	Nitrous oxide emission in agricultural fields
[49]	loamy sand, sandy clay loam	Cyst nematode population density
[50]	sandy loam, silty clay loam	Yield response to salinity
[51]	loamy textures, sandy textures, clay textures	Effect of root zone temperatures on interorganismal signal molecules
[52]	loamy textures, clay textures, sandy textures	Carbon and nitrogen mineralization
[53]	clay, sandy clay loam	Canopy dry matter production
[54]	sandy loam, silt loam, loam, sandy clay loam, clay loam, silty clay loam, clay	Incidence of brown stem rot in conservation-till fields

The most suitable soil texture classes per study were bolded.

. All evaluated studies noted a higher level of suitability for loamy textures compared to other classes present in the study area. Loamy soil textures with higher silt or clay content were regarded as more suitable compared to higher sand content, such as sandy loam or loamy sand. Recent studies based on GIS-based multicriteria land suitability for soybean cultivation containing soil texture as an influential criterion are shown in

Table 2. Soil texture suitability for soybean cultivation from previous GIS-based multicriteria analysis studies.

Source	Suitability Level for Soybean Cultivation
	S1	S2	S3	N1	N2
[55]	loam	sandy loam, clay loam, silt loam	sandy clay, silty clay	other classes	sand
[56]	loam, clay loam, silty clay loam, silt loam	sandy loam, clay, silty clay, silty clay loam	loamy sand, silt	sand
[57]	silt loam, silty clay loam, loam	sandy clay loam	sandy loam	/	/
[58]	silty clay loam, sandy clay loam, clay loam	loam, silt, silty loam, sandy clay, silty clay	sandy loam, clay	loamy sand	sand
[59]	clay loam	sandy clay loam	loamy sand	sandy loam	/

. The suitability of soil texture classes was natively represented according to Food and Agriculture Organization (FAO) specifications [47] or adjusted to this specification from original classes. FAO specifications for land suitability consist of five suitability zones: highly suitable (S1), moderately suitable (S2), marginally suitable (S3), marginally not suitable (N1) and permanently not suitable (N2). Similar to agronomic studies, loamy textures were considered as the most suitable, with clay loam and loam being regarded as highly suitable. Suitability levels generally decreased with the increase of clay, silt or sand soil contents. The final classification of soil texture suitability zones for soybean cultivation was performed according to FAO specifications and is presented in Figure 4.

3. Results

Descriptive statistics used for the estimation of training sample sets’ normality and stationarity are shown in

Table 3. Descriptive statistics of random training sample sets.

Set No.	Soil Parameter	Mean (%)	CV	SK	KT
	Clay	30.997	0.368	0.678	0.228
Set 1	Silt	58.141	0.223	−0.760	0.833
	Sand	10.862	1.207	2.322	6.515
	Clay	31.355	0.332	0.584	0.336
Set 2	Silt	57.945	0.211	−0.599	0.391
	Sand	10.701	1.150	2.065	5.093
	Clay	29.889	0.343	0.581	0.159
Set 3	Silt	57.958	0.225	−0.800	0.874
	Sand	12.153	1.141	2.060	4.861
	Clay	29.899	0.354	0.633	0.334
Set 4	Silt	58.276	0.225	−0.805	0.860
	Sand	11.826	1.173	2.116	5.000
	Clay	30.795	0.353	0.594	0.241
Set 5	Silt	57.305	0.228	−0.688	0.674
	Sand	11.900	1.151	2.082	5.050
	Clay	30.587	0.350	0.614	0.260
Mean	Silt	57.925	0.222	−0.730	0.726
	Sand	11.488	1.164	2.129	5.304

. Silt soil content was observed as the dominant soil texture parameter in the study area, followed by clay. Clay and silt soil content data resulted in similar CV values through all five repetitions, representing low-variance data. Sand soil content produced CV values higher than 1.100 in all repetitions, indicating a high-variance data. SK remained consistent for all three soil parameters through repetitions. Clay and silt soil content resulted in right-tailed and left-tailed moderate skewness, respectively. Sand soil content consistently produces the highest SK values of all soil parameters, indicating a highly skewed data. KT resulted in the most variability of the calculated descriptive statistics for all soil parameters. However, values in all repetitions indicated a considerable platykurtic distribution for clay and silt soil contents and a considerable leptokurtic distribution for sand soil content. Sand soil content data values indicate a highly skewed distribution and moderate data stationarity.

The Moran scatterplot per training set is shown in Figure 5. Positive spatial autocorrelation was observed for all three soil parameters, with mean Moran’s I resulting in 0.396, 0.259 and 0.333 for clay, silt and sand, respectively. Interpolation ranges for IDW and SI per training sample set are presented in Appendix B. Mean distances until negative autocorrelation values were 81,600 m for clay, 29,920 m for silt and 54,400 m for sand. Variograms and interpolation parameters using OK for five training sample sets are presented in Appendix B. Square root and Gaussian theoretical mathematical models allowed the best fitting on most occasions, being applied 6 times each.

OK resulted in the highest mean RMSE and R² values for all soil parameters and was determined as the most accurate interpolation method in this study (

Table 4. Accuracy assessment of the evaluated interpolation methods per soil parameter.

Set No.	Stat. Value	OK			IDW			SI
		Clay	Silt	Sand	Clay	Silt	Sand	Clay	Silt	Sand
Set 1	RMSE	1.85	3.04	1.79	3.13	2.91	2.51	3.50	4.80	3.79
	R2	0.786	0.632	0.740	0.456	0.560	0.645	0.490	0.328	0.444
Set 2	RMSE	2.01	1.80	1.52	2.51	2.06	2.68	3.19	4.15	2.90
	R2	0.733	0.727	0.846	0.720	0.703	0.680	0.475	0.437	0.597
Set 3	RMSE	2.70	3.59	1.60	3.00	3.28	3.28	5.50	3.85	3.77
	R2	0.655	0.469	0.842	0.598	0.422	0.537	0.313	0.385	0.497
Set 4	RMSE	2.78	3.08	2.45	2.75	2.58	2.62	3.24	3.27	2.61
	R2	0.538	0.630	0.673	0.552	0.632	0.637	0.441	0.472	0.528
Set 5	RMSE	2.13	3.98	1.70	1.81	3.14	2.42	4.03	3.66	2.40
	R2	0.788	0.503	0.820	0.809	0.528	0.710	0.446	0.392	0.580
Mean	RMSE	2.29	3.10	1.81	2.64	2.79	2.70	3.89	3.94	3.09
	R2	0.700	0.592	0.784	0.627	0.569	0.642	0.433	0.403	0.529

). It was followed by IDW as the second most accurate method according to both statistical values, while SI was the least accurate interpolation method. The interpolation results within the same soil parameter and interpolation methods generally produced high variability in accuracy values. Most notably, the difference between maximum and minimum values for RMSE was 2.26 for clay interpolation using SI, while R² was 0.353 for the interpolation of clay using IDW. All three interpolation methods produced the highest accuracy for the interpolation of sand, followed by clay and silt soil contents.

Autocorrelogram and OK interpolation parameters of full sample sets are shown in Figure 6. Deterministic interpolation methods were performed with maximum ranges of 81,600 m for clay, 27,200 m for silt and 68,000 m for sand, matching maximum distances with positive autocorrelation. The specific local variabilities of the interpolation results were analysed according to Figure 7. OK produced interpolation results with low extreme values, compared to IDW and SI results. Consequentially, local variability in all soil fractions was minimized in interpolation results. Generalization of the local variabilities for OK is particularly noticeable in the south-west part of the study area for sand soil content interpolation. IDW results displayed the effects of all local variabilities to some degree but restricted their impact on interpolated results on local surroundings of the training samples. SI was the most balanced interpolation method considering the effects of local variabilities. The effects of local variations were retained in the interpolation results, having higher impacts on the surroundings of the training samples than IDW.

Histograms of interpolation results enabled quantification of the effect of local variabilities in the interpolation results (Figure 8). IDW resulted in the narrowest value interval of the tested interpolation method while having the top frequency value for every soil parameter. In contrast, SI produced the most variability regarding the interpolation values interval. Standard deviations of sand interpolation resulted in the greatest variability, with 7.91 for OK, 9.90 for IDW and 12.80 for SI.

The correlation matrix of tested interpolation methods for the respective soil fractions is displayed in

Table 5. Correlation matrix of interpolation results using Pearson’s coefficient of correlation.

	Clay			Silt			Sand
	OK	IDW	SI	OK	IDW	SI	OK	IDW	SI
OK	1.000			1.000			1.000
IDW	0.869	1.000		0.826	1.000		0.818	1.000
SI	0.775	0.868	1.000	0.803	0.914	1.000	0.726	0.869	1.000

. The maximum correlation was achieved for silt interpolation between IDW and SI. The lowest correlation was observed between OK and SI for both clay and sand interpolation. The visualization of classified soil texture in soil texture triangles is performed in Figure 9. OK produced the lowest dispersion of soil texture classes, dominantly covering silty clay loam and silt loam. IDW and SI covered a slightly wider range of soil texture classes. This observation particularly refers to loam, which is present in ground truth data (255 soil samples), but under-classified by OK. Both deterministic methods showed a notable presence of results in silty clay and clay loam classes, however, these classes had a low presence in ground truth data.

The interpolation of final clay, silt and sand soil fractions was performed by OK, as the most accurate interpolation method, with all 255 input samples. The results of the soil texture classification using Python script are displayed by coverage per soil texture class (

Table 6. Coverage per soil texture class in the study area.

Coverage	Silty Clay	Silty Clay Loam	Clay Loam	Loam	Sandy Loam	Silt Loam
Area (%)	6.04	53.26	0.99	4.74	0.01	34.96
Area (km2)	1861	16,412	306	1460	2	10,773

). Loam and clay loam, as the two most suitable soil texture zones for soybean cultivation, were detected in the northern part of the study area, mostly covering the border area near the river Drava. Silt loam and silty clay loam, the two moderately suitable zones, covered 88.22% of the study area. Silty clay, as the only marginally suitable zone, was mostly present in the south-east on the study area. The soybean suitability thematic map according to soil texture classes is shown in Figure 10.

4. Discussion

The first objective of the study regarding the selection of the most accurate interpolation method was fulfilled by accuracy assessment based on five independent random split samples. Clay and silt possessed values moderately deviated from a normal distribution and moderate data stationarity, while sand possessed highly skewed values. Previous studies [60,61,62,63] suggested that these cases are common in analyses of soil properties, so the transformation of skewed data was determined as a necessary procedure in the algorithm. The application of logarithmic transformation in the pre-processing to OK minimized the effect of skewed distribution, as suggested by [64]. The descriptive statistics of clay, silt and sand in this study suggest that the thresholds for performing logarithmic transformation could be CV > 0.500 as the first condition, and the second condition that SK and KT deviate 0.500 or more from values of the normal distribution sample of zero and three, respectively. However, it is necessary to test this assumption on the different soil samples to evaluate its reliability. Similar studies imply the applicability of statistical tests like the Shapiro–Wilk [65] or Kolmogorov–Smirnov test [66] for the evaluation of data normality and could present an upgrade to the soil classification algorithm in future research. Sparse sampling density of one sample per 118 km² used in this study is expected to have an impact on overall interpolation accuracy results. Denser soil sampling is expected to result in the lower CV of soil parameter values [67] and more accurate interpolation results [68,69]. However, the mapping of soil properties on a large scale was successfully performed using low-density samples, using one sample per 199 km² in [70] and even one sample per 443 km² in [71]. Due to the structure of the workflow and automatic soil texture classification algorithm in a GIS environment, a similar performance is expected by adding more samples in the existing area or using the same number of samples on a smaller area. Slightly more time-expensive computation in these cases is the only notable difference expected, as more input data should be processed.

The highest fitting R² values for OK were achieved on shorter distances with full sample sets (71,400 m, 54,400 m and 44,200 m, compared to 84,320 m, 69,360 m and 119,680 m) in contrast to mean values of five training sample sets, which correspond better to observed spatial autocorrelation. OK was determined as the most accurate interpolation method for all three soil fractions, with R² ranging from 0.592 to 0.784 and RMSE from 1.81 to 3.10. IDW was the second most accurate method, resulting in the most uniform accuracies through the interpolation repetitions. Moderate interpolation accuracies imply the necessity for denser soil samples in soil texture classification on a national level [72]. Sampling using regular grids proved in most cases to be the optimal suitable field sampling method [73]. The evaluation of interpolation methods displayed the importance of testing multiple methods. Mean R² differences between OK as the most accurate and SI as the least accurate methods were 0.267 for clay, 0.189 for silt and 0.255 for sand. Mean RMSE differences showed the same trend, amounting to 1.60 for clay, 0.84 for silt and 1.28 for sand. The correlation matrix supported the previous statement, as the lowest correlation was observed between OK and deterministic interpolation methods, especially SI. OK also produced the most accurate interpolation results in a similar study, compared to deterministic interpolation methods [31]. However, a study by [34] proved that in the case of low sample data normality and stationarity, deterministic methods like IDW produce more accurate results than OK. The accuracy assessment in five independent repetitions was proven to be effective, as it allowed testing of interpolation performances in variable conditions. The interpolation results of IDW and SI were suitable for maintaining the effects of local soil variabilities, however the same can be achieved by the modifications of OK parameters if necessary [74].

The fulfilment of the second objective of the study regarding the automation of soil texture classification enabled the operator to avoid manual reclassification of classes. Automation of the classification process naturally resulted in shorter processing time than by manual tools execution [75]. Lower human interference in the data processing also provided a uniform classification result regardless of the input data, as well as its universal applicability. Recent advances in machine learning algorithms development allowed their implementation and automation in soil texture mapping, also presenting a new possibility in future research [76]. The current limitation of this study is the full automation of the interpolation and classification algorithm, which could be performed by R or Python programming languages.

The applied soil texture classification algorithm allowed the delineation of soil texture suitability zones, which is often used as a criterion in various GIS-based multicriteria analysis studies. Soil texture was an impactful criterion in many natural resource management studies: grassland conservation [77], irrigation management [78], soil salinity analyses [79] and drought risk assessment [80]. Future directions of suitability analyses using soil texture zones are planned for vegetable suitability multicriteria analyses, which are highly beneficial in sustainable planning [81].

5. Conclusions

The automation of soil texture classification enabled time-efficient delineation of suitability zones for soybean cultivation. Clay, silt and sand soil contents were evaluated using a random split-sample method in five repetitions, which ensured the objective determination of training data normality and stationarity. The descriptive statistics of evaluated training data indicated moderate data normality and stationarity for clay and silt, while sand soil content resulted in highly skewed distribution and moderate stationarity. The application of logarithmic transformation as a pre-processing to OK reduced the effect of non-normal distribution, particularly for the sand soil fraction. Tested deterministic methods, IDW and SI, produced lower interpolation accuracy values than OK, but were considered as viable methods in the case of low data normality and stationarity. The authors recommend the evaluation of interpolation methods for suitability calculations, as their performance is variable and depends on the properties of input samples. It was also observed that these methods were effective in retaining the local variation data in the interpolation result. Loam and clay loam were determined as the most suitable soil texture zones for soybean cultivation, covering the northern part of the study area near the river Drava. Classified soil texture classes with the highest presence in the study area were moderately suitable silty clay loam (53.26%) and silt loam (34.96%). According to the proposed methods, the study area is dominantly moderately suitable for soybean cultivation, based on the soil texture data. The integration of Sentinel-2 multispectral satellite images and GIS-based multicriteria analysis with current results presents the foundation of future research.