
Research Article


Spatial Prediction of Surface Soil Properties Using Terrain and Remote Sensing Data 

Jafar Yasrebi,
Mahboub Saffari,
Hamed Fathi,
Mostafa Emadi,
Majid Baghernejad,
AbdolMajid Ronaghi
and
Mehdi Emadi



ABSTRACT

The main objective of this research is to enhance prediction of soil properties such as Electrical Conductivity (EC_{e}), Exchangeable Sodium Percentage (ESP), available Phosphorus (P), Organic Matter (OM), Total Nitrogen (TN) and pH by making use of the ancillary variables as covariates. Methods that was used for this purpose may be divided into two groups: (i) those that use only a single variable in the prediction process Simple Linear Regression (SLR), Ordinary Kriging (OK)) and (ii) another that make use of additional variables as a part of prediction Simple Kriging with a Locally Varying Mean (SKLVM)). LISSIII data from Indian remoter sensing satellite (IRSP6) were used as secondary data with SKLVM method. Mean Square Error (MSE) was used to evaluate the performance of the map prediction quality. It was concluded that SKLVM method provided the most accurate predictions based on the summary statistics of prediction errors from crossvalidation for mapping OM, pH and EC_{e}. Maps from these kriged estimates showed that a combination of geostatistical techniques and digital data from LISSIII receiver could improve the prediction of quality soil management zones, which is the first step for sitespecific soil management.





