Bulletin 609 February 1991
Alabama Agricultural Experiment Station Auburn University Lowell T. Frobish, Director Auburn University, Alabama

O*-

0

0. 0

A Simplified Theory of Point Kriging and Its Extension to Co-kriging and Sampling Optimization

CONTENTS
INTRODUCTION ...... ............. ............. KRIGING AS A LEAST SQUARES PROCEDURE .........
THE KRIGING EQUATIONS......................................4 EXTENSION TO CO-KRIGING...................................6

Page 3 4

ESTIMATION AND VALIDATION....
MOVING NEIGHBORHOOD KRIGING AND

...........

6

CROSS-VALIDATION...........................................6 ESTIMATION OF THE COVARIANCE STRUCTURE................7

DISCUSSION AND CONCLUSIONS... APPENDIX I...................................
DEFINITIONS AND NOTATIONS.................................9

............ ............

8 9

APPENDIX II......................................13
COMPUTER PROGRAMS............................................. AND GENERAL DESCRIPTION

13

INTRODUCTION MODEL

.................... 13
14

LIMITATIONS...............................................

PROGRAM DESCRIPTIONS.......................................... 14

Example ..................................................... 14 Developing and Using a Kriging Model for Water Content...............................................14 Developing and Using a Kriging Model for
Clay

Content.

..............................................
1991

28

Co-kriging Water with Clay ................................ 32
Finding an Optimal Sampling

Strategy....... ............ 35

FIRST PRINTING 3M, FEBRUARY

Information contained herein is available to all persons without regard to race, color, sex, or national origin.

A Simplified Theory Of Point Kriging And Its Extension To Co-Kriging And Sampling Optimization'
G.P.Y. Clarke and J.H. Dane
2 3

INTRODUCTION SIMPLIFIED explanation of the linear prediction technique of kriging was derived from a statistical point of view. The kriging technique allows values of a given, spatially dependent, variate to be predicted at points where no measurements were made. It is then possible to construct a contour map for that variate. Based on the theoretically developed equations, computer programs were written to carry out the predictions. Besides assisting the computer user, the aim of this research was also to point out the similarities between kriging and standard least squares, of which it is indeed a special case. Six computer programs were developed: (1) instead of determining a semi-variogram to specify the spatial interdependency of a given variate, a more general, cross-validation method was developed to determine the range parameter as needed in the kriging equations; (2) the selected model was then more thoroughly examined, e.g., for outliers, by yet another cross validation program; (3) in addition to simple kriging, equations and computer programs were developed for universal kriging, i.e. kriging under a lack of stationarity; (4) these equations and programs were subsequently extended to co-kriging to improve the estimation process by using information on auxiliary variables; (5) a program was written to validate the co-kriged model; (6) a procedure was developed to optimize sampling schemes with the use of kriging and co-kriging, i.e., the best sampling locations in the field and their minimum number according to preset criteria were determined.

A

'The authors wish to acknowledge the following organizations for their financial support culminating in this report: (1) The National Aeronautics and Space Administration for supporting the project "Spatial and Temporal Variability of Soil Temperature, Moisture and Surface Soil Properties" under agreement NAGW-758; (2) The Soil Conservation Service of the United States Department of Agriculture (SCS-A-4101-12-86); and (3) The Cooperative State Research Service of the United States Department of Agriculture through the Alabama Agricultural Experiment Station for their support of Regional Project S-185 "Spatial and Temporal Variability of Soil Characteristics and Material Fluxes in Field Soils." -Professor, Department of Statistics and Biometry, University of Natal, Pietermaritzburg,
Republic of South Africa.

'Professor, Department of Agronomy and Soils, Auburn University, Auburn, Alabama
36849-5412.

4

ALABAMA AGRICULTURAL EXPERIMENT STATION

The six developed computer programs were applied to data collected in the following experiment. A 50- by 100-m field was sampled on two occasions at a total of 60 locations and at five to 11 depths to study the spatial variability of pH, Ca, Mg, K, P, organic matter, texture, bulk density, soil water retention, and saturated hydraulic conductivity. The programs were written in such a manner that data sets obtained at two different times could be combined into one for prediction purposes. The equations developed in the theoretical sections are the same as those used in the computer programs presented in the Appendix, with most of the notation remaining consistent throughout.

KRIGING AS A LEAST SQUARES PROCEDURE The Kriging Equations Let Z denote a vector of n random variables with mean (column vector with n elements) and variance-covariance matrix V (size n x n). This can be written as Z
(,V).

I

It should be noted that Z denotes a vector with n observations if only one random variable is considered, as is often done with kriging. Now consider Z in partitioned form with vector Z, of dimension (n-i) and vector Z2 of dimension 1. Then

A Z1 1v172]
2 + v,12
and

where v,2 is a covariance row and v1 a covariance column vector. 2 The prediction problem is to predict Z2, denoted as Z2, given a set of observations ZI = z,. Assuming a linear model relationship among the means, the conditional mean of Z2, given Z, = z,, is:

E (Z2 I Z1 - z,) -

V11 (?1 - L1) ,

(1)

where z, is a column vector containing the observed values and E (Z2 Z, = z,) required Best Linear Unbiased Estimate (BLUE).

I

Now suppose

=

Dg

[1]

=

[d,]

i

(2)

KRIGING

5

where the elements of the n x m matrix D are known and g is a m- dimensional parameter vector. The matrix D is called the design matrix. For example, a particular case might be D = 1, a single column of ones, in which case we have stationarity in the mean. Other models for the mean reflect what is called "trend". The row vector, d', of dimension m, is that row of the design matrix which corresponds to the random variable Z2, which depends on the location at which an estimate will be made. The submatrix D1 is that part of the design matrix corresponding to the column vector of the observed values z1. It should be borne in mind that any linear mean model can be used. Assuming V to be known, the BLUE derived from equations 1 and 2 is the least squares estimate, namely: Z2 = d + v 2 V1 (z1 - D,) , (3A) where: = (D' V11 D1)-'D' V11 z, (3B)

1

_

is the estimate for 0. Equations 3A and 3B are the kriging equations. When the mean model is not stationary, a drift or trend needs to be taken into account. The process is now called "universal kriging." Note that, if a prediction, 22, is made at some point which coincides exactly with an existing observation, zj, then, the vectors of covariances for Z2 and Zj are equal. In other words, the column vector v12 belongs to the column space of V,,. Consequently:
v,':V =

is a vector of zeros except for 1 in place j, and
z2= d2 + (z
-

d2j)=

The above state of affairs, usually referred to as the prediction "honoring" the observation may not always be considered desirable. For example, it may produce rough looking contour maps. This can be avoided by postulating that v12 does not belong to V, 1 , which is equivalent to defining a non-zero variance, called the "nugget variance," for the difference between two observations on Z made at the same site. The classic reference to kriging is the book by Journel and Huijbregts, 4 where it is explained in a geological setting. The basis of their derivation of the kriging equations is to find a function of the observations which is unbiased in estimating the conditional mean and has
'Journel, A. and C. Huijbregts. 1978. Mining Geostatistics. Academic Press, New York.

6

ALABAMA AGRICULTURAL EXPERIMENT STATION

minimum variance. For that reason, their "kriging equations" may appear at first sight to differ from equation 3, but it is easy to prove their equivalence. Extension To Co-kriging Let

Z = ZA

where ZA, of dimension nA, denotes a vector of variates of primary of interest and Z,,, dimension n,, a vector of co-variates. Now suppose:
Z
V VB
1

