Comparison of Spatial Interpolation Methods for Mapping the Qualitative Properties of Soil
Mohammadali Nikpey ^{1*}, Mahdi Sedighkia ^{2}, Mohammad Bagher Nateghi ^{1}, Javad Robatjazi ^{3}
Comparison of Spatial Interpolation Methods for Mapping the Qualitative Properties of Soil
Mohammadali Nikpey ^{1*}, Mahdi Sedighkia ^{2}, Mohammad Bagher Nateghi ^{1}, Javad Robatjazi ^{3}
ABSTRACT
The high proliferation rate of population and their over growing demands for foods that rely on environment and natural resources resulted in evaluation and recognition of more inputs, in particular soil and water resources. In this regard, mapping the soil properties is considered as welladopted approach to provide fundamental information related to land resources. As a branch of applied statistics, geostatistics uses the collected information sampled location to provide a broad range of estimates concerning the land properties in unsampled locations. In the present study, we evaluated the performance of estimates kriging, inverse distance weighting, and Cokriging to map some of the soil quality properties in chat fields of Golestan province located in northern part of Iran. We measured percentage of clay, silt, sand, calcium carbonate, organic carbon, and concentrations of micronutrients such as iron, copper, zinc, manganese, and concentrations of macronutrients including nitrogen, potassium, and phosphorus as main parameters affecting the soil quality. Sampling was performed on grids 250 × 250 meters in 35 points of the depth of 20 cm. After data normalization, experimental semivariograms was drawn. For Kriging and Cokriging estimation we utilized spherical, exponential, circular, and Gaussian models, and to estimate variables using Inverse Distance Weighting, the parameters of 15 power was used. Kriging method was found to the two other methods for estimation of soil properties under study. For the estimation of electrical conductivity and soil’ iron content, the Kriging and Cokriging showed the best estimation of results. However, to estimate the other parameters of soil (sand, organic carbon, pH and lime) inverse distance weighting method with power one was found to have the least error leading to the best interpolation. Finally, considering the best method of interpolation, zoning maps of soil quality properties was created in the GIS.
Keywords: Geostatistics, Geographic Information System, Soil Quality
cultivated land is sensed in large part due to rapid population growth along with development of industrial and residential zones, limitation of land’s fertility, and a serious global crisis (Ayobi et al., 2006; Khaledian et al., 2016). A sustainable agriculture can be achieved upon classifying the existing lands resources based on their suitability for various agricultural purposes in addition to establishing strategies for better exploitation of resources in long term period (Sys et al., 1993). In this regard, soil map is among the most fundamental land resources (Webster and Oliver 2007). In the past, soil mapping and spatial heterogeneity were often accomplished by the classical statistical methods with the assumption that any changes in the properties of the soil map units was considered to be random, while nowadays it is clear that the accuracy of geostatistical methods to estimate the spatial distribution of the measured data is significantly higher than the classical methods due to consideration of spatial continuity (Uyan, 2016).
The spatial statistics based on the theory of regional variable socalled “Geostatistics” aim to model the regional variable within the framework of the theory of probability (Uyan and Cay, 2013). Using variogram as a a tool for quantitative spatial variability of phenomena, the geostatistics have been able to evaluate the spatial continuity of the observations, and upon spatial correlation of the observation, the unpredictable variables at the sampling points is turned into predictable values to provide estimates maps (Acerbi Júnior et al., 2015). In the other words, the geostatistical estimation consists of two steps: (i) understanding and modeling the spatial structure of the variable that can be investigated by variogram analysis, and (ii) the estimation of desirable variable using geostistical approaches, such as Kriging, and Cokriging (Davis, 1987). The use of Geostatistical and GIS techniques have captured significant attention in the management and decisionmaking of agricultural research. For instance; zoning crop pests (Bouwnmeester et al., 2016; Pizzato, 2014) prediction of crop phenology (Gerstmann et al., 2016), evaluation of physical, chemical and biological properties of soil (Bitencourt et al., 2016; Jeelani et al., 2017; Wang and Shao, 2013) and so on. In light of these efforts, the use of theoretical methods such as GIS and geostatistics has potential to significantly improve the efficiency of traditional cropping systems through providing high accuracy estimation of the soil properties of the target area (da Silva et al., 2015). In this work, we have adopted GIS method based on kriging, Cokriging and inverse distance weighting interpolation schemes to study and evaluate the soil quality properties of Chat land in Gonbad province located in northern part of Iran.
Materials and Methods
Study area
The study was conducted in farmlands area of Gonbad province located in northern part of Iran with a total area over 207 hectares and within longitude of 55° 57´ E, latitude 37° 54´ N (Figure 1). Alluvial plains are the main types of topography in that area. It is characterized by dry semiarid climate according to Domarten and coupons dry classification method and based on Ambrose, with a mean annual temperature of 18.6◦C, a mean annual relative humidity of 65%, and a mean annual precipitation of 500 mm represent the ranges of means for different portions of the study area. According to America soil classification, the main soil of that area types are Aridisol (Figure 1).
Sampling and physicalchemical analysis of soil
The area under study was gridded (250 m × 250 m) and with regard to soil and land units separated on map, a total of 35 samples (025 cm) were collected. Laboratory analysis were conducted on the collected samples after air drying and sieving to pass a 2mm mesh (F. Dai et al., Ecological Indicator 2014). (Black et al., 1965). Organic carbon (Walkley and Black) (Jou, 1978), soil pattern (Hydrometer method) (Bouyoucos, 1962), carbonate calcium (Calcimeter) and soil acidity were analyzed using pH meters after diluting the sample with 1:1 ratio of double distilled to soil. Iron, zinc, manganese, copper (DTPA extraction and atomic absorption) and total nitrogen content were quantified using Microdigital technique. We measured P through Olsen method, potassium content through ammonium acetate method and Electrical conductivity with an electrical conductivity meter (Sparks and Bartels, 1996).
Geostatistical Analysis
Variogram was used for spatial analysis to evaluate the reliability of assumptions of truth. Variogram is utilized to assess and evaluate the structural features of the regional variable and to explain its variation. If variogram reaches to a certain limit that results in certain effective range, the spatial structure and conditions for inherent assumption can be achieved. Given that the calculation of the semivariance is not possible for entire society under study, one can use a function (1) in order to estimate the semivariance in a specific separation distance (Oliver and Webster 2014);
Where, Y_{i} (h) is the experimental semivariogram; and Z(x), and Z (x_{i} + h) are the observed value of the variable in x_{i} and x_{i }+ h points, respectively.
We used the C_{0} / C + C_{0} ratio to obtain the power of spatial structural of the variable in which C_{0} is the nugget and C_{0}+ C is the threshold of semivariogram (sill). The C_{0} / C + C_{0} ratio, in the other words, represents a portion of overall variation in the system. A ratio of below ≤ %25 indicates a strong spatial dependence, while a ratio between % 25 to %75 is indication of median spatial dependence and a ratio of above %75 shows a weak spatial dependence. (Qu et al, 2014)
One of the ways to assess the validity of the geostatistical methods is the method of Cross validation that is a statistical method for testing the integrity of the covariance function and the neighborhood design. The procedure utilizes the kriging algorithm and allows for comparison of the estimated values at the measured value locations, just as one computes residuals between predicted and observed values in regression or analysis of variance (Deutsch and Journel, 1998; Hohn, 1999; Isaaks and Srivastava, 1989; Olea, 1994). Here, our calculation is established based on the root mean square error, mean absolute error and mean error deviation.
Where, N is the number of observation, P_{i} is the Parameters estimated by the model, Q_{i} is the Parameters measured in the laboratory, and P̄, and Q̄ are the mean values for P_{i} and Q_{i.}
Exploratory data analysis
We calculated the primary statistical information of the collected samples, including the distribution of population and statistical index such as mean, median, variance, skewness and kurtosis. For investigating the distribution and normality of the population at 95% confidence level, the KolmogorovSmirnov test was used. In the case of nonnormal distribution of original data, a normal distribution of population was obtained through appropriate logarithmic transformations or the second root of the variable. Descriptive statistical analysis was performed using Statistical Package for the Social Sciences (SPSS) software (version, company name and license number is required). Front shifting pattern was drawn using Geographic Information System (GIS) software for selecting the best pattern.
Prior to applying interpolation methods of spatial statistics the Outliers, trends and isotropic tests were performed in order to achieve a proper data set. We investigated the existence of the trend through applying linear and nonlinear polynomial models on variable values relative to the interval of sampling distance in the X (EastWest) and Y (northsouth) directions.
Interpolation methods
For interpolation, we used inverse distance weighting (IDW) and kriging methods. IDW is an intuitive and efficient methods that is based on the distance between a certain point and a target point that is subject of estimation. A higher weight contribution to the final interpolated values is allocated to the neighboring values than that of contribution from distant observations. In the other words, the influence of a specific data point in interpolated values is inversely proportional to its distance from an estimated target location. IDW works best with evenly distributed points (Azpurua and Dos Ramos., 2010). The weighting factor is calculated using the equation (2):
In which, ʎi is the i station’s weight, Di is the distance between i station and unknown point, and α is a weigh power.
Similar to IDW method, Kriging also uses a weighting approach, which ascribes more weight to the points those are closer to the unknown value for interpolation. Yet, Kriging extends the weighting approach of IDW to include random components where exact point location is not known by the
function. Although it requires substantially more processing time, Kriging is able to incorporate the interdependence of variables and variable error surface output (Legendre and Legendre 1998).
It is said to be the best linear unbiased estimator and is defined as equation (3)
In which, the estimated Z * is the weight or importance of the quantity dependent on the sample i and Z (x_{i}) is the value of the measured variable. Given the linear combination of n, the condition of using this estimate is to have a normal distribution of variable Z. Otherwise, either nonlinear kriging should be used or a normal distribution of variable should be obtained using various transform function (Hassanipak, 1998). The important part of Kriging algorithm is to determine the statistical weights ʎ_{i}. For an unbiased estimation, these weights should be determined such that their sum is equal to one ( ). In addition, to use the Kriging algorithm, the estimated variance should be calculated to minimize the resulting function (see equation 4). (Davis, 1987; Hassanipak, 1998).
Var [z*(x)] = [(z*(x)z(x)) 2] = Min
The multivariate version of Kriging estimation is called cokriging algorithm that is based on correlation between variable of interests in a set of variable (Odeh et al., 1995). It works best where the primary variable of interest is less densely sampled than the others. The Cokriging is formulated as equation (5) (Davis, 1987).
Where, is the estimated value for the point , is the weight of the variable z, is the weight of the auxiliary variable y, Z (xi) is the observed value of the main variable, and Y (xi) is the observed value of the auxiliary variable (Mohammadi, 2006). To estimate the weights, we first need to calculate the variogram according to equation (6):
Where z ^ * (x_i) is the estimated value of the given variable, is the measured value of the given variable (observation value) and N is the number of observations.
Results and discussions
Table 1 summarizes the results of statistical analysis for obtained from surface soil samples (030 cm) of agricultural land in the chat plain. The skewness of data represents their abnormal distribution, and is directly related to the degree of nonuniformity in variance and existence of outliers (Behera and Shukla, 2015). The KolmogorovSmirnov test revealed abnormal distribution of data for phosphorus content in Table.1 The logarithm method was used to normalize data point for phosphorous content that resulted in desirable skewness of about 0.5. The calculate value of %68 for coefficient of variation (CV) of electrical conductivity is also notably larger than that of other parameters. This can be ascribed to the relatively large variability of electrical conductivity of the region under study (Warrick, 1980) caused by natural processes and level of management. In the other words, the changes in degree of drainage and agronomy management along with nonuniform topography of the land can directly influence the electrical conductivity of the soil in different regions (Shakouri et al., 2011). Analysis shows the trends in two main directions EastWest and NorthSouth for all measured parameters.Using variogram analysis of normal data, we investigated the existence of spatial structure among the data that is the first step in employing kriging and Cokriging methods. Due to the symmetry of the semivarogram, all measured parameters were found to be isotropic. This means that we should be able to use Omni directional variogram for next step in our calculations.
Experimental variogram models with models fitted to that in Figure 3 with the validation indices in Table 2 are presented for all variables.
Table 1. The results of SPSS analysis of physicochemical properties of surface soil samples (030 cm) of agricultural land in the chat plain

