Exploring the Effect of Sampling Density on Spatial Prediction with Spatial Interpolation of Multiple Soil Nutrients at a Regional Scale

Abstract

Essential soil nutrients are dynamic in nature and require timely management in farmers’ fields. Accurate prediction of the spatial distribution of soil nutrients using a suitable sampling density is a prerequisite for improving the practical utility of spatial soil fertility maps. However, practical research is required to address the challenge of selecting an optimal sampling density that is both cost-effective and accurate for preparing digital soil nutrient maps across regional extents. This study examines the impact of sampling density on spatial prediction accuracy for a range of soil fertility parameters over a regional extent of 8303 km² located in eastern India. Surface soil samples were collected from 1024 sample points. The performance of six levels of sampling densities for spatial prediction of 14 soil properties was compared using ordinary kriging. From the sample points, randomization was used to select 224 points for validation and the remaining 800 for calibration. Goodness-of-fit for the semi-variograms was evaluated by R² of model fit. Lin’s concordance correlation coefficient (CCC) and root mean square error (RMSE) were evaluated through independent validation as spatial prediction accuracy parameters. Results show that the impact of sampling density on prediction accuracy was unique for each soil property. As a common trend, R² of model fit and CCC scores improved, and RMSE values declined with the increasing sampling density for all soil properties. On the other hand, the rate of gain in the accuracy metrics with each increment in the sampling density gradually decreased and ultimately plateaued. This indicates that there exists a sampling density threshold beyond which the extra effort on additional sampling adds less to the spatial prediction accuracy. The findings of this study provide a valuable reference for optimizing soil nutrient mapping across regional extents.

1. Introduction

In the present era of precision agriculture, site-specific nutrient assessment is essential for correcting nutrient deficiencies and optimizing crop yield by adopting variable-rate fertilizer applications [1]. However, adopting site-specific soil management is much more challenging in small land-holding settings, where farmers regularly lack the resources to test their farm soils. For instance, Indian agriculture is dominated by small and marginal land holdings, with 86% of individual land sizes being less than 2 ha and an average of 1.08 ha [2]. Moreover, NPK fertilizer consumption in India has reached 30.6 million Metric Tons during 2023-24 [3]. Correspondingly, the subsidy bill on fertilizers in India has reached 1.9 trillion rupees (22.5 billion USD) during 2023-24 [4]. In these scenarios, various stakeholders, including farmers, researchers, and policymakers, ask for soil nutrient information in digital formats to make informed decisions about nutrient application and subsidy management [5].

Soil is heterogeneous at all scales due to dynamic interactions among natural and human-induced factors [6]. Measuring soil properties everywhere and all the time is impractical. Over the past decades, digital soil mapping (DSM) techniques have been extensively used to prepare spatial soil maps based on limited soil observations [7]. Among different mapping approaches, kriging [8] is a long-established geostatistical method for spatial interpolation, which produces unbiased estimates at unsampled locations with minimum error variances. Spatial interpolation adheres to the fundamental principle of geography, i.e., spatial autocorrelation encompassing ‘stationarity’ and ‘linearity’ [9,10]. It considers that all things are related, but near things are more related than farther apart [11]. Among different forms of kriging, ordinary kriging is regularly employed for soil mapping, especially when maps need to be prepared without any auxiliary information or covariates [12,13].

Because kriging works on the principle of spatial autocorrelation, sample locations must fully cover the map area. In addition, predictions perform best overall when the distances among sample points are minimized, providing sufficient support for the semi-variogram [14,15]. Ordinary kriging interpolation is especially sensitive to the number of sample points or sampling density [16,17]. A suitable sampling density is a critical prerequisite for accurate prediction of the spatial distribution of soil properties [18,19]. In general, more samples would be expected to capture the variation of soil properties more accurately [20]. However, operationally, the quantity of samples cannot be increased infinitely as scientists and policymakers must work with limited budgets [21,22,23]. Therefore, practical research is required to address the challenge of selecting an optimal sampling density that balances both economic efficiency and map accuracy [24,25,26].

Several studies have reported their observed threshold sample sizes at various mapping extents. While predicting soil organic carbon (SOC) at a small spatial extent of 0.24 km², Vašát et al. observed that accuracy increased up to an optimal sample size of 37 samples (156,118 samples per 1000 km²) [27]. With a mapping area nearly twice the size (0.44 km²), Wang et al. concluded that 40 samples (90,909 samples per 1000 km²) were optimal [28]. Although somewhat different from SOC concentration, Sherpa et al. found 160 samples (6897 samples per 1000 km²) to be optimal for mapping SOC stock for a 23 km² mapping area [29]. Zhang et al. while predicting SOC stock in a geographical extent of 72 km², observed that the prediction accuracy initially decreased and then increased with the number of samples up to 357 (4958 samples per 1000 km²) [30]. Similarly, Yu et al. observed the accuracy of SOC prediction in a 927 km² area improving with each increment of sample size to the maximum of 525 (566 samples per 1000 km²) [31]. Continuing this trend, Ye et al. found 669 samples (41 samples per 1000 km²) to be optimal for mapping SOC over a 16,400 km² area [32]. Finally, in a study with a much larger mapping extent (180,000 km²), Sun et al. concluded that a sample size of 800 samples (4.4 samples per 1000 km²) was optimal for mapping SOC [33]. Trends suggest that optimal sampling densities are situation-specific, and the ranges vary significantly across the mapping extents. However, most studies so far on the impact of sampling density on spatial prediction accuracy are around smaller mapping extents ranging from a few ha to less than 100 km². And, we have less certainty that these patterns continue to larger extents. Therefore, there is a need to unfold these results for regional extents to facilitate large-scale soil management for small land holdings.

The sampling requirement for one soil property might not be the same for another [27,34,35,36]. Information on the impact of sampling density on spatial interpolation accuracy has been reported only for a few soil properties, viz., SOC [27,28,31,32,33,37]; SOC stock [29,30], soil organic matter (SOM) [26,35,38,39,40], soil depth [41], clay [38], texture [16], granularity [35], pH [27,35], cation exchange capacity [26], available P, K [38,42], and available micronutrients [43]. This indicates that the literature so far is limited to only a few target soil properties, majorly SOC and SOM. Information on the impact of sampling density on spatial prediction accuracy for soil fertility parameters, especially soil nutrients, is very limited. These warrants experimenting on multiple soil properties, especially soil fertility parameters, which need to be managed through the application of external inputs.