(V4

It is apparent, that by proceeding exactly as for kriging, but with the vector A and matrix V defined as in equation 4, there is no extra difficulty in finding the least squares estimate of the conditional mean of any element of Z, using equation 3. In the simplest case, where there is stationarity of the mean, D will have two columns with each element 0 or 1 depending on whether that observed value of Z is a primary variate or a co-variate. It is worth mentioning that the formulation given in equation 4 allows observations to be made at any place on either the primary variate or the co-variate; they need not be measured together or equally frequently. It also can be seen that an extension to any number of co-variables is simple in theory. ESTIMATION AND VALIDATION Moving Neighborhood Kriging and Cross-validation In practice, when prediction is required at some point P, only those observations within a given distance of P are used in the prediction. This so-called "moving neighborhood," therefore, contains a specified subset of all the observations. It can be defined by requiring a minimum number of observations (Dubrule5 recommends about 15) or by spanning a given distance from P. The first advantage of such a strategy is that solution of the kriging equations involves a matrix of reasonably small dimensions. The second advantage is that pragmatic and sensible models can be fitted to the data. For example, it is often perfectly reasonable to fit a simple stationary mean model
'Dubrule, O. 1983. Two Methods with Different Objectives: Mathematical Geology Vol. 15(2):245-257. Splines and Kriging.

KRIGING

7

knowing that this in no way implies statiopnarity over the whole region. Perhaps one of the most important aspects of any analysis is the validation of the model used. We recommend cross-validation, which is a sequential procedure whereby each datum point in turn is removed from the set of observations and its value is then predicted by the model based on those observations remaining. A number of criteria can now be used to decide on the "goodness of fit" between predicted and observed values. However, the special spatial nature of the problem should not be ignored and it is often useful to portray the deviations graphically. This may, e.g., illustrate poor predictions on the edges of the map. Another useful aspect of cross-validation is its ability as an estimation technique in its own right. This matter is pursued in the next section. Estimation of the Covariance Structure The basic assumption usually made at this stage is that vj, the covariance between Zi and Z, depends only on their distance apart, d;j, and possibly on their orientation with respect to each other (e.g., Clark6 ; Journel and Huijbregts). A very simple model, and the only one pursued here, is the spherical model which specifies that the correlation between Z and Zj is 1 3 Pj = 1 - 2 (dj/a) + 2 (d;,/a) 3 , d, < a
PJ= 0 , d; >a

where the parameter a is called the "range." Points farther apart than the range are uncorrelated. The literature on spatial statistics invariably recommends that the nature of the spatial correlation should be investigated by means of a "semi-variogram." Suppose observations z; and z; have been made at points a distance d;; apart. Then, if all pairs of observations that are a distance d, apart are collected, and if there are n;; such pairs, it follows that:
-

1
2n,

i

where the summation is only over those pairs that are a distance d, apart. Now, if a value -j is calculated for every possible distance d;; in the data and j is plotted versus di, such a plot is called a semivariogram. Various models may be fitted through the generated data points and their parameters estimated.
'Clark, 1. 1979. Practical Geostatistics. Academic Science Publishers, London.

8

ALABAMA AGRICULTURAL EXPERIMENT STATION

The authors' experience with chemical and physical measurements on soil samples suggests that typically the semi-variogram exhibits poor structure. In addition, any effects of "drift" must be accounted for. These considerations have motivated the authors to first consider only the simplest model, namely, the spherical model, and second to investigate alternative methods for the estimation of the single covariance parameter a. Consider now the dual problem of estimating c in each given moving neighborhood and then combining these estimates to form one global value. Initially there is the additional problem of defining the collection of moving neighborhoods for the entire procedure. One elementary definition, and the only one pursued in this publication, is to have one moving neighborhood centered on each observation. Given an assumed value for a and a specific neighborhood with n, observations z, possible cross-validation estimators will be considered. Let z; = observed value of Zj at site j and 2 = predicted value of Z at site j for a given value of a. Note that this predicted value must be determined from the data set which excludes the actual observation at site j, so it is a cross-validation estimate. For a given value of a Ej = Zj -Zj for j = 1 to n. can be defined. The choice of an estimator for a, based on the crossvalidated residuals, is difficult to justify on the grounds of a single criterion in all situations. Among the desirable criteria to consider are the following: (1) minimum E , (2) minimum inter-quartile range (E), (3) minimum inter-decile range (E), and (4) minimum rank mean (E). For example, to define (4), suppose that for each site j E; is evaluated for Lvalues of a (a = a, a2 , ..... , a_), i.e., E;(a,),.......,E,(aL). Rank the absolute values of the E (a) and then calculate the means of these ranks. DISCUSSION AND CONCLUSIONS There has been a considerable upsurge of interest in spatial statistics in soil science over the past decade. This reflects a growing need by the scientific community to take into account the non-random

KRIGING

9

nature of the distributions of variables measured at different locations. One tool which can meet many of these needs is kriging. It is imperative, however, that the entire process of model building should be seen as a whole. One aspect of the statistical model is the mean structure and the other the covariance structure. In the opinion of the authors there has been an over-emphasis on estimation of the latter. This has led to a number of studies which dealt solely with semi-variograms. In this bulletin an attempt has been made to focus attention on the predictive part of the process. How to fit a model, how to validate it, and how to portray the results. Another aspect of the publication deals with the theoretical explanation from a least-squares point of view, because this is familiar to the authors. APPENDIX I DEFINITIONS AND NOTATION
DRIFT MODEL

The mean value of a spatial variate Z at some location identifed by coordinates (x,y) may be written as mean(Z) = j(x,y). This is the drift model. Typical examples are: (stationary model), Mean(Z) = o30 (linear drift model), Mean(Z) = Rio + /, x + 32 y Mean(Z) = Xo + ,3x + 32y + 112 + 322y2 + 31xy x 2 (quadratic drift model) where o, O,, /32, ...are unknown parameters, denoted by the vector 3.
DESIGN MATRIX

This matrix, denoted by D, contains the coefficients of the drift model for all sample points. If the drift model is stationary, D consists of a single column of ones (m = 1). If the drift model is linear, the ih row of D would contain the following 3 (m= 3) elements: { 1,x,y,}. If the drift model is quadratic, the ith row of D would contain the following 6 (m =6) elements: {l,xi,yi,xi, Yxi Yi }" The mean vector of Z is thus denoted as = D_.
MOVING NEIGHBORHOOD

In the kriging procedure, when a prediction is required at coordinates (x*,y*), this prediction is made (using kriging or co-kriging), based on the nearest ng sample points to the given point. This subset of ng points defines the neighborhood of the given point. As kriging pro-

10

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

ceeds to predict from one point to another, the neighborhood "moves." This procedure is analogous to the "windows" used in calculating moving averages in one dimensional time series.
NEIGHBORHOOD SIZE

Given some point with coordinates (x*,y*) we define a "neighborhood" of that point as the smallest subset of all observations Z,,i = 1 to n, such that: (1) it contains at least the closest ng observations, and (2) the distance from (x*,y*) to the most distant member is at least dg. In most applications it is safe to set dg = 0, so that neighborhoods depend on numbers of members only. A reasonable value for ng is 12 to 15.
PREI)ICTEDI) VALUE

