Soil Fertility Mapping: Comparison of Three Spatial Interpolation Techniques

Soil Fertility Mapping: Comparison of Three Spatial Interpolation Techniques

(Ordinary Kriging, Inverse Distance Weighted, and Spline)

Ismaili Samira, Douaik Ahmed

Dept. Environment and Conservation of Natural Resources National Institute of Agricultural Research (INRA)

Rabat, Morocco

Moughli Lhoussaine

Dept. Natural Resources and Environment

Hassan II Institute of Agronomy and Veterinary Medicine Rabat, Morocco

Abstract: – Soil fertility mapping is essential when planning land use and developing crop fertilization strategies. Interpolation techniques are widely used for the mapping processes in varied fields of soil sciences to estimate the soil property values at unsampled sites. The aim of this work was to evaluate and compare the performances of three spatial interpolation methods (inverse distance weighting IDW; ordinary kriging OK; and spline) using the statistical criterion of root mean square error (RMSE) for cross validation and then generate a set of accurate soil property maps (pH, organic matter, phosphorus, and potassium). The study covers an area of 45 000 ha in the Loukkous irrigated district, Northwest of Morocco and includes 934 soil samples. These samples were collected from irregular cross-line nodes grids and were analyzed in the laboratory. Exploratory data analyses were first adopted to identify and remove all spatial outliers and to validate the normal distribution required for geostatistical analyses. In all cases, the distributions were found strongly skewed and needed to be transformed. Box- Cox transformations were used; they performed well for all soil properties. Experimental variograms were fitted with the exponential, spherical and Gaussian models. With the use of the lowest RMSE approach, ordinary kriging model was selected as the best method compared to IDW and spline for interpolating with the exponential variogram for pH, organic matter and potassium, and the spherical model for phosphorus. A map was generated for each soil property. The maps indicated the low level of potassium and organic matter soil content throughout the study area. These maps could be used for optimizing crop fertilizing considering the different soil fertility levels.

Keywords-: Soil Properties; geostatistics; spatial interpolation; soil maps; kiging; inverse distance weighting; spline,; prediction accuracy.

1. INTRODUCTION

   The inappropriate agricultural practices followed by the farmers and the lack of technical know-how associated with conventional and intensive rainfed agriculture are responsible for most degradation of land resources [1].

   In Morocco, the traditional agriculture caused rapid soil structure degradation with loss of soil fertility rendering conventional agriculture unsustainable. Ait Kadi and Benoit (2012) [2] assert that if the current production practices are continuing, Morocco will face serious food shortages in the very near future, given also the rapid population growth.

   Therefore, the Moroccan agricultural policy has been promoting environment-friendly agriculture concept through the introduction of good farming strategies. The minimization of fertilization inputs is one of the key strategies of this concept. For this purpose, developing skills on the spatial distribution of physical and chemical soil properties is strongly needed for an accurate estimation of plant fertilization requirements in different areas.

   Soil mapping is a process in which the spatial variance of topsoil properties are estimated and exposed in a way that can be understood and analyzed by a wide range of users [3, 4]. A standard method for generating these maps is to sample the area of interest using a grid sampling design and then interpolate the measured variable values of the samples using one of the existing interpolation methods. The spatial interpolation methods make available a tool for estimating the values of soil variable at unsampled points using data from point observations [5, 6].

   There exist a large number of deterministic and stochastic spatial interpolation methods [7]. Spatial interpolation allows converting irregular discrete point data into regular continuous data that can be processed in a Geographic Information System GIS for a better decision making.

   In the literature, there were many published research works that compared the prediction performance of two or more spatial interpolation methods in varied fields of research like climatology [8,9], environmental studies [10,11], oceanography [12], etc. In particular for soil sciences, studies were done for erosion [13, 14], salinity [15, 16], etc. More specifically for soil fertility, comparisons were done for pH [17, 18], organic matter [19], and macronutrients [20, 21].

   However, the results of comparing spatial interpolation methods were not concordant. In some published works, the stochastic geostatistical method of kriging was found better than the deterministic method of the inverse distance weigthed (IDW) [17, 22, 15, 20, 23], whereas other researchers found that this was not true [24, 25, 26]. Therefore, the question about which method is most appropriate to specific conditions is still unsolved. Li and Heap (2011) [27] have compared the key features of the frequent used methods. Unfortunately, at

   present there is no compromise regarding most robust methods, the published studies are not clear and obvious as to which interpolation approach is most accurate to map soil properties [18].

   The principal aim of this work is to evaluate the consistency of maps of topsoil fertility in an irrigated area in Morocco, the Loukkous district, made by kriging, inverse distance weighted, and splines interpolation techniques. More specific objectives are:

   • To produce maps of soil fertility: pH, organic matter, phosphorus and potassium content;

   • To assess and compare the performance and accuracy of the topsoil maps created by kriging, IDW and splines methods,

   • To recommend the most appropriate technique for the soil properties studied in the specific site,

   where N (h) is the number of data pairs within a given class of distance and direction. The semivariances are smaller at shorter distances and then they stabilize at some distance. This can be explained as follows: the values of a target variable are more similar at shorter distances, up to a certain distance where the differences between the pairs are more or less equal to the global variance [6].

   B. Inverse Distance Weighting

   This is one of the simplest and most available methods. Inverse distance weighting directly implements the assumption that a value of an attribute at an unsampled location can be approximated as a weighted average of values at points within a certain cut-off distance, or from a given number m of the closest points (typically 10 to 30). Weights are usually inversely proportional to a power of distance [30, 31]. The formula of this exact interpolator is [7]:

   ( )

   ( )

   (3)

   o To develop soil data systems by improving the prediction accuracy of soil properties.