In this context, this study was conducted to examine the effect of varying sampling densities on (a) goodness-of-fit for semi-variograms and (b) spatial prediction accuracy of ordinary kriging interpolation at a regional scale. A range of relatively static to dynamic soil properties were considered as target soil properties to identify an optimal sampling density that would suffice to generate reasonably accurate and economically efficient spatial soil maps across the parameters. The evaluation was performed in two stages. First, the goodness-of-fit for the resulting semi-variograms was measured quantitatively i.e., R² of semi-variogram model fit. Then, the prediction accuracies in terms of Lin’s concordance correlation coefficient (CCC) and root mean square error (RMSE) of the resulting maps were evaluated using a consistent independent validation data set. Briefly, this study is unique in terms of (1) its larger mapping extent, (2) consideration of multiple target soil properties, especially soil nutrients, (3) reporting the results in terms of sampling density instead of sample size, (4) exploring quantitative relationships between goodness-of-fit for semi-variograms and prediction accuracies of the resulting maps, and (5) performing a brief literature review of the previous works to identify a general trend of optimal sampling density with mapping extent.

2. Materials and Methods

A methodology was developed using evidence synthesis to support the development of a conceptual framework that would assess the relationship between sampling density and prediction accuracy of spatial soil maps for a suite of 14 soil properties at a regional extent. A brief description of the study area and conceptual framework are outlined below.

2.1. Study Area

The study area is Keonjhar district of Odisha state, located in eastern India, extending from 21.01° N to 22.16° N latitude and 85.18° E to 86.37° E longitude with a geographical extent of 8303 km² (Figure 1). The region experiences a hot and moist sub-humid tropical to sub-tropical climate with hot summers and cool winters. May is the warmest month of the year, with mean minimum and maximum temperatures of 25 °C and 46 °C, respectively. December is the coolest month of the year, with mean minimum and maximum temperatures of 12 °C and 25 °C, respectively. The mean annual rainfall is 1488 mm, more than half of which is received in June, July, and August during the movement of the SouthWest monsoon [44]. Common soil types of this region include red, yellow, alluvial, and lateritic with sandy loam to silty clay loam soil textures [45]. The parent material of the study area is predominately granite-gneiss. Soils belong to Alfisols, Inceptisols, and Entisols soil orders as per the USDA Soil Taxonomy [46].

2.2. Conceptual Framework

A conceptual framework was developed using literature review and field visits. A diagrammatic representation of the conceptual framework is presented in Figure 2. The components of this framework include soil sampling, laboratory analysis, sampling strategy for modeling, geostatistical modeling, soil mapping, and evaluation of spatial prediction accuracy. Individual components of this workflow are described below.

2.2.1. Soil Sampling

Topsoil samples (0.30 m depth) were collected randomly from 1024 sample points during the year 2019–2020 (Figure 1). Inaccessible areas were avoided for sampling. Local differentiations in soil colour and texture were given importance while choosing the sample points to capture the local variabilities. Approximately 1 kg of soil was collected from each point. The geographical coordinates (latitude and longitude) of the sample points were recorded using a GPS instrument (Garmin 76MAPCSx). The collected soil samples, after processing, were stored within properly levelled polyvinyl bottles for further laboratory analysis.

2.2.2. Laboratory Analysis

Soil samples were analyzed for 14 soil properties viz., SOC, pH, electrical conductivity (EC), available nitrogen (N), phosphorus (P), potassium (K), calcium (Ca), magnesium (Mg), sulfur (S), iron (Fe), manganese (Mn), copper (Cu), zinc (Zn), and boron (B). SOC was determined by the wet digestion method using the Ferroin indicator, and the carbon oxidized was multiplied by the efficiency factor of the experiment, i.e., 77% [47]. Soil pH at 1:2.5 w/v was determined using a glass electrode digital pH meter, and EC was estimated at 1:2.5 w/v by an electronic conductivity meter [48]. Available N was determined by the alkaline potassium permanganate method using an automatic nitrogen auto analyzer [49]. Available P was estimated by Bray’s No. 1 method using a spectrophotometer at 660 nm [50]. Available K was determined using the neutral normal ammonium acetate extraction method using a flame photometer [51]. Exchangeable Ca and Mg were determined by extraction with neutral normal ammonium acetate followed by titration against 0.01N EDTA (Ethylenediamine tetraacetic acid) using Eriochrome Black T and Calcon indicators [52]. Available S was determined using 0.15% CaCl₂.2H₂O extractant and BaCl₂ reagent followed by measuring the turbidity at 410 nm using a spectrophotometer [53]. Available Fe, Mn, Cu, and Zn were determined by the DTPA (diethylenetriaminepentaacetic acid) extraction method using an atomic absorption spectrophotometer [54]. Available B was estimated by the hot water extraction method using the Azomethine-H indicator [55].

Individual soil properties were grouped into three categories of soil properties, i.e., (i) physico-chemical soil properties including SOC, pH, EC, (ii) macronutrients including N, P, K, Ca, Mg, S, and (iii) micronutrients including Fe, Mn, Cu, Zn, and B. Descriptive statistics for all the soil properties at different sampling densities is presented in

Table 1. Descriptive statistics of the soil properties at different sampling densities.

Soil Property	Sampling Density (Number of Samples/1000 km²)	Minimum	Maximum	Mean	SD	CoV (%)	Skewness
SOC (%)	3	0.30	1.38	0.76	0.28	37.40	0.17
	6	0.10	1.38	0.72	0.29	40.82	0.26
	12	0.10	1.38	0.72	0.29	40.30	0.02
	24	0.09	1.38	0.72	0.29	39.83	0.03
	48	0.09	1.38	0.72	0.29	39.67	−0.03
	96	0.09	1.38	0.73	0.29	39.25	−0.09
	Validation dataset (N = 224)	0.10	1.32	0.75	0.29	37.92	−0.08
pH	3	4.15	6.90	5.58	0.86	15.48	0.15
	6	4.15	7.35	5.65	0.88	15.63	0.09
	12	4.15	7.35	5.70	0.88	15.47	0.12
	24	4.15	7.37	5.67	0.89	15.62	0.15
	48	4.15	7.37	5.63	0.89	15.76	0.23
	96	4.15	7.38	5.63	0.88	15.58	0.26
	Validation dataset (N = 224)	4.18	7.35	5.57	0.88	15.77	0.33
EC (dS/m)	3	0.10	0.50	0.30	0.11	35.80	−0.16
	6	0.03	0.50	0.27	0.12	43.86	−0.20
	12	0.03	0.50	0.27	0.11	40.71	−0.26
	24	0.03	0.50	0.26	0.10	39.85	−0.18
	48	0.03	0.50	0.27	0.11	39.82	−0.18
	96	0.03	0.50	0.27	0.11	39.56	−0.11
	Validation dataset (N = 224)	0.04	0.50	0.27	0.10	38.34	0.03