The BLUE of Z, namely the "kriged" estimate is given in equation 3 (page 5) which is repeated here: 22 = d 2' j + v,' VI (z, - D,f3) S= (D,' V,2 D,) 'D,' V,1 z,
PREDIICTION STANDARD ERROR

If Z denotes the predicted values of Z at location i (either predicted by kriging or co-kriging), and if Z, is the true unknown value of Z at this point, then the prediction SE is our estimate of: (Variance (Z 2))':. The prediction standard error is a measure of the reliability of the prediction and can be used in tests of significance and probability statements. It is estimated by the following procedure: Since Z; is a linear function of the observations z, we may write it as:

and it may be verified that
c' = ((d,' - v,', VI D,) HD,' + v,',) V,

where
H = (D,' V,, D,) '.

Hence it follows that:
Var (Z, Zi) =
2

(1 -

2v,'2c + c' V, c).

To get the prediction SE, a- must be estimated and the best quadratic estimate is given by:

KRIGING

KRIGING

11 11

a2 = e'Viie/(ng - m)
where: e = Z - Dg.
NEIGHBORHOOD VARIANCE

This quantity, referred to as a2 above, measures the conditional variance at any point about the local mean. Within a given neighborhood we expect it to be constant at all points. However, variation from one neighborhood to another is expected.
LOCAL MEAN AND LOCAL SD

In a given neighborhood, the mean value of Z at the point (x*,y*), as defined by the drift model, is called the local mean. Within the entire given neighborhood, a constant conditional variance is assumed, being that variance of the variate Z about its local mean. The local SD is the square root of this variance. Naturally, the local mean and local SD may vary from one moving neighborhood to another. In other words, at any point (x*,y*) there is a corresponding row d*'of the design matrix. For example, if (x*,y*) = (4,5) and we fit a linear drift model, then: d*' = {1 4 5} and the local mean at this point is given by = d* ' p. Thus, e.g., if the mean model exhibits linear drift, all means in the neighborhood lie on the same plane. In the absence of spatial correlation, the local mean estimate and the kriged estimate coincide.

I

COVARIANCE MODEL

Given a random variable Zi at location (xi,y) and another Z;, at location (x3,y;), and if dijis the Euclidean distance between these two locations, we might propose a model to describe the covariance of Zi and Z; as a function of di only. This is a simple covariance model (said to be anisotropic).
SPHERICAL COVARIANCE MODEL

If cov(Zi,Z) = 1 - 0.5(di;/a) + 1.5(di/a)3 for di; < a and zero otherwise, this is the spherical covariance model and the parameter a is called the "range."
PROBABILITY STATEMENTS

On the assumption of normality, the probability that a value of Z at (x*,y*) should exceed some given quantity G, may be estimated by: w ((G
-

d ) /

),

where 1 denotes the cumulative standardized normal distribution.

12

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

12 x z

ALABAMA AGRICULTRLEPRMN

TTO

rhe three columns of (x,y)-coordinates and observed values z of the variable of primary interest (each of length n).
n

The number of observed points for the variable of primary interest.
d

The design matrix for the variable of primary interest (of dimension nx m).
m

The number of columns of the design matrix D.
ng

The number of points in each moving neighborhood.
k

The smoothing parameter for kriging. When k = 0 the predictions exactly honor the observations. When k = 1, maximum smoothing takes place and a predicted value equals the local mean. Note 0 < k < 1.
a

The value of the range parameter in the spherical covariance model.
x* y*

The (x,y) coordinates for the prediction points, _.
z

The column of kriged or co-kriged predicted values with ith member
z!.

D*

The design matrix for the prediction points.
as,ad,an

When running program AX, as = starting value for the range parameter, ad = increment by which the value of the range parameter is increased in the sequence; an = total number of values of the range parameter to be tried in sequence. The column of cross-validated predicted values. Thus is the predicted value at location i when the observed value zi is itself deleted from the observations.
E

The n x an matrix of cross-validated residuals from program AX. Thus if the current trial value of the range parameter is a; and this results in a predicted value _Ci location i during cross-validations, for
then ei
e
=

Zij

-

ci

·

The single column of cross-validated residuals produced by fitting a

KRIGING

13

single specified model by kriging or co-kriging. Note that if program AX was used to fit a single value of the range parameter only, then the matrix E reduces to this single column vector.
Xbybzb

The (x,y)-coordinates and observed values of the covariate used in co-kriging.
db

The design matrix for the covariate points.
ah

The value of the range parameter for the covariate.
mh

The number of columns in the covariate design matrix.
ng b

The number of points in each moving neighborhood for the covariate.
U
2

The average local variance for the primary variate. The average local variance for the covariate. P The correlation between primary variate and covariate. APPENDIX II COMPUTER PROGRAMS'
INTRODUCTION AND GENERAL DESCRIPTION

The suite of computer programs presented in this APPENDIX has been written to carry out point kriging, incorporating nonstationary means (universal kriging), and co-kriging. In addition, there is a program SAM, which searches for the optimal locations of sample sites. The language used is GAUSS. All programs, and files as used in the examples, are stored on the diskette. Model fitting is done by cross-validation. The programs are not intended to be used as prescriptions, but require interaction and understanding from the user. It is to be expected that a user will interface analyses available here with software from other sources. For example, to use the kriging program KRIG, the user must supply details of his/her chosen drift model and covariance parameter. Although these details can be found using program AX, some people might prefer to use a more conventional procedure based on variograms.
'Copies of the computer programs and the data sets as used in this bulletin will be provided on a floppy disk upon request from J.H. Dane, Department of Agronomy and Soils, Auburn University, Auburn, Ala. 36849-5412.

14

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

Likewise, once the program KRIG has produced a grid of predicted values, the graphing of contours can be accomplished with any of the user's favorite graphical packages. No graphics procedures are given here. Perhaps a user may be interested in testing a model produced from some other source. The program CROSSK can be used to validate such a model. Although users of these programs need not have knowledge of GAUSS, such knowledge makes the file manipulation and rearrangement of data infinitely easier and is strongly recommended.
MODEL LIMITATIONS

Any drift model which is linear in its unknown parameters can be fitted. Examples are quadratic drift and multiple regression on any co-variates whose values are known at the prediction points. Where the covariate values have themselves to be estimated at the prediction points, it becomes a co-kriging problem. The only covariance model which is fitted is the spherical model. This model is expanded to handle co-kriging without introducing extra parameters.
PROGRAM DESCRIPTIONS

There are six main programs, which together with supporting procedures, are stored in a library called SPATIAL. Descriptions of each program are given in Appendix tables 1 and 2. The programs all require input to have been stored previously in files with designated names and output is placed in named files.
Example

Water and clay content were recorded at 60 locations in a rectangular field of 50 m X 100 m. Appendix figure 1 shows the layout of the sample points and Appendix table 3 gives the raw data with their (x,y)-coordinates. The purpose of the analysis is to produce a contour map of predicted water content values. Both kriging and cokriging will be used. The sequence of analytic steps is shown in Appendix figure 2. A detailed description now follows. Note that all the input files needed to run the examples are stored on diskette. Appendix table 4 lists their file names. Developing and Using a Kriging Model for Water Content
Using Program AX to Find the Best Water Model

When fitting a model to spatial data there are two aspects to consider jointly. First, what model fits the drift or mean and second, what model fits the covariance structure. We adopted the simplifying stand

Appendix Tabl le 1. Description of computer programs for kriging. Name: Purpose: AX To evaluate a sequence of values for the range parameter, a, given a specific drift model. Cross-validation is used to define various "goodness-of-fit" statistics for each value of a. CROSSK KRIG To produce cross-validated predicted To calculate predicted values by universal kriging, given the range parameter a and a values based on universal kriging. specific drift model. Moving neighborhoods may be used. 1. STORE a, ng, m, k in the file TEMPPAR. 2. Store [x y z D] in the file TEMPAIN. 3. Store [x* y* D*] in the file TEMPCOR. Execution: Output: Run AX. On the file TEMP_OUT:
1. The matrix E = {e}, i = 1, n andj

Preparatory 1. Store as, ad, an, ng, m in the file TEMP_PAR. steps: 2. Store [x y z D] in the file TEMPAIN.

1. Store a, ng, m, k in the file TEMPPAR. 2. Store [x y z D] in the file TEMP_AIN.

Run KRIG. On the file TEMPOUT: For each prediction point (x*, y*), the predicted value, Z, the prediction SE, the local conditional SD, and the local mean.

Run CROSSK. On the file TEMP-OUT: For each observation point (x,y), the crossvalidation predicted value Z a, the prediction SE, the local SD, the local mean.

= 1, an, consists of cross-validation residuals. 2. The "goodness-of-fit" statistics: rank means, mean sum of squared deviations, mean absolute deviations, inter-quartile range, inter-decile range.

Name: Purpose:

Appendix Table 2. Description of computer program for co-kriging and optimal sampling. COK SAM CROSSC To find an optimal subset of the sample To calculate predicted values by coTo produce cross-validated predicted kriging. Any form of drift models may points from a given larger set of potential values based on co-kriging. sample points. A co-kriging model is used be specified for the primary and coand therefore a set of covariates is revariates. The two variates need not have quired. the same drift or co-variance models.
j 2 1. Store a, ab,P , , a,,m,mb,ng,ngb, k in the file TEMPPAR.

Preparatory 1. Store a, ab, P, 2, ob,m,mb,ng,ngb, k in steps: the file TEMP-PAR. 2. Store [x y z D] in the file TEMP_AIN.
3. Store [xb yb Zb Db] in the file

1. STORE a, ab,P,a , Ub,m,mb,ng,ngb,

2

maxout, in the file TEMP_PAR. 2. Store [x y D] where, these entries refer to potential sample points, TEMP_AIN.
3. Store [xb Yb Zb Db] in the file

2. Store [x y z D] in the file TEMPAIN.
3. Store [xb Yb Zb Db] in the file

TEMP-BIN. 4. Store [x* y* D*] in the file TEMP_COR. Execution: Run COK.

TEMP_BIN.

TEMP_BIN. Store [x* y* D*], where the entries refer to the prediction points in file TEMPCOR.

Run CROSSC.

Run SAM. Note that this puts the program into an inter-active mode and you will now have to respond to the question asking for a value for parameter ns.

Output:

On the file TEMPOUT: For each prediction point (x*, y*), the predicted value (co-kriged), z, the prediction SE, the local conditional SD (for the primary variate), the local mean (for the primary variate). 1In

The sequence of points added to and subtracted from the optimal subset. The coorFor each observation point (x,y), the cross- dinates of the best subset with their validation predicted value zc, the predicaverage prediction SE's and those points tion SE, the local conditional SD, the local discarded and also the prediction SE's for mean. the prediction points. On the file TEMPOUT:

KRIGING

KRIGING
Appendix Table 3. Water and clay content at 80 cm depth at 60 locations. Water
Xcontent 3 3

17

1

m 0.00 6.25 12.50 18.75 25.00 37.50 10.93 17.18 23.43 29.68 6.25 34.37 17.18 12.50 25.00 7.81 3.12 4.60 0.00 12.50 18.75
25.00

m 0.00 0.00 0.00 0.00 0.00 0.00 3.12 3.12 3.12 3.12 6.25 6.25 9.37 12.50 12.50 15.62 18.75 21.87 25.00 25.00 25.00
25.00

cm /cm 0.220 0.235 0.243 0.225 0.243 0.241 0.203 0.232 0.229 0.240 0.230 0.258 0.246 0.249 0.262 0.243 0.243 0.265 0.255 0.270 0.274
**

31.25 37.50 4.68 10.93 3.12 9.37 12.50 25.00 3.12 31.25 7.81 14.06 20.31 0.00 12.50 25.00 31.25 37.50 7.81 14.06 3.12 18.75 21.87 12.50 15.62 10.93 3.12 18.75 21.87 34.37 1.56 7.81 14.06 29.68 0.00 12.50 25.00 37.50

25.00 25.00 28.12 28.12 31.25 37.50 27.50 37.50 43.75 43.75 46.87 46.87 46.87 50.00 50.00 50.00 50.00 50.00 53.12 53.12 56.25 56.25 56.25 62.50 62.50 65.62 68.75 68.75 68.75 68.75 71.87 71.87 71.87 71.87 75.00 5.00 75.00 75.00

"missinfg

I II

0.261 0.256 0.236 0.276 0.280 0.301 0.295 0.278 0.296 0.242 0.294 0.287 0.290 0.266 0.306 0.268 0.274 0.238 0.302 0.313 0.299 0.269 0.307 0.291 0.293 0.267 0.292 0.302 0.289 0.234 0.286 0.303 0.299 0.285 0.295 0.305 0.280 0.255

Clay content % 15.4 19.2 17.6 17.7 17.4 17.7 13.4 18.4 17.3 17.1 18.1 19.4 21.5 20.7 25.0 19.5 21.5 19.4 19.2 22.0 24.0 29.0 23.0 22.5 15.3 21.4 26.6 29.9 28.0 24.8 28.9
*

24.7 27.0 27.6 24.2 32.0 26.0 20.7 24.4 25.5 26.8 26.0 23.6 32.9 22.8 24.3 22.0 28.2 25.6 22.7 16.0 20.7 23.5 26.5 29.5 26.9 25.7 20.1 19.8

data

18
v

ALABAMA AGRICULTURAL EXPERIMENT STATION

that the spatial correlations may be modelled by the spherical model with a single parameter, a~, the range. Program AX will be executed three times; first for a stationary drift, drift and finally for a quadratic drift model. In each then for a case, a range of values for a~ between. 4 m and 13 m will be used, because this represents slightly more than the minimum. distance between sample sites and the distance needed to ensure spanning a neighborhood of 20 points respectively.

linear
.5

Fitting the Stationary Model

(1) Store 4

19 20 0 on file

Note that the first 3 numbers, as
.--

TEMP-PAR.
=

4, ad

=

.5,
1

and an
--

=

19,

-7--------,

\--------%

-

,

1~

5

%

;I1----- -'---%%% %

-----%% %

%-I1
/ %

%4--- -- D

-----------

----

V

%

%

I

%

/t1%

~/
50X

-PEDXFGR

tigus tetw

ies that

i f smlelctin .Dstribution /ampled hildwas

h fed.Tecicesadsqaesds

KRIGING

KRIGING1

19

Footnote : Shaded boxes indicate steps which do not involve SOILKRIG programs.

[Clear

boxes involve SOILKRIG programs.

APPENDIX FIGURE 2. A schematic diagram to illustrate the steps in kriging water content and

clay content, respectively, and co-kriging water with clay. Shaded boxes indicate steps which do not involive SOILKRIG programs. Clear boxes involve SOILKRIG programs.

Appendix Table 4. The location of input and output computer files used in or produced by the example. Input files TEMP Page 18 21 21 24 25 Type of Analysis Using AX to fit a stationary drift model through the water content data. Using AX to fit a linear drift model through the water content data. Using AX to fit a quadratic drift model through the water content data. Using Crossk to validate the linear drift model for the water data. Using KRIG to predict water content. PAR
TEMPAIN TEMPCOR

TEMPOUT
WATOUT.1 WATOUT.2 WATOUT.3 WAT._OUT.4

Output files

WATPAR.1 WATPAR.2 WATPAR.3 WATPAR.4 WATPAR.4

Location of Files on Diskette WATAIN.1 WATAIN.2 WAT AIN.3 WATAIN.2 WATAIN.2 WATCOR. 1

WATOUT.5

28 28 28 32 32 33 3535

Using AX to fit a stationary model through the clay content data. Using AX to fit a linear drift model through the clay content data. Using AX to fit a quadratic drift model through the clay content data. Using CROSSK to validate the linear drift model for the clay data. Using KRIG to predict clay content. Using COK to predict water content. Using CROSSC to validate the co-kriging model. Using SAM to optimize sampling.

WAT

PAR.1

CLAAIN.1 CLAAIN.2 CLAAIN.3 CLAAIN.2

CLAOUT. 1 CLA-OUT.2 CLAOUT.3 CLAOUT.4 CLAOUT.5 WCOUT.2 WCUOUT. 1 SAMOUT. 1

WATPAR.2 WATPAR.3 CLAPAR.4

V

WAThCOR.1 CLAAIN.2 CLAPAR.4 Required input files Location WATCOR. 1 TEMPCOR TEMPYPAR WC-PAR.1 TEMPAIN WATAIN.2 CLAAIN.2 TEMPBIN TEMPYPAR SAM-PAR. I TEMP AIN SAMAIN.l TEMP-BIN CLA-AIN.2 TEMPIUCOR SAMCOR. 1

KRIGING

21

define the range of values for a, and the next two numbers, ng = 20 and m = 0, specify the neighborhood size (a good default value) and a stationary model, respectively. (2) Store the matrix of values [x y z], where z represents water content, on file TEMPAIN. There is no need to include a design matrix since the model is stationary. (3) Run the program by entering: RUN AX. When the execution is complete, examine the output on file TEMP OUT. We postpone discussion of the output until all three drift models have been fitted.
Fitting the Linear Drift Model

(1) Amend the parameter values in file TEMP PAR by changing the value of m from 0 to 3, because the design matrix will now contain 3 columns. (2) Store the array [x y z D] on file TEMPAIN. Appendix table 5 contains a partial listing of this required input. Note that in the GAUSS language we would simply calculate. S = x y z ones (n,l)~ x~y and store S in TEMPAIN. (3) Run the program by entering: RUN AX. The output is filed on the data file TEMP OUT.
Fitting the Quadratic Drift Model

(1) Change the value of m to 6 in TEMP PAR. (2) Calculate the elements of the quadratic design matrix D and store [x y z D] on file TEMP__AIN. Note that the GAUSS language command is S = S" (x.*x) ~ (y.*y) "" (x.*y) assuming the old array S is retrieved first. This is now stored on TEMPAIN. (3) Run the program by entering: RUN AX. The output is filed on data file TEMP OUT. Interpreting the Output There are two parts to the output from program AX. The first part comprises the matrix of cross-validated residuals E
=

{ei3}, i = 1 to 59,

j = 1 to 19

where e = the residual value for location i using the jth a-value.

22

22

ALABAMA

ALABAMIA AGRICULTUA

AGRICULTURAL

EXPERIMENT

XEIMN

STATION

TTO

This matrix may be useful for extra analysis; for example in graphical representation. The second part gives measures of spread for the residual values, calculated separately for each a-value. The program AX calculates the following specific measures of goodness of fit for each of the 19 values for a considered: (1) Rank means. For each row of B, allocate a rank to each column element, thereby generating a n x p matrix R of ranks. Calculate the column means.
Appendix Table 5. The input file TEMPAIN required for fitting a linear drift model to water data. water Y x content Y X m m cm 3/cm 3 m r 000.000 000.000 1.000 0.221 0.000 0.000 000.000 12.500 1.000 0.235 0.000 6.250 000.000 12.500 1.000 0.243 0.000 12.500 000.000 18.750 1.000 0.225 0.000 18.750 000.000 25.000 1.000 0.244 0.000 25.000 000.000 37.500 1.000 0.241 0.000 37.500 3.125 10.938 1.000 0.204 2.125 10.938 3.125 17.188 1.000 0.233 3.125 17.188 3.125 23.438 1.000 0.229 3.125 23.438 3.125 29.688 1.000 0.240 3.125 29.688 6.250 6.250 1.000 0.230 6.250 6.250 6.250 34.375 1.000 0.259 6.250 34.375 9.375 17.188 1.000 0.246 9.375 17.188 12.500 12.500 1.000 0.250 12.500 12.500 12.500 25.000 1.000 0.262 12.500 25.000 15.625 7.813 1.000 0.243 15.625 7.813 18.750 3.120 1.000 0.243 18.750 3.125 21.875 4.688 1.000 0.266 21.875 4.688 25.000 0.000 1.000 0.256 25.000 0.000

14.063 3.125 18.750 21.875 15.625 10.938 3.125 18.750 21.875 34.375 1.563
7.813

12.500

53.125 56.250 56.250 56.250
62.500

0.314 0.299 0.269 0.308
0.291

1.000 1.000 1.000 1.000
1.000

14.063 3.125 18.750 21.875
12.500

53.125 56.250 56.250 56.250
62.500

62.500 65.625 68.750 68.750 68.750 68.750 71.875
71.785

0.294 0.267 0.292 0.302 0.289 0.234 0.287
0.303

1.000 1.000 1.000 1.000 1.000 1.000 1.000
1.000

15.625 10.938 3.125 18.750 21.875 34.375 1.563 14.063 29.688
0.000

14.063 29.688
0.000

71.875 71.875
75.000

0.300 0.286

12.500 25.000 37.500

75.000 75.000 75.000

0.2960.305 0.280 0.255

1.000 1.000
1.000

7.813

62.500 65.625 68.750 68.750 68.750 68.750 71.875
71.875

71.875 71.875
75.000

1.000 1.000 1.000

12.500 25.000 37.500

75.000 75.000 75.000

KRIGING

KRIGING

23

2

The remaining measures are standard, and each is calculated on every column of B separately: (2) Mean sum of squares. (3) Mean absolute value. (4) Inter-quartile range (IQR). This is the difference between the 25th and 75th percentile. (5) Inter-decile range (IDR). This is the difference between the 90th and 10th percentile. Output data for the different measures of goodness of fit (second part of program AX) are presented in Appendix table 6 when fitting a linear model. Observe from this table that an a-value between 6
« m 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5
9.0

Appendix Table 6. Output of program AX when fitting to water content. Rank Mean Mean means SS(resid) Abs(resid) 3 3 (cm /cm)2 cm 3/cm 3 8.20339 0.00020 0.01163 8.25424 0.00020 0.01160 8.37288 0.00021 0.01161 8.37288 0.00021 0.01164 8.57627 0.00021 0.01165 8.71186 0.00021 0.01166 8.44068 0.00022 0.01169 8.40678 0.00022 0.01178 8.50847 0.00023 0.01192 8.57627 0.00024 0.01215 8.98305 0.00026 0.01249 9.42373 0.00027 0.01286 10.11864 0.00028 0.01325 10.98305 0.00029 0.01362 11.91525 0.00031 0.01401 12.79661 0.00032 0.01438 13.52542 0.00032 0.01464 13.89831 0.00033 0.01477 13.93220 0.00033 0.01480 Percentage of worst value 58.88 61.63 78.57 59.25 61.91 78.41 60.10 62.40 78.47 60.10 63.07 78.62 61.56 63.89 78.73 62.53 64.94 78.82 60.58 66.35 79.02 60.34 68.24 79.60 61.07 70.76 80.56 61.56 74.00 82.11
64.84 77.86 84.36

a linear drift model I.Q.R. cm /cm 0.01514 0.01439 0.01436 0.01450 0.01415 0.01439 0.01504 0.01597 0.01665 0.01729 0.01732 0.01783 0.01891 0.01992 0.02076 .0.02125 0.02182 0.02320 0.02454 61.69 58.63 58.52 59.11 57.66 58.65 61.29 65.08 67.87 70.46
3 3

I.D.R. cm 3/cm 3 0.03749 0.03741 0.03779 0.03747 0.03748 0.03739 0.03697 0.03666 0.03670 0.03698 0.3760 0.03972 0.04163 0.04312 0.04420 0.04491 0.04527 0.04484 0.04421 82.81 82.63 83.49 82.76 82.80: 82.59 81.66 80.98 81.07 81.69
83.06

9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0

67.64 72.63 78.83 85.52 91.85 97.08 99.76 100.00

81.91 85.85 89.53 93.01 96.26 98.61 99.85 100.00

86.87 89.53 92.05 94.70 97.18 98.91 99.82 100.00

72.67 77.07 81.19 84.62 86.61 88.84 94.54 100.00

70.59'

87.75 91.96 95.26 97.63 99.21 100.00 99.05 97.66

24

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

and 8.5 will optimize the IQR and IDR, which are preferred measures of spread, since they are relatively robust to outliers. In order to compare all three drift models simultaneously, the three outputs from program AX must be examined. Appendix table 7 summarizes the results for the mean sum of squared deviations and interquartile range. The linear drift model is clearly superior to the other two. Given the linear drift model, a reasonable compromise value for a is 6.0. Appendix figure 3 gives a graphical description of variations in measures of spread for different a-values. This confirms a value for a of 6.0 as being reasonable. It is noteworthy that no single choice of parameter values is best by every criterion.
Using Program CROSSK to Validate the Given Water Model

When a tentative model to describe spatial interdependency has been adopted (e.g. by using the program AX), it can be more thoroughly examined by the use of cross-validation. The process of cross-validation for validating any given model is no different from the procedure used by program AX in the search for a best model. The only real advantage from running program CROSSK in this example is to obtain a more detailed output. Program CROSSK uses a specified model and

APPENDIX FIGURE 3. Measures of deviation (as

o of worst value) of cross-validated residuals using a linear drift model on soil water content.

KRIGING

KRIGING25 Appendix Table 7. Measures of deviation of cross-validated residuals for three different drift models on water content. Sum of squared deviations Inter-quartile range (cm 3/cm 3 )2 (cm 3/cm3) S LQS L 0.0176 0.0120 0.0132 0.0219 0.0151 0.0175 0.0120 0.0135 0.0207 0.0144 0.0173 0.0121 0.0137 0.0203 0.0144 0.0172 0.0122 0.0141 0.0207 0.0145 0.0171 0.0124 0.0145 0.0211 0.0141 0.0169 0,0126 0.0149 0.0220 0.0144 0.0169 0.0155 0.0215 0.0150 0.0169 0.0132 0.0160 0.0197 0.0160 0.0171 0.0137 0.0167 0.0192 0.0167 0.0175 0.0143 0.0175 0.0191 0.0173 0.0181 0.0151 0.0184 0.0192 0.0173 0.0187 0.0159 0.0194 0.0189 0.0178 0.0194 0.0166 0.0202 0.0195 0.0189

25

a m 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0

-0.0129
0.0173

Q 0.0182 0.0183 0.0184 0.0185 0.0184 0.0181
0.0175

0.0178 0.0192 0.0199 0.0206 0.0217 0.0226

10.5

0.0199

0.0210

0.0201

0.0199

0.0233

11.0 0.0205 0.0181 0.0218 0.0211 0.0208 0.0236 11.5 0.0211 0.0186 0.0225 0.0227 0.0212 0.0241 12.0 0.0214 0.0191 0.0230 0.0242 0.0218 0.0245 12.5 0.0216 0.0194 0.0232 0.0251 0.0232 0.0247 13.0 0.0215 0.0194 0.0231 0.0256 0.0245 0.0242 S = stationary drift model, L = linear drift model, Q = quadratic drift model.

makes cross-validation predictions at each observed point. Output of the observed values (Z), the predicted values (2), the prediction standard errors, the local means and the local standard deviations may be useful in the following ways: (1) To test for suitability of the model and detection of outliers. The ratio of (Z - 2) /(prediction SE), at any point, should follow the t-distribution. Any values in excess of about 3 can be seriously considered as outliers. (2) To find an estimate of the average local variances (actually conditional variances) to be used in a subsequent co-kriging program. Assuming a, linear drift model with a = 6, we proceed as follows: (1) Store 6 20 3 0 in file TEMP .PAR. Note that this sequence of values assigns a = 6, ng = 20 (the neighborhood size adopted in all previous analyses), m = 3 (because we are fitting a drift model), and k = 0 (because we require predictions to "honor" the observations). (2) Store [x y z D] in TEMP_AIN. Note that these are exactly the same values as those used in program AX with the linear drift model. (3.) Run program CROSSK and examine file TEMP OUT for the output. Part of this output is shown in Appendix table 8. The crossvalidated residuals should be stored for later use in co-kriging.

linear

Using Progratm KRIG to Predict Water Content Values

Given a set of spatially oriented data, categorized by coordinates

26

26

ALABAMA

ALABAMA AGRICULTRLEP

AGRICULTURAL

EXPERIMENT

IMN

STATION

SAIO

Appendix Table 8. Output for cross-validation of linear drift model for water content, with parameter values a = 6, ng =20, m =3, and k =0. prediction local SD local z z y x mean SE 3 3 cm /cm m m 0.221 0.0128 0.0155 0.216 0.216 0.000 0.000 0.220 0.0123 0.235 0.0137 0.220 0.000 6.250 0.224 0.0093 0.0094 0.243 0.219 0.000 12.500 0.232 0.0115 0.0117 0.232 "00 0.225 18.750 0.0118 0.0120 0.234 0.244 0.232 00 25.000 0.0122 0.0144 0.239 0.239 10 0.241 37.500 0.233 0.0076 0.0076 0.236 0.204 x.1 25 10.938 0.235 0.0109 0.0109 0.234 0.233 3.125 17. 188 0.239 0.0117 0.0117 0.241 0.229 23.438 3.125 0.242 0.0121 0.0129 0.242 0.240 3.125 29.688 0.0128 0.0136 0.232 0.232 0.230 6.250 6.250 0.243 0.01 17 0.0129 0.243 0.259 6.250 34.375 0.245 0.0121 0.0118 0.245 0.246 9.375 17. 188 0.245 0.0128 0.0131 0.245 0.250 12.500 12.500

14.063 3.125 18.750 21.875 12.500 15 .625 10.938 3.125 18.750 21.875 34.375 1.563 7.813 14.063 29.688 0.000 12.500 25.000 37.500 ------

53. 125 56.250 56.250 56.250 62.500 62.500 65.625 68.750 68.750 68.750 68.750 71 .875 71.875 71.875 7 1.875 75.000 75.000 75.000 75.000 -----

0.314 0.299 0.269 0.308 0.291 0.294 0.267 0.292 0.302 0.289 0.234 0.287 0.303 0.300 0.286 0.296 0.305 0.280 0.255 r-

0.294 0.290 0.293 0.273 0.288 0.293 0.294 0.290 0.286 0.285 0.263 0.293 0.290 0.291 0.264 0.289 0.289 0.274 0.256 ___

0.0149 0.0150 0.0144 0.0149 0.0133 0.0138 0.0134 0.0154 0.0162 0.0166 0.0176 0.0147 0.0158 0.0166 0.0180 0.0167 0.0168 0.0184 0.0220

0.0150 0.0141 0.0148 0.0153 0.0142 0.0142 0.0135 0.0149 0.0166 0.0170 0.0159 0.0149 0.0147 0.0166 0.0166 0.0149 0.0163 0.0170 0.0179

0.290 0.290 0.286 0.277 0.293 0.292 0.294 0.292 0.283 0.280 0.263 0.293 0.290 0.287 0.264 0.290 0.286 0.274 0.256

KRIGING

27

x, y and variate values z, values of Z will be predicted at points not necessarily covered by the original data. Best linear unbiased estimates are obtained by the procedure known as kriging. The background to the theory is given in an earlier section (page 4,5). The program implements this theory for any given mean model and co-variance parameter a. In addition, it requires specification of neighborhood size. Note the following necessary inequalities. n ng m+1.

The program works by choosing a subset of the data, namely points in the neighborhood of the prediction point (x*,y*), and then making predictions from this neighborhood. Each time the prediction point changes, a new neighborhood is defined and the process repeated. It may happen that two neighborhoods coincide, but the program does not take advantage of this possibility. It should be noted that, when a prediction point coincides with an observed point, the kriged prediction value equals the observed value. If smoother plots are required, which can be justified on the grounds of the presence of a "nugget" variance, the smoothing parameter, k, can be set non-zero, but now the kriged and observed values differ. Assume that prediction of water content values at 400 regularly spaced points on a grid covering the entire field is desired. This can be accomplished as follows: (1) Store 6 20 3 0 in file TEMP PAR. Note that a = 6, ng= 20, m = 3, and k = 0 are exactly the same values as those used in the earlier program CROSSK. (2) Store [x y z D] in file TEMP_AIN. Note that these also are exactly the same values as those stored in this file for both earlier programs. (3) Store the array [x* y* D*], made up of the 400 prediction coordinates and corresponding design matrix, on the file TEMP COR. 4) Run the program KRIG and examine the file TEMP OUT for the output. Appendix table 9 lists part of this output. One would normally use this output to create a contour plot, for example by using a package such as PLOT88 (Plotworks, Inc., La Jolla, CA 92037-0635). Appendix figure 4 gives such a graphical display for the predicted values, while Appendix figure 5 gives a similar display of the prediction standard errors.

28

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

Developing and Using a Kriging Model for Clay Content The reason for developing a kriging model for clay is that such a model is needed later in the co-kriging of water content. The percentage clay at 80 cm was chosen as the variable of interest since it is closely correlated with water content and is a physical property of the soil which hardly changes over time.
Using Program AX to Find the Best Clay Model

The sequence of models, namely the stationary, linear and quadratic drift model, with a varying in value from 4 m to 13 m, was fitted to determine the best value for a. Appendix table 10 shows part of the output from the fitting of a linear drift model. The same results are depicted graphically in Appendix figure 6. From these results, depending heavily on the inter-quartile range, we chose a value for oT

E

)

0. O0

7.50

15.00

22.50

30.00

37.50

X,

m

APPENDIX FIGURE 4. Kriging estimates of soil water content using a linear drift model.

KRIGING

29

KRIGING Appendix Table 9. Kriging estimates of water content using a linear drift model, with parameter values az=6, ng =20, and k =0.

2

m=3,

X m 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

Y m 0.000 3.215 6.250 9.375 12.500 15.625 18.750 21 .875 25.000 28. 125 0 .221 0 .223 0 .227 0..232 0..237 0..242 0..245 .252 0..256 0. .258

z

prediction SE 0.00000 0.01404 0.01528 0.01478 0.01466 0.01482
0.01389

0.01349 0.00000 0.01269

local conditional SD cm 3/cm 3 0.01271 0.01366 0.01366 0.01346 0.01346 0.01370 0.01365 0.01365 0.01452 0.01281

local mean 0.219 0.222 0.227 0.232 0.237 0.242 0.248 0.253 0.253 0.256

37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 ~r \IIII

15.625 18.750 21.875 25.000 28.125 31.250 34.375 37.500 40.625 43.750 46.875 50.000 53.125 56.250 59.375 62.500 65.625 68.750 71.875 75.000 1\11111

0.252 0.255 0.249 0.256 0.251 0.251 0.251 0.250 0.247 0.247 0.245 0.239 0.246 0.247 0.247 0.250 0.248 0.245 0.251 0.255

0.01696 0.01697 0.01480 0.00000 0.01403 0.01537 0.01569 0.01549 0.01719 0.01608 0.01466 0.00000 0.01493 0.01605 0.01601 0.01885 0.01849 0.01759 0.01723 0.00000 111111111111

0.01549 0.01549 0.01481 0.01363 0.01405 0.01412 0.01433 0.01402 0.01566 0.01469 0.01463 0.01488 0.01488 0.01453 0.01453 0.01706 0.01699 0.01706 0.01706 0.01706

0.252 0.255 0.247 0.248 0.249 0.25 1 0.251 0.250 0.247 0.247 0.248 0.249 0.249 0.247 0.247 0.250 0.250 0.252 0.253 0.253

30

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

Appendix Table 10. Output program AX when fitting a linear drift model to clay content. a m 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 Rank means 9.12 9.17 9.20 9.25 9.56 9.92 10.10 9.98 9.53 9.05 8.44 8.22 8.76 9.61 10.39 11.22 12.24 12.85 13.39 68.10 68.48 68.73 69.11 71.39 74.05 75.44 74.56 71.14 67.59 63.04 61.39 65.44 71.77 77.59 83.80 91.39 95.95 100.00 Mean Mean SS(resid) Abs(resid) 2 (pct.) pct. 10.92 2.43 11.04 2.43 11.22 2.43 11.44 2.44 11.72 2.47 12.07 2.50 12.42 2.53 12.74 2.55 13.00 2.56 13.22 2.56 13.43 2.57 13.66 2.59 13.92 2.63 14.25 2.67 14.66 2.72 15.17 2.80 15.68 2.88 16.15 2.95 16.59 3.00 Percentage of worst value 65.83 80.93 66.55 80.93 67.62 80.95 82.34 68.96 70.62 82.17 72.74 83.33 74.86 84.30 76.77 84.93 78.35 85.20 79.71 85.22 85.53 80.98 86.24 82.31 87.45 83.89 85.88 88.90 88.38 90.73 91.45 93.21 94.50 95.80 97.36 98.06 100.00 100.00 I.Q.R. pct. 3.43 3.41 3.26 3.17 3.11 3.08 3.13 3.20 3.12 3.92 2.79 2.64 2.69 2.90 3.08 3.26 3.46 3.77 4.40 84.77 84.44 80.68 78.40 77.05 76.08 77.39 79.16 77.23 72.26 68.92 65.26 66.47 71.70 76.18 80.65 85.49 93.31 100.00 I.D.R. pct. 8.19 8.50 9.01 9.12 9.14 9.16 9.20 9.22 9.28 9.58 9.88 10.15 10.40 10.45 10.51 10.66 10.76 10.87 10.96 74.77 77.60 82.80 83.27 83.41 83.57 83.95 84.16 84.70 87.41 90.16 92.64 94.92 95.36 95.91 97.28 98.22 99.15 100.00

KRIGING
0

31

o
0

~I
0.016
O pC

o

M

CD
d

O"

00.00

7.50

15.00

22.50

30.00

37.50

APPENDIX FIGURE 5. Prediction standard errors of soil water content values

based on a linear drift kriging model.

0)90