2. SPATIAL INTERPOLATION TECHNIQUES TESTED

   A. Kriging

   Kriging is a geostatistical interpolation technique which considers both the distance and the degree of variation between known data points when estimating values in unknown areas. It is a weighted linear combination of the known sample values around the point to be estimated [28].

   Following Goovaerts (1997) [29], all kriging estimators are variants of the basic linear regression estimator Z*(u) defned as:

   ( ) ( ) ( ) , ( ) ( )- (1)

   where x0 is the estimation point and xi are the data points within a chosen neighborhood. The weights (r) are related to distance by dij, which is the distance between the estimation point and the data points. The disadvantage of the IDW interpolation technique is that it handles all sample points that fall within the search radius the same way.

   C. Splines

   A spline function is based on a set of interpolating polynomials and an ascending array of domain knot points, determining the intervals over which the spline function is defined by the constituent polynomials [32].

   The prediction value by spline or radial basis functions can be expressed as the sum of two components [33]:

   with

   u,u: location vectors for estimation point and one of the

   ( )

   ( )

   ( ) (4)

   neighboring data points, indexed by ;

   n(u): number of data points in local neighborhood used for estimation of Z*(u);

   m(u), m(u): expected values (means) of Z (u) and Z(u); (u): kriging weight assigned to datum z(u) for estimation at location u; the same datum will receive different weights for

   where w(dj) shows the radial basis functions and dj the distance from sample site to prediction point x, fi(x) is a trend function, a member of a basis for the space of polynomials of degree <m. The coefficients and bj are obtained by solving the system:

   different estimation locations.

   ( )

   ( )

   ( ) for k = 1, 2,.,n (5)

   Kriging uses a property called the semivariance to identify the degree of relationship between points on a surface. The semivariance is simply half the variance of the differences

   ( ) for k =1, 2,, m. (6)

   between all possible points spaced a constant distance apart. The value of the experimental variograms for a separation distance of h (referred to as the lag) is half the average squared difference between the value z(u) and the value z(u + h) separated by a distance and a direction h [30]:

   Splines generate good results with softly varying surfaces,

   and are often not appropriate when there are considerable changes in the surface values within a short horizontal distance [34].