INTRODUCTION
Efficient tools to measure withinfield spatial variation in soil are important when establishing agricultural field trials and in precision farming. Spatial prediction of soil properties has become a common topic in soil science research. Variability in soil properties can present management challenges to producers (Goovaerts, 2000). This is enhanced by the advancement of technology that enabled collection of remotelysensed imagery for using in precision agriculture and digital soil mapping (LopezGranados et al., 2005). Sitespecific management or precision agriculture seeks to identify, analyze and manage spatial and temporal variability within fields in order to optimize profitability, sustainability and environmental protection (Robert et al., 1996; Duffera et al., 2007).
Chaplot et al. (2000) demonstrated that quality of soil hydromorphy prediction was highly improved by cokriging 60 pedological data points with a topographical regression model. Simbahan et al. (2006) reported that for reducing uncertainties, it was recommended that use independently measured, multivariate secondary information in regression kriging approaches for mapping of soil organic carbon. They indicated that geostatistical methods that utilized spatially correlated secondary information increased the quality of maps of soil organic carbon stock as compared to OK. Also, showed that regression kriging with EC_{e} performed better than OK, kriging with an external drift or cokriging. Minasny and McBratney (2007) concluded that improvement in the prediction of soil properties does not rely on more sophisticated statistical methods, but rather on gathering more useful and higher quality data. Due to high cost and timeconsuming nature of soil sampling and their analysis, research in developing methods for the creation of soil maps from sparse soil data is becoming increasingly important. In the past 20 years, the development of prediction methods that use cheap secondary information to spatially extend sparse and expensive soil measurements has been a sharpening focus of research (Odeh and McBratney, 2000; Minasny and McBratney, 2007; McBratney et al., 2003; LopezGranados et al., 2005).
Examples of secondary information include remote sensing imagery, elevation data and crop yield data. Furthermore, a number of proximal soil sensors are becoming more available; examples are the Soil Doctor Colburn (1999) and the VERIS conductivity cart (Lund et al., 1999). Consequently, the potential for using the secondary information to aid soil mapping withinfield extent is greater than ever before.
The introduction of ancillary, exhaustive spatial information linked to salinity might improve the mapping of this attribute. Many researchers have shown that remote sensing, particularly within the visible spectral range, yields spatial information strongly correlated with salinity (Mougenot, 1993; Rahman et al., 1994; Khan et al., 2001; Metternicht and Zinck, 2003) or soil surface features with additional microwave remotely sensed data (Metternicht and Zinck, 1998; Sumfleth and Duttmann, 2007). Moreover, several researchers have demonstrated the advantage of combining data from remote sensing with pinpoint information observed on the ground (Bishop and McBratney, 2001; Carré and Girard, 2002). Several studies have explored the potential benefits of using secondary information to map different variables. For example, Goovaerts (2000) compared Thiessen polygons, Inverse Distance Weighting (IDW), SKLVM, kriging with an external drift and collocated cokriging for soil mapping properties. With the latter three approaches elevation was used as the secondary variable and these techniques resulted in smaller prediction errors than the univariate ones.
Remotely sensed data can be useful for improving existing coarsescale soil survey information at a regional scale. Thus, Odeh and McBratney (2000), Bishop and McBratney (2001) and Kerry and Oliver (2003) demonstrated that AVHRR (Advanced Very High Resolution Radiometer) data from the NOAA (National Oceanic and Atmospheric Administration) series of satellites, bare soil LANDSAT TM (Thematic Mapper) imagery and bare soil aerial colour photograph have been useful for the fieldextent creation of different soil property maps using different prediction models (statistical and geostatistical techniques). They used the soil spectral variation for spatial prediction of soil attributes at a regional scale. The major objectives of our study was to compare different prediction methods such as Simple Linear Regression (SLR), Ordinary Kriging (OK) and method that make use of secondary or ancillary variables (LISSIII receiver data from IRSP6) as a part of prediction such as Simple Kriging with a Locally Varying Mean (SKLVM) in order to determine the most appropriate approach to spatially transfer data from a limited number of sampling points to unsampled locations.
MATERIALS AND METHODS
Study site: The study area is found in Badjgah plain that located in Fars province, Southern Iran at geographical coordinates of 29°42` to 29°46` N latitude and 53°10` to 53°17` E longitude. Soils of the region were developed over the limestone parent material. The mean annual precipitation, evaporation and temperature are 333.4 and 919.1 mm and 15.2°C, respectively. Soil moisture and temperature regime are xeric and thermic, respectively. According to the USDA Soil Taxonomy (Soil Survey Staff, 2006), the soil at the study region was classified as fine, mixed, mesic, Fluventic Calcixerepts.
Soil sampling, laboratory analysis and remotely sensed data: Eightyfive soil samples, on 08 September 2006, from the top soil (depth of 030 cm) were collected and georeferenced using GPS receiver (accuracy of ±5 m), analyzed for EC_{e}, P, OM, pH, TN and ESP. ESP was determined using the ammonium acetate (NH_{4}OAc) method (Thomas, 1982); soil pH was measured with a glass electrode pH meter (McLean, 1982). Soluble salts were calculated by measurement of EC_{e} in the soil extraction by the use of a conductivity meter (Rhoades, 1982). OM was determined dichromate oxidation procedure (Allison, 1965). Total Nitrogen (TN) was determined with the Kjeldahl method (McGill and Figueiredo, 1993), available phosphorus (P) was measured by the Olsen method (Olsen et al., 1954).
Remote sensing data is now considered as an appropriate tool for deriving information in spatial and temporal domains by providing multispectral reflectance data at regular intervals in a synoptic mode. The satellite data used in this research is IRSP6 scene, dated 08 September 2006. Geometric corrections were carried out in conjunction with the atmospheric correction. The imaging sensors on IRSP6 that was used is a multispectral Linear Imaging SelfScanner (LISSIII) in visible (0.520.59 μm, band 1; 0.620.68 μm, band 2), nearIR spectral bands (0.770.86 μm, band 3) with spatial resolution of around 23 meters and a Short Wave IR (SWIR) band (1.551.75 μm, band 4) with a resolution of around 70 m.
Prediction methods: A brief description of the prediction methods used is given below.
Simple Linear Regression (SLR): Every sampled soil point was located
in the satellite image and its corresponding digital value in four bands (1,
2, 3 and 4) was extracted. It was verified that all variables (i.e., soil properties
and spectral values in visual range) were normally distributed. Pearson linear
correlations were determined between soil variables and spectral values in four
bands, accepting a confidence level of 95%. Regression equations were calculated
for those soil variables that showed higher significant correlations with digital
data. It should be imply that band combinations and principal component analysis
obtained from these four main bands had not more accuracy; therefore, the results
presented of the four main bands in Table 1.
Ordinary Kriging (OK): OK with a global variogram was used as a basis
of comparison with other methods, as predictions may only be derived from ground
measurements
where, m is the number of neighbours considered and wi (si) are the weights derived from variogram fitting (Goovaerts, 1997).
Simple Kriging with Locally Varying Mean (SKLVM): Simple Kriging (SK)
is the most basic form of kriging. With SK, the mean is assumed to be constant
and known. If we can estimate the mean at locations in the domain of interest
then this locally varying mean can be used to inform prediction. SKLVM prediction
is defined as:
where, simple
kriging is a known locally varying mean. The locally varying mean can be estimated
in various different ways. One approach is to use regression (presented in SLR
approach) to predict at all observation locations and all locations where SKLVM
predictions will be made. Then, the semivariogram of the residuals was computed,
modelled, crossvalidated and simple kriging on the residuals was carried out.
The final estimate of every soil property was obtained by adding the trend estimate
to the simple kriged estimate of the residuals (Goovaerts, 1997; Vanderlinden,
2001). This method was applied to the soil variables showing significant correlations
with digital values in four bands at p<=0.01, i.e., OM, pH and EC_{e}
with band 1 (Table 1).
Comparison between the different methods: For the purpose of comparison,
several comparison indices can be used as a measurement of the prediction quality,
however, the most common of which is the Mean Square Error (MSE) which measures
the average square difference between the actual soil variable z(xi) and its
estimate z*(xi)where, n = soil variable data set (Goovaerts, 2000). The comparative
performance of the prediction models was measured by using MSE of OK as the
standard, which did not take into account the digital numbers (Bishop and McBratney,
2001). The MSE of OK was calculated as reported in LopezGranados et al.
(2005).
Table 1: 
Pearson linear correlations between soil parameters and spectral
values for four bands 