S80

S70
0

IJ

4

6

8
Range,

10

12

14

M

APPENDIX FIGURE 6. Measures of deviation (as 01oof worst value) of cross-validated

residuals using a linear drift model on clay content.

32

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

a of 9. Analysis of all three models reveals that the linear drift model is to be preferred to either a stationary or quadratic model. Appendix table 11 summarizes some of these results.
Using Program CROSSK to Validate the Given Clay Model

The given model is that of linear drift with a = 9. The procedure of validation is an exact repetition of that for the water model. Part of the output is shown in Appendix table 12. It is important to store the cross-validated residuals as these are needed in the co-kriging analysis to come.
Using Program KRIG to Predict Clay Values

The input requirements reflect an exact repetition of those for the water model. Part of the output is shown in Appendix table 13. Co-kriging Water with Clay To run the co-kriging program, the best kriging models for clay and water individually must be known. These models were determined in sections "Developing and Using a Kriging Model for Water Content" and "Development and Using a Kriging Model for Clay Content." In addition, the variances of the cross-validated residuals for clay and for water and the correlations between these two sets of residuals are required. Appendix table 14 lists these residuals for the 58
Appendix Table 11. Measures of deviation of cross-validated residuals for three different drift models on clay content. Sum of squared deviations Inter-quartile range pct. (pct.)= Q S L S L 3.43 644.31 722.67 3.64 695.01 3.56 3.41 651.37 736.38 695.38 3.61 3.26 661.89 753.65 698.57 3.17 772.92 3.63 704.55 674.96 3.11 795.31 3.58 691.20 713.39
724.80 734.85 711.97 732.74 823.15 850.63 3.60 3.38 3.08 3.13