3. MATERIALS AND APPROACH

   ( )

   ( )

   ( ), ( ) ( )-

   (2)

   1. Study area, sampling design and laboratory analysis

      The study area is located in the middle of Kenitra and Larache Provinces, in the North West of Morocco and covering 45 000 Hectares (Fig. 1). The climate is typical

      Mediterranean, with average air temperature between 11°C in winter and 25 °C in summer and mean annual rainfall of 700 mm distributed between October and April. To ensure high added value of agricultural production all this area had been converted to irrigated lands.

      The local topography is divided into three different geomorphologic categories: plains, plateaus and hills. The main soil types are the alluvial soils extending in the plains and the sandy soils dominating the plateaus.

      A total of 934 surface soil samples (030 cm) were collected from the study area (Fig. 1). Soil samples were gathered from irregular cross-line nodes of approximately average distance 700 m * 700 m grids. Each sampling point was geo-referenced using a GPS receiver.

      The soil samples were taken to the laboratory, air-dried and passed through a 2 mm sieve. Available phosphorus (P) was determined by Olsen method (1982) [35], available potassium (K) was determined by extraction with ammonium acetate [36], pH was measured in a 1:2.5 soil/water suspension with the pH meter method [37], and organic matter (OM) was determined by using oxidation method of Walkley-Black [38].

   2. Exploratory data analysis

   The spatial prediction and comparative assessment of the soil properties begin with basic summary statistics, including mean, median, variance and skewness. Summary statistics is the first step in data treatment previous to the application of any data analysis techniques [39].

   transformation and its application in geostatistics is at the origin of the log-normal kriging [42]. Log-normal variables seem to be best performed when outcomes are influenced by various independent elements [43].

   Other transformation functions are also used to achieve the normality such as the Box-Cox transformation [44]. Box-Cox represents a potential best practice to normalize data and offers a range of power transformations that incorporate and extend the traditional options (square root, logarithm, and inverse) to help researchers easily identify the optimal transformation techniques [43].

   D. Interpolation and comparison

   The interpolations are performed by using the Geostatistical Analyst extension of the Arc GIS software [45]. To validate the accuracy of those predictive models, the cross- validation statistical method is applied [46]. Cross-validation involves consecutively eliminating a data point, estimating the value from the remaining observations and comparing the predicted value with the measured one [21]. The cross- validation technique is used generally to choose the best variogram model among proposed models for kriging and also the best parameters from those tested for IDW and splines [47].

   To compare different interpolation methods, the mean error (ME) and the root mean square error (RMSE) calculated from the measured and interpolated values at each sample location are used:

   * ( ) ( )+

   (7)

   Other statistics tools such as histograms, box-plots and

   normal probability plot are available to identify the outliers.

   * ( ) ( )+

   (8)

   Outliers affect the performance of spatial interpolation

   techniques. The variogram is sensitive to outliers and to extreme values because the exceptionally big or small values will distort the average of semivariance [5].

   Outliers should be removed if they are thought to not belong to the population; the elimination of outliers can improve considerably the performance of spatial interpolation methods [40].

   C. Data transformation

   Geostatistical analyses, like many statistical procedures, make the assumption that the data distribution is normal. There are several ways to verify if the variables are normally distributed. These tools range from simple evaluation of the skewness coefficient (ideally closer to 0) and kurtosis (closer to 3) to the evaluation of P-P plots (plotted percentages have to be close to the diagonal line to indicate normality) and inferential tests of normality such as the Kolmorogov- Smirnov, Shapiro-Wilkor or Cramer-von Mises tests (a probability value > 0.05 indicates the normality) [41].

   When non-normality is revealed, the strongly skewed distribution needs to be transformed to make it approximately Gaussian, or at least symetric [5]. The transformation widely applied in different fields of science is logarithmic

   where z(xi) is the observed value at location i, z*(xi) is the interpolated value at location i, and n is the sample size. Ideally ME should tend to zero and RMSE should be as small as possible to indicate less error and more accurate spatial interpolator, respectively.

   Fig.1. Study area location and map of the distribution of soil samples.