N (kg/ha)	3	95.28	531.42	295.49	93.12	31.52	0.28
	6	53.72	531.42	272.83	100.38	36.79	0.19
	12	50.99	531.42	278.69	92.97	33.36	−0.05
	24	50.99	531.42	275.62	90.93	32.99	0.07
	48	50.99	531.42	277.85	90.22	32.47	−0.08
	96	50.99	533.29	278.92	88.29	31.65	0.09
	Validation dataset (N = 224)	77.00	501.19	277.19	82.41	29.73	0.35
P (kg/ha)	3	8.53	44.48	22.64	9.38	41.41	0.31
	6	5.87	44.48	21.60	9.69	44.85	0.42
	12	5.87	44.48	21.99	8.62	39.20	0.22
	24	5.76	44.48	22.12	8.73	39.48	0.19
	48	5.76	44.48	21.93	8.87	40.44	0.20
	96	5.68	44.48	22.11	8.81	39.86	0.17
	Validation dataset (N = 224)	5.87	41.15	23.44	9.16	39.08	−0.04
K (kg/ha)	3	241.87	657.67	418.32	95.75	22.89	0.40
	6	173.22	724.68	428.84	124.42	29.01	0.51
	12	92.00	724.68	425.92	127.61	29.96	0.08
	24	76.90	728.24	425.91	130.02	30.53	0.01
	48	76.90	728.24	421.35	132.05	31.34	0.11
	96	76.90	728.24	419.37	130.17	31.04	0.14
	Validation dataset (N = 224)	98.24	722.22	411.48	125.06	30.39	−0.01
Ca (cmol(+)/kg)	3	0.49	5.24	2.90	1.48	50.95	0.05
	6	0.49	6.02	3.12	1.51	48.33	0.00
	12	0.49	6.09	3.23	1.41	43.72	−0.01
	24	0.49	6.39	3.27	1.39	42.67	−0.01
	48	0.49	6.60	3.24	1.40	43.38	0.12
	96	0.41	6.60	3.25	1.40	43.01	0.16
	Validation dataset (N = 224)	0.48	6.55	3.14	1.39	44.29	0.19
Mg (cmol(+)/kg)	3	0.13	3.23	1.71	0.91	53.07	0.05
	6	0.13	3.43	1.81	0.91	50.14	−0.04
	12	0.13	3.62	1.89	0.87	45.89	−0.01
	24	0.13	3.62	1.90	0.86	45.25	−0.07
	48	0.06	3.68	1.86	0.87	47.02	0.03
	96	0.06	3.68	1.88	0.87	46.13	0.03
	Validation dataset (N = 224)	0.07	3.68	1.82	0.87	47.72	0.07
S (mg/kg)	3	1.73	29.60	15.53	7.34	47.25	−0.05
	6	1.35	29.60	13.90	7.67	55.19	0.18
	12	1.35	29.60	13.96	6.96	49.86	0.12
	24	1.10	29.60	13.63	7.05	51.72	0.24
	48	1.06	29.60	13.93	6.90	49.57	0.07
	96	1.06	29.71	14.01	6.93	49.44	0.14
	Validation dataset (N = 224)	1.21	29.48	13.84	6.68	48.28	0.24
Fe (mg/kg)	3	35.55	180.60	81.59	35.09	43.01	0.85
	6	10.31	180.60	78.00	35.25	45.20	0.93
	12	10.31	180.60	76.09	31.60	41.53	0.56
	24	10.31	180.60	76.58	30.72	40.12	0.45
	48	10.05	180.60	76.44	32.95	43.10	0.37
	96	5.61	180.60	77.34	33.44	43.24	0.33
	Validation dataset (N = 224)	16.57	180.07	79.74	34.60	43.39	0.45
Mn (mg/kg)	3	17.49	85.99	44.92	17.45	38.85	0.42
	6	0.52	85.99	43.95	18.22	41.45	0.45
	12	0.52	85.99	43.46	17.69	40.71	0.07
	24	0.52	91.28	43.62	17.07	39.15	0.15
	48	0.52	91.28	42.65	16.75	39.28	0.17
	96	0.52	91.28	43.20	16.96	39.26	0.06
	Validation dataset (N = 224)	0.59	84.52	43.91	16.97	38.64	0.18
Cu (mg/kg)	3	1.29	5.51	3.21	1.16	36.10	0.09
	6	0.03	5.51	2.99	1.13	37.91	0.14
	12	0.03	5.86	2.90	1.23	42.35	0.29
	24	0.03	5.86	2.95	1.20	40.47	0.18
	48	0.03	5.86	3.01	1.19	39.42	0.05
	96	0.03	5.90	3.01	1.17	39.06	0.02
	Validation dataset (N = 224)	0.21	5.85	3.08	1.14	36.91	0.04
Zn (mg/kg)	3	0.30	2.32	1.00	0.52	51.97	0.57
	6	0.02	2.42	0.92	0.52	57.11	1.01
	12	0.02	2.42	0.89	0.52	58.65	1.04
	24	0.02	2.42	0.91	0.51	55.86	0.89
	48	0.02	2.49	0.93	0.50	54.45	0.74
	96	0.01	2.49	0.92	0.49	53.48	0.72
	Validation dataset (N = 224)	0.06	2.33	0.94	0.48	50.47	0.58
B (mg/kg)	3	0.25	1.63	0.78	0.40	51.44	0.57
	6	0.08	2.39	0.78	0.47	59.71	1.16
	12	0.07	2.43	0.77	0.46	60.08	1.13
	24	0.07	2.63	0.78	0.47	59.95	1.27
	48	0.06	2.63	0.83	0.51	61.11	1.15
	96	0.04	2.63	0.82	0.49	59.15	1.12
	Validation dataset (N = 224)	0.06	2.48	0.84	0.48	57.98	1.07

. The mean, standard deviation (SD), coefficient of variation (CoV), and skewness values of the soil properties did not differ remarkably across the sample subsets. Overall, CoV for different soil properties of the whole dataset varied between 16% (pH) and 59% (B). The CoVs for the physico-chemical soil properties, macronutrients, and micronutrients varied in the range of 16–44%, 31–49%, and 39–59%, respectively.

2.2.3. Sampling Strategy for Modeling

