Australian Cotton Cooperative Research Centre
Department of Agricultural Chemistry & Soil Science
The University of Sydney, Sydney, NSW
Phone: (02) 9351 4178; Fax: (02) 9351 5108
Email (Corresponding author): firstname.lastname@example.org
A common misgiving by traditional pedologists is that quantitative methods, especially geostatistical interpolation techniques, are not suitable for regional soil survey and land evaluation. A major reason is that these methods require too many field samples for them to accurately predict land features. This notion is not always plausible. With ubiquitous auxiliary data describing the soilscape, especially remote sensing imageries, it is possible to apply the quantitative methods that can produce good results. To illustrate this several examples of soil inventories at the regional extent (≥ 50 km x 50 km), using regression models, the geostatistical technique of kriging and the hybrid regression-kriging, are provided. The examples prove that soil classes and soil attributes, can be predicted to very fine resolutions at which remotely sensed data were measured or observed. The resulting inventories can be gauged by their relative accuracy, which is not always provided by conventionally produced land feature maps or inventories.
Traditional soil survey methods have in the past been criticised, perhaps justifiably, for being too qualitative in character. In response to these criticisms, quantitative models have been developed, especially within the last thirty years or so, which are being used to describe, classify and study the spatial distribution patterns of soil in a more objective way. These quantitative methods enable precise statements about the soil to be made. The methods are collectively categorised in the emerging field of soil science known as pedometrics. On the other hand, the new and emerging pedometric methods have also been critisised for being two expensive to adopt, as they require too many samples for them to be worthwhile, especially at the regional extent. This criticism is becoming untenable, especially with increasing availability of ancillary data, quantitative methods have proved to be equally efficient, if not better than, the conventional methods of soil inventory, even at the regional scale. This paper is aimed at providing examples to illustrate the hypothesis that “pedometric soil inventory methods are applicable at the regional extent.
The task of a soil survey is to provide soil information for either general purposes or for a specific use, with the latter task being dominant. In the past, soil scientists based their approach on the qualitative analysis of the landscape either by physiographic analysis or by aerial photographic interpretation or both. These methods were all used to enrich the soil information through the use of ancillary data. Due to increasing awareness of environmental pollution and associated problems by the general community, quantitative soil information is now required to enable more precise statements on the status of the environment to be made, even at the regional scale. Pedometric techniques have been developed to meet these requirements. Basically, these techniques stem from the classical approach when the soil scientist generally would study the climate, geology, geomorphology, vegetation, and land use and land use history, prior to any soil survey. Lately, technical advances have provided us with a wealth of new environmental data sources, which are summarised in Table 1.
Table 1. Examples of sources of exogenous data for soil inventory
spectral imageries (albedo)
DEM, crop growth, vegetation
Multi-spectral thematic mapper, etc
radiance energy and albedo
High resolution visible
radiance energy and visible spectral imageries (albedo)
DEM, landscape, vegetation and moisture
Advanced very high resolution radiometers(AVHRR)
surface reflectance, surface temperature
soil, vegetation, moisture (drought)
Proximal Sensors and Scanners:
Time domain reflectrometer (TDR)
apparent dielectric constant
existing ground-truth data and maps, etc
e.g., topographic, geologic, vegetation
Air-borne and space sensors are now generating terabytes of data on a daily basis. The air-borne sensors vary from the traditional aerial photography to air-borne videography (Table 2). The space sensors generally consist of series of satellites, many of which were originally designed and launched for general studies of the various earth resources. Increasingly, proximal sensors are now being used for the purpose of research and development for soil-specific crop management (McBratney and Pringle, 1997).
There are a variety of in situ soil-measuring systems, and they may determine variability internally or externally from the soil. For example, internal detectors include the neutron probe and time domain reflectometry (TDR), both providing an estimate of soil moisture. The soil can be analysed from an external, or remote position, and gamma-radiometrics (e.g Cook et al., 1996), radiation, or electromagnetic induction (EM) (Lesch et al., 1995) provide another quantitative measurement for use in pedometrics. There has been an attempt in recent years to obtain the quantitative data in a continuous fashion across the landscape, and at the catchment and continental extent, this may be provided by satellite imagery (Curran, 1998). At the local extent newly developed sensors attached to mobile machinery are providing continuous measurements of a field. Examples of this include the conductivity cart, which measures the EC of soil, and the ground-penetration radar. We now examine pedometric techniques that are used for the analysis of the sensed data for more efficient and inexpensive, quantitative soil inventories at the regional scale.
There are a variety of quantitative techniques available for analysing the spatial distribution of soil, with the most common methods in use at present being geostatistics, classical statistics and the combination of the two.
Generic pedometric techniques
Focussing on the basic types of pedometric techniques used for spatial prediction of soil, and hence in general sense soil survey, two generic techniques come to mind. The first, the classical approach, collectively referred to here as the CLORPT methods, and the second, geostatistical methods. To these we add the third group, the hybrid methods, which are some combinations of techniques from the two generic methods, above, in order to optimize prediction of soil properties. The techniques, and their styles are summarised in Fig. 1.
Figure 1. Pedometrics techniques
The CLORPT methods are based on the empirical-deterministic models that originated from Jenny’s (1941) “Factors of Soil Formation” (Fig. 3). Jenny’s (1941) state-factor equation is expressed as:
S = f(CL, O, R , P, (2)
where S is some soil properties as a function of the state factors: CL as climate, O as organisms, R as relief, P as parent material, and T as time. Soil spatial variability is therefore considered as being causative realisation of the complex combinations of soil-forming processes as influenced by the soil forming factors. The CLORPT function (Eq. 2) earlier in the 19th century stimulated numerous studies. Much of the earlier studies, and indeed some recent ones, were based on general and bivariate-simple linear regression (e.g. Furley, 1971; Moore et al., 1993) although multiple polynomial regression models were applied by Ruhe and Walker (1968). The realisation that many of these studies do not accommodate the non-linearity in the relations has led to recent application of the more robust methods such as generalised linear models (GLM), generalised additive models (GAM) and regression trees (RT) (Odeh et al., 1994, Gessler et al., 1995). Another developing method is the artificial neural network model, which is a non-parametric modelling technique, which mimic the neural networks of the brain (Venables and Ripley, 1994). The networks are composed of processing units, or neurons, which are organised into layers, i.e. input, hidden and output layers. Neural networks are “trained” by the data used to create them. Inputs are distributed from the input units into the hidden layer where the inputs are weighted and summed, a bias added and a fixed function of the result taken. The neural network fits iteratively, minimizing the differences between predictions and actual values during training (Venables and Ripley, 1994). But the question to ask is while the classical models or the more robust methods may take care of the deterministic relations, do they account for spatial autocorrelations of the soil properties, especially at the local level? To answer this question, the pioneer pedometricians initiated the application of geostatistics (which was primarily developed for the mining industry), the compendium of which is given below.
Geostatistical methods are based on the theory of regionalised variables (Matheron, 1965), which allows us to consider spatial variability of a soil property as a realisation of a random function represented by a stochastic model. The geostatistical method of spatial interpolation is termed kriging. The first major applications of ordinary kriging in soil studies in the early 1980’s (e.g., Burgess and Webster, 1980) was mainly at the field scale. Since then it has been extended to larger areas, upscaling to the regional extent (Odeh et al., 2000). Major limitations of the univariate geostatistics technique of kriging are due to the assumptions of stationarity which are not often met by the field-sampled data sets and, of course, the often cited requirement of large amount of data to define the spatial autocorrelation. However, as stated above and with increasing availability of ancillary information, the lack of adequate samples has been partially solved. The univariate usage of kriging is also limiting in situations of complex terrain where the soil-forming processes are themselves complex. In such situations there is the need to model both the structured and the spatially dependent components of the soil variable. Also there are economic and logistic reasons for including the ancillary variables influencing soil variability, especially if the latter are more readily and cheaply available. As both the soil and the exogenous factors are multivariate, the most obvious choices are appropriate combinations of multivariate/univariate analysis using the CLORPT factors and the geostatistical methods. These combinations constitute the hybrid techniques.
The hybrid (quantitative) techniques for soil survey are based on various combinations of the geostatistical and multivariate/univariate CLORPT methods. Let us suppose that a data vector describing a soil property is a random variable Z, determined at locations in a region, X = x1 ..., xN, and consisting of three components as
Z = m + Z1 + ε (3)
where m is the local mean for the region, Z1 is the spatially dependent component and ε the residual error term, spatially independent. Now there may be situations where m is varying and dependent on some exogenous factors such as the CLORPT factors. In other words it is deterministically related to some causative factor and, in geostatistics parlance, when the causative factor is related spatially, the variable is said to exhibit a trend. Wherever trend exists, the ordinary univariate kriging is inappropriate. Several methods have been designed to accommodate the trend.
(a) Universal kriging
Universal kriging (Matheron, 1969) has been the commonly used method to accommodate the trend or the “changing drift”, as it is sometimes known, in a soil variable. The universal kriging is a combination of the standard model of multiple-linear regression and the geostatistical method of ordinary kriging (Webster, 1994), which is also analogous to combining CLORPT methods with the univariate kriging using the geographical coordinates as determining the drift. More recently, a more advanced approach, the Intrinsic Random Function of Order k (IRF-k), has been used to accommodate the varying nature of the trend in a regionalised soil variable (McBratney et al., 1991). The term k represents the order of polynomial trends- k = 0 means constant drift, and the IRF-k is equivalent to ordinary kriging system of equations; if k = 1, we have linear drift; k = 2 yields quadratic drift. But where there is not trend but deterministic relationships are with some known or readily-available and inexpensive covariates (CLORPT factors) or other easy-to-measure soil variables, cokriging has played a major role in efficiently predicting the target soil variable (Stein et al., 1989; Odeh et al., 1995). Universal cokriging is also possible when considering the trend and covariation with one or more secondary variable (Stein et al., 1988).
Cokriging is the multivariate extension of kriging that allows the inclusion of more readily available and inexpensive attributes in the prediction process. There are many instances in soil survey where the CLORPT factors such as topography, time, variable parent material, etc, are easily discernible or are either readily available and/or are cheap to obtain. The most efficient way to predict the expensive-to-measure target soil variable, the variation of which is affected by the CLORPT factors, is to use the factors in cokriging the target soil variable- sampled at fewer locations, into dense grid nodes. This is termed as heterotopic cokriging (Wackernagel, 1994) in comparison with isotopic cokriging which requires that data on both the target variable and co-variables be available at all sample locations. A variant of both is the generalised cokriging (Myers, 1982) that involves simultaneously prediction of all the correlated variables into more dense locations. Heterotopy can either be complete or partial (Wackernagel, 1995). The complete case is the case where the covariates and the target variable do not share any common locations. A third is termed collocated cokriging, whereby covariates are available at all interpolation nodes, even though the target variable is available at only a few locations. This is often the case with using exogenous variables especially landform attributes derived from DEM (Odeh et al., 1995) and satellite imageries for predicting the target soil variable. The partial heterotopy involves cases where there are some coincidence of the locations of the target variable and the ancillary variables. The latter is often the case when other soil covariates are used (McBratney and Webster, 1983).
(c) Regression kriging
Regression-kriging (RK) is another hybrid method that combines either a simple or multiple-linear regression model (or a variant of GLM, GAM and regression trees) with ordinary, or simple, kriging of the regression residuals (Odeh et al., 1995; Goovaerts, 1997). The assumption here is that the deterministic component (m in equation 3) of the target (soil) variable is accounted for by the regression model, while the model residuals represent the spatially varying but dependent component Z1 in equation 3). If the exogenous variables used in the regression equation are available at more dense locations than the target variable, the equation can then be used to predict the m onto those locations. The Z1 can also be predicted to the same locations by ordinary kriging system of equations, and then added to the m to obtain Z*. Odeh et al. (1995) and Odeh and McBratney (2000) have demonstrated the superiority of RK to other prediction methods such as ordinary kriging, universal kriging, multiple-linear regression and cokriging. A variant of RK is kriging with uncertainty (Ahmed and DeMarsily, 1987). The kriging with uncertainty introduces the regression residuals (as representing model uncertainty) into the kriging system, which is then used to predict the target soil variable (Knotters et al., 1995). This reduces the extrema of the target variable and therefore produces smoother function of the predicted values. While Odeh et al. (1995) found kriging with uncertainty as not good as regression-kriging, nevertheless they reported it to be better than ordinary kriging or cokriging alone.
(d) Kriging with external drift
Kriging with external drift (KED) is a somewhat different hybrid technique, which integrates the universality conditions into the kriging system using one or more of the ancillary drift variables (Wackernagel, 1995). It is similar to universal kriging, but using an ancillary variable to represent the trend (Goovaerts, 1997). These variables could be digitised covariate derived from digital elevation model (DEM), or rainfall data or scanned images. The universality conditions need to be known not only at the sampled locations but also at the prediction locations. KED has not been widely used in soil science, but as remotely sensed data become more readily available, it may well be the method to be used along with regression-kriging (Odeh and McBratney, 2000), especially for soil inventory at the regional/catchment (200 m –2 km) resolution.
All of the quantitative methods described above are, by nature, data-hungry. This is why they are often quoted as being too inexpensive and inefficient to implement at the regional scale. However, are all of the generic and hybrid quantitative methods useful for regional soil inventory? We attempt to answer this question by providing a few examples and case studies.
Until the late 1980s when the Edgeroi area was intensively soil surveyed, there was lack of quantitative soil data any where in the cotton-growing region of N NSW. The survey was carried out using an equilateral triangular grid with a spacing of 2.8 km so that site placement was random with respect to landscape or landform (McGarry et al., 1989). The main purpose of the survey was to provide more accurate soil information for extension purpose, but has now found a more useful purpose in environmental modelling and for total catchment management. At each site the soil was characterised fully to a depth of 2.6m (where possible).
The soil characteristic information was first used to group the soil individuals at the sample locations into various classes at the suborder level of The Australian Soil Classification (Isbell, 1996). The AVHRR coverage of the area recorded on May 22, 1996, during the non-growing period was obtained to maximise the registration of soil exposure to the satellite radiometer, and thus provide optimal data for soil classification and prediction. The AVHRR data is at an approximate spatial resolution of 1 km. Before using the satellite images the data were kriged to a target grid of 200-m spacing (Odeh and McBratney, 2000). In addition, physiographical map of the area was digitised from the geologic and geomorphic map produced by William Ward (pers. comm..) onto the same grid matrix. The features used to predict the soil classes were NDVI, mid infrared (MIR), second thermal infrared band (TIR2), elevation (also interpolated to 200 m grid), profile curvature, plan curvature and the physiographic data. These features were used in a classification tree model to rapidly allocate the soil individuals onto the 200-m grid, based on soil classes obtained at 209 sites.
The resulting soil class map is shown in Fig. 2. As shown on the map, two soil types, the Black Vertosol and the Grey Vertosol, are predominant in the area occurring mainly in the central to western parts of the sub-catchment. To the east are mainly the Rudosols and the Kandosols, interspersed with the some Vertosols, Red Dermosols and Red Chromosols. Thus soil variability is much higher in this part than to the west of it. Physiologically there are two main landforms in the study area that are shaping the soil distribution patterns in the sub-catchment: the plains to the west, and the hilly region in the east. Looking at the map in Fig. 2 almost the whole of the plains are classified as Vertosols with few patches of Dermosols and Kandosols and Kurosols occurring in slightly elevated levees and prior stream formations (Triantafilis, 1997). A variable mantle related to the baseline geology underlies the hilly and undulating land to the east. The soil types in this part are a function of geology and topography, hence the high variability of soil types. Generally, the pattern of soil distribution is almost well aligned to the physiographical setting of the hills and valleys in the area. Even though, misclassification by our model was shown to be low (15%), the extent to which the classification tree model rapidly classify the soil into classes based on coincident of soil classes with physiographic unit, topography and AVHRR attributes, will need to be corroborated further.
Figure 2. The soil classes of the Edgeroi County of NW NSW, in accordance with Ibell (1997)
This section illustrates the fact that even without ancillary data, geostatistical interpolation techniques can be used at the regional level, especially if there is sufficient data in some localities within the region of interest. Due to high cost of getting accurate and quantitative information, much can be relied upon the available data at different scales. There is paucity of formal documentation of reliable, repeatable salinity assessment and spatial and temporal transfer of information through the scale continuum. This section focuses on the preliminary work done prior to and during the acquisition of quantitative information on soil salinity with the purpose of determining broad processes of salinisation in the region.
Figure 3. The variograms of a) Edgeroi subsoil (30-40 cm) % clay, and b) of soil salinity at 0.7-0.8 m and 1.2-1.3 m depth in the Macintyre region.
Upscaling (county - regional): The need to extend the survey to the greater cotton-growing region necessitates the use of the Edgeroi data for optimal upscaling. Percent clay content at 0.30-0.40 m depth for the Edgeroi area was selected as the primary variable for our sampling design at the regional scale. We started by fitting a spherical model (γ(h) =109.25 + 0.231h(km) up to the sill) to the experimental variogram from which we derived sample observation density and the kriging error for the primary variable using the resulting variogram (Fig. 3a) parameters. We determined the relationship between the density of sampling and kriging error for the % clay. Up to a certain level, decreasing the prediction error leads to improved accuracy of interpolated soil information as the sampling intensity increases. Obviously the accuracy of soil information improves with increasing density of observations that generate the information, and this could only be to a degree determined by the accuracy of the prior information used for the analysis. After taking into consideration the time, cost and logistic constraints we determined the sampling density of 0.354 observations per km2 for the lower Macintyre valley, much lower than the optimal (about 2.4 observations per. km2). A total of 120 sites were visited and sampled to 2.0 m depth. However, only ECe data at two depths (0.70-0.80 and 1.20-1.30 m) were used for this study.
Variography analysis at different scales: We first determined the variograms of the ECe data at both depths. We modeled the variogram based on the spherical models. The comparative variogram models are shown in Figure 3b. There is larger variability of subsoil ECe at 1.2-1.3 m depth than the shallower depths because the measured ECe at the regional scale encompasses more of the population variance than those at the county or the field levels. Evidently, there is a weak trend in the regional variability of subsoil (1.2-1.3 m depth) salinity. As shown in Figure 3, the variograms model for the ECe in the subsoil (1.2-1.3 m) consistently exhibit the larger spatial ranges and spatial variances than the upper subsoil (0.60-0.70 m). This is consistent with greater variation in salinity with increasing depth.
Figure 4. ECe map of the lower Macintyre region, as derived by geostatistical method of kriging
The nugget effect, representing the uncertainty of the kriging system, also varies at different depths and scales. The nugget generally increases with depth and increases as we upscale. Obviously, the results of the lower Macintyre translated into sub-optimal predictions of ECe at unsampled locations and at undersampled areas.
In Figure 4, there is a general increasing trend in salinity from east to west, with few patches of areas of high salinity, around 6-10 dSm-1, south-west of Goondiwindi and just north-east of Boomi. The most salt-affected areas are found around Mungindi at the southwest corner of Figure 4. The main cause of salinity distribution patterns in the lower Macintyre appears to be climate, with the salinity trend following the patterns of rainfall and temperature. In the east where the average annual rainfall is large, the soluble salts are probably flushed below 2 m depth. In the west where there is greater proportion of evaporation than rainfall, there is upward flux of soluble salts, which accumulate within 2 m depth.
This case study is based on the combined area of the Macintyre, Gwydir and Namoi valleys. The total area covered is approximately 45,600 km2. This region has been sampled extensively for a major study of irrigated cotton-growing regions of eastern Australia. The objective was to obtain a soil database so that a quantitative statement on the status of soil could be made. For this purpose, a total of 734 locations were visited and sampled. The AVHRR data combined with the elevation were used to predict the CEC over the entire region using the MLR, GAM and RK. First the AVHRR and elevation data were interpolated onto a grid of 500 m spacing for higher resolution target grid matrix for our modelling.
Figure 5 shows the map of CEC predicted by the RK model. The distribution patterns highlight the influence of physiography on the soil types – and hence the CEC. Clearly the light-textured soil on the levees are identified by relatively low CEC (CEC<150 mmol kg-1). These areas are illustrated by light grey patches on the map (Fig. 5). The dark grey patches represent areas of relatively large CEC values – mainly low depressions that are components of the clay plains mostly consisting of cracking clays (Stannard and Kelly, 1977).
Figure 5. The CEC map of the cotton-growing region of N NSW.
The predicted CEC values in some areas of the map are spurious to say the least (these are enclosed by the dotted lines). This is especially true for the southeastern corner of the map, the northern strip, and the areas between the valleys (Fig. 5). The uncertainty of prediction for the northern section of the map is less than that for the southeastern portion. This is due to the fact that the physiography of northern portion is more similar to the sampled areas along the valleys, and so extending the prediction model, using the AVHRR data, produced more realistic values. On the other hand the eastern section is the western fringe of the Great Dividing Range, which is undulating and hilly. The landforms are less similar to the clay plains. This accounts for low precision of prediction for the section, especially as the section was scantly sampled, or not sampled at all. Nevertheless, the map generally depicts a reasonably acceptable distribution of CEC for the region – adequate for regional planning and elucidation of areas of low CEC – that needs further investigation.
Application of each quantitative techniques depends on the purpose and the scale of the survey as the ultimate use of soil survey information determines the accuracy required. Many of methods are applicable at the regional extent. There are many examples in literature– e.g., De Gruijter et al. (1994), Collins et al. (2000), and Odeh and McBratney (2001). From all the examples given and new methods being generated, it appears that:
- Quantitative pedology (pedometrics) and the traditional pedology are growing closer together;
- There is a range of new data sources, which makes more plausible to extend pedometric techniques to larger area extents;
- Hybrid methods of CLORPT and geostatistics offer powerful spatial prediction methods at the regional extent.
Burgess, T. M. Webster, R. (1980) Optimal interpolation and isarithmic mapping of soil properties. I. The semivariogram and punctual kriging. Journal of Soil Science 31, 315-331.
Collins, M. McBratney, A. and Voltz, M. (2000). Developments in quantitative soil resource assessment (Pedometrics ’98). Geoderma 97: 147424.
Cook, S.E., Corner, R.J., Groves, P.R., Grealish,G.J. (1996) Use of gamma radiometric data for soil mapping. Australian Journal of Soil Research 34, 183-194.
Curran, P.J. (199). Remote sensing technologies for global ecological applications: an extended abstract. 9th Australasian Remote Sensing and Photogrammetric Conference, July 20-24 1998, University of New South Wales, Sydney, Australia.
de Gruijter, J.J., , R. Webster, K., and D.E Myers (Eds) (1994) Pedometrics-92: Developments in spatial statistics for soil science. Geoderma 62, 1-326.
Furley, P. A. (1971) Relationships between slope form and soil properties developed over chalk parent materials. In: Brusden D. (Ed), Institute of British Geograpghers Special Publication 3, 141-164.
Gessler P.E., Moore, I.D., McKenzie, N.J., Ryan, P.J. (1995) Soil landscape modelling and spatial prediction of soil attributes. International Journal of Geographical Information Systems 9, 421-432.
Goovaerts, P. (1997) Geostatistics for Natural Resources Evaluation. Oxford University Press, New York.
Isbell, R.F. (1996). The Australian Soil Classification. Australian Soil and Land Survey Handbook. CSIRO Publishing Melbourne.
Jenny, H. (1941) Factors of Soil Formation, a System of Quantitative Pedology. McGraw-Hill, New York, 281pp.
Knotters, M., Brus, D.J. and Oude Vishaar, J.H. (1995) A comparison of kriging combined with regression for spatial interpolation of horizon depth with censored observations. Geoderma, 63 :227-246.
Lesch.S.M., Strauss, D.J., Rhoades, J.D. (1995) Spatial prediction of soil salinity using electromagnetic induction techniques. 1. Statistical prediction models: a comparison of multiple linear regression and cokriging. Water Resources Research 31, 373-386.
Matheron, G. (1965) Les Variable Regionalisees et leur Estimation. Masson, Paris.
Matheron, G. (1969) Le krigeage universel. Cahiers du Centre de Morphologie Mathematique, no 1, Ecole des Mines de Paris, Fontainebleau.
McBratney, A.B. (1998) Some considerations on methods for spatially aggregating and disaggregating soil information. Nutrient Cycling in Agroecosystems 50, 51-62.
McBratney, A.B., Pringle, M.J. (1997) Spatial variability in soil- implications for precision agriculture. In: Stafford, J.V. (Ed.) Precision Agriculture 1997. BIOS Scientific Publ. Ltd. Oxford UK, pp 3-33.
McBratney A.B., Hart, G.A., McGarry, D. (1991) The use of region partitioning to improve the representation of geostatistically mapped soil attributes. Journal Soil Science 42, 513-532.
McGarry, D., Ward, W.T., McBratney, A.B. (1989) Soil Studies in the lower Namoi: Methods and Data. 1. The Edgeroi Data Set. CSIRO, Australia.
Moore, I.D., Gessler, P.E., Nielsen, G.A., Peterson, G.A. (1993) Soil attribute prediction using terrain analysis. Soil Science Society of America Journal 57, 443-452.
Myers, D.E. (1982) Matrix formulation of cokriging. Math Geol., 14: 249-257.
Odeh, I.O.A., and McBratney A.B (Eds.). (2001) Estimating Uncertainty in Soil Models (pedometrics ’99). Geoderma (in press).
Odeh, I.O.A., and McBratney A.B. (2000) Using AVHRR images for spatial prediction of clay content in the lower Namoi valley of eastern Australia. Geoderma 237-254..
Odeh, I. O. A., McBratney, A. B., Chittleborough, D. J. (1995) Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging. Geoderma 67, 215-225.
Odeh, I.O.A., Gessler, P.E., McKenzie, N.J. McBratney A.B. (1994) Using attributes derived from digital elevation models for spatial prediction of soil properties. Resource Technology’94, Melbourne Australia. pp 451-463
Ruhe, R. V., Walker, P. H. (1968) Hillslope models and soil formation: I. Open systems. Transactions of the 9th Congress, International Society Soil Science, Adelaide 4, 551-560.
Stein, A., van Dooremolen, W., Bouma, J., Bregt, A. K. (1988) Co-kriging point data on moisture deficit. Soil Science Society of America Journal 52, 1418-1423.
Stein. A, Bouma, J., Mulders, M.A. and Weterings, M.H.W. (1989) Using cokriging in variability studies to predict physical land qualities of a level terrace. Soil Techn. 2: 385-402.
Triantifilis, J. and McBratney, A. B. (1993) Application of continuous methods of soil classification and land suitability assessment in the lower Namoi valley. CSIRO Div. Soils Divisional Report No. 122. Canberra Australia.
Venables, W.N , Ripley, B.D. (1994) Modern Applied Statistics with S-Plus, Verlag, New York.
Wackernagel, H. (1988) Geostatistical techniques for interpreting multivariate spatial information. In: C.F. Chung et al. (Eds.). Quantitative Analysis of Mineral and Energy Resources. Reidel, Dordrecht, pp 393-409.
Wackernagel, H. (1994) Multivariate spatial statistics. Geoderma 62, 83-92.
Wackernagel, H. (1995) Multivariate Geostatistics. Sringer, Berlin, pp 256.
Ward, W.T., Soils and landscapes near Narrabri and Edgeroi, New South Wales, with data analysis using fuzzy k-means. CSIRO Division of Soils Divisional Report, (pers. comm.)
Webster, R. (1994) The development of pedometrics. Geoderma 62, 1-15.