4. RESULTS AND DISCUSSION

   1. Summary statistics

      The exploratory analysis revealed outliers in the different soil properties data as shown in Fig. 2. The pH, organic matter, and phosphorus variables showed 2, 6, and 14 outliers, respectively. The exploratory analysis suggested that there were no potential outliers for potassium. In total, the 22 outliers were removed from the original dataset, thus finally we kept 912 samples for the spatial interpolation.

      A statistical summary of the new trimmed data is presented in Table I. Phosphorus and potassium have high coefficient of variation (CV) fluctuating between 64 and 75%. The valuesof their coefficients of skewness and kurtosis do not fit the normality standards, indicating that the variables are highly skewed to the right with a high peak and thin tails.

      However, these descriptive statistics do not provide conclusive information about normality. Further statistical tests (Shapiro- Wilk, Kolmogorov-Smirnov, etc) were applied to examine normality and all of the variables showed a non-normal distribution.

      The Box-Cox transformation was applied to the data. Different parameters for this transformation were found, for pH, OM and phosphorus variables whereas for potassium a particular parameter was found which corresponds to the logarithmic transformation. Those transformations resulted in the best fit of a normal distribution, as the skewness and kurtosis values were near zero.

      All the data processing later, from the variogram computation and the validation tests to the spatial prediction were carried out with the transformed data.

   2. Interpolation

      For every soil property the experimental variogram was calculated. Among the exponential, spherical and Gaussian models the best fitted model to these experimental variograms were chosen using the lowest RMSE as presented in Table II. Fig. 3 illustrates the structure of the fitted variogram models obtained for all soil properties. The pH, OM, and K semivariograms were fitted to an exponential model and the P semivariogram was fitted to a spherical model.

      In both cases of IDW and splines, the best weighting parameters were found using the optimizer parameter tools of the Geostatistical Analyst extension of the Arc GIS software. For IDW the optimal power value was found to be one in all cases.

      On the other hand, the precision of all these interpolation methods is strongly affected by the number of the closest neighbors used for estimation. Therefore, kriging, IDW and splines were implemented using the same neighborhood structure. This latter was divided into eight sectors including a maximum of three and a minimum of two neighbors per sector.

      TABLE I: Summary statistics for pH, organic matter, potassium and phosphorus

      Vol. 3 Issue 11, November-2014

      Soil properties	Statistical analysis
      	N	Min	Max	Mean	Median	Range	Variance	CV (%)	Skewness	Kurtosis
      pH	912	5,58	8,93	7,66	7,77	3,35	0,39	8,1	-0,64	0,04
      OM	912	0,15	3,17	1,26	1,2	3,02	0,33	45,4	0,63	0,25
      K	912	14	822,24	139,02	106,05	808,24	11247	76,3	1,36	2,68
      P	912	2,43	92	26,74	23,79	89,57	287,95	63,5	1,07	1,14

      Fig. 1. Box plots of soil properties; pH, organic matter content, phosphorus content and potassium content.

      Fig. 3. Variograms and fitted models for soil properties: (a) soil pH, (b) organic matter content, (c) phosphorus content and (d) potassium content.

      Table II shows the RMSE based on the cross-validation of all the spatial interpolation techniques tested. The method with the lowest RMSE was chosen and applied to create the soil property maps presented in Fig. 4.

      This figure shows the interpolation of soil pH using the ordinary kriging technique. It reveals that 80 % of the enclosure has a pH higher than 7.5, which means that the soils are alkaline. Barrow (1984) [48] proved that the availability of the phosphorus for plants is optimal in the soil pH range between 5.5 and 7.5. For that reason, it would be beneficial to consider this soil pH map while advising fertilization strategies.

      According to the interpolation map of organic matter using the ordinary kriging, the entire area presents a clear shortage of organic matter with all values below 2.5 %. These results illustrate perfectly the inappropriate agricultural practices followed by the farmers. Recently, it is quite common to see that no crop residues are left on the fields for improving soil fertility.

      Ordinary kriging proved to be the best method for interpolating soil phosphorus as shown in Fig. 4. Since a value greater than 14 ppm suggests good nutrient storage capacity, it can be seen that almost all the soils have sufficient phosphorus. This region has adopted an efficient phosphorus fertilization programs which explain the high values presented in the interpolated map.

      The ordinary kriging is again the best performed method to create the maps of soil potassium as presented in Fig. 4. The map brings out a very low level of soil potassium in this area; over 65 % of the values are below 100 ppm. The lowest potassium is located along the coastal side of the area, which join the results obtained from soils study of the Atlantic coast of Morocco [49, 50]. On the other hand, the highest potassium is located over clayey soils. Mutscher (1978) [51] asserts that the low level of soil potassium is mostly due to a low bioaccumulation and a small degree of alteration of primary minerals in the deep horizons.

      TABLE II: Summary of ME and RMSE based on the cross-validation of all the models tested.

      Soil property	Cross validation results
      	Interpolation method	ME	RMSE
      pH	Spherical*	-0,3367	30,39
      	Exponential*	-0,2789	30,36
      	Gaussian*	-0,3765	30,49
      	IDW	0,028	30,96
      	Spline	-0,1992	30,60
      OM	Spherical*	0,0035	0,4365
      	Exponential*	0,0022	0,4353
      	Gaussian*	0,0044	0,4384
      	IDW	-0,0039	0,436
      	Spline	0,0029	0,4433
      P	Spherical*	-0,0051	1,448
      	Exponential*	-0,0044	1,451
      	Gaussian*	-0,0055	1,448
      	IDW	0,0153	1,472
      	Spline	-0,0084	1,489
      K	Spherical*	-0,0023	0,644
      	Exponential*	-0,0016	0,643
      	Gaussian*	-0,0028	0,646
      	IDW	0,0046	0,652
      	Spline	-0,0004	0,647
      *: kriging using spherical, exponential or Gaussian model for the variogram.

      IJERTV3IS111043

      www.ijert.org

      1640

      Fig. 4. Soil fertility maps created using the best spatial interpolation method and the best variogram model: pH (upper left); OM (upper right); phosphorus (lower left); and potassium (lower right).

      IJERTV3IS111043

      www.ijert.org

      1641