a m 4.0 4.5 5.0 5.5 6.0
6.5 7.0

Q 3.88 3.94 4.06 4.07 4.18
4.24 4.08

7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0
S = stationary

742.81 748.43 753.00 758.83 767.98 781.32 799.33 822.97 852.65 882.10 909.47 934.47
model, L

751.81 766.91 780.21 792.65 805.65 821.13 840.61 865.11 895.14 924.99 952.97 978.80 :liea

876.29 899.74 921.54 942.50 963.53 986.58 1013.37 1045.05 1082.17 1118.31 1151.79 1182.76
Q qadati

3.04 2.87 2.78 2.87 3.00 3.04 3.01 2.99 3.14 3.36 3.59 3.76
diftmoel

3.20 3.12 2.92 2.79 2.64 2.69 2.90 3.08 3.26 3.46 3.77 4.40

4.06 3.98 3.85 3.83 3.87 3.89 4.00 4.13 4.31 4.52 4.75 4.96

drftmoelan

KRIGING

KRIGING

33

3

locations where both water and clay were recorded. An analysis reveals the following: variance(water residuals) = a2 = 0.0002116, variance(clay residuals) = 2 = 13.36, and the
correlation coefficient = p = 0.6273.
Using Program COK to Predict Water Values by Co-Kriging