Figure 3 Illustrates that the multidimensional variogram can be obtained in terms of the quality of soil properties depending on the crossvalidation. After testing different models using validation indices and taking the spatial structure of each of the models into consideration, the suitable one was selected for fitting on the experimental Variogram. Table 2 shows the parameters of the fitting model on experimental variogram and measures of Kriging credit control. The ratio of nugget effect to the Sill can be used to assess the spatial structure of the data. if the ratio was ≤25%, the variable was considered to be strongly spatially dependent; if the ratio was between 25%and 75%, the variable was considered to be moderately spatially dependent; and if the ratio was >75%, the variable was considered to be weakly spatially dependent (Dai et al., 2014; Robinson and Metternicht, 2006).
Given the fact that the ratio of nugget effect to sill for some of the soil properties such as Organic carbon, zinc, Calcium carbonate is less than 0.25, such properties have Strong spatial continuity. In addition, the properties such as Manganese, copper, iron, silt, potassium, pH, and electrical conductivity have moderate spatial structure. The best theoretical model fitted to the experimental variogram is exponential model. The best fitted model for parameters such as: sand, silt and clay, based on cross validation method was spherical. The range parameter in this table is the average of two diameters of an ellipse around sampling points. As the result of the variography is clear, effective range of electrical conductivity, calcium carbonate,
Table 2. The best fitted model to the variogram and its related parameters
Table 3. The correlation between the given variable and dependent variable
**Correlation significance at p<0.01؛* Correlation significance at p<0.05 
potassium, sand, silt, zinc, copper and manganese above 1000 meters and effective range of nitrogen, organic carbon, soil pH, calcium carbonate, phosphate, clay and iron less than 1000. It can be concluded that for each soil sample gathered for measuring the amount of potassium, the radius can be applied to 1108 meters. Thus, if in the studied region, the amount of potassium in soil is required, soil sampling interval can be considered 1108 meters. Relatively short range of influence of nitrogen (503m), is due to its ion mobility. The previous study (Cahn et al., 1994) showed The nitrate had the shortest range of impact and variables such as phosphorus and potassium had average correlation which was because of Ion mobility and management factors such as fertilization and irrigation. The impact range of soil properties is a function of the studied area’s scale, Portrait of sampling distance and location of land. Obviously, the larger the range of impact, the wider spatial structure and in fact the more spatial continuity in the values of variable.
Also, utilizing the results from table 3, auxiliary variables for each of the main variables can be observed. The selection of the appropriate auxiliary variable was based on the coefficient correlation. Based on existing experience, this ratio should be over 0.5, but sometimes was much less in the study. The results showed that despite the fact that the correlation coefficient is less than 0.5, but auxiliary variable helps increase the accuracy of the main variable interpolation. It should be noted that in some cases an element may show a significant correlation with more than just one dependent variable. In this case, the dependent variable is assumed to have the highest correlation with the independent variable. This is reflected the correlation of electrical conductivity, manganese, copper, potassium and nitrogen. Cokriging method can be operating in conditions in which there is a high correlation between the given variable and the dependent variable.
Also a physical relationship between independent and dependent variables is very important and necessary. In this study, 11 variables were considered as independent variables such as (EC, nitrogen, soil pH, phosphorus, potassium, sand, silt, iron, zinc, copper, manganese) and according to investigation, past studies, and existing data, 3 variables were considered as dependent ones such as (Percentage of organic carbon, clay content, and the percentage of calcium carbonate).
After variogram modeling, Kriging, IDW and Cokriging were used to predict soil spatial variability. Then, using the technique of cross validation and the use of validation indices (RMSE, MBE, MAE) the most suitable interpolation method was selected. The results of the evaluation statistic showed that Kriging estimator for parameters of nitrogen, phosphorus, potassium, clay, silt, zinc, copper, manganese is of less error with root mean square error and mean error and mean absolute error compared to the Cokriging and IDW because this method is the best unbiased estimator and shows better results than other interpolation methods. (Schelin and Luna, 2014). The results also indicate that the parameters of the variogram model of nitrogen, phosphorus, potassium, clay, silt, zinc, copper, manganese were well selected. The accuracy of the Cokriging method for interpolating the electrical conductivity and iron parameters in this study was more than the IDW and kriging. Because the accuracy indicators for the Cokriging method were less. These results are consistent with some of the research results conducted by (Nourzadeh et al., 2012), which in their study identified the Cokriging as the best interpolation method for iron, zinc and cobalt. Terms of using Cokriging is when the number of available samples are low. In places, which the number of available sample is low, the estimation is done with help of secondary variables and using crosscorrelation between the primary and secondary variais consibles. According to table 4, in order to estimate some of the soil properties (Iron, and The electrical conductivity) the Cokriging method was of superior than to the other two methods because it can be used when sample size is small and there is lack of data (Odeh et al., 1995). Theoretically, if the initial and secondary variables are nearly equal in size or if the shape and type of the variogram model are identical, the results of the two methods of kriging and cocrigging will be the same. These results were obtained in studies conducted by Voltz and Goulard (1994). Although the Cokriging method is very strong and justifiable in terms of theoretical foundations, this study did not show a significant advantage over the Kriging method. Therefore, although the Cokriging method has a slight advantage over the Kriging method, but due to the computational complexity of the Cokriging method, as well as the difficulty of the multidimensional variogram modeling, the Kriging method is prioritized. The weighted average estimator with a power parity of 1 to 5 for the measured characteristics was evaluated Table (4). The accuracy of the weighted average method depends on the number of near samples used to estimate which can range from 5 to 30 samples for different strengths. The results showed that for all parameters (except for phosphorus, potassium and zinc), the power of one has the least amount of estimation error compared to other powers. For phosphorus, potassium and zinc parameters, power of 5 has the least estimated error. To estimate the parameters of soil (sand, organic carbon, pH and calcium carbonate) inverse distance method with power of one which is of the least error was selected as the best method of interpolation. One of the main reasons for this can be attributed to the semivariogram model with a large segmental variance (C_{0}). It appears that the estimated Kriging methods are stronger than the inverse distance and the Cokriging, but when there is a weak spatial dependency, or there is a limitation in the number of data, the weighted average method is a suitable alternative for different methods of Kriging
Table 4. Accuracy assessment indices of predicted soil properties using different interpolation methods from validation sites