5. CONCLUSION

The elaboration of soil property maps is an important step in implementing the good farming concept. These maps will be used as basis to measure spatial variance, in order to enhance crop productivity and optimize fertilization strategies. In this direction, this study intended to map soil properties, using three spatial interpolation techniques, and reveal among these techniques the most suitable method for each soil attribute.

According to the results of our case study, the ordinary kriging (using either exponential or spherical models) is more accurate for predicting the spatial patterns of the four soil properties (pH, OM, P, and K) than the two other methods (IDW and splines). Actually, the RMSE values do not present sharp fluctuations, but the differences obtained between the methods are sufficient to prove the efficiency of ordinary kriging.

For our dataset, the variables showed a non-normal and strongly skewed distribution which required introducing different transformation methods. The best transformation application was found by applying the Box-Cox method.

The structure of the four variograms points out the weak spatial dependence existing in this specific landscape. Overall, it has been proved that the topography, land use and farming practices affect considerably the spatial dependence of soil properties. These results lead to future research to assess if additional kriging techniques, such as stratified kriging, regression kriging and co-kriging should be applied to include this auxiliary information and improve the accuracy of the spatial predictions and then the quality of soil fertility maps.

Finally, the maps generated highlight the low level of both potassium and organic matter soil content which is clearly explained by the soil typology and the agricultural practices. These conclusions could be useful while implementing new farming management or fertilization orientations.

ACKNOWLEDGMENTS

Special thanks are due to the Regional Center of Agricultural Research-Kenitra (CRRAK), for supporting the field survey. The authors would like to thank the staff of the Laboratory of Soil Analysis of the Agronomy and Veterinary Institute Hassan II-Rabat (IAV).

REFERENCES

1. M. Falkenmark and J. Rockström, Building resilience to drought in desertification-prone savannas in Sub-Saharan Africa: the water perspective, Natural Resources Forum 32, pp. 93-102, 2008.

2. M. Ait Kadi and G. Benoit, Agriculture 2030: a future for Morocco. The futures of agriculture, Brief No. 41, English. Rome: Global Forum on Agricultural Research (GFAR), 2012.

3. P.H.T. Beckett, Soil survey, Agric. Prog. 51, pp. 33-49, 1976.

4. D. Dent and A. Young, Soil Survey and Land Evaluation, George Allen &Unwin, London, England, vol. xiii, 1981.

5. R. Webster and M.A. Oliver, Geostatistics for environmental scientists, Wiley : Chichester, UK, 2001.

6. T. Hengl, A practical guide to geostatistical mapping of environmental variables, JRC scientific and technical report, 143p, 2007

7. P.A. Burrough and R.A. McDonnell, Creating continuous surfaces from point data, In: Burrough, P.A., Goodchild, M.F., McDonnell, R.A., Switzer, P., Worboys, M. (Eds.), Principles of geographic information systems. Oxford University Press, Oxford, UK, 1998.

8. A.D. Hartkamp, K. De Beurs, A. Stein, and J.W. White, Interpolation techniques for climate variables, Natural Resources Group- Geographic Information Systems Series 99-01, CIMMYT, Mexico, p.34, 1999.