A soil dataset of 14 target soil properties for the 1024 sample points was generated by laboratory analysis. At first, 224 points were randomly selected out of the 1024 sample points as the independent validation dataset. The remaining 800 sample points were considered as the calibration dataset. The selection of samples was random using the =RAND () function in Microsoft Excel Version 2408, which returned evenly distributed random real numbers from the dataset. Using this function, new random real numbers were returned every time the worksheet was calculated. During each random selection of sample points, the function was repeated for 10 iterations to avoid the possible instability caused by a single random sample and to ensure good coverage of sample points for the entire study area. A similar procedure was also followed by Lai et al. [40]. The distribution of calibration and validation sample points has been presented in Figure 1.

Secondly, the calibration dataset was segregated into six levels of calibration sampling densities (Figure 2). The calibration sample sizes have been expressed in terms of sampling densities since sampling density also represents the mapping extent. The original calibration dataset of 800 samples (per 8303 km²) corresponds to a sampling density of 96 samples per 1000 km². This original sampling density was subjected to several consecutive random selections to generate six levels of calibration sampling densities. The original sampling density of 96 samples per 1000 km² (800 points) was considered as the first calibration dataset. For generating the second calibration subset, only 50% of the samples of the first calibration dataset were selected randomly to generate a sampling density of 48 samples per 1000 km² (400 points) using the =RAND () function in Microsoft Excel (Version 2408) with 10 iterations. This process of randomly selecting 50% of the sample points was carried out subsequently four more times to generate a total of six overlapping sampling densities viz., 96, 48, 24, 12, 6, and 3 samples per 1000 km² with corresponding sample sizes of 800, 400, 200, 100, 50, and 25, respectively.

2.2.4. Geostatistical Modeling and Soil Mapping

Geostatistical analyses were performed in the ‘gstat’ package of R version 4.1.0 [56]. Considering the values of a soil property z at two different points viz., x and x + h, which are separated at a distance (lag) of h from each other, the semi-variance (γ) between the pair of points was calculated following Equation (1). While building the empirical variogram, the default lag provided by the R software (version 4.1.0) was used.

$$\mathsf{\gamma }(\mathrm{h})={\displaystyle \frac{1}{2}}\{\mathrm{z}(\mathrm{x})-\mathrm{z}{(\mathrm{x}+\mathrm{h})\}}^2$$

(1)

where γ (h) is the semi-variance at lag h, h is the vector denoting the distance and direction between the two points x and x + h; z (x) and z (x + h) are the observed values of the soil property z at locations x, x + h, respectively. Here, x denotes the coordinates of the sample points in two dimensions.

Similarly, considering the soil property z being measured at numerous places in the region, and there are m pairs of observed values (considering all neighbouring pairs of observations), the mean semi-variance at any particular lag (h) was calculated following Equation (2) [14,57].

$$\mathsf{\gamma }(\mathrm{h})={\displaystyle \frac{1}{2\mathrm{m}}}{\sum }_{\mathrm{i}=1}^{\mathrm{m}}\{\mathrm{z}(\mathrm{x})-\mathrm{z}{(\mathrm{x}+\mathrm{h})\}}^2$$

(2)

The functions relating the γ and h were expressed as sample semi-variograms (also known as experimental semi-variograms). Different mathematical models viz., circular, spherical, exponential, gaussian, and bessel were tested on each sample semi-variogram, to fit continuous smooth curves, i.e., model semi-variograms (also known as theoretical semi-variograms) to the respective sample semi-variograms. The best-fitted mathematical models for each soil property and each sampling density were chosen based on whichever provided the lowest mean standard error. Semi-variogram parameters, namely nugget (C₀; uncorrelated variance), partial sill (C; spatially correlated variance), sill (C₀ + C), and range (r; limit of spatial correlation) were calculated. Since skewness values of all the soil properties were close to zero, the datasets were normally distributed, and thereby, data transformation was not performed [39,58]. Nevertheless, preserving the original values was also useful to avoid complications of data back-transformation while estimating error components at the validation stage [15].

Predictions for the unsampled locations were carried out using the ordinary kriging interpolation method (Equation (3)) [59,60]. Respective prediction maps were prepared using the ‘Krige ()’ function of the ‘gstat’ package of R (version 4.1.0) at a 100 m spatial resolution with a projected coordinate system of ‘WGS 1984 UTM Zone 45N’.