We will make predictions at the same 400 regularly spaced grid

points used earlier in kriging. (1) Store the following parameter values on file TEMP.PAR: 6.0 9.0 0.6273 0.0002116 13.36 3 3 20 20 0.0.

X m 0.000 6.250 12.500 18.750 25.000 37.500

,oi:

Appendix Table 12. Output for cross-validation of linear drift model

-9, ng - 20, for clay content, with parameter values co Y Z 15.4 19.2 17.6 17.7 17.4 17.7 13.4 18.4 17.3 17.1 18.1 21.5 20.7 z 16.66 15.72 15.26 17.71 16.92 18.49 18.25 17.97 18.33 18.65 16.81 19.45 20.13 20.02 prediction SE pct. 3.097 2.479 1.898 2.139 2.198 2.932 1.484 2.081 2.077 2.401 2.611 2.581 2.390 2.595

m =3,

and k

=

0.

local SD

local mean

10.938 17.188 23.438
29.688 6.250 34.375 17.188 12.500

m 0.000 0.000 0.000 0.000 0.000 0.000 3.125 3.125 3.125 3.125 6.250 6.250 9.375 12.500

19.4

2.549 2.425 2.201 2.389 2.422 2.420 1.707 2.386 2.340 2.396 2.543 2.548 2.401 2.644

16.29 16.12 16.81 17.30 17.61 18.51 17.94 18.03 18.55 18.89 17.47 19.81 19.86 19.89

14.063 3.125 18.750 21 .875 12.500 15 .625 10.938 3.125 18.750 21 .875 34.275 1.563

7.813 14.063
29.688 0.000 12.500 25.000 37.500

53.125 56.250. 56.250 56.250 62.500 62.500 65.625 68.750 68.750 68.750 68.750 71.875 7 1.875 71.875 7 1.875 75 .000 75 .000 75.000 75.000

26.8 26.0 23.6 32.9 22.8 24.3 22.0 28.2 25.6 22.7 16.0 20.7 23.5 26.5 29.5 26.9 25.7 20.1 19.8

27.59 25.29 28.83 23.99 23.78 24.01 23.95 21.76 24.10 23.79 24.09 27.03 24.81 25.12 19.42 19.26 24.33 24.39 18.68

3.110
3.962 2.758 3.155 2.870 3.318 3.575 3.292 3.879 3.968 4.407 2.781 4.089 4.009 3.441 3.780 4.230 4.664 5.139

3.570 3.842 3.273 3.625 3.640 3.826

3.966
3.614 4.589 4.558 3.828 3.503 4.119 4.597 3.414 3.710 4.599 4.446 4.181

26.01 2576 25.79 24.24 24.85 25.56 25.06 23.78 23.95 23.18 22.88 24.92 24.55 23.89 20.23 20.61 23.35 23.00 19.08

34

34

ALABAMA

ALABAMA AGRICULTRLEP

AGRICULTURAL

EXPERIMENT

IMN

STATION

SAIO

Appendix Table 13 Kriging estimates of clay content using a linear drift model,
with parameter values

a=9,

n

=

20,

m =3, and

k = 0.

x m 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

yz m 0.000 3.125 6.250 9.375 12.500 15 .625 18.750 21.875 25.000 28.125 15.4 16.2 17.0 17.6 18.3 19.9 20.8 20.3 19.2 20.9

prediction SE 0.00 2.86 3.42 3.51 3.48 3.28 2.91 2.74 0.00 2.79

local conditional SD pct 2.51 3.18 3.18 3.17 3.17 3.22 3.25 3.34 3.40 3.39

local mean 16.19 16.25 16.91 17.48 18.14 19.32 20.37 20.69 20.37 21.15

37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500 37.500

15.625 18.750 21.875 25.000 28. 125 31.250 34.375 37.500 40.625 43.750 46.875 50.000 53.125 56.250 59.375 62.500 65.625 68.750 7 1.875 75.000

22.3 22.8 21.6 22.5 22.4 22.7 23.1 24.1 22.4 22.5 23.2 24.4 22.8 22.4 21.8 20.8 18.9 17.3 18.4 19.8

2.65 2.58 2.30 0.00 2.37 2.97 3.95 4.10 4.61 4.45 3.59 0.00 3.79 4.65 4.79 4.64 4.34 3.76 3.46 0.00

2.37 2.37 2.57 2.55 2.65 2.78 3.55 3.64 4.14 4.11 4.00 4.23 4.23 4.25 4.25 4.18 4.22 4.17 4.17 4.08

22.26 22.97 21.24 21.67 21.93 22.60 23.10 24.14 22.40 22.10 21.95 21.66 21.47 22.15 21.81 21.34 21.20 20.72 20.18 19.98