9. W. Luo, M.C. Taylor, and S.R. Parker, A comparison of spatial interpolation methods to estimate continuous wind speed surfaces using irregularly distributed data from England and Wales. Int. Jour. Climato. 28, pp. 947959, 2008.

10. D.W. Wong, L. Yuan, and S.A. Perlin, Comparison of spatial interpolation methods for the estimation of air quality data, Journal of Exposure Analysis and Environmental Epidemiology.14, pp. 404415, 2004.

11. J. Hooyberghs, C. Mensink, G. Dumont, and F. Fierens, Spatial interpolation of ambient ozone concentrations from sparse monitoring points in Belgium,. Jour. Environ. Monit. , pp. 1129-1135,2006.

12. V.M. Merwade and D.R. Maidment, J.A. Goff, . Anisotropic considerations while interpolating river channel bathymetry, Jour.

    Hydrology. 331, pp.731-741, 2006

13. C. S. Renschler, C. Mannaerts, and B. Diekkrüger, Evaluating spatial and temporal variability in soil erosion risk – rainfall erosivity and soil loss ratios in Andalusia, Spain, J. Catena. 34, 209-225, 1999.

14. M. Angulo-Martinez, M. Lopez-Vicente, S.M. Vicente-Serrano, and

    S. Beguera, « Mapping rainfall erosivity at a regional scale: a comparison of interpolation methods in the Ebro Basin (NE Spain),J. Hydrol. Earth Syst. Sci. 13, pp.1907-1920, 2009.

15. E. Hosseini, J. Gallichand, and D. Marcotte, Theoretical and experimental performance of spatial interpolation methods for soil salinity analysis, J. Transactions of American Society Of Agricultural Engineers. ASAE.37, pp. 1799-1807, 1994.

16. A. Douaik, R. Moussadek, K. Barhmi, and H. Bouksirat, Spatial interpolation methods for soil salinity mapping, In Hamdi A (ed.). Proceeding of the International Conference on Water, Land and Food Security in Arid and Semi-arid Regions.Mediterranean Agronomic Institute, Bari (Italy), 6-11 September, p.11p, 2005.

17. G.M. Laslett, A.B. McBratney, P.J. Pahl, and M.F. Hutchinson,

    Comparison of several spatial prediction methods for soil pH,J. Soil Sci. 38, pp. 325341, 1987.

18. W. Shi, J. Liu, Z. Du, Y. Song, C. Chen, and T. Yue, Surface modelling of soil pH,J. Geoderma.150, pp. 113119, 2009.

19. L. Mabit, and C. Bernard, Spatial distribution and content of soil organic matter in an agricultural field in eastern Canada, as estimated from geostatistical tools,J. Earth Surf. Process. Landforms. 35, pp. 278283,2010.

20. A. Kravchenko, and D.G. Bullock, A comparative study of interpolation methods for mapping soil properties, Agron. Jour. 9, pp. 393-400, 1999.

21. T.G. Mueller, N.B. Pusuluri, K.K. Mathias, P.L. Cornelius, R.I. Barnhisel,and S.A. Shearer, Map quality for ordinary kriging and inverse distance weighted interpolation, J. Soil Sci. Soc. Am. Jour. 68, pp. 20422047, 2004.

22. H. Leenaers, J.P. Okx, and P.A. Burrough, Comparison of spatial prediction methods for mapping floodplain soil pollution, J. Catena. 17, pp. 535-550, 1990.

23. T. Panagopoulos, J. Jesus, M.D.C. Antunes, and J. Beltrao, « Analysis of spatial interpolation for optimising management of a salinized field cultivated with lettuce, J. Europ. Jour. Agronomy. 24, pp 1-10, 2006.

24. D. Weber, and E. Englund, Evaluation and comparison of spatial interpolators,J. Math. Geol. 24, pp. 381-389, 1992.

25. N.C. Wollenhaupt, R.P. Wolkowski, and M.K. Clayton, Mapping soil test phosphorus and potassium for variable-rate fertilizer application,

    J. Production Agric. 7, pp. 441-448, 1994.

26. C.A. Gotway, R.B. Ferguson, G.W. Hergert, and T.A. Peterson,

    Comparison of kriging and inverse-distance methods for mapping soil parameters, J. Soil Sci. Soc. Am. Jour. 60, pp.1237-1247, 1996.