$$\widehat{\mathrm{z}} \left({\mathrm{X}}_0\right)={\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\mathrm{z}\left({\mathrm{X}}_{\mathrm{i}}\right)$$

(3)

where, $\widehat{\mathrm{z}}\left({\mathrm{X}}_0\right)$ is the predicted value of soil property ‘z’ at the unsampled location ${\mathrm{X}}_0$, z (${\mathrm{X}}_{\mathrm{i}}$) is the observed value of the same soil property ‘z’ at a neighbouring sampled location around the unsampled location, n is the number of neighbouring points, and ${\mathsf{\lambda }}_{\mathrm{i}}$ are the respective weights assigned to each measured value at the sampled locations.

2.2.5. Evaluation of Spatial Prediction Accuracy

Goodness-of-fit for the semi-variogram models was measured using the coefficient of determination (R²; Equation (4)). Subsequently, the prepared prediction maps were compared against the 224 independent validation points. Prediction accuracy of the maps were evaluated by determining CCC (Equation (5)) and RMSE (Equation (6)).

$${\mathrm{R}}^2=1-{\displaystyle \frac{{\mathrm{S}\mathrm{S}}_{\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}}{{\mathrm{S}\mathrm{S}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}}$$

(4)

$$\mathrm{C}\mathrm{C}\mathrm{C}={\displaystyle \frac{2\mathrm{r}{\mathsf{\sigma }}_{\mathrm{o}\mathrm{b}\mathrm{s}}{\mathsf{\sigma }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}}{{\left({\overline{\mathrm{Z}}}_{\mathrm{o}\mathrm{b}\mathrm{s}}-{\overline{\mathrm{Z}}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\right)}^2+{\mathsf{\sigma }}_{\mathrm{o}\mathrm{b}\mathrm{s}}^2+{\mathsf{\sigma }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^2}}$$

(5)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left[\mathrm{z}\left({\mathrm{S}}_{\mathrm{i}}\right)-\widehat{\mathrm{z}}\left({\mathrm{S}}_{\mathrm{i}}\right)\right]}^2}$$

(6)

where SS_error is the sum of the squares of the errors, SS_total is the total sum of the squares, $\mathrm{r}$ is the Pearson’s correlation coefficient, $\overline{\mathrm{Z}}$ is the mean of the variable, ${\mathsf{\sigma }}^2$ is the variance, $\mathrm{z}\left({\mathrm{S}}_{\mathrm{i}}\right)$ is the observed value of variable z at the location ${\mathrm{S}}_{\mathrm{i}}$, $\widehat{\mathrm{z}}\left({\mathrm{S}}_{\mathrm{i}}\right)$ is the predicted value of variable z at the location $\left({\mathrm{S}}_{\mathrm{i}}\right)$, and n is the number of sample points.

3. Results

In Section 3.1., goodness-of-fits for semi-variograms in terms of R² of model fit are presented. In Section 3.2., the prediction accuracies of the spatial soil maps in terms of independently validated CCC and RMSE are presented. In each section, trends are presented in three subsections viz., (a) general trends across all soil properties, (b) variations across categories of soil properties, and (c) variations across individual soil properties. Trends for RMSE values are presented across individual soil properties only since the units of RMSEs are specific to respective soil properties. Variations in the visual representation of the prediction maps are presented in Section 3.3.

3.1. Goodness-of-Fit for Semi-Variograms (R² of Model Fit)

3.1.1. General Trends across All Soil Properties

The sample semi-variograms obtained with smaller sampling densities viz., 3, 6, or 12 samples per 1000 km², were observed to be noisier and more erratic compared to those obtained with higher sampling densities, viz., 24, 48, and 96 samples per 1000 km² (Figure 3). Correspondingly, the goodness-of-fit for the model semi-variograms to the sample semi-variograms was also poorer for the smaller sampling densities than those obtained with the higher sampling densities. Semi-variograms and semi-variogram parameters, including the type of semi-variogram model, nugget, sill, and range for all the soil properties at different sampling densities, have been presented in the Supplementary Materials.

Irrespective of the soil property, the R² of the semi-variogram model fit initially increased with the increasing sampling density and then reached a plateau (Figure 4). The range for R² of model fit across soil properties gradually became narrower with the increasing sampling density, indicating that as the sampling density increased, the reliability of goodness-of-fit across all soil properties improved. In other words, the probability of obtaining a higher R² of model fit, regardless of soil property, increased with increased sampling densities. The rate of gain in the mean R² of model fit decreased gradually from 20% for increasing the sampling density from 3 to 6 samples per 1000 km² to 0% for increasing the sampling density from 48 to 96 samples per 1000 km².

3.1.2. Variations across Categories of Soil Properties

All three categories of soil properties, viz., physico-chemical soil properties, macronutrients, and micronutrients, followed similar trends as the general trend for R² of model fit (Figure 5). For instance, for all three categories of soil properties, the minimum, maximum, and mean values of R² of model fit initially showed increasing trends, which later increased relatively slowly and finally remained constant with the increasing sampling density. For all three categories of soil properties, sampling densities of 48 and 96 samples per 1000 km² were able to generate smooth semi-variograms with R² of model fits varying between 0.92 and 0.99. However, at any sampling density, the mean R² values for model fit per category consistently decreased from physico-chemical soil properties (R² = 0.92) > macronutrients (R² = 0.90) > micronutrients (R² = 0.79).

The rate of gain in the mean R² of model fit values also followed similar trends as the general trend. Irrespective of the soil property category, the rate of gain in the mean R² of model fit reduced gradually and ultimately became 0% with the increasing sampling density. Nevertheless, the goodness-of-fit for micronutrients category responded more to the rise in sampling density compared to other categories of soil properties.

3.1.3. Variations across Individual Soil Properties

All soil properties showed increasing trends for R² of model fit with the increasing sampling density (Figure 6). For each soil property, the lowest R² of model fit was observed at the sampling density of 3 samples per 1000 km². However, the highest R² of model fit was observed at varying sampling densities for each soil property. For instance, SOC, pH, K, Ca, Mg, Mn, Cu, and Zn reached their highest R² of model fits at the sampling density of 48 samples per 1000 km², whereas for P, Fe, and B, the highest R² of model fits were obtained at the sampling density of 96 samples per 1000 km². For a few soil properties viz., EC and N, the sampling density of 24 samples per 1000 km² was enough for generating the highest R² of model fit that could be obtained with even higher sampling densities. For the soil property S, even the sampling density of 12 samples per 1000 km² was enough to obtain the highest R² of model fit, which stayed constant thereafter, even with further increments in the sampling densities. Nevertheless, each soil property except P could achieve a R² of at least 0.95 or greater, usually at the higher side of sampling densities.

As an exception, the R² of the model fits decreased in a few instances of increasing sampling density. For example, the R² of model fit for S and Mn decreased by 1–3.2% while increasing the sampling densities from 48 to 96 samples per 1000 km². Also, for K, the R² of the model fit decreased by 6.3% while increasing the sampling density from 6 to 12 samples per 1000 km². In terms of patterns for R² of model fit with increasing sampling density, SOC, pH, EC, N, Ca, Mg, and S had similar trends with the maximum rates of gain in R² of model fits while increasing the sampling density from 3 to 6 samples per 1000 km², which thereafter showed an almost constant level with further increments in the sampling density. Soil property K also possessed a similar trend as above, except that there was also a decline in the R² of model fit while increasing the sampling density from 6 to 12 per 1000 km². Soil properties Cu and Zn possessed trends that showed exceptional improvements in the R² of the model fit while increasing the sampling densities from 6 to 12 samples per 1000 km². The trends of R² of model fits were distinguishably different for P, Fe, Mn, and B, which had their highest rates of gains in R² of model fit while increasing the sampling density from 12 to 24 samples per 1000 km².

3.2. Prediction Accuracy (CCC and RMSE)

3.2.1. General Trends across All Soil Properties

The mean prediction accuracy, as measured by CCC, increased with increasing sampling density. The minimum and maximum CCCs obtained across soil properties increased with increasing sampling density. The rate of accuracy improvement with increasing sample density was initially higher but gradually slowed down and finally remained almost constant (Figure 4). The range of prediction accuracy across soil properties was widest at the sampling density of 6 samples per 1000 km². It was narrowest at the sampling density of 48 samples per 1000 km². Considering the range of CCC as an indicator of the reliability of prediction across soil properties, the sampling densities of 6 and 48 samples per 1000 km² were observed to be the most unreliable and most reliable, respectively.

The rate of gain in the mean CCC while increasing the sampling density from 3 to 6 samples per 1000 km² was observed to be 20%. However, this rate of gain decreased for each increment in the sampling density, and finally, the rate of gain in the mean CCC was only 1%, while the sampling density was increased from 48 to 96 samples per 1000 km². Thus, even though the mean CCC gradually increased with increasing sampling density, there was a simultaneous decline in the rate of gain in the mean CCC with each sequential increase in the sampling density. Moreover, the R² of model fits, and CCCs were positively correlated (r = 0.80).

3.2.2. Variations across Categories of Soil Properties

Similar to the general trend, all three categories of soil properties showed improvements in the independent validation CCC with increasing sampling density (Figure 5). The mean CCCs at different sampling densities for different categories of soil properties indicated that mean CCCs initially increased rapidly, then slowly, and finally remained almost constant with increasing sampling density. In other words, the rate of improvement in the mean CCCs ultimately plateaued with the increasing sampling density for all three categories of soil properties.

The level of prediction accuracy decreased in the order of physico-chemical properties > macronutrients > micronutrients with mean CCCs of 0.74, 0.67, and 0.53 for physico-chemical soil properties, macronutrients, and micronutrients, respectively. Standard deviation values of CCCs for physico-chemical soil properties, macronutrients, and micronutrients were 0.09, 0.15, and 0.13, respectively, indicating that CCCs for the macronutrients and micronutrients were more variable across their constituent soil properties compared to that for the physico-chemical soil properties. However, the variability of independent validation results across micronutrients was much higher at lower sampling densities than at higher sampling densities.

For all three categories of soil properties, CCCs were observed to be positively correlated with the R² of model fits for semi-variograms. The respective correlation coefficients for the relationship between semi-variogram model fit and prediction performance for physico-chemical soil properties, macronutrients, and micronutrients were 0.89, 0.77, and 0.84, respectively. Although prediction accuracies were successively lower for macro and micronutrients, the generalization of sampling density being connected to the semi-variogram quality and prediction accuracy—with diminishing returns—remained true also for these soil property categories.

3.2.3. Variations across Individual Soil Properties

The CCCs of each of the soil properties increased with increasing sampling density (Figure 6). Each soil property possessed its lowest CCC at the sampling density of 3 samples per 1000 km². For SOC, N, K, Mg, S, Zn, and B, the highest CCCs were obtained at the highest sampling density of 96 samples per 1000 km², whereas the remaining soil properties reached their highest CCCs with 48 samples per 1000 km².

The low rate of improvement for CCC while increasing the sampling density from 48 to 96 samples per 1000 km² is noteworthy because of the relatively high effort required to produce that increase in sample density. This increase in sample density for the present study represents efforts of sampling an additional 400 sample points along with the associated costs for laboratory analysis. For SOC, N, K, Mg, S, Zn, and B, the rates of gain in CCC varied in the range from 1.3% to 7.3% while increasing the sampling density from 48 to 96 samples per 1000 km². However, in the same situation, pH, EC, Ca, and Cu did not show any improvement in CCC. In contrast, for P, Fe, and Mn, the CCCs decreased in the range of 1.6–5.8% while increasing the sampling density from 48 to 96 samples per 1000 km².

In terms of patterns of prediction accuracy, soil properties viz., EC, N, K, Ca, and Mg possessed similar trends having the maximum rates of gains in CCCs while the sampling density was increased from 3 to 6 samples per 1000 km², which thereafter decreased with the further increments in the sampling density. Boron possessed a similar trend, except that the rate of gain in CCC became negative (−1.6%) while increasing the sampling density from 12 to 24 samples per 1000 km², which again increased and later decreased with further increments in the sampling densities. Soil properties SOC, pH, and S possessed considerably higher rates of gains in CCCs at the initial two stages of increments in the sampling density, i.e., while increasing the sampling density from 3 to 6 and from 6 to 12 samples per 1000 km². Zinc possessed the highest rate of gain in CCC while increasing the sampling density from 6 to 12 per 1000 km², which further decreased, increased, and ultimately decreased with the further increments in the sampling density. Micronutrients Fe, Mn, and Cu possessed similar patterns of CCCs having higher rates of gains in CCCs at the second and third increments in sampling density viz., while increasing the sampling density from 6 to 12 and 12 to 24 samples per 1000 km², which then decreased thereafter with further increments in the sampling density. The trend of P suggests that the rate of gain in CCC was considerably higher while increasing the sampling densities from 12 to 24 and 24 to 48 samples per 1000 km², which then decreased with further increments in the sampling density.

In terms of levels of prediction accuracy achieved, maximum CCCs greater than 0.80 were obtained for soil properties SOC, pH, N, Mg, and S, whereas Ca and EC had maximum CCCs of 0.79 and 0.76, respectively. On the other hand, the maximum CCCs obtained for P, K, Fe, Mn, Cu, Zn, and B were observed to be varying between 0.52 and 0.69, which are poorer CCCs than that for other soil properties. In the context of the soil property categories, these results indicate that EC was the more difficult property to predict in the category of physico-chemical properties, while P was the more difficult property to predict in the macronutrient category. The high prediction accuracy for N was intriguing, considering how dynamic this soil property tends to be, both spatially and temporally.

Independent validation RMSEs of individual soil properties have been presented in Figure 7. All soil properties showed declining trends of RMSEs with increasing sampling density. Soil properties N, K, Ca, Zn, and B had sharp declines in RMSEs between 3 and 6 samples per 1000 km². SOC, pH, Mg, S, and Cu, in contrast, had gradual declining trends of RMSEs at each increment of sampling density. For EC and P, the RMSEs initially decreased with the increments in sampling density up to 12 samples per 1000 km² and thereafter remained constant with further increases in the sampling density. Although still following the general trend, RMSEs for Fe and Mn were more erratic.

3.3. Prediction Maps

Soil properties SOC, EC, N, P, S, Fe, Mn, Cu, Zn, and B were observed to have the highest concentrations in the North-West part of the study area, whereas the remaining soil properties were observed to have higher concentrations towards the South-East part of the study area (Figure 8a,b). These patterns suggest that many of the soil properties covary with each other, while the other soil properties (i.e., pH, K, Ca, Mg) appear to have an inverse relationship with those properties.

For any soil property, the prediction maps generated similar patterns of spatial distribution at all the sampling densities. However, the prediction maps obtained from lower sampling densities were comparatively smoother across the study area compared to the ones produced with larger sampling densities. The map smoothness could be caused by either a lower proportion of the map being outside the range of spatial autocorrelation identified by the respective semi-variograms or the semi-variograms identifying longer ranges of spatial autocorrelation with increased sampling densities. An intriguing pattern in the semi-variogram ranges with different sampling densities was that the shortest ranges were usually around 12 to 24 samples per 1000 km², with lower and higher sampling densities tending to have semi-variograms with longer ranges (Supplementary Materials). This pattern suggests that a consistent increase in map detail with increasing sampling density is more related to the proportion of the map area observed rather than differences in the semi-variograms generated at different sampling densities. In other words, once a plateau for improving the quality of the semi-variogram is reached, additional samples distributed across the mapping area can still improve the level of detail represented in the map.

4. Discussion

4.1. Effect of Sampling Density on Goodness-of-Fit for Semi-Variograms

For an effective kriging interpolation, the goodness-of-fit between the sample (experimental) and model (theoretical) semi-variograms should be satisfactory [12,61]. A reliable semi-variogram requires sufficient sample points for measuring and describing the characteristics of spatial autocorrelation [62]. Previous recommendations for sampling in support of spatial interpolation by kriging have tended to focus on minimum sample quantities. Past studies have recommended minimum quantities of samples ranging from 50 to 225 [13,29,63,64,65]. The improved model fits in this experimental dataset at 200 sample points, suggests a minimum quantity of points towards the higher end of this range for generating reliable semi-variograms at the regional scale.

Recommendations for minimum sample size have sometimes been linked to methods of semi-variogram estimation but likely also correspond with the sample set’s ability to capture the pattern of spatial autocorrelation in the soil landscape of interest. This spatial pattern could be produced by physical processes in the landscape or moving average processes [66]. In any case, the chances of capturing that full pattern with the sample set decrease as sample locations become further apart. As an extreme example, a sample set of 100 data points distributed over a 100 km² study area would not be expected to capture the same patterns as they would if they were distributed over a 100,000 km² area. For this reason, each of the sample sets in the present study has been expressed in terms of sampling density. However, this concept does not exclude the possibility that there is also a minimum quantity of samples needed to construct a reliable semi-variogram, regardless of the study area size.

In this experiment, the sample semi-variograms obtained with sampling densities below 24 samples per 1000 km² (sample size below 200 points) were observed to be noisy and erratic, whereas those obtained at the higher sampling densities were smooth and stable (Figure 3). Reflecting this pattern, the goodness-of-fit between the sample and model semi-variograms (R² of model fit) was also observed to improve with increasing sampling density (Figure 4). Thus, higher sampling densities or higher sample sizes favoured generating semi-variograms with less noise and improved the goodness-of-fit. As an example of an analysis to identify sufficient sampling for supporting a reliable semi-variogram, Oliver and Webster reanalyzed exchangeable potassium data from the Broom’s Barn farm [57,67]. They observed that the sample semi-variogram obtained from a smaller sample size (87) was much more erratic and noisier compared to the one obtained from a larger sample size (434). Accordingly, the goodness-of-fit for semi-variograms was also poorer in the case where a smaller sample size was used. Ye et al. while predicting SOC over a spatial extent of 16,400 km², also observed that the goodness-of-fit for semi-variograms improved with an increasing sampling density [32]. The reason behind the poor semi-variogram fittings at lower sampling densities could be that too few sample points might not sufficiently capture the complex spatial structures of soil properties. In such situations, semi-variogram estimates generated from very few sample points are usually unreliable [57,68].

4.2. Effect of Sampling Density on Spatial Prediction Accuracy

Trends of increasing CCCs and decreasing RMSEs for independent validation with the increasing sampling density indicate that irrespective of the soil property, prediction accuracy improved with sampling density. These results are in line with many other studies, both at local [28,41] and regional scales [16,31,33,35]. Improvements in prediction accuracy with the increase in sampling density can be attributed, at least in part, to the increased amount of representativeness of the spatial variation of the soil properties. In other words, more of the map area had been directly observed. The mean distance between samples decreased with the increased sampling density, which resulted in a reduced overall standard error of kriging interpolation. Additionally, increased sample density could have better captured the spatial autocorrelation structures occurring at smaller phenomenon scales [39]. Decreased distances between sample locations associated with the increased sampling density might have also added useful spatial information, leading to improved prediction accuracies at the higher sampling densities.

The decrease in the rate of gain in prediction accuracy with the increasing sampling density could be because the sample points became increasingly spatially clustered, causing the information to reach a saturation level [39]. This is in line with McBratney and Webster’s observation that spatial autocorrelation reduces the quantity of samples needed to achieve lower standard errors [69]. The relatively consistent trend observed in the present study of increasing prediction accuracy with increasing sample density up to the point of diminishing improvement largely agrees with trends reported in the literature. Nonetheless, some deviations from that trend have been observed. For example, similar to the results in the present study, Long et al. observed that the rate of improvement in the prediction accuracy of SOM was gradual in the first stage of increasing sampling density, which was then followed by the rate of improvement being stagnated in the second stage [39]. However, in the third stage of increasing sampling density, the rate of improvement in prediction accuracy increased again and at a faster rate. In their case, the mapping extent was sufficiently larger (geographical area of 124,000 km²), and the calibration sample sizes were also sufficiently larger (varied between 9365 and 188,247 soil samples). The findings of Li are also interestingly different as the validation RMSE of ordinary kriging interpolation initially decreased and then again started increasing linearly with the increasing sampling density while predicting SOM content over a geographical area of 400 km² in China [70]. Recently, Sun et al. observed that even though prediction accuracy in general improved with the sampling density, the increasing trend was abrupted at certain levels of increasing sampling density [33]. These exceptions suggest that plateaus of prediction accuracy may occur in certain situations.

4.3. Effect of Sampling Density on Spatial Patterns Observed in the Prediction Maps

For any soil property, the prediction maps prepared from different sampling densities followed similar patterns of spatial distribution. Long et al. also observed that the spatially interpolated prediction maps of SOM obtained at 20 different sampling densities were visually similar [39]. Similarly, Domenech et al. observed that the prediction maps for soil depth obtained with six different sample sizes possessed similar spatial patterns [41].

Although the prediction maps across the sampling densities possessed similar spatial patterns, the ones obtained from smaller sampling densities were smoother. The more generalized maps could roughly represent the spatial patterns but lacked details about the local variations. Another problem with the smoothed maps is that they might increase the low values and reduce the high values of soil properties [71]. Even though the smoother maps might be useful for certain broad-level applications, they may be less suitable for the more precise management recommendations for soil properties at finer scales [72]. In contrast, the maps obtained with higher sampling densities were more fragmented, which could represent more details about the local variations. Similarly, Sun et al. observed that spatial heterogeneity increased with increasing sampling density, which could capture the spatial variations of soil properties caused by anthropogenic activities at smaller distances [37]. Thus, even though the lower sampling densities can produce the overall pattern of spatial distributions, the local variations are less likely to be represented [32]. Overall, results from these and the present study agree with Loiseau et al., who observed that with an increasing sampling density, the stabilization of spatial prediction increased, the smoothing effect decreased, and spatial structures were described in finer detail with the kriged prediction maps [16].

In contrast, John et al. observed that SOM maps prepared with ordinary kriging of 150 samples differed from those prepared with 250, 450, 600, and 1200 sample points [17]. However, those differences were related to semi-variogram ranges that were relatively short to sample spacing, leading to different representations of ‘hot spots’ surrounding point locations. Lai et al. also observed that the spatial distribution pattern of prediction maps derived from ordinary kriging interpolation with lower sampling densities differed from those obtained with larger sampling densities [40]. In this case, the map of SOM became smoother with the use of additional samples. The patterns in the resulting map suggest that the additional samples modified the modeling of the semi-variogram such that the range was likely increasing with the additional samples, which converted the local ‘hot spots’ around a few points to more regional groupings with the added samples.

Irrespective of the similarity of the overall spatial pattern of the prediction maps across all the sampling densities, the spatial patterns obtained with higher sampling densities (viz., 24, 48, and 96 samples per 1000 km²) were observed to be visually more identical to each other. Visual similarities among prediction maps prepared with different sampling densities could be related to the reuse of sample locations across the sampling densities in the present study’s experimental design [40]. However, the experimental sampling densities were increased exponentially to keep the proportion of new sample locations at 50% for each increment. Visual comparison of the whole study area will tend to emphasize the major trend. The increase in CCC is likely a better indicator of how well the map represents more local variations.

4.4. Which Sampling Density Should Be Considered for Optimal Interpolation of Soil Properties?

The utility of soil maps and the requirements of the stakeholders are the ultimate driving forces in deciding the accuracy requirements. For example, if the stakeholders are policy makers or decision makers, the coarse resolution maps prepared from a small sampling density with a relatively lower accuracy might be sufficient, but not for farmers for their field management. While lower levels of accuracy may be acceptable for lowering costs, this study focused on the inflection point for the rate of gain in improved model performance. At the sampling density of 3 samples per 1000 km², all the soil properties possessed their lowest R² of model fits and CCCs for independent validation. By these metrics, model quality improved—with diminishing returns—for increasing sample density until the improvement was negligible, moving from 48 to 96 samples per 1000 km². Across the soil properties examined in this study, a sampling density of 48 per 1000 km² was identified as that inflection point where there was no substantial improvement in model performance for the next increment of increased sampling density. A possible exception to that was the lack of a plateau for independent validation of the physico-chemical category of soil properties. Although the goodness-of-fit for semi-variogram models did not change, the prediction accuracy of the resulting maps may have still benefited from additional increases in sampling density.

A brief literature review was conducted to compare the results of this study to those of previous studies. For this purpose, the results of the previous studies, either reported an optimal sampling density or sample size, were converted into sampling density for quantitatively identifying the relationship between sampling density and mapping extent. When considered in the context of previous research along with the size of the mapping extents, the results of this study fit a pattern of optimal sampling density decreasing as a function of increasing mapping extent (Figure 9). A decrease in optimal sampling density for larger mapping extents suggests potential explanations. One possibility is that more samples are providing less additional information for spatial modeling due to spatial autocorrelation. McBratney and Webster demonstrated how spatial autocorrelation reduces the number of samples needed to achieve a desired level of standard error, which may translate to lower densities of samples needed to produce maps with optimized prediction performance with kriging [69]. Another possibility is that the statistics focus on larger-scale trends with larger mapping extents, with more local variation being treated as noise. However, the maps produced in this study with lower sampling densities were the most generalized and became more spatially detailed up to optimal sample density. The degree to which spatial detail was added to the resulting maps corresponded to the decrease in the rates of gain for semi-variogram goodness-of-fit and independent validation CCC.

In the literature review of optimal sampling density for kriging, a few studies did not fit the trend as well as the others (Figure 9). The one study that mapped soil depth was an outlier with the identification of around 323 samples per km² for fields ranging from 0.2 to 0.7 square km [41], which was a much higher sample density than would be expected for those mapping extents. Nanni et al. identified different optimal sampling densities for mapping clay content, SOM, P, K, and base saturation for an area of 1.84 square km [38]. The optimal sampling densities for clay content and SOM aligned with the observed trend with 35 and 48 samples per km², respectively. However, the authors concluded that a sampling density greater than 100 samples per km² was needed to map P, K, and base saturation in the area of highly managed sugar cane.

The optimal sampling densities identified in the present study for macro and micronutrients were consistent with the general trends, but the area being mapped has likely been less intensely managed than the fields in [38]. This difference suggests that management may induce an additional level of spatial variability that requires a higher sampling density. Moreover, landform complexity and land use patterns are also supposed to affect the sampling density requirement and prediction accuracy of soil mapping. For example, the findings of Sun et al. demonstrate that for the same set of sample quantities, accuracy metrics were better in plain than in the hill mountainous landform [33]. Similarly, Safaee et al. observed that prediction performance was more consistent and accurate in uniformly managed lands than in intensively managed ones [26]. Therefore, the sampling density-prediction accuracy relationships could be rather more complex and unique (site-specific) to an area depending on the local factors such as existing spatial autocorrelations, climate, topography, landform complexity, and management practices. Thus, future studies may include these additional local factors in addition to the experience gained in this study to explore the impact of sapling density on spatial prediction accuracy further.

5. Conclusions

This study highlights the critical role of sampling density in generating reliable predictions of soil properties through semi-variogram models. Sampling densities below 24 samples per 1000 km² produced noisy and unstable results, but as densities increased, the semi-variograms stabilized, and model fit improved, with higher R² values and better prediction accuracy (increased CCCs and decreased RMSEs). However, the rate of improvement diminished as sampling density increased, indicating diminishing returns beyond a certain point. For instance, increasing sampling density from 48 to 96 samples per 1000 km² had minimal effect on the model’s goodness-of-fit and prediction accuracy, suggesting a saturation threshold.

A strong correlation was observed between model fit (R²) and prediction performance (CCC), especially for physico-chemical properties, which showed higher accuracy compared to macronutrients and micronutrients. The reliability of prediction accuracy or probability of obtaining a higher R² of model fit or CCC, regardless of soil property, increased with increased sampling densities. Prediction accuracy across all the soil properties was most reliable at the sampling density of 48 samples per 1000 km². The mean CCC values decreased in the order of physico-chemical properties (0.74) > macronutrients (0.67) > micronutrients (0.53). This indicates that different soil properties respond differently to increases in sampling density, with physico-chemical properties being the most predictable. Nevertheless, considering results in the context of previous work on mapping with kriging, this study identified a general trend of optimal sampling density decreasing with exponential increases in mapping extent. The study emphasizes that regional soil maps, despite lower prediction accuracy for nutrients, provide valuable spatial information that can guide policymakers in soil management. These maps can inform planners where nutrient management is less intensive, offering a practical tool for broad-scale agricultural recommendations when site-specific data are unavailable or too costly to obtain.