The meaning of each component of this parameter vector is explained in Appendix table 2 (under program COK on page 16). (2) Store the array [x y z (water) D] on file TEMPAIN. (3) Store the array [x y z (clay) D] on file TEMPBIN. Note that these file contents are exactly the same as those originally stored on file TEMPAIN for the earlier kriging analyses. (4) Store the array [x y* D*I, made up of the 400 prediction coordinates and corresponding design matrix on file TEMPCOR. Note that this is identical to the contents stored in the same file earlier when using program KRIG. (5) Run the program COK and examine the output on file TEMPOUT. Appendix table 15 lists part of this output. Appendix figures 7 and 8 give the contour plots for predicted

KRIGING Appendix Table 14. Cross-validated residuals from the kriging models on water and clay respectively. water clay water cm 3/cm 3 pct. cm 3/cm 3 0.0041 -1.2574 0.0170 0.0157 3.4809 0.0090 0.0245 2.3384 -0.0015 -0.0065 -0.0069 0.0125 0.0112 0.4754 -0.0199 0.0023 -0.7857 0.0116 -0.0325 -4.8478 -0.0035 -0.0013 0.4267 0.0149 -0.0119 -1.0305 -0.0117 -0.0019 -1.5540 0.0105 -0.0015 1.2948 0.0201 0.0153 -0.0457 0.0089 0.0014 1.3742 -0.0246 0.0047 0.6846 0.0348 0.0117 4.5687 0.0034 -0.0059 -0.9914 0.0008 -0.0094 1.8468 -0.0265 0.0104 -1.1814 0.0018 0.0031 0.0068 0.0162 0.0047 -0.5285 0.0038 0.0102 -0.0624 -0.0291 0.0079 -0.1707 -0.0066 0.0089 0.9369 0.0128 -0.0321 -8.0718 0.0086 0.0064 -0.4702 0.0213 0.0210 7.3490 0.0069 0.0186 3.8222 0.0163 0.0127 0.4496 0.0060 0.0113 0.0092 -0.0002 The 58 observations relate to those 58 places where both variables were recorded.

35

clay pct. 5.5313 -2.3833 -1.8158 2.3710 0.0764 5.1165 0.2382 -3.5636 3.5411 -1.7545 -0.7928 0.7074 -5.2342 8.9130 -0.9759 0.2948 -1.9522 6.4434 1.4985 -1.0939 -8.0895 -6.3272 -1.3082 1.3795 10.0780 7.6365 1.3736 -4.2853 1.1182

water and prediction standard errors, respectively, producef by

the package PLOTT88.
Using Program CROSSC to Validate The CO-KRIGED Model.

(1) The contents of files TEMP PAR , TEMP__AIN and TEMP BIN remain exactly the same as those used in the previous program, COK. (2) Run program CROSSC and examine file TEMP OUT for the output. Appendix table 16 lists part of that output. It is immediately apparent, both from the smaller residuals and smaller prediction standard errors, that the co-kriged predictions are considerably more precise than those based on kriging alone. Finding an Optimal Sampling Strategy The program SAM is used to answer the following question. Requiring minimal loss of information, which locations should be selected to measure water content values at only 16 sites in the future? There are 60 locations that were used in the past; call these the future potential sample points. Now define 18 prediction points, none of which should coincide with the potential sample points. The idea is to choose a subset

36 36

ALABAMA

ALABAMA AGRICULTURAL EXPERIMENT STATION
AGRICULTURAL EXPERIMENT STATION

E

>-

X,
water and clay content.

M

APPENDIX FIGURE 7. Co-kriging estimates of soil water content using linear drift models for

KRIGING
o o

37

o

o
0
COO

=:ID

0

0

o0.01
0

O0

0

C0

0.00

7.50

15.00
X,

22.50 m

30.00

37.50

APPENDIX FIGURE 8. Prediction standard errors of soil water content values based on the co-kriging model.

38

38

ALABAMA

ALABAMA AGRICULTRLEPIMNSAIO

AGRICULTURAL

EXPERIMENT

STATION

of 16 points from the full set of 60 potential sample points in such a way that the maximum prediction standard error of all 18 prediction points is minimized. The actual technique used in the program is to calculate the prediction SE's for the closest 6 prediction points to any potential sample point and hence the average of these 6 points. This average SE is used as the criterion for retaining a potential sample point; those 16 points having largest average SE's are retained. Appendix figure 9 shows the layout of the potential sample points, the prediction points, and the finally selected optimal 16 points. In
using program SAM we make use of our knowledge of the best co-

kriging model for water and clay.
Appendix Table 15. Output of co-kriging water and clay with linear drift models. Parameter values are a~(water)= 6.0, as(clay) = 9.0, p = 0.63, u2(water)= 0.0002116, u2(clay) = 13.40, m(water)= 3, m(clay) = 3, ng (water) = 20, ng(clay) - 20, and k = 0.0 x U311 y water m 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 m 0.0000 3. 1250 6.2500 9.3750 12.5000 15 .6250 18.7500 21.8750 25.0000 28. 1250

1

prediction SE water

local conditional SD 3VI 3 ~~ VUI U3~water cm /cm 0.0127 0.0137 0.0037 0.0135 0.0135 0.0137 0.0136 0.0137 0.0146 0.0128

local mean water

0.2205 0.2230 0.2273 0.2323 0.2373 0.2434 0.2471 0.25 17 0.2558 0.2573

0.0000 0.0140 0.0153 0.0148 0.0147 0.0148 0.0138 0.0135 0.0000 0.0127

0.2 189 0.2225 0.2273 0.2323 0.2373 0.2427 0.2481 0.2522 0.25 18 0.2557

37.5000 37 .5000 37 .5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000 37.5000

15 .6250 18.7500 21.8750 25.0000 28. 1250 31.2500 34.3750 37.5000 40.6250 43.7500 46.8750 50.0000 53. 1250 56.2500 59.3750 62.5000 65 .6250 68.7500 71.8750 75.0000

0.2536 0.2564 0.2496 0.2560 0.2522 0.2523 0.2520 0.25 18 0.2485 0.2480 0.2477 0.2389 0.2480 0.2473 0.248 1 0.2491 0.2473 0.2423 0.2485 0.2554

0.0170 0.0170 0.0148 0.0000 0.0140 0.0153 0.0157 0.0155 0.0172 0.0161 0.0146 0.0000 0.0149 0.0161 0.0160 0.0189 0.0185 0.0176 0.0172 0.0000

0.0156 0.0156 0.0148 0.0137 0.0141 0.0142 0.0144 0.0141 0.0157 0.0147 0.0146 0.0149 0.0149 0.0146 0.0146 0.017 1 0.0170 0.017 1 0.017 1 0.017 1

0.2536 0.2564 0.2476 0.2490 0.2498 0.2523 0.2520 0.25 18 0.2485 0.2480 0.2484 0.2480 0.2485 0.2473 0.248 1 0.2491 0.2507 0.2504 0.2507 0.2524

KRIGING

KRIGING

39

3

Appendix Table 16. Output of cross-validation of co-kriged model of water and clay. Parameter values are a (water) = 6.0, a (clay) = 9.0, P = 0.63, a2 (water) = 0.0002116, 2 v (clay) = 13.40, m (water) = 3, m (clay) = 3, ng (water) = 20, ng (clay) = 20, and k = 0.0 xx y Z water water prediction local SE conditional water SD water cm 3/cm 3 0.0123 0.0128 0.0102 0.0124 0.0070 0.0093 0.0089 0.0116 0.0092 0.0119 0.0115 0.0123 0.0057 0.0076 0.0081 0.0109 0.0090 0.0118 0.0099 0.0121 0.0105 0.0128 0.0100 0.0118 0.0094 0.0118 0.0101 0.0128 local mean water

m 0.0000 6.2500 12.5000 18.7500 25.0000 37.5000 10.9375 17.1875 23.4375 29.6875 6.2500 34.3750 17.1875 12.5000

m 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3. 1250 3.1250 3.1250 3.1250 6.2500 6.2500 9.3750 12.5000

0.2205 0.2352 0.2430 0.225 1 0.2436 0.2410 0.2038 0.2328 0.2290 0.2400 0.2301 0.2580 0.2462 0.2499

0.2139 0.2308 0.2243 0.2340 0.2323 0.2383 0.2219 0.2330 0.2375 0.2394 0.2326 0.2450 0.2500 0.2477

0.2157 0.2219 0.2248 0.2323 0.2341 0.2387 0.2323 0.2352 0.2393 0.2417 0.23 19 0.2443 0.2457 0.2454

14.0625 3. 1250 18.7500 21.8750 12.5000 15 .6250 10.9375 3. 1250 18.7500 21.8750 34.3750 1.5625 7.8125 14.0625 29.6875 0.0000 12.5000 25.0000 37.5000 IL)~~nL~

53. 1250 56.2500 56.2500 56.2500 62.5000 62.5000 65.6250 68.7500 68.7500 68.7500 68.7500 71.8750 71.8750 71.8750 71.8750 75.0000 75.0000 75.0000 75.0000

0.3136 0.2991 0.2691 0.3077 0.2912 0.2936 0.2672 0.2922 0.3023 0.2893 0.2341 0.2868 0.3033 0.2999 0.2856 0.2957 0.3053 0.2800 0.2554 ~n ~I~I Ir

0.2912 0.2917 0.2828 0.2957 0.2840 0.2907 0.2878 0.3082 0.2915 0.2839 0.2413 0.2749 0.2882 0.2979 0.2923 0.3054 0.2944 0.2639 0.2554

0.0110 0.0114 0.0106 0.01 14 0.0101 0.0104 0.0103 0.01 17 0.0120 0.0126 0.0134 0.0110 0.01 16 0.0123 0.0135 0.0129 0.0128 0.0141 0.0172

0.0150 0.0141 0.0148 0.0153 0.0142 0.0142 0.0135 0.0150 0.0166 0.0170 0.0160 0.0149 0.0147 0.0166 0.0167 0.0150 0.0163 0.0170 0.0179

0.2894 0.2906 0.2853 0.2776 0.2926 0.2915 0.2941 0.2949 0.2836 0.2800 0.2588 0.2909 0.2906 0.2879 0.2686 0.2952 0.2869 0.2725 0.2554

(1) Store the following parameter vector (as a row vector) in file TEMP PAR: ca(water) 6.0 9.0 a(clay) 0.6273 cr(water) 0.0002116 13.36 u2(clay) 3 m (water) 3 m(clay)

40

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