27. J. Li, and A.D. Heap, A review of spatial interpolation methods for environmental scientists: performance and impact factors, J. Ecological Informatics. 6, pp. 228-241, 2011.

28. E.H. Isaaks, and R.M. Srivastava, An introduction to applied geostatistics,Oxford University Press, New York, p. 542, 1989.

29. P. Goovaerts, Geostatistics for Natural ResourcesEvaluation.Applied Geostatistics, Oxford University Press, New York, p. 483, 1997.

30. P. A. Burrough, .Principles of geographical information systems forland resources assessment, Oxford, Clarendon Press, 1986.

31. D. F. Watson, Contouring: a guide to the analysis and display of spatial data, Oxford, Pergamon, 1992.

32. APACHE, http://www.apache.org/. Consulted on 26-9-2014, 2014.

33. H. Mitasova, and L. Mitas, Interpolation by regularized spline with tension: I. Theory and implementation, J. Math. Geol. 25, pp. 641 655, 1993

34. T.P. Robinson, and G. Metternicht, Testing the performance of spatial interpolation techniques for mapping soil properties, J.

    Computers and Electronics in Agriculture.50, pp. 97-108.2006

35. S.R. Olsen, and L.E. Sommers, Phosphorus, In: Page AL, et al (eds), Methods of soil analysis, Part 2, 2nd edn, Agron Monogr 9. ASA and ASSA, Madison WI, pp. 403430, 1982.

36. L.A. Richards, Diagnosis and improvement of saline and alkaline soils, USDA Hand b. No. 60. U.S. Gov. Print. Office, Washington DC, p. 160, 1954.

37. A. L. Page, R. H. Miller and D. R. Keeney, Methods of soil analysis, chemical and microbiological properties, The Series Agronomy, Part 2, Soil Science Society of America, Madison, 1982.

38. C.A. Black, Methods of soil analysis, chemical and microbiological properties, The Series Agronomy, Part 2, 9, pp.771-1572, 1965.

39. C. De Fouquet, From exploratory data analysis to geostatistical estimation: examples from the analysis of soil pollutants, J. European Journal of Soil Science. 62, pp. 454 466, 2011.

40. G.M. Laslett, and A.B. McBratney, Further comparison of spatial methods for predicting soil pH, J. Soil Sci. Soc. Am. Jour. 54, pp.1553-1558, 1990.

41. J. W. Osborne, Best practices in data transformation: the overlooked effect of minimum values, In J. W. Osborne (Ed.), Best practices in quantitative methods. Thousand Oaks, CA: Sage Publishing, 2008.

42. N.A. Cressie, Statistics for spatial data (revised edition). John Wiley & Sons, Inc., New York.1993.

43. J.W. Osborne, Improving your data transformations: applying the Box- Cox transformation, J. Practical Assessment, Research & Evaluation,. 15 (12), pp. 1-9, 2010.

44. G. E. P., Box, and D. R. Cox, An analysis of transformations , J. Journal of the Royal Statistical Society, B. 26, pp. 211-234, 1964.

45. ESRI., ArcGIS geostatistical analyst: statistical tools for data exploration, modeling and advanced surface generation. An ESRI white paper, 2001.

46. M. Voltz, and R. Webster, A comparison of kriging, cubic splines and classification for predicting soil properties from sample information, J. Soil Sci. 41, pp. 473490, 1990.

47. M. Tomczak, Spatial interpolation and its uncertainty using automated anisotropic inverse distance weighting (IDW) crossvalidation/jackknife approach, J. Geographic Inf. Decis. Anal. 2, pp. 1830, 1998.

48. N.J. Barrow, Modelling the effects of pH on phosphate sorption by soils, J. Journal of Soil Science. 35, pp. 28329, 1984.

49. SASMA., Caractérisation des sols dans le secteur tramé de Sidi- Bennour, Etude réalisée au compte de la sucrerie des Doukkala, p. 48, l987.

50. A. Aït Houssa, Etude du potassium dans divers types de sols et systèmes de culture au Maroc, Thèse Dr. Sc. Agro., I .A.V Hassan Il, Rabat, 1989.

51. M. Mutscfier, Recherches sur le régime du K dans les sols typiques du Nord de I'Algérie. Lévaluation globale de l'état d'approvisionnement en K des sols, Revue de la potasse, sect. n°4, 1978.