*Significant at 0.05 level, **Significant at 0.01 level, ^{†}Not
significant at 0.05 level 
RESULTS AND DISCUSSION
Simple Linear Regression (SLR): Pearson linear correlations between
soil parameters and spectral values (Table 1) revealed that
OM, pH and EC_{e} showed highest correlations by using the spectral
data in the band 1, although these correlations were relatively moderate (0.50.6).
Negative correlations meant that small digital numbers for the band 1 corresponded
to high values of pH. On the other hands, P, TN and ESP were the soil properties
having the lowest correlation coefficients with remotely sensed data. Regression
equations for the highest correlation coefficients between soil properties and
spectral data are shown in Table 2. In all cases, the band
1 from LISSIII receiver data was used for fitting the regression equations
because its correlation coefficients were higher and significant. Maps of OM,
pH and EC_{e} could be easily illustrated in remote sensing software
such as ILWIS using presented equations in Table 2. High negative
correlation of soil reaction (pH) and digital number of remotely sensed data
that was demonstrated by LopezGranados et al. (2005) could be attributed
to indirect effect of soil pH on soil physical and chemical properties.
Ordinary Kriging (OK): Table 3 shows mean, Coefficient
of Variation (CV), standard deviation and skewness of the soil parameters. Skewness
is the most common form of departure from normality. If a variable has positive
skewness, the confidence limits on the variogram are wider than they would otherwise
be and consequently, the variances are less reliable. A logarithmic transformation
is considered where the coefficient of skewness is greater than 1 and a squareroot
transformation applied if it is between 0.5 and 1 (Webster and Oliver, 2001).
Therefore, a logarithmic transformation performed for EC_{e}, pH and
ESP parameters because, their skewness was greater than 1. The CV of soil properties
except pH and TN were fairly high, indicating that soil properties were generally
heterogeneous (Table 3). The highest CV value was for ESP,
while the CV value for pH was the lowest. Anisotropic semivariograms did not
show any differences in spatial dependence based on direction and therefore
isotropic semivariograms were chosen. The geostatistical analysis indicated
two spatial distribution models and spatial dependence levels for the soil parameters.
Exponential and spherical models were used to define soil properties (Table
4). Nugget effect was higher for ESP, TN and P compared to pH, OM and EC_{e}.
This indicated that these soil properties had spatial variability in small distances.
The large nugget semivariance suggest that an additional sampling of these variables
at smaller lag distances and in larger numbers is needed to detect spatial dependence.
When the distribution of soil properties is strongly or moderately spatially
correlated, the mean extent of these patches is given by the range of the semivariogram.
A larger range indicates that observed values of the soil variable are influenced
by other values of this variable over greater distances than soil variables
which have smaller ranges (LopezGranados et al., 2005). Range value
varied from 1811 m (for pH) to 4924 m (for OM). Generally, range values of EC_{e}
and pH were smaller than that of the other soil properties. The low nugget variance/total
variance ratio and small range values for some soil properties exhibited patchy
distribution pattern.
Table 2: 
Regression equations between soil parameters and spectral
values for the first band 