ng(water) 12 (must be less than 18) ng(clay) 20 maxout 20 (maximum number of iterations) (2) Store [x y D], the coordinates and design matrix for the potential sample points, in file TEMPAIN. Note that the design matrix D corresponds to the linear drift model. Also note that no knowledge of the observed values for water content at these sites is needed. (3) Store [x y z(clay) D] in the file TEMP BIN. Note that these stored values are identical to those used in the earlier co-kriging analysis of water with clay. (4) Store [x* y* D*], the coordinates and design matrix for the 18

0 ,,,Y,,L~ LY\I
0i2

"
9
ml

0iXEl

9

2 X By

X E 0 a X

E
e

X a2
.1

2

H
01 2

E
0 O

2
a1iX

0
X
2

0
2

2

X 2

X 2
O1S1

2
S1

X 21 D1

8
. 0
2

0
D

X

91 @1

XE:p pJ2

2
@1

D

50 m
APPENDIX FIGURE 9. Lay-out of the potential sample points (solid circles and squares), the prediction points (crosses), the selected optimal 16 points after one run (1), and after a second run (2).

KRIGING

41

prediction points, in file TEMPCOR. These values are shown in Appendix table 17. (5) Run the program by entering: RUN SAM. We are now in an interactive mode and the command "Enter the required sample size, ns. Note ns = 0 stops the program," will be displayed on the screen. By entering 16, the program will begin the optimization process, searching for an optimal 16 points. It starts by simply taking the first 16 potential sample points and then runs through 20 iterative cycles (20 being the specified value of the parameter "maxout''), in each cycle replacing one of the 16 points by a better point from the remaining 44. At the end of each cycle it prints out the details of the points added and dropped. When the 20 cycles are complete (or earlier if no improvement can be found), the earlier question is repeated on the screen. If a non-zero value for ns is entered the whole procedure is repeated but with the important difference that it starts with the 16 points from the final stage of the earlier run. Appendix tables 18, 19, and 20 give the output, while Appendix figure 9 illustrates the positions of the best 16 points. Note that after cycle number 14 the program keeps repeating the same exchange between points (34.375, 68.75) and (3.125, 43.75). At this stage the current run is unable to effect any improvement. The final arrangement arrived at can never be guaranteed to be best and a good idea in practice is to run an extra sequence of cycles with a larger ns value and then revert to the real target sample size. This is easily
Appendix Table 17. The input file TEMPCOR required by program SAM to optimize water sampling x m 6.250 18.750 31.250 6.250 18.750 31.250 6.250 18.750 31.250 6.250 18.750 31.250 6.250 18.750 31.250 6.250 18.750
31.250

y m 6.250 6.250 6.250 18.750 18.750 18.750 31.250 31.250 31.250 43.750 43.750 43.750 56.250 56.250 56.250 68.750 68.750
68.750

1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
1.000

x m 6.250 18.750 31.250 6.250 18.750 31.250 6.250 18.750 31.250 6.250 18.750 31.250 6.250 18.750 31.250 6.250 18.750
31.250

y m 6.250 6.250 6.250 18.750 18.750 18.750 31.250 31.250 31.250 43.750 43.750 43.750 56.250 56.250 56.250 68.750 68.750
68.750

42

42

ALABAMA

ALABAMA AGRICULTRLEPIMNSAIO

AGRICULTURAL

EXPERIMENT

STATION

achieved by entering the required non-zero ns value at each stage when the computer screen requests such input. In fact, Apendix figure 9 shows two sets of "best" 16 points. The program was run once as previously reported and then re-run with inproved starting values to eventually yield an even better set of points. Note also that the program is very computer-intensive and it is more efficent to run several sequences of cycles than one long sequence with a high value of the parameter maxout. Finally, note that a different set of prediction points will yield a different (often markedly so) set of optimal sampling sites.
Appendix Table 18. Optimization of sample locations for water content determination. Output from program 5AM - part 1. Parameter values are a (water)= 6.0, a (clay)= 9.0, p = 0.63, a2 (water) = 0.0002116, Q2 (clay)= 13.36, m (water)= 3, m (clay)= 3, ng (water)= 12, ng (clay)= 20, maxout= 20. Cycle number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 sample site added y x m m UVIIIL3 68.7500 34.3750 31.2500 25.0000 3.1250 68.7500 31.2500 3.1250 25.0000 37.5000 37.5000 25.0000 31.2500 43.7500 43.7500 3.1250 .50.0000 37.5000 25.0000 25.0000 43.7500 3.1250 50.0000 31.2500 43.7500 3.1250 68.7500 34.3750 43.7500 3.1250 68.7500 34.3750 43.7500 3.1250 68.7500 34.3750 43.7500 3.1250 68.7500 -~~'--34.3750 sample site dropped y x m m 0.0000 6.2500 3.1250 10.9380 0.0000 0.0000 0.0000 12.5000 6.2500 6.2500 3.1250 17.1880 29.6880 3.1250 43.7500 31.2500 0.0000 18.7500 43.7500 3.1250 6.2500 34.3750 43.7500 3.1250 68.7500 34.3750 43.7500 3.1250 68.7500 34.3750 43.7500 3.1250 68.7500 34.3750 43.7500 3.1250 68.7500 34.3750 43.7500 3.1250

KRIGING

43

Appendix Table 19. Optimization of sample locations for water content determination. Output from program SAM - part 2. x m 3.1250 34.3750 3.1250 37.5000 25.0000 37.5000 31.2500 25.0000 23.4380 25.0000 37.5000 31.2500 17.1880 12.5000 25.0000 7.8130 3.1250 4.6880 0.0000 12.5000 18.7500 34.3750 10.9380 6.2500 4.6880 10.9380 12.5000 9.3750 12.5000 17.1880 31.2500 29.6880 7.8130 14.0630 20.3130 0.0000 12.5000 25.0000 3.1250 18.7500 7.8130 14.0630 3.1250 18.7500 21.8750 12.5000 15.6250 10.9380 0.0000 18.7500 21.8750 6.2500 1.5630 7.8130 14.0630 Continued y m The best subset of sample points with their average SE's 68.7500 68.7500 31.2500 50.0000 0.0000 0.0000 25.0000 37.5000 3.1250 25.0000 25.0000 50.0000 9.3750 12.5000 12.5000 15.6250 The subset of unused sample points with their average SE's 18.7500 21.8750 25.0000 25.0000 25.0000 6.2500 3.1250 6.2500 28.1250 28.1250 0.0000 37.5000 37.5000 3.1250 43.7500 3.1250 46.8750 46.8750 46.8750 50.0000 50.0000 50.0000 43.7500 0.0000 53.1250 53.1250 56.2500 56.2500 56.2500 62.5000 62.5000 65.6250 0.0000 68.7500 68.7500 0.0000 71.8750 71.8750 71.8750 Average SE 3 3 cm /cm 0.0150 0.0163 0.0146 0.0157 0.0147 0.0149 0.0154 0.0153 0.0147 0.0152 0.0154 0.0157 0.0146 0.0146 0.0148 0.0145 0.0145 0.0146 0.0145 0.0148 0.0151 0.0149 0.0145 0.0145 0.0146 0.0147 0.0146 0.0147 0.0147 0.0146 0.0153 0.0148 0.0145 0.0147 0.0149 0.0146 0.0146 0.0152 0.0145 0.0146 0.0148 0.0145 0.0149 0.0150 0.0146 0.0149 0.0148 0.0150 0.0144 0.0151 0.0150 0.0145 0.0149 0.0151 0.0150

44

ALABAMA

AGRICULTURAL

EXPERIMENT

STATION

Appendix Table 19 (continued). Optimization of sample locations for water content determination. Output from program SAM - part 2. x m 29.6880 0.0000 12.5000 25.0000 37.5000 y m The best subset of sample points with their average SE's 71.8750 75.0000 75.0000 75.0000 75.0000 Average SE 3 cm 3/cm 0.0161 0.0149 0.0151 0.0154 0.0158

Appendix Table 20. Optimization of sample locations for water content determination. Output from program SAM - part 3. The prediction SE's. x m 6.2500 18.7500 31.2500 6.2500 18.7500 31.2500 6.2500 18.7500 31.2500 6.2500 18.7500 31.2500 6.2500 18.7500 31.2500 6.2500 18.7500 31.2500 y m 6.2500 6.2500 6.2500 18.7500 18.7500 18.7500 31.2500 31.2500 31.2500 43.7500 43.7500 43.7500 56.2500 56.2500 56.2500 68.7500 68.7500 68.7500 Initial SE 3 cm /cm' 0.0000 0.0139 0.0139 0.0163 0.0177 0.0195 0.0244 0.0262 0.0295 0.0348 0.0375 0.0424 0.0456 0.0475 0.0540 0.0556 0.0588 0.0675 Final SE cm 3/cm' 0.0140 0.0144 0.0148 0.0145 0.0151 0.0155 0.0144 0.0151 0.0154 0.0142 0.0149 0.0158 0.0155 0.0131 0.0166 0.0151 0.0143 0.0173

With an agricultural research unit in every major soil area, Auburn University serves the needs of field crop, livestock, forestry, and horticultural producers in each region in Alabama. Every citizen of the State has a stake in this research program, since any advantage from new and more economical ways of producing and handling farm products directly benefits the consuming public.

0

ww. O

0

00Co

1a

is

® Main Agricultural Experiment Station, Auburn.
i E. V. Smith Research Center, Shorter.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Tennessee Valley Substation. Belle Mina Sand Mountain Substation. Crossv ile North Alabama Horticulture Substation. Cullman Upper Coasta Pain Substation. Winfield Forestry Unit. Fayette County Cn iton Area Horticulture Substat on Clanton Forestry Unit. Coosa County Piedmont Substation. Camp Hill Plant Breeding Unit. Tallassee Forestry Unit. Autauga County Prattvi le Experiment Field. Prattville Black Belt Substation Marion Junction The Turnipseed-Ikenberry Pace. Union Springs Lower Coasta Plain Substation. Camden Forestry Ur t. Barbour County Monroevi le Exper ment Field. Monroeville Wiregrass Substation. Headland Brewton Experiment Feld. Brewton Solon Dixon Forestry Education Center. Covington and Escambia counties Ornamenta Horticulture Substation. Spring Hill Gulf Coast Substation Fairhope