estimation. (Lu and Wong, 2008). Another point to note is that the range of soil quality changes in the estimation with ordinary kriging is less than the range of changes in the IDW estimation. The main reason beyond this difference is the higher Kriging power in softening the values. That is, in Kriging the large values of the parameter are estimated less than the real value and the small ones are estimated more than the real value causing the estimation range to become more limited.
The trend of changes in accuracy assessment indicators (RMSE & MSE) has been the same for all variables. In other words, in all methods and for all variables, increasing the amount of MAE, the RMSE rate also increased. This associated trend between MAE and RMSE was reported in studies conducted by (Robinson and Metternicht, 2006; Yang et al., 2009). Finally, taking into account the best interpolation method, soil zoning maps in “GIS” environment was prepared. These maps can be used for the separation of farms into homogenous smaller pieces for the same management because based on that, plowing, irrigation and fertilization can be done taking into account the ratio of changes and plant requirements. Figure 4 shows the equivalence maps of soil quality variables using three estimation methods (kriging, Cokrigging and weighted average.
Conclusion
The high proliferation rate of population and their demands for foods as well as human beings growing need for environment and natural resources resulted in evaluation and recognition of more inputs, in particular soil and water resources. In this regard, one of the fundamental information regarding land resources is soil properties map. Therefore, with the help of new methods, including geostatistical approaches it could be possible to prepare the soil map with minimal data available. Overall, the survey results showed that geostatistical methods are appropriate methods for estimating soil properties. Parameters measured included percentage of clay, silt, sand, calcium carbonate, organic carbon and concentrations of micronutrients such as iron, copper, zinc, manganese, and concentrations of macronutrients including nitrogen, potassium, and phosphorus. After data normalization, experimental semivariograms was drawn. Estimation by kriging and Cokriging was done utilizing spherical, exponential, circular, and Gaussian models. To estimate variables using Inverse Distance Weighting, the parameters of 15 power was used. To compare the accuracy of estimates crossvalidation techniques and validation indices of the estimated was used. Based on results, in order to estimate some soil properties (nitrogen, phosphorus, potassium, Zinc, copper, manganese, soil electrical conductivity, iron, clay and silt) Kriging method was superior to the two other methods. To estimate the other parameters of soil (sand, organic carbon, pH and lime) inverse distance weighting method with power of one has the least error and was selected as the best method of interpolation. Finally, considering the best method of interpolation, zoning maps of soil quality properties was created in the GIS.
References
Acerbi Júnior, F.W., Silveira, E.M.d.O., Mello, J.M.d., Mello, C.R.d. and Scolforo, J.R.S. 2015. Change detection in Brazilian savannas using semivariograms derived from NDVI images. Ciência e Agrotecnologia, 39(2): 103109.
Azpurua, M.A. and Ramos, K.D. 2010. A comparison of spatial interpolation methods for estimation of average electromagnetic field magnitude. Progress In Electromagnetics Research M, 14: 135145.
Behera, S.K. and Shukla, A.K. 2015. Spatial distribution of surface soil acidity, electrical conductivity, soil organic carbon content and exchangeable potassium, calcium and magnesium in some cropped acid soils of India. Land Degradation & Development, 26(1): 7179.
Black, C.A., Evans, D. and Dinauer, R. 1965. Methods of soil analysis (Vol. 9): American Society of Agronomy Madison, WI.
Bouyoucos, G.J. 1962Hydrometer method improved for making particle size analyses of soils. Agronomy journal, 54(5): 464465
Cahn, M., Hummel, J. and Brouer, B. 1994. Spatial analysis of soil fertility for sitespecific crop management. Soil Science Society of America Journal, 58(4): 12401248.
da Silva, A.F., Barbosa, A.P., Zimback, C.R.L., Landim, P.M.B. and Soares, A. 2015. Estimation of croplands using indicator kriging and fuzzy classification. Computers and Electronics in Agriculture, 111: 111.
Dai, F., Zhou, Q., Lv, Z., Wang, X. and Liu, G. 2014. Spatial prediction of soil organic matter content integrating artificial neural network and ordinary kriging in Tibetan Plateau. Ecological Indicators, 45: 184194.
Davis, B.M. 1987. Uses and abuses of crossvalidation in geostatistics. Mathematical geology, 19(3): 241248.
Deutsch, C.V. and Journel, A.G. 1998. Geostatistical software library and user’s guide. Oxford University Press, New York.
Gerstmann, H., Doktor, D., Gläßer, C. and Möller, M. 2016. PHASE: A geostatistical model for the Krigingbased spatial prediction of crop phenology using public phenological and climatological observations. Computers and Electronics in Agriculture, 127: 726738.
Hassanipak, A. 1998. Geostatistics. Tehran: University of Tehran: Persian
Hohn, M. 2013. Geostatistics and petroleum geology: Springer Science & Business Media.
Isaaks, E.H. Srivastrava RM (1989) An introduction to applied geostatistics: Oxford University Press NY.
Jou, A. 1978. Selected methods of soil analysis. Int. Inst. Trop. Agricult., Ibadan, Nigeria.
Klute, A. and Page, A.L. 1996. Methods of Soil Analysis: Chemical methods (Vol. 3): American Society of Agronomy.
Legendre, P. and Legendre, L.F. 2012. Numerical ecology (Vol. 24): Elsevier.
Lu, G.Y. and Wong, D.W. 2008. An adaptive inversedistance weighting spatial interpolation technique. Computers & Geosciences, 34(9): 10441055.
Mubarak, I., Mailhol, J.C., AnguloJaramillo, R., Ruelle, P., Boivin, P. and Khaledian, M. 2009. Temporal variability in soil hydraulic properties under drip irrigation. Geoderma, 150(1): 158165.
Nourzadeh, M., Mahdian, M.H., Malakouti, M.J. and Khavazi, K. 2012. Investigation and prediction spatial variability in chemical properties of agricultural soil using geostatistics. Archives of agronomy and soil science, 58(5): 461475.
Odeh, I.O., McBratney, A. and Chittleborough, D. 1995. Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regressionkriging. Geoderma, 67(34): 215226.
Oliver, M. and Webster, R. 2014. A tutorial guide to geostatistics: Computing and modelling variograms and kriging. Catena, 113: 5669.
Pizzato, J.A., Araújo, D.V., Galvanin, E.A., Júnior, J.R., Matos, Â.N., Vecchi, M. and Zavislak, F.D. 2014. Geostatistics as a methodology for studying the spatiotemporal dynamics of Ramularia areola in cotton crops. American Journal of Plant Sciences, 5(15): 2472.
Qu, M., Li, W. and Zhang, C. 2014. Spatial distribution and uncertainty assessment of potential ecological risks of heavy metals in soil using sequential Gaussian simulation. Human and Ecological Risk Assessment: An International Journal, 20(3): 764778.
Robinson, T. and Metternicht, G. 2006. Testing the performance of spatial interpolation techniques for mapping soil properties. Computers and Electronics in Agriculture, 50(2): 97108.
Schelin, L. and Sjöstedtde Luna, S. 2014. Spatial prediction in the presence of leftcensoring. Computational Statistics & Data Analysis, 74: 125141.
Shakouri, K.M., Shabanpour, M., Asadi, H., Davatgar, N. and Babazadeh, S. 2011. Evaluation efficiency spatial interpolation techniques in mapping organic carbon and bulk density paddy soils of Guilan.
Sys, I.C., Van Ranst, E., Debaveye, I.J. and Beenaert, F. 1993. Land evaluation (Part IIII). Crop Requirements. 199p. General Administration for Development Cooperation, Brussels. Belgium.
Uyan, M. 2016. Determination of agricultural soil index using geostatistical analysis and GIS on land consolidation projects: A case study in Konya/Turkey. Computers and Electronics in Agriculture, 123: 402409.
Uyan, M. and Cay, T. 2013. Spatial analyses of groundwater level differences using geostatistical modeling. Environmental and ecological statistics, 20(4): 633646.
Voltz, M. and Goulard, M. 1994. Spatial interpolation of soil moisture retention curves. Geoderma, 62(13): 109123.
Voltz, M. and Goulard, M. 1994. Spatial interpolation of soil moisture retention curves. Geoderma, 62(13): 109123.
Wang, Y. and Shao, M. 2013. Spatial variability of soil physical properties in a region of the Loess Plateau of PR China subject to wind and water erosion. Land Degradation & Development, 24(3): 296304.
Warrick, A. 1980. Spatial variability of soil physical properties in the field. Applications of soil physics: 319344.
Yang, P., Mao, R., Shao, H. and Gao, Y. 2009. The spatial variability of heavy metal distribution in the suburban farmland of Taihang Piedmont Plain, China. Comptes rendus biologies, 332(6): 558566.