Table 3: 
Descriptive statistics for studied soil properties 

^{a}: Standard deviation; ^{b}: Coefficient
of variations 
Table 4: 
Semivariogram models and models parameters for studied soil
properties 

(+): Spatial distribution (Sstrong spatial dependence (<=25%);
MModerate spatial dependence (2675%); Wweak spatial dependence (>75%);
Pure Nugget no spatial correlation (100%) and their spatial distribution
model, *Residual sum of squares (often the model with the lowest RSS chooses
as optimal) 
Simple Kriging with Locally Varying Mean (SKLVM): This method was applied
to the soil variables showing significant correlations with digital values in
four main bands at p<=0.01, i.e., O M, pH and EC_{e} with band 1 (Table
1). This kriging method is an interpolation that incorporates secondary
information into the kriging system. It uses the ancillary (or secondary) information
to characterize the spatial trend of the primary (target) variable and performs
simple kriging on the residuals (Goovaerts, 1997). Table 5
shows the semivariogram of residuals for OM, pH and EC_{e} with the
fitted model. Figure 1 shows the maps of OM, pH and EC_{e}
estimates obtained by SKLVM. The Nugget effect, sill, semivariogram model and
range of the residuals semivariograms were approximately similar of raw semivariogram.
Table 5: 
Semivariogram models of the residuals 

Sill semivariance of OM was 2.247 and 2.127, for EC_{e} it was 1.1806
and 1.1906, for pH it was 1.2148 and 1.2811 for the raw soil parameter semivariogram
and the residual semivariogram, respectively, indicating the lag distance between
measurements at which one value for a variable does not influence neighboring
values (Table 4, 5). Goovaerts (1999) found
a similar trend when he incorporated a digital elevation model into the mapping
of annual erosivity values using the same kind of kriging. The residual semivariogram
model of pH was pure nugget, which was similar to that of the raw semivariogram
and means that pH and EC_{e} were considered strongspatially correlated.
There is some similarity in the map pattern of OM, pH and EC_{e} as
produced by OK and SKLVM methods (Fig. 1). However, OK oversmoothed
the spatial variability of OM, pH and EC_{e}. Comparatively, it seems
that SKLVM reflects local variation more than OK, but it is necessary to compare
the MSE to evaluate this.
Comparison between different prediction methods: MSE for different methods
for pH, EC_{e} and OM (Table 6) shows that the generic
geostatistical technique, such as OK, exhibited the highest MSE because it does
not take into account the secondary information and only uses the primary soil
variable. Comparing the other prediction methods, higher predictions errors
were obtained when only SLM was considered for ESP, P and TN. SKLVM was clearly
the best method for the prediction of OM, pH and K showing the lowest MSE values.
Also, the best prediction method for mapping ESP, P and TN was obtained from
SLM.

Fig. 1: 
Soil maps obtained using OK method: (a) pH, (c) ECe (dS m¯1)
and (e) OM (%) and soil maps obtained using SKLVM: (b) pH, (d) ECe (dS m¯1)
and (f) OM (%) 
Table 6: 
Mean square errors for the compared prediction methods 

Bourennane et al. (2000) compared linear regression with kriging with
an external drift (This method also uses secondary variable) for mapping soil
horizon thickness, using slope gradient as the secondary variable. They found
that kriging with an external drift provided more accurate predictions than
linear regression. PardoIgu´zquiza and Dowd (2002) compared OK, SKLVM, kriging
with an external drift, cokriging and Bayesian integration for prediction using
wirelinelogs of acousticimpedance recorded at nine boreholes, with acousticimpedance
values from a 3D seismic survey as secondary data. On the basis of the Mean
Absolute Relative Error (MARE) cokriging provided the most accurate predictions,
followed by SKLVM. However, on the basis of the Mean Squared Error (MSE) and
other statistics the authors considered SKLVM the preferred method in their
application.
In general, estimation method using spectral data had more favorable MSE results than prediction methods using only soil variables, indicating that the correlation of soil variables with spectral data is more important for mapping soil variables than the spatial correlation of available soil measurements. Thus, the least accurate estimation for OM, pH and EC_{e}, with highest MSE values, was OK because the spectral data were ignored and only the spatial component of soil variables was considered. On the other hand, SLR method resulted in the poorest prediction because low correlations between soil attributes and spectral values were obtained, probably due to the spatial component being ignored. Despite this, when secondary information is available it should be incorporated into map soil attributes because the MSE for simple linear regression is lower than the MSE for OK. Bishop and McBratney (2001) found that when kriging of the residuals was incorporated to different prediction models studied, the RMSE (root mean square error) was lower. They indicated that kriging with an external drift, an interpolation method very similar SKLVM which also performs kriging on the residuals, was the best prediction method for mapping soil properties using bare soil aerial photograph. They also concluded that when secondary information is available, it should be used because generic geostatistical techniques that only use the primary variable, such as OK, do not obtain the prediction performance of the methods that incorporate that secondary information.
CONCLUSIONS
The results here indicate that when secondary information especially remote sensed data is available, it should be used to model the deterministic trend in the variation of a soil surface attribute. Also, this study has demonstrated that sparse and expensive soil measurements combined with secondary information, such as remotely sensed (spectral) data from IRSP6 and geostatistical techniques were adequate to map soil properties accurately. The relationship between spectral (digital) data in the band 1 from LISSIII receiver data and OM, pH and EC_{e} developed in this research might be applied to other fields in southern Iran, especially in all Badjgah plain soil that have the same qualities. For variables presenting a high or moderate correlation with spectral data, as secondary information, kriging with varying local means results accurate estimates because it uses spectral data to derive the local mean or trend of any soil property. Therefore, precision farming could be benefited by such enhanced technique where the data of remote sensing or other cheap valuable secondary variable of soil parameter are available.

REFERENCES 
1: Allison, L.E., 1965. Organic Carbon. In: Methods of Soil Analysis, Black, C.A. (Ed.). American Society of Agronomy, Madison, WI., pp: 13671382.
2: Bishop, T.F.A. and A.B. McBratney, 2001. A comparison of prediction methods for the creation of fieldextent soil property maps. Geoderma, 103: 149160. CrossRef  Direct Link 
3: Bourennane, H., D. King and A. Couturier, 2000. Comparison of kriging with external drift and simple linear regression for predicting soil horizon thickness with different sample densities. Geoderma, 97: 255271. CrossRef  Direct Link 
4: Carre, F. and M.C. Girard, 2002. Quantitative mapping of soil types based on regressionkriging of taxonomic distances with landform and landcover attributes. Geoderma, 110: 241263. Direct Link 
5: Chaplot, V., C. Walter and P. Curmi, 2000. Improving soil hydromorphy prediction according to DEM resolution and available pedological data. Geoderma, 97: 405422. CrossRef  Direct Link 
6: Colburn, J.W., 1999. Soil doctor multiparameter real time soil sensor and concurrent input control system. Proceedings of the 4th International Conference on Precision Agriculture, July 1922, 1999, ASA, CSSA and SSSA, Madison, WI., pp: 10111022.
7: Duffera, M., G.W. Jeffrey and R. Weisz, 2007. Spatial variability of southeastern U.S. Coastal plain soil physical properties: Implications for sitespecific management. Geoderma, 137: 327339. Direct Link 
8: Goovaerts, P., 1997. Local Estimation: Accounting for Secondary Information: Geostatistics for Natural Resources Evaluation. Oxford University Press, UK., pp: 185197.
9: Goovaerts, P., 1999. Using elevation to aid the geostatistical mapping of rainfall erosivity. Catena, 34: 227242. CrossRef  Direct Link 
10: Goovaerts, P., 2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. J. Hydrol., 228: 113129. CrossRef 
11: Kerry, R. and M. Oliver, 2003. Variograms of ancillary data to aid sampling for soil surveys. Precision Agric., 4: 261278. CrossRef 
12: Khan, N.M., V.V. Rastoskuev, E. Shalina and Y. Sato, 2001. Remote sensing indicators. A simple approach with the use of GIS Idrissi. 22nd Asian Conference on Mapping salt affected soil using Remote Sensing, Singapore.
13: LopezGranados, F., M. JuradoExpósito, J.M. PenaBarragan and L. Garc2aTorres, 2005. Using geostatistical and remote sensing approaches for mapping soil properties. Eur. J. Agron., 23: 279289. Direct Link 
14: Lund, E.D., P.E. Colin, D. Christy and P.E. Drummond, 1999. Applying soil conductivity technology to precision agriculture. Proceedings of the 4th International Conference on Precision Agriculture, July 1922, 1999, ASACSSASSSA, Madison, WI., pp: 10891100.
15: McBratney, A.B., M.L. MendonçaSantos and B. Minasny, 2003. On digital soil mapping. Geoderma, 117: 352. Direct Link 
16: McGill, W.B. and C.T. Figueiredo, 1993. Total Nitrogen. In: Soil Sampling and Methods of Analysis, Carter, M.R. (Ed.). Lewis Publishing, Canadian Society of Soil Science/Lewis Publishers, Boca Raton, pp: 201211.
17: McLean, E.O., 1982. Soil pH and Lime Requirement. In: Methods of Soil Analysis, Part 2: Chemical and Microbiological Properties, Page, A.L., R.H. Miller and D.R. Keeney (Eds.). 2nd Edn., ASA and SSSA, New York, USA., pp: 199224.
18: Metternicht, G.I. and J.A. Zinck, 1998. Evaluating the information content of JERS1 SAR and landsat TM data for discrimination of soil erosion features. ISPRS. J. Photogramm. Remote Sens., 53: 143153. CrossRef  Direct Link 
19: Metternicht, G.I. and J.A. Zinck, 2003. Remote sensing of soil salinity: Potentials and constraints. Remote Sens. Environ., 85: 120. Direct Link 
20: Minasny, B. and A.B. McBratney, 2007. Spatial prediction of soil properties using EBLUP with the matérn covariance function. Geoderma, 140: 324336. Direct Link 
21: Mougenot, B., 1993. Effets des sels sur la réflectance et télédétection des sols salés. Cah. ORSTOM, Ser. Pédol., 28: 4554.
22: Odeh, I.O.A. and A.B. McBratney, 2000. Using AVHRR images for spatial prediction of clay content in the lower Naomi Valley of eastern Australia. Geoderma, 97: 237254. Direct Link 
23: Olsen, S.R., C.V. Cole, F.S. Watanabe and L.A. Dean, 1954. Estimation of available phosphorus in soils by extraction with sodium bicarbonate. USDA Circular No. 939, United States Department of Agriculture, Washington, DC., USA., pp: 118.
24: PardoIgúzquiza, E. and P.A. Dowd, 2002. Geostatistical Integration of Primary and Secondary Data for Generating More Realistic Stochastic Petrophysical Models. GeostatsUK, Book of Abstracts.
25: Rahman, S., G.F. Vance and L. Munn, 1994. Detecting salinity and soil nutrient deficiencies using spot satellite data. Soil Sci., 158: 3139. Direct Link 
26: Rhoades, J.D., 1982. Soluble Salts. In: Methods of Soil Analysis, Page, A.L., R.H. Miller and D.R. Keeney (Eds.). American Society of Agronomy, Madison, WI., ISBN: 0891180729, pp: 167185.
27: Robert, P.C., R.H. Rust and W.E. Larson, 1996. Precision Agriculture. ASA, CSSA, SSSA, Madison, WI.
28: Simbahan, G.C., A. Dobermann, P. Goovaerts, J. Ping and M.L. Haddix, 2006. Fineresolution mapping of soil organic carbon based on multivariate secondary data. Geoderma, 132: 471489. CrossRef  Direct Link 
29: Soil Survey Staff, 2006. Kets to Soil Taxonomy, Natural Resources Conservation Service. US Department of Agriculture.
30: Sumfleth, K. and R. Duttmann, 2007. Prediction of soil property distribution in paddy soil landscapes using terrain data and satellite information as indicators. Ecol. Indicators (In Press).
31: Thomas, G.W., 1982. Exchangeable Cation. In: Methods of Soil Analysis, Part 2, Page, A.L., R.H. Miller and D.R. Keeney (Eds.)., 2nd Edn. American Society of Agronomy, Madison, WI., ISBN: 0891180729, pp: 159164.
32: Vanderlinden, K., 2001. An´alisis de procesos hidrol´ogicos a diferentes escalas espaciotemporales. Ph.D Thesis, Universidad de C´ordoba, Spain.
33: Webster, R. and M. Oliver, 2001. Local Estimation or Prediction: Kriging, Geostatistics for Environmental Scientists. John Wiley and Sons, England, pp: 149191.



