in .,.. _,,. .
. .3... aw, hm...
11m? ..

.i 5.35.555

3:;
fit
than».
.Eogfifiu
.
u... vauxhbasfi.
.. .

 

 

l.

I: 2 . .3 »n. 4
{Jrz «53.3%..

~ .

 

THESIS
4.

{\fl) 3‘ i
\ ' I

illllllllllllllllllll lllllllllLllllll

3129

This is to certify that the

dissertation entitled

ACCURACY OF SOIL PROPERTY MAPS
FOR SITE-SPECIFIC MANAGEMENT

presented by
Thomas G. Mueller

has been accepted towards fulfillment
of the requirements for

Ph . D . degree in So 11 Management

 

 

myé? Mam/Le

Major professor

Date December 16, 1998

 

MS U is an Affirmative Action/Equal Opportunity Institution 0- 12771

 

 

LIBRARY
Michigan State
University

 

 

PLACE IN RETURN BOX
to remove this checkout from your record.
TO AVOID FINES return on or before date due.

 

MTE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

use WWI-m.“

 

ACCURACY OF SOIL PROPERTY MAPS
FOR SITE-SPECIFIC MANAGEMENT

By

Thomas G. Mueller

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Crop and Soil Science

1998

ABSTRACT

ACCURACY OF SOIL PROPERTY MAPS
FOR SITE-SPECIFIC MANAGEMENT

By

Thomas G. Mueller

The accuracy of soil property maps for site-specific management may be
inadequate at sampling intensities recommended by commercial agriculture. Since the
success of site-specific fertilizer applications depends on the quality of soil property
maps, it is critical for Michigan farmers who are adopting these practices to have an
understanding of the accuracy associated with different soil sampling strategies, soil
sampling intensities, and interpolation techniques. This thesis evaluates how grid
sampling and interpolation schemes affected map accuracy based on measures of map
error. The second objective was to evaluate how different interpolation techniques that
incorporate terrain attributes affects the spatial predictions of soil properties and whether
relative performance of these techniques is affected by the scale of soil sampling. In
addition to the soil samples used for spatial interpolation, samples were collected to
assess the quality of the predictions. Grid point sampling at the industry standard
intensity (100 m regular grid), grid cell sampling (100 m grid cells) and directed
sampling based on soil type were not adequate to produce accurate nutrient condition
maps for this field even though most of the variables were spatially structured. Prediction
efficiencies were 0.5 to 10.5 % greater for inverse distance weighted interpolation than

for kriging using a distance exponent of 1.5 at the 30-m grid sampling intensity. At high

sampling densities (30 m regular grid), interpolation methods that utilized terrain
attributes had similar prediction errors to interpolation methods that did not utilize terrain
attributes. At a lower sampling intensity (61 m regular grid), methods that utilized terrain
attributes, especially multiple regression, were more accurate than methods that did not.
At the 61 m grid, the RMSE for multiple regression was lower (3.3 g kg") than the
RMSE for ordinary kriging (4.1 g kg’l). A 100 m grid was not of sufficient intensity to
be used to create accurate maps of soil properties for site-specific nutrient management,
enhancing spatial estimates with terrain attributes can reduced the number of samples
required to create an accurate map. It was not necessary to use complex, time consuming
geostatistical techniques to use terrain attributes because multiple regression was

sufficient.

Copyright by
Thomas G. Mueller
1 998

This dissertation is dedicated

to my wife, Angela,
for her love and perseverance;

to our parents, Philip, Camille, Peter, and Gabriela,
for their support and encouragement,

and to our children, Daniel and Christina,
for the joy they have given us,

ACKNOWLEDGMENTS

I would like to thank Dr. Francis J. Pierce believing in me and for his friendship,
guidance, encouragement, and financial support. I would also like to express my
gratitude to Dr. Darryl Wamcke for serving as the co-chair of my graduate committee
and for his helpful advice. I am very grateful to both Dr. G. Philip Robertson and Dr.
Irvin E. Widders for their advice and for serving on my committee. Again, I would like
to thank my entire committee for their precious time.

A special thanks goes to the farmers, John Anibal and Patrick Feldpausch, for
there generosity. Also thanks to Jim Walseth at Springhill Engineering for mapping
elevation at one of the research locations. Also thanks to Steve Law at the NRCS office
in Clinton County for his help soil sampling.

I am grateful to Dr. Oliver Schabenberger for many hours of consultation. I
appreciate the invaluable advice of Dr. Pierre Goovaerts from the University of
Michigan. I am especially grateful to Dr. Jose E. Cora for his friendship and
collaboration. I am also grateful to for the technical assistance of Brian Long. I also
would like to thank Cal Bricker for his advice and friendship. Thanks to Brian Long, Cal
Bricker, and Gary Zehr for their help in the field near Kalamazoo.

I appreciate the friendship and technical assistance of Brian Baer, Sven Bohm, Dr.
Dave Harris, Dr. Samira Daroub, Brian Graff, and Elaine Parker. Thanks also goes to
Dan Brown for teaching two very useful classes in GIS and for being so helpful. A
special thanks goes to Demetrios Gatziolis for help with GIS and his friendship. I also
appreciate the help of Jon Dahl and Vickie Smith from the Michigan State University soil

testing laboratory for always being so helpful with research and teaching soil samples.

vi

Thanks for all of the hard work' to the students I have supervised: Curtis Beard,
Bill Bower, Frank Boyd, Raulie Casteel, Jolene Friedland, Jeanette Makries, and Edwin
Toledo. I appreciate their friendship and their hard work.

Thanks to Dr. Daniel Rasse, Dr. Joseph Masabni, Mohamed Elwade, James
Sallee, and Djail Santos, Laurent Gilet, Carlos Paglis and Dr. Helvecio De-Polli, and Dr.

Tom Willson, and Jeff Smeenk for their friendship.

vii

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................. ix

LIST OF FIGURES ............................................................................................................ x

INTRODUCTION .............................................................................................................. 1

Prediction Errors and Efficiencies .......................................................................... 2

Spatial Interpolation ................................................................................................ 3
CHAPTER 1 ASSESSING MAP ACCURACY FOR SITE-SPECIFIC

FERTILITY MANAGEMENT .................................................................. 7

INTRODUCTION .............................................................................................................. 7

MATERIALS AND METHODS ...................................................................................... 10

Site description ...................................................................................................... 10

Soil Sampling Design and Laboratory Analysis ................................................... 12

Data Analysis ........................................................................................................ 14

Phase 1: FULL data set .............................................................................. 14

Phase 11: Accuracy of SSF M sampling and interpolation ......................... 15

RESULTS AND DISCUSSION: ...................................................................................... 16

Phase 1. FULL Data Set. ....................................................................................... 16

Phase 11: Accuracy of SSFM sampling and interpolation ..................................... 26

Grid Sampling ........................................................................................... 26

CELL and directed sampling .................................................................... 36

CONCLUSIONS ............................................................................................................... 38

CHAPTER 2 COMPARISON OF TECHNIQUES TO OPTIMIZE
SPATIAL ESTIMATES OF SOIL PROPERTIES USING

TERRAIN ATTRIBUTES ........................................................................ 40
INTRODUCTION ............................................................................................................ 40
MATERIALS AND METHODS ...................................................................................... 45

Site description ...................................................................................................... 45
Soil Sampling Design and Laboratory Analysis ................................................... 47
DEM Creation and Terrain Analysis .................................................................... 48
Data Analysis ........................................................................................................ 48
RESULTS AND DISCUSSION ....................................................................................... 50
Nature of the Data ................................................................................................. 50
Evaluation of Interpolation Methods .................................................................... 58
CONCLUSIONS ............................................................................................................... 65
SUMMARY AND RECOMENDATIONS ...................................................................... 65
BIBLIOGRAPHY ............................................................................................................. 68

viii

LIST OF TABLES

Table 1.1 Map symbols, soil series or complex name, NRCS soil taxonomic

Table 1.2.
Table 1.3.
Table 1.4.

Table 2.1.

Table 2.2.

Table 2.3.

Table 2.4.
Table 2.5.

description (Pregitzer, 197 8), and area occupied in the field ........... 10
Data set descriptions .......................................................................... 13
Summary statistics for the FULL data set. ........................................ 16
Directional and omnidirectional semivariogram model parameters

for the FULL data set ....................................................................... 22

Applications of regression, kriging with an external drift, and
cokriging for enhancing spatial estimates of soil properties using
auxiliary information. ....................................................................... 41

Directional (G30) and omnidirectional (G30 and G61)
semivariogram model parameters for carbon to be used for

ordinary kriging ................................................................................ 55
Correlations between total carbon and terrain attributes and

explained variability (G30 data set). .................................................. 56
Linear model of coregionalization (G30) used for cokriging. . ........... 57

Semivariograrn model parameters fo the direction of maximum
spatial continuity (N 90 E) for kriging with a trend model and
kriging with an external drift. ........................................................... 57

ix

LIST OF FIGURES

Figure 1.1. Modeled empirical semivariogram. .................................................... 6
Figure 1.2. A digital orthophotograph overlain by the 30 m (black square),
the 100 111 (blue cross), the random validation data set (green
circle) and NRCS soil types (CaA = Capac loam with 0 to 4 %
slope; MeA = Metamora-Capac sandy loam with 0 to 4% slope;
MoB = Morley Loam with 2 to 6% slope; and WbA = Wasepi

sandy loam with 2 to 6 % slope)... .................................................. 11
Figure 1.3. Normal probability (Q-Q) plots for variables and log

transformations with 95% confidence intervals. .............................. 18
Figure 1.4. Directional variogram contour maps for the FULL data set. ........... 20
Figure 1.5. Directional and omnidirectional semivariograms for the original

and transformed FULL data sets. ..................................................... 21

Figure 1.6. Contour maps of kriged soil variables and fertilizer
recommendations using the FULL data sets overlain by the soil
type boundaries ................................................................................ 24
Figure 1.7. Omnidirectional experimental and modeled (exponential)
semivariogram for the FULL, Gwmb, G30, G61, and G100 data sets ...27
Figure 1.8. The RMSE for several estimation approaches: IDW interpolation
for distance exponents ranging from 0.1 to 5.0 in 0.1 increments
(line plots on the left), kriging, whole field, and grid cell

sampling; and zone sampling. .......................................................... 28
Figure 1.9. The relationship between the optimal distance exponent and the

CV. ................................................................................................... 30
Figure 1.10. Predicted vs. measured for kriging with the Gcomb, G30, and G100

data sets: ........................................................................................... 31

Figure 1.11. The relationship between prediction efficiency for kriging the
G30 data set, the range of spatial correlation of the FULL data
set, and the RSV of the FULL data set. ............................................ 33
Figure 1.12. The relationship between prediction efficiency for the kriged G30
data set and the range of spatial correlation for the FULL data set

for variables with RSV values greater than 70% .............................. 34
Figure 1.13. Prediction efficiency for kriging versus the prediction efficiency
for IDW.’ ........................................................................................... 35

Figure 1.14. Average semivariograms (large solid circles) and
semivariograms for individual pairs (hallow circles) for P using

the Gm“, data set. ............................................................................. 37
Figure 2.1. Study location: (a) areal photograph and (b) soil and kinematic

GPS measurement locations overlain by soil type boundaries ......... 46
Figure 2.2. Surface maps (a) of elevation and (b) slope and a contour map

(c) of elevation overlain with arrows indicating aspect .................... 51
Figure 2.3. Normal probability (Q-Q) plots for variables and log

transformation with 95% confidence intervals. ................................ 52
Figure 2.4. Sernivariogram surfaces for total carbon (G30) and elevation (n =

1000). ................................................................................................ 53

Figure 2.5. Total carbon (g kg") contour map created with ordinary kriging

using the G30 data set overlain with elevation contours (lines). ....... 55
Figure 2.6. Prediction efficiencies and RMSE's. ................................................ 59
Figure 2.7. Predicted versus measured values for several prediction methods

at two scales of measurement. ......................................................... 60
Figure 2.8. The relationship between prediction efficiency and sample point

density for multiple regression interpolation .................................... 62
Figure 2.9. Predicted total carbon (g kg!) using kriging with an external drift

(G30) and multiple regression (G30 and G61) ..................................... 64

xi

INTRODUCTION

Precision agriculture is about managing soil and crop variability in space and time
in order to improve crop performance and environmental quality (Pierce and Nowak,
1999). Soil chemical and physical properties and crop yield are spatially quite variable.
While some attributes are stable over time (e.g. total carbon, texture), others exhibit a
great deal of temporal variability (e. g. soil N, soil water content water, grain yield, and
pest infestation). The area of precision agriculture related to nutrient management is
often referred to as site-specific fertility management (SSFM). In this dissertation, I am
specifically concerned with the SSF M of soil nutrients that that are spatially variable but
relatively temporally stable, specifically, total carbon, pH, P, K, Ca, and Mg.

Several conditions are essential for successful SSF M, the most basic of which is
that variation in soil properties is adequately known (Pierce and Nowak, 1999 and
Sawyer, 1994). Some of the studies that have tested the validity of this premise have
shown that it is not always valid (Wollenhaupt et a1, 1994; Gotway et al., 1996). Some
SSFM agronomic and economic studies that have had mixed or negative results (Snyder
et al., 1997; Wibawa et al., 1993; Wollenhaupt and Bucholz, 1993) might be explained by
poor prediction accuracy of soil properties. Fortunately, there are good methods for
measuring prediction and map accuracy.

The kriging variance can not be used to estimate the map accuracy because it is
independent of the data values (Deutsch and Journel, 1998) and only dependent on the
covariance model and the data configuration (Goovaerts, 1997). Cross-validation is a
technique where sample data points are sequentially dropped from the prediction data set,

estimated from the neighboring point, and then replaced (Deutsch and Joumel, 1998).

The measured values are subtracted from the predicted values to calculate the residuals
which can be used to assess the accuracy of the predictions. For a regular gridded data
set, cross-validation tends to over estimate prediction errors. A better approach would be
to jack-knifing with an independent validation data set. This approach was used in this
dissertation. The calculation of these errors is discussed in the first section of this

introduction.

The first objective of this dissertation was to evaluate how grid sampling and
interpolation schemes affected map accuracy and prediction efficiency. The next
objective was to evaluate different analytical techniques that using terrain attributes affect
the spatial prediction of soil properties and whether relative performance of these
techniques was affected by the scale of soil sampling. This dissertation requires some
understanding of geostatistics, ordinary kriging, and inverse distance interpolation. The

theory is presented in the second section of this introduction.

Prediction Errors and Efficiencies

In this study, measures of map error included MSE, RMSE (root mean squared
error), and bias. Let vi denote the difference between predicted value and observed value
at location 8;, i=1, ..., nv, where nV is the number of values in the validation data set. A

map correct on average should have EN] = 0. The bias of the map is estimated as

11

V

l ""
Bias = ——Zvi
i-l

and the MSE of the map as

MSE = Bias2 +—l—1—Zv(v,.—v)2

1 "V
——zvs
11v —1i=1

The RMSE is the square root of the MSE, RMSE-:JMS . The mean square
error combines accuracy (biasz) with precision, the variance of the residuals. Prediction

efficiency referred to as goodness by Gotway et a1. (1996) is calculated as
Prediction efficiency = 100% * (MSE .30 — MSEMM XMSE .30 )“. [1]

Positive prediction efficiencies can be interpreted as the percent reduction in MSE as
compared with the field average approach (A30). A prediction efficiency of 15% for
kriging can be interpreted as "compared to the field average approach, kriging reduced
the MSE by 15%." The optimal distance exponent for IDW interpolation was determined

as the one yielding the lowest RMSE.

Spatial Interpolation

A spatial estimator at an unobserved location 3 of an attribute Z is defined as a
weighted average of the observed values of the attribute at spatial locations 5:. The
weights Wsm may be restricted to non-zero values in some neighborhood N“) of the target
location 5. If n(s) is the number of observed values in the neighborhood, the predicted

value is calculated as

,. n(s)
2(s) = Zwsizc.)
i=1
The essential difference between kriging and IDW estimation lies in the
determination of the weights. The weights w for IDW are based on the distance between

the point to be estimated and each of the n sample data points d(sa, 3) within the search

neighborhood N(s).

n(s)
w. = [dew]: Eiders)?

where e is the user defined distance exponent (Gotway, 1996). Smoothness of
interpolated values decreases with the magnitude of the exponent. The kriging weights

are calculated by solving the kriging system presented here in an expanded matrix form

 

 

 

 

 

_ _ _ 1-1 .

W1 CovH . - - Covln 1 Cov” 1
w" Covn1 . Covm 1 COVm
M- _ 1 1 0_ . 1 .

 

where )i. is a Lagrange parameter needed to satisfy certain constraints. Cov“
through Covm are the covariances among the sample data points, and C1, to Cns are the
covariances between each sample data point and the unobserved location 5 (Goovaerts,
1997). For details of constrained minimization through Lagrange multipliers in this
context see Isaaks and Srivastava (1989). Under weak stationarity conditions (Cressie,

1993) the covariances between two data points depends not on their actual coordinate, but

matters. Since the kriging weights are functions of the covariances which in turn depend
on spatial separation of data points, kriging weights are distance related weights too.
The metric in which distances are assessed is not, however, Euclidean distance alone, but
depends on the degree of spatial dependencies. Under second order stationarity, the
covariance model are related to semivariogram models through the relationship COVij =
C(O) - y(hi,~); where C(O) is the variance and hij is the Euclidean distance between
locations 3; and Si the distance between sample locations. The semivariogram models are
fit to empirical semivariograms 11(h) which is the average sample variance of points
separated by distance h but computationally defined as V2 of the average squared

difference of points separated by distance h,

n(h)

)=n[2 1(h)]Z[S (ti)— S(ai+h)] 2

i=1

where n(h) is the number of pairs at lag distance h or in some neighborhood of b. Figure
1.1 illustrates an isotropic semivariogram modeled with an exponential function. The
plateau the variogram reaches is called the sill. The discontinuity at the origin is referred
to as the nugget variance and is attributable to the additive effects of a white noise
process and measurement error (Cressie, 1993). The sill less the nugget variance is
referred to as the structural variance. The separation distance (lag) at which the
variogram reaches the plateau (spherical models) or 95% of the sill (exponential and

Gaussian models) is called the range or range of spatial correlation.

 

12

Variance at 95 % of sill
\
_ _ _5. °

 

0

——f———e—~————e .H__u

    

sm————->

   

Semivariance
O)

 

. / Range 1
[if _ _Nugget variance

0 I.
0 75 150

lag (m)

Figure 1.1. Modeled empirical semivariogram.

l
l
.(_ _ _Structural variance l
I
l
I
I
l

 

Relative structural variability (RSV) was used as a normalized measure of spatial

dependence (Robertson et al., 1993; 1997) and is defined as

RSV = structurallyanance = 1 _ RNE [2]
$1

 

and is related to what is commonly referred to as the relative nugget effect (RN E).

CHAPTER 1

Assessing Map Accuracy for Site-Specific Fertility Management

INTRODUCTION

Several conditions are essential for successful site-specific fertility management
(SSFM). The most basic condition is that variation in soil properties is adequately known
(Pierce and Nowak, 1999 and Sawyer, 1994). Soil property predictions across landscapes
are affected by soil sampling, laboratory analysis, prediction, and cartographic errors.
Poor map quality may explain why results of some SSF M agronomic and economic
studies have had mixed or negative results (Snyder et al., 1997; Wibawa et al., 1993;
Wollenhaupt and Bucholz, 1993). Some have suggested alternatives to grid sampling,
including directed sampling (Pocknee et al., 1996). In general, condition and
management maps are rarely examined for quality, which is unfortunate because methods

exist to assess map accuracy.

Measures of accuracy and goodness have been used in SSFM research to assess
quality of maps and soil properties predictions, mostly by studying the impact of grid.
sampling intensity (Wollenhaupt et al., 1994; Franzen and Peck, 1995; Gotway et al.,
1996, Mohamed et al., 1996) and interpolation techniques (Wollenhaupt et al., 1994;
Gotway et al., 1996, Mohamed et al., 1996) on map or prediction error. Wollenhaupt et
a1. (1994) and Mohamed et al. (1996) considered map accuracy to be the percentage of
areal overlap of mapping categories between maps in question and maps considered

representative of the true spatial distribution of a soil property in space. These truth maps

were arbitrarily defined to be contour maps created either with Delaunay triangulation of
soil properties sampled at a 32-m grid (Wollenhaupt et al., 1994) or with the kriging of

soil properties sampled on a 20 x 40-m grid (Mohamed et al., 1996).

The correlation between predicted and observed data sets has been used as a
measure of prediction accuracy (Franzen and Peck, 1995; Goovaerts, 1997). A
correlation approach alone is problematic. While predicted and observed data sets may
be highly correlated, they may deviate greatly from a 1:1 relationship and therefore be of
low predictive value in a mapping context. Franzen and Peck (1995) assessed map error
by assigning both measured values and their associated predictions to classes and then
determining the percentage of measured values that were assigned to the same class as

their predictions.

Gotway et al. (1996) used mean square error (MSE) as a measure of prediction
accuracy calculated using independent validation data sets. Map goodness was calculated
by comparing the MSE for interpolation to the MSE for the field average approach. To a
large extent, the appropriate grid spacing depends on the spatial structure of soil
properties (Sadler et al., 1998) and the range of spatial correlation (Mohamed et al.,
1996). Furthermore, the appropriate interpolation method may depend on specific
coefficients of the interpolation procedure used. For example, the optimal distance
exponent for the inverse distance weighted (IDW) interpolation procedure may depend

upon the coefficient of variation (CV; Gotway et al., 1996).

While grid sampling and interpolation approaches have limitations, a shifi from
grid sampling to directed sampling schemes or grid cell based approaches may not

improve map accuracy. Directed sampling is a technique wherein samples are collected

and composited from specified areas within a field. Pocknee et al. (1996) suggest that the
areas be delineated based on established differences within a field. Differentiating
criteria might include soil map delineations (Bell et al., 1995, Moore et al., 1993,
Windawa et al. 1993), management history, yield potential maps, aerial imagery
(McCann et al., 1996; Pocknee et al., 1996), or electromagnetic induction (Jaynes, 1996).
Directed sampling techniques require prior knowledge about the factors that regulate crop
yield and nutrient availability. In general, map or prediction errors associated with
directed sampling schemes have not been reported. Instead of sampling on grid points,
composite samples can be taken from the area between grid points, a practice that is
referred to as grid cell sampling. Grid cell sampling schemes have not found extensive
use in SSFM. This may be related to earlier reports that grid point sampling more
accurately described soil properties than did grid cell sampling, such as that reported by

Wollenhaupt et al. (1994) for two fields in central Wisconsin.

The purpose of this study was to evaluate how grid sampling and interpolation
schemes affect map accuracy based on measures of bias, precision, and prediction
efficiency, or what Gotway et al. (1996) termed goodness. For this study, a field was
subjected to several soil sampling strategies including grid point sampling at several grid
spacings, grid cell sampling, and directed sampling based on soil type. Soil fertility and
fertilizer recommendation data were interpolated using kriging and a range of IDW
coefficients for the various sampling schemes. Each resulting prediction map was tested

against a random validation set to evaluate map accuracy.

MATERIALS AND METHODS

Site description

This study was conducted within a 20.4-ha field (42° 57' 54" N, 84° 43' 38 W) in
Clinton County, Michigan, 6-km south of Fowler. The field has been in a corn (Zea
mays L.)-soybean (Glycine max L. (Mcrr.)) rotation for 22 years and the southeastern
portion of the field has been sub-irrigated since 1988. The field was selected because it
contained multiple soil map units and exhibits a range of terrain features, both of which
would be conducive to SSFM. Pregitzer (1978) described and mapped the soils in
Clinton County, MI at a scale of 1 :15,840. The great group taxonomic classifications of
the soils in this field were either Ochraqualfs (Capac and Metamora) or Hapludalfs
(Morley and Wasepi; Table 1.1). The soils were somewhat poorly drained except for the
Morley, which was moderately well drained. All of the soils formed in glacial till except

for the Wasepi, which formed in loamy glaciofluvial deposits.

Table 1.1 Map symbols, soil series or complex name, NRCS soil taxonomic description
(Pregitzer, 1978), and area occupied in the field

 

 

Symbol Name and slope Taxonomic Family and Subgroup Area
(ha)

MoB Morley loam Fine, illitic, mesic 1.3 ‘
(2 to 6 % slope) Typic Hapludalfs

CaA Capac loam Fine-loamy, mixed, mesic 7.2
(0 to 4 % slope) Aeric Ochraqualfs

MeA Metamora-Capac sandy F ine-loamy, mixed, mesic 8.9
loams Udollic and Aerie Ochraqualfs
(0 to 4 % slope)

WbA Wasepi sandy loam Coarse-loamy, mixed, nonacidic, 3.0
(0 to 3 % slope) mesic Aquollic Hapludalfs

 

10

70

600

500

Northing (m)
A
O
O

300 If

 

Easting (m)

Figure 1.2. A digital orthophotograph overlain by the 30 m (black square), the 100 m
(blue cross), the random validation data set (green circle) and NRCS soil
types (CaA = Capac loam with 0 to 4 % slope; MeA = Metamora-Capac
sandy loam with 0 to 4% slope; MoB = Morley Loam with 2 to 6% slope;
and WbA = Wasepi sandy loam with 2 to 6 % slope).i

T red solid and dashed lines indicate boundaries of directed sampling zones (DS 1-
11) and dashed yellow lines indicate the boundaries of the CELL based sampling.

11

Soil Sampling Design and Laboratory Analysis

Soil samples were obtained from the field using four sampling designs: a 30.5-m
regular grid, a 100-m grid cell, a 200-m unaligned grid, and random sample design.
These samples were used as data sets or to calculate data sets as described in Table 1.2.
The sample point locations were flagged using a DGPS system with a base station for on-
the-go differential correction. At each grid and randomly distributed point, 5 sub-
samples (1 at the grid point and 4 within a 1.5-m radius) were obtained to a depth of 20-
cm using a 2.5-cm diameter core and composited. The grid cell samples were taken by
compositing 9 individual cores of the same diameter and depth taken at regular intervals
within each 100-m grid cell (Figure 1.2). Soils were dried under forced air at 35° C for 3
days and ground to pass a 2-mm sieve. Standard soil analyses were conducted by the
Michigan State University Soil and Plant Nutrient Laboratory using the recommended
chemical soil test procedures for the North Central Region (Brown, 1998). Analyses
included pH (1 :1 soil water mixture), BpH (SMP buffer), P (Bray P-l extractable) K, Ca,
and Mg (lM/L NILOAc extractable). Cation exchange capacity (CEC) was calculated
by summation and lime (an), P (Pm), and K (Km) fertilizer recommendations were
calculated using the tri-state fertilizer recommendations (V itosh et al., 1995) for corn

with a uniform yield goal of 11.3 Mg ha'l (180 bu acre").

12

 

Table 1.2. Data set descriptions

 

 

 

 

 

 

 

 

Name N Description
Prediction Data Sets
G30 215 Soil samples taken on a 30-m regular grid.
G100 24 Soil samples taken on a 100-m regular grid.
CELL 15 Soil samples taken in IOO-m grid cells
Gcomb 239 Created by combining G30 and G100 grid data sets.
Gel-a 54 Created by separating the 30-m grid data set into four separate
Gal-b 54 6l-m grids. Each grid was analyzed separately, but the results are
G6i-c 54 given as the average results from the four grids.
Gm-d 53
DIR 11 Created by overlaying the soil map units onto the 30-m grid and
determining the mean fertility value for each management unit.
A30 1 Field average value for the entire field calculated from the means
of the 30-m grid.
A100 1 Field average value for the entire field calculated from the means
of the IOO-m grid.
Validation Data Set
VAL 62 Independent validation created by combining a 200-m unaligned
grid with additional random samples.
Data set for geostatistical analyses
FULL 301 Created by combining G30, G100, and VAL.

 

13

Data Analysis

Data analysis was performed in two phases. Phase I involved a quantitative
analysis of the FULL data set, which consisted of the combination of the G30, the G100,
and the random (VAL) data sets (Table 1.2). The purpose of Phase I was to assess the
extent to which the spatial variability of soil fertility in this field lends itself to SSF M.
The steps in phase I included an assessment of normality, calculation of descriptive
statistics, and analysis of spatial variability. Phase II was an assessment of the accuracy
of sampling designs and interpolation procedures using quantitative measures of map

quality.

Phase I: FULL data set

The FULL data set consisted of the 301 sample locations corresponding to an
average sampling intensity of 14.8 samples ha'l (Table 1.2). Normal probability (Q-Q)
with 95 % confidence intervals (Friendly, 1991) were used to assess normality of the
FULL data set. When the Q-Q normal probability plots indicated large deviations from
normality, the natural log of these variables was calculated and the resultant tested for
normality. If the log transformations were also non-normal, then the original variables
were power transformed and again tested for normality. The power transformations were
not successful in inducing normality and will not be discussed further. Finally, normal
score transformations (Deutsch and Joumel, 1998; Goovaerts, 1997) were used to
normalize the remaining variables that could not be normalized with log or power
transformations. Contour maps of variogram surfaces were created for each original

variable to determine the direction of the anisotropic axes if anisotropy existed

14

(Goovaerts, 1997; Isaaks and Srivastava, 1989). For directional (anisotropic)
semivariograms, an angular tolerance of i 40° was used because it allowed the
variograrns to be well defined while still preserving their essential features (Isaaks and
Srivastava, 1989). Nested semivariogram models (combinations of spherical,
exponential, and/or Gaussian models) were chosen based on their fit to the empirical
variograms, if warranted. Using the modeled semivariograms, GSLIB (Deutsch and
Joumel, 1998) was used to create kriged 4 by 4-m grids at search radii equal to the
distances to which the semivariograms were modeled. Contour maps for each kriged grid

were created with Surfer® (Golden Software, Golden, CO).

Phase 11: Accuracy of SSFM sampling and interpolation

Phase II was concerned with how different sampling schemes and estimation
procedures affected map accuracy. For the FULL, G30, 661.4, (3100, and Goon“, grid data
sets, empirical semivariograms were calculated using Variowin (Pannatier, 1996) and an
omnidirectional exponential model was fit to the empirical semivariograms. Surfer was
used to interpolate 4 x 4 m grids by kriging each data set with the modeled
semivariogram and using IDW for distance exponents from 0.1 to 5.0 in 0.1 increments.
For comparison purposes, omnidirectional semivariograms based on an exponential

model were also developed for the FULL data set.

To quantify the error of prediction for each attribute, sampling scheme, and
interpolation method, the difference between the predicted surface and the validation
points were estimated. Since the VAL points did not always coincide with the 4 x 4 m
predicted grid locations, bilinear interpolation in Surfer was used to estimate the

predicted value at each VAL grid location. Because only one value is assigned to each

15

soil management unit, each cell, and to the entire field, residuals for the DIR, CELL, A30,
and A100 predictions were calculated as the distance between the VAL points and the
measured value for the area containing that point. Prediction errors and efficiencies were

calculated as described in Chapter 1.

RESULTS AND DISCUSSION:

Phase 1. FULL Data Set.

While the average soil tests for this field indicated that soil fertility was adequate,
there was considerable range in each parameter (Table 1.3), indicating that some portion
of the field may respond to SSFM. For SSFM to be applicable to this field, the variation
of soil fertility must be spatially structured, of sufficient magnitude, and within the
manageable range. Semivariance analysis was conducted to quantify the extent to which

these conditions were met.

Table 1.3. Summary statistics for the FULL data set.

 

 

Variate mean Median MinT MaxT CV (%)
pH 6.1 6.0 5.2 8.0 9
P (mg kg") 26 25 8 86 44
K (mg kg“) 167 165 95 291 23
Ca (mg kg") 1482 1450 650 3053 22
Mg (mg kg") 274 267 80 474 23
CEC (cmole kg“) 13.0 13.0 7.7 20.5 17
L...c (Mg ha") 4 5 0 13 76
Free (kg ha") 63 75 0 114 42
Kmc (kg ha") 25 o 0 100 136

 

1 Min = minimum, Max = maximum.

16

While normality is not a requirement for developing a semivariogram or kriging,
the classical linear kriging predictor does not retain optimal properties if the underlying
spatial process is not Gaussian (Cressie, 1993). Only Mg and CEC were normally
distributed, while P, K and Ca were log-normally distributed (Figure 1.3, Table 1.3). Soil
pH, Lm, Pm, and Krec could not be transformed to normality with log or power
transformations but could be through a normal score transformations (Figure 1.3). The
reason Lm, Pm, and Km deviated so drastically from normality is because the tri-state
lime and fertilizer recommendations (V itosh et al., 1995) combine stair stepping, nested
functions. Because of the nature of these calculations, however, the back transformations
for the normal score transformations of LM, Pm, and Km were problematic and therefore
should not be used for SSFM of this field. It may be that native Mg and CEC levels were
distributed normally and remained so because lime applications have been primarily
calcitic and CEC has minimally been affected by management practices. In another geo-
spatial study in Michigan, soil properties in an uncultivated landscape appeared to be
distributed more normally than those in an adjacent cultivated field (Robertson et al.,

1993).

Anisotropy occurs when the semivariogram depends not only on separation
distance, but also on the angular relationships between data points. The presence and
direction of anisotropy must be known to create anisotropic models and can easily be
detected and measured using contour maps of semivariogram surfaces (Goovaerts, 1997).

Only CBC and Lrec were considered to be isotropic (Figure 1.4). Directional or

17

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I 7 I 'J//4 In
a 6 . 2 0 q
5 .1 -3 .4) ‘J .
at; 80 _‘ fl. 4 ...,: 4 d /-
g 40 - “/4 a‘ 2 3 4
0 - g n
T: 300 - /J a: 6 “ ,x
x E 150 I Q ‘2 /.
, 3;“ a
V 0 ‘ i? -
“A 1 5':- 4 g /
3 '3. 3000 - “.3; 7:. /-/e
g 1500 4/ E g 8.0 1
V v a ‘
o .1 1 E E. L I 1
Je- 500 4 3:6
2 g 250 d/
V o q u
o; 21 - _4
8 0 14 4
i 7 .1
V 1
v-A 15 .1 2?: 3 -4 I“ i
§ ‘2 0 d/ :8; 0 /
" 2" “2 d
v -15 _ -3 _’ .-
e" 80 - _r' g 3 r
g8 '2 40 4/ If. 0 ‘ /
(L g ..- g,
V 0 4,...__.' -3 g,
w: 200 .. 47/18 3 -1 . ~ 1
g f, 0 F/ x 0 '1 /
x 2
‘“ -200 . -3 .. "
-3-2-1012 3 -3-2-1012 3
Normal Quantiles Normal Quantiles

Figure 1.3. Normal probability (Q-Q) plots for variables and log transformations with
95% confidence intervals.

1 ns indicates a normal score transformation; pH already log transformed

18

omnidirectional semivariograms were calculated and modeled for each variable and their
normal transformations (Figure 1.5). The semivariogram models accounted for geometric
(range changes with direction, sill does not) and zonal (sill changes with direction; range
does not) anisotropy and mixtures of both (Isaaks, and Srivastava, 1989). The directional
and non directional variograms were described with one, two, or three nested spherical,
exponential, or Gaussian transitional structures in addition to a nugget structure (Figure
1.5 and Table 1.4). Overall, anisotropy was not severe and, removable with
transformations (pH and Ca). It would not be cost effective for anisotropy to be modeled
for this field on a commercial basis because the anisotropy is not strong, its modeling
time and resource expensive, and the coarseness of commercially accepted SSFM would

not allow the short range anisotropy to be resolved.

The kriging system of equations requires that a second order stationary model be
fit to the empirical semivariograms. Second order stationarity can be inferred the
semivariogram reaches a plateau as occurred for most variables (Figure 1.5). When
semivariograms do not reach a plateau, intrinsic stationarity is ofien assumed (Pannatier,
1996). But, by limiting the search radius, a stationary model can still be applied. The
semivariograms for pH in the N 54° E direction and and Ca in the N 49° E direction
exhibited intrinsic stationarity in these directions. This behavior was attributable to a
"hot spot" in the southwestern comer of the field where lime had been stored (Figure 1.6)
rather than a gradual trend in the data. Normalizing the pH and Ca data removed the
directional, intrinsic behavior. Second order and intrinsic stationarity was assumed for
the FULL data set and issues of non-stationarity are irrelevant for SSFM management of

this field.

19

pH

N$°\N

 

North-south direction (m)

 

-200'

 

Mg cec

N91°E

   

-200 0 200 -200 0 200
East-west direction (111)

Figure 1.4. Directional variogram contour maps for the FULL data set. I

T white and black represent the minimum and maximum semivariogram,
respectively of the FULL data set (not transformed).

20

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

"5(pH) .
0.7 -
/ Omnidirectional
0.0 - ‘
ln(P)
0.15 "‘ . N1go‘wo
0 00 y . N 71 ° E
. ln(K)
0.04 -fiflz‘:: N 00 W i
. N 90° E
0 . , 0.00 1
Ca ln(Ca) . o 0.-
80000 ~ ‘ N . E 0.03 g
8 N w Omnidirectional
5 0.00 .
'c
(U
.2
E
<13
0 CEO
0 1
3.5 ‘f—‘fifiu—f‘
Omnidirectional
0.0 m
Lroc ns(Lm) .
V
7 - 0.7 4 ’ '
Omnidirectional / Omnidirectional
0 . 0.0 1
PM , . ns(Pm)
mufflw. 1‘ . 0°w
- N 71 ° E - N 90° E
0 L 0 1
Km 4 a Ans(KnJ
n ‘ ‘ ‘
8501M ‘ 1“: ::.o.o O
M ‘ N0°W no ". N09‘w..
' N 90° E . o
o I 0 IN 90 E
0 175 350 0 175 350
L39 (m) Lag (m)

Figure 1.5. Directional and omnidirectional semivariograms for the original and
transformed FULL data sets. I

1' ns indicates a normal score transform.

21

How—o Tr 980302: Ea 03393283— moamfiiomama 580— $385688 men 9a mCrr 38 m9

 

 

 

 

 

man—8:6 _ mag—2:3 m
5533 2.me Zeno—4 a: 9832.4 25mm _wm<H 2.0%— mE E398: Rama
1%

z 3.. ’2 we 2 mac 2 38

U: Cbhm m Obccoc .2 mg. m m¢ mm 0 Ta .2 who m _Omc
3%: obc m 93 O .8 S O c. _ u 0 m8
241m 3 Ziom as
_u wk» m wg 2 _co 2 _o 0N m was 2 $6 ’2 _we
zdom 3 zflom 3c
5 3 ob. m o._ _m 2 3.. $4 q 3 m Pom 2 So It mac
z 8.. m .3 z 8.. m «8
_A ”Ewe m mama 2 co m 3 Q 3 _m z 9. new
z 8.. m mo Z 8.. m m3
_= :3 982 o 0.23 2 co 3 A: 0 ob; 2 co _:

 

ca m .1. mvronomr m n 9683:? Q n Guzman? O n ogmmmaomoawr am Samoa—Sm m nonsu— moono ”35.333398:

H ”m5 magaos N. 8.9—.83 comma—£920? 55m :5 a: men 90 mam" 3:088 on?

22

edge ..A 00:33:09

 

 

 

 

 

madnER _ mango N @2883 will
$53.6 _Emmo. Zen—6.4 mm: finance—.4 Rama _wm<H Zona— a: cacao: Rama Zena. .2: E398: Sumo
axe
226$ mu 251$ 3.3
On 38 m 288 2 3.. m S. we 0 $988 2 3o m 53
533 o m 9?... O 8. _oo Q 92: 0 3c
2 3.. m no 2 8.. m Zoo
in o m 33o 2 co 3 _oo 0 $3 2 co :3
0mm T: m u. 73 O 3 3
~42 Nmam m 939 O .3 d
airs 93 m cum O 3 ac
Ziom mm Ziom moo
Pa _3 Q 23 2 3.. é _o oo o N: 2 Go 5. NS
28%. c 280m :N 286 So
are... 93 m 8 z o. 3 a o 93 z c.- t o 93 z 9. wow
z 8.. m d Z 8.. m muo
Kan So Q So 2 co _ w u. D Nmm 2 co .3

 

23

Northing (m)

 

1:
high MeA & w N®
slightly high CM

“WW

optimal
slightly low
low

 

 

 

very high
high
optimal

slightly low
44
very high
170
slightly high

optimal
120

3000
2000 T
1400

adewate

   
  
  
  
  

 

‘7oo

mar inal
350 9

464

ads ate
63 _L_

 

marginal

 

East‘ng (m) Easting (m)

Figure 1.6. Contour maps of kriged soil variables and fertilizer recommendations using

the FULL data sets overlain by the soil type boundaries *.

‘i' CaA = capac loam with 0 to 4% slope; MeA = Metamora-Capac sandy loam
with 0 to 4% slope; MOB = Morley Loam with 2 to 6 % slope; and WbA =
Wasepi sandy loam with 2 to 6 % slope.

1: indicates the soil types symbol designations

24

The RSVs for the first structure were greater than 50% with exceptions of K,
ln(K) and Km so most of the variables could be described as being spatially structured
(Figure 1.5; Table 1.4). The directional semivariogram models had ranges of spatial
correlation of 34-m or greater in one of the two directions and the omnidirectional
semivariogram models had ranges of at least 43-m. The degree to which these

parameters affect the accuracy of SSFM predictions is of great interest but unknown.

The interpolated maps show that considerable areas of the field were low or
slightly low in pH and slightly low in P. However, the entire field had greater than
optimal K levels and greater than adequate Ca and Mg levels (Figure 1.6). As evident in
these maps, the variables did not relate well with soil type. The apparent discrepancy
between the ranges in Table 1.3 and Figure 1.5 is due to the smoothing effect of kriging
that occurs with a non zero nugget effect (Goovaerts, 1997; Figure 1.5). Only when the

nugget effect is zero does kriging behave as an exact interpolator.

Soil pH, P, and K had significant variability, the variability was spatially
structured, and was in the manageable range. The variables were either normal or could
be transformed to normality with log or normal score transformation. For some
variables, the transformations also removed anisotropic features including directional,
intrinsic stationarity. Unfortunately, normal score back transformations were problematic
because of the nested, stair stepping features of the fertilizer recommendations. Because
Lac, Pm, and Km, cannot be back transformed and deviate severely from normality, their
SSFM may be difficult. Based on the presence of spatial structure and the fact that the
variability was within the manageable range, there is potential for SSFM in this field but

to determine if variability is adequately known, interpolations must be evaluated.

25

Phase II: Accuracy of SSFM sampling and interpolation

Grid Sampling

Grid sample design had an effect on semivariograms for all parameters (Figure
1.7). Anisotropy was not modeled for the reasons listed in the Phase I analysis. The
semivariograms for the FULL (average grid size of 26.0-m) and the Gcomb (average grid
size of 29.2-m) data sets were similar. Because of the similarity, the test of the Geomb set
against the VAL set discussed later should be indicative of the performance of the FULL
data set. As sampling density decreases, the semivariograms deviate from those of the
FULL data set. The semivariogram for the G30 data set has a higher nugget variance and
range but a similar sill. The semivariograms for the four Ga and the 6100 data sets

deviated greatly from the FULL data set and from each other.

For kriging, the RMSE of the residuals of the predicted grid versus measured
(VAL data set) generally increased as sampling intensity decreased with Goomb < G30 <
G61< 6100 (Figure 1.8). A spatial sampling scheme with utility should improve
prediction over the whole field sample approach. The A30 field sampling scheme had the
lowest RMSE for K and Km, which may relate to the fact that soil test K levels were high
(low Km) and both K and Krec had low RSV (Figure 1.7) or high nugget effects. It
should also be noted that K was the only variable with substantial bias. The kriged 6100
data sets had similar RMSE values to the field average approach for the remaining
variables. The fact that the A100 field average had higher RMSE values for some
variables than the A30, illustrates the importance of obtaining a good estimate of the field

mean, particularly in variables that are poorly spatially structured like K and Km.

26

 

 

 

 

Gcomb

 

FULL

 

 

 

 

 

 

 

/

855283

175
ponential) semivariogram for

175

175
Lag distance (m)

27

1i5

 

. -iV|.
u_

 

 

I:
wllilli ii.iiiii Vi. l. .

1%5

 

the FULL, Gm, G30, G6,, and Gm data sets. T
t The semivariogram models for 661 are not shown.

Figure 1.7. Omnidirectional experimental and modeled (ex

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mmim

g gggggg mafia.
in“ WE: _. .- w "3:2 E mmfihwzjfifi - SE... jwo
E E SE SE SE, mam, 552:2: 53.. 8..
flgagggfl a...
E _HDDHD EH 3.13444 ages...
a, a a a VII, VIII 7! ..o as.
E QM RN NNM E S, RN .0 so...
a a mm mm a. a, @503...
a H __ H_ _. _. _ L L a
._ ._ . __ .._ ._. ._ . _ ..4
._ ._ _ ._ ._. . _ . at
._ H. . ._ L. ..p .. ._ .3 mm
__ .. ._ L c ,1. a L . mm
I. ..-
k _. . .\ . x _ .. . _ -

_. L . . _ . L . ._. 0
Nu. W m 0 “v 0 . M 2 0 w 0 w 0

A29. 95 A29. 95 ...9. as. $9. 95 if .23 A22 9: A22 9; A22 9:

_ IDW cs30
_ - IDW (3,,
low Gm

Figure 1.8. The RMSE for several estimation approaches: IDW interpolation for distance

exponents ranging from 0.1 to 5.0 in 0.1 increments (line plots on the left),

kriging, whole field, and grid cell sampling; and zone sampling.

28

 

Although kriging is often considered to be the most accurate interpolator,
commercial SSFM applications commonly use IDW interpolation in lieu of kriging
because it does not require modeling the semivariogram. Like kriging, the RMSE for
IDW interpolation was affected by grid increment, but was also affected by the choice of
the IDW distance exponent G’igure 1.8). Generally, RMSE decreased with increasing
sampling intensity G30 < G61 < Gm with the exceptions of CBC and L“... The RMSE
decreased markedly for P and Prec with an increase in sampling intensity from G61 to G30
while the RMSE for K and Krec was reduced minimally. The RMSE decreased modestly
for all variables with an increase in sampling intensity from G 100 to G51. The optimal
distance exponent depended on the variable and sampling intensity. Consistent with the
findings of Gotway et al. (1996), the value of each distance exponent was inversely

related to the CV but only for the G30 and G61 data sets (Figure 1.9).

Measures of map error and prediction are indicators of map or prediction
goodness but can not be completely understood without graphically comparing predicted
and measured values along a 1:1 line (Figures 1.10). At the 6100 sampling intensity, this
relationship was generally weak for most variables. At the G30 sampling intensity, the
relationship between predicted and measured values was moderate for P, and Free and
marginal for the others. Even though, for Goomb, the relationship between predicted and
measured markedly improved over this relationship for G30, their regression lines were

not parallel with the 1:1 line.

29

 

 

 

 

 

 

pH
G30
\ Mg
\Ca Free
2 ..
CEC \ lime
p
\
1 - K \
\
\
y=2.3-0.011 * CV K’e‘
...- r2=0.48
E, 0 l l l A
a G.
x
0)
§ 2‘
0
§ \‘01 o
‘5 Q \\
3 1 ‘ .‘\\
e 0 \..
'5 y = 1.55 + .0055 * CV
,g r==o.41
@- 0 1 1 l i
G100
2 ..
. I
o. iz”
1 f i '

O
y=1.18+.0035*CV
r’=0.07

0 I I I l
0 50 100 150
CV (%)

Figure 1.9. The relationship between the optimal distance exponent and the CV.

30

 

Gcomb

 

Predicted

 

 

~K' LCEC

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

       

 

 

 

 

 

G30
'0
O
72"
E
a
G100
8 :. .
E o
E Z 1330.10 ’/ 330.00 / 330.17
Kri-ODO/ CEC .-/ . /
- .I . . .. ' . '
i ' // €20.02 '7 R'e'ooo
Measured

Figure 1.10. Predicted vs. measured for kriging with the Gwmb, G30, and 0100 data sets."r

‘l' Dashed lines are the 1:1 lines. Solid lines are the regression lines.

31

Sadler et al. (1998) and Mohamed et a1. (1996) suggest that the performance of
SSFM at different grid intensities is a function of the spatial structure and range of spatial
correlation of the spatial processes. The relationship between prediction efficiency, the
range of spatial correlation, the RSV, and grid sampling increment are described in
Figure 1.11. At each grid intensity, when the RSV values were less than 60 %, e.g. for K
and Km, prediction efficiencies are low regardless of the fact that their ranges of spatial
correlation were relatively high (greater than 90-m). For the G30 data set, with the
exception of Mg there was a quadratic relationship between prediction efficiency and the
range of spatial correlation when RSV values exceeded 75% (Figure 1.12). While this
relationship held for the Ger data set, it does not exist for the 0100 data set. Therefore,
inferences between geo-spatial studies and the expected performance of SSFM can be

drawn but only when grid sampling intensities are adequate. _

Using the optimal distance exponent, IDW was equal to or superior to kriging for
most of the variables sampled on the 30, 61, and IOO-m grid excluding log-transformed
attributes (Figure 1.13). Unfortunately, there is no a priori knowledge of the optimal
distance exponent. Visual inspection of Figures 1.13 suggests that a distance exponent of
1.5 would be a reasonable choice for this field. When the regression equations in Figure
1.12 were used to calculate the optimal distance exponent, the prediction efficiencies for
the IDW approach were nearly the same as or slightly better than when the distance
exponent was 1.5 with the exceptions of K and Krec (Figure 1.13). How the relationship
in Figure 1.12 holds for other soils or locations is not known. For both kriging and IDW
interpolation, prediction efficiencies were lower for ln(P), and ln(Ca) than for P, and Ca.

Although normalizing K improved the prediction efficiencies, they were

32

 

 

 

 

 

 

 

 

 

 

 

 

l
. l
l i
/ : l I
\.\ l l
40 / // . :5 i
/ ' :
Ex; /\e~\ f
:1 / a \i\_
8 2° / ; l
.:._ /\l\ l l
g // \i\“ ‘
9. 0 ,/ \ . \-
5 Z / 7\.-7L\4
E “’K\ \ / 100
-20 / \ 80

{\3
. / 7\ .
120 /

Figure 1.11. The relationship between prediction efficiency for kriging the G30 data set,
the range of spatial correlation of the FULL data set, and the RSV of the
FULL data set.

33

 

 

 

 

 

 

63° /o’ ‘\
40 4 / .
0/
/ outlier
20 ~ ..9/
/ y=-75.l +2.41 x
/ - 0.0119 x2
0 - o r2=0.97
.3 GG‘ 0 //“\\
> / .
° / .
QC) 10 _, / outlrer
'§ 0
f) / Mg
c /
3% g/ = -34.9 + 1.06 x
‘6 o . , - 0.00534 x2
93 ' r2=0.74
o- 1 1 l
(3100
20 2
o
o o
0 ‘ o
l
-20 ‘ .
o
40 ~
0
'60 f I T
40 80 120

Range of spatial correlation

Figure 1.12. The relationship between prediction efficiency for the kriged G30 data set
and the range of spatial correlation for the FULL data set for variables with
RSV values greater than 70%.

34

 

 

 

 

 

 

 

 

30 ‘ /"i
/’o’
/
0 - // to
/ to
to
I//
yfi
30 .1
.E’ :z
F” /
i2 0 - ./ o
r. / 0
£2 0
E /
.‘2 / O
E 30 ‘ e).
9’ ,/d'
C /
23 O<- / o
.2 /°
8 ' °
5 L l
/
/o¢
30 ~ .9”
,/,o
/
0 — / O
or
/
/’e'
30 ~ I,”
/<V°
0 ~ / 0
/'<0
0 30

 

 

 

 

 

 

 

 

0 ‘ 3'.
l'{// .
~15 - /
15 - ’/
e,
0 ‘ /dfio
/ O
‘15 " /”
./
15 r /
0 ‘ 2’.
-15 n //
15 4 9'
.> .15
0 ‘ oi/g‘
° ,/
-15 - /
15 - j.
0 . °/
0 / O
’/
-15 -i /
/
-15 0 15

 

 

 

 

 

 

 

 

6100
J 0”
o 4
///"
l/
-30 ~ / .°
. 3
2‘
0 q
/ 0
4
O/
-30 - /
. . O
.2:
0 .i
/ 0
4
q 0/
-30 . /
. . o
’9
o _
A.
4b
or
-30 - /
1 1 O
‘V
o - 6.
/
/ O
/
o
e
'30 q /
I r .
-3o 0

Prediction efficiency for IDW

Optimum
IDVV
Exponent

IDW
Exponent
= 1.0

IDW

Exponent
=1.5

IDW
Exponent
= 2.0

IDW
Exponent
cakuuahxi

flonr

regnmunons
in
Hgme111

Figure 1.13. Prediction efficiency for kriging versus the prediction efficiency for IDW.'*

1' Prediction efficiencies are for the G30 data sets.

35

negative for K and ln(K), indicating that the field average approach was superior.

Normalizing the data had an insignificant or negative effect on the accuracy of the data.

The poor performance of kriging may be explained by the large variability of the
semivariograms for individual pairs (Figure 1.14). A fitted semivariogram model to the
average semivariograms would be accurate for a small fraction of the pairs. So for a
kriged estimate, two sample points the same distance from the prediction point will be
weighted the same even though the variances between the prediction points and the point
to be estimated (if it were known) would vary widely. It may be that the scatter of the
semivariogram cloud may be one of the best indicators of the spatial predictability of a

given variable.

CELL and directed sampling

Alternatives to grid sampling include cell sampling and directed sampling. The
CELL sampling approach had the highest RMSE for most variables (Figure 1.8) and the
lowest prediction efficiencies, confirming the results of Wollenhaupt et al. (1994). The
directed sampling approach had RMSEs that were similar to the A30 field approach
(prediction efficiencies near 0), indicating no advantage of directed sampling over the use

of a mean value for this field (Figure 1.8).

36

Semivariogram

14000

 

12000 ~

10000 a

8000 q

6000 -

4000 2

2000 e

 

 

O
o
o o
o o
O
o o
oo
o
o
o o
o
o o
o o
o
o
o o
oo
o o
o
o
o
o
o
oo
0 0
oo
o o
o o
o
°8
o o
o o o
o 8
o o 0 0°
0 8 g
o o o
o o o
o
o o
a 03
<2. °°§
o
o°8 °§°
o 808 g °o 8
8 o a °8§
8
o o 88 0°
0 3 8 8
0 0 g 5 38°
0
oo 0 3°88)
fig 3

 

Separation distance

Figure 1.14. Average semivariograms (large solid circles) and semivariograms for

individual pairs (hallow circles) for P using the Goomb data set.

37

CONCLUSIONS

Analysis of the FULL data set indicated that most soil fertility variables were
spatially structured. Nevertheless, the presence of spatial structure alone did not prove
sufficient for producing accurate yield maps, as evidenced by the plots of measured vs.
predicted in Figures 1.10. Sampling at lower intensities increasingly diminished the
delectability of spatial structure and generally increased the error of prediction as
measured by RMSE. Where spatial structure was poor, particularly for K and Km,
accurately sampling the field average was sufficient for nutrient management because
SSFM for these variables was not appropriate. These data suggest that grid sampling at
coarse grids and directed sampling were not adequate to produce accurate nutrient
condition maps for this field. Cell sampling at least at the course 100-m grid intensity

was also inadequate.

These data suggest that grid point sampling at the industry standard 100-m
intensity was inadequate. Sampling at greater intensities only modestly improved
prediction accuracy, likely not enough to justify the geometric increase in sampling costs.
In the second chapter of this dissertation, I will examine methodology for incorporating

secondary landscape information into spatial estimates of a soil property at several scales.

In this study, the accuracy of IDW interpolation with a distance exponent of 1.5
generally equaled or exceeded the accuracy of kriging at each scale of measurement. If
the data had been strongly anisotropic or was not second order stationary, kriging may

have been superior to IDW interpolation because the semivariogram model could account

38

for these peculiarities. The poor performance of kriging may be explained by the large
variability of the semivariograms for individual pairs (Figure 1.14). Some measure of the

scatter of the semivariogram cloud may be an indicator of the predictability in space.

39

CHAPTER 2

COMPARISON OF TECHNIQUES TO OPTIMIZE SPATIAL ESTIMATES OF
SOIL PROPERTIES USING TERRAIN ATTRIBUTES

INTRODUCTION

Traditional survey methods and the more recent use of grid sampling and
interpolation methods have not produced maps of soil properties with the accuracy
needed for soil surveys (Bell et al., 1995; Moore et al., 1993) and precision farming
(Pierce and Nowak, 1999; Robert, 1993). New analytical approaches are being used to
utilize geometric properties of a landscape (slope, aspect, and curvature), collectively
referred to as terrain attributes, to improve spatial estimates of soil properties. Terrain
attributes are predictive of soil properties because topography is a soil forming factor. A
high resolution, digital elevation model (DEM) is needed to calculate terrain attributes.
Only recently has the technology become generally available to map elevation at the
needed resolution, achieved through advances in high resolution global positioning
system (GPS) now universally available. Several methods exist to use terrain attributes
in spatial estimates of soil properties, ranging in complexity from simple regression to
geostatistical methods. However, there is little consensus regarding which terrain
attributes are most useful or which analytical method is most appropriate for a given soil

property (Table 2.1).

40

Table 2.1. Applications of regression, kriging with an external drift, and cokriging for
enhancing spatial estimates of soil properties using auxiliary information.

 

 

 

Primary Secondary Measure of Reported results 1'
accuracy and or
goodness 1
Regression
Moore et al., A horizon depth 4 terrain attributes Visual Regression
1993 (R2 = 0.51) comparison and approach was
soil P (R2 = 0.48) the regression R2 considered good
soil pH (R2 = 0.41) values because terrain
particle size (R2 = 0.64) attributes could
explain substantial
43 points hal 43 points ha'l variability
Bell et al., A horizon depth 3 terrain attributes VDS; plots of A- horizon and
1995 (R2=0.51) predicted versus depths to carbonate
depth to free carbonates measured predicted within 20
(18:04:) cm for 70% of
validation samples
7.25 points ha’1 75 points he"
Thomson and Total carbon 2 terrain attributes None Unclear
Robert, 1995 (R2 = 0.66 to 0.69) and photographic
tone
2.7 points hal 30 points ha’l
Gessler et al., A horizon depth 2 terrain attributes Prediction error Regression reduced
1995 Solum depth (not specific) deviance by 63 and
68%

scale not given scale not given

 

Bourennane et
al., 1996

Gotway and
Harford, 1996

Goovaerts,
1997

Thickness of silty-clay- one terrain attribute

loam pedological

horizon

0.62 samples ha'l 4.8 samples ha'I
residual soil NO;l corn grain yield
11 samples ha" 66 samples ha"
Soil Cd, Cu, Pb, and Zn blocked

Co estimates???
(259 samples per field (359 samples per
area) ' field area)

Kriging with an external drift (KED)

 

VDS; ME and
RMSE

improvement over
OK and KT

Cross validation, COK increased the

compared MSE MSE by 7% over
for OK and KT OK

with the MSE

for KED

VDS; rank Correlations
correlations for between predicted
predicted and and measured for

measured and % COK explain 16 to
misclassifrcation 35% more of the
for both SK, OK variability than for
and KED OK

 

41

Table 2.1, continued.

 

 

 

Primary Secondary Measure of Reported Results
accuracy and or
goodness
J‘vokriging (COK)
Zhang et al., 1992 N03" and Ca electrical Cross validation, cokriging reduced
(COK with pseudo- conductivity compared MSE the MSE by 78%
crossvariogram) for OK with MSE
0.8 samples ha" 1.3 samples ha'l for COK
Vaughan et al., Water content and surface electrical Visual inspection improvement over
1995 soil salinity conductivity ordinary kriging

0.12 points ha'l 0.15 samples ha'l

Rosenbaum and Soil Cd Soil Zn Independent Correlation
SOderstrOm, 1996 validation data set; between predicted
(standardized 0.00034 points ha‘l 0.0010 samples ha“n correlation of and measured was
ordinary cokriging) predicted and greater for COK
measured (p = 0.85) than
OK (p = 0.68)
Gotway and residual soil N03 corn grain yield Cross validation, COK reduced the
Harford, 1996 compared MSE MSE by 2%
10.8 samples ha'1 66 samples ha’l for ordinary
kriging with MSE
for cokriging
Zhang et al., 1997 soil solute soil solute Cross validation, COK reduced the
concentrations concentrations compared MSE MSE between 30
measured at for OK with MSE and 60 %

shallower depth for COK
1.3 - 1.8 points ha'l 1.3 - 1.8 samples ha'l

Goovaerts, 1997 Soil Cd, Cu, Pb, Four combinations of VDS; correlations Correlation

(isotropic and and Co (259 Ni, Zn, Pb, and or between predicted greater and errors
anisotropic samples per field Cu. and measured and lower for COK
standardized area) ??? ME for both OK compared with
ordinary cokriging) and COK OK
Juang and Lee, Soil Cd and Zn Soil Cd and Zn at Compared OK and COK improved
1998 same or lower depth COK predictions the r2 values by 6
(scale = 2 points to 60% over
2 points ha" 5.5 points ha'I ha") with OK COK.
predictions (scale
= 5.5 points ha")

 

1' VDS = Validation data set; MSE = Mean squared error; RMSE = Root mean squared error;
ME = mean absolute error; SK = Simple Kriging; OK = Ordinary kriging; KT = Kriging with a trend
model; COK = cokriging.

42

The simplest approach has been the use of simple or multiple regression in which
a soil property is regressed on a single or multiple terrain attributes and the regression
equation used to predict the soil property at unsampled locations within the field where
terrain attributes are mapped. Success has been measured by the magnitude of the
regression coefficient of determination (R2), which ranged from 0.41 to 0.69 for the
studies reported in Table 2.1. In some cases, the accuracy of regression prediction was
evaluated using a validation data set (Bell et al., 1995), which is a more robust accuracy

measure .

The regression approach relies solely upon the relationship between the soil
property of interest and the selected terrain attributes. Geostatistical prediction
approaches utilize a statistical model of the spatial variability either using distance alone
(ordinary kriging) or in conjunction with other measured variables (e.g., co-kriging). For
ordinary kriging, a search radius is defined for each point that is to be assigned an
estimate. A mean attribute value is calculated from sample data points within this radius
and subtracted from each sample data point value. A weighted average of the residuals is
calculated. The weights are based on an empirical, statistical model of the relationship
between the separation distance and sample variance (semivariogram model). The
neighborhood mean is added to the weighted average of the residuals to calculate an
ordinary kriging point estimate. Estimation may be improved by incorporating secondary
information into kriged estimates by substituting the mean term with a smoothly
changing, rescaled variable (e.g. terrain attribute) that is linearly related to the variable
being predicted, a procedure known as kriging with an external drift (Goovaerts, 1997;

Deutch and Joumel, 1998). Multiple secondary variables can be incorporated in this

43

fashion with a procedure known as random field analysis. However, this procedure has
traditionally been used to remove spatial correlation from an analysis of variance (Stroup
et a1, 1994). Standardized ordinary co-kriging does not use the mean to incorporate
secondary information. Rather, the prediction is the sum of the weighted averages of the

primary variable and each of the secondary variables.

The performance of the various kriging approaches is mixed (Table 2.1).
Goovaerts (1997) reports that, while correlated, kriging with an external drift and
cokriging performed better than simple or ordinary kriging. In addition, cokriging
performed better than kriging with an external drifi for three of the four variables and
anisotropic cokriging performed better than isotropic cokriging. However, Gotway and
Hartford (1996) found that cokriging and kriging were respectively worse or only slightly
better than ordinary kriging. The fact that correlations between residual nitrate and yield

were not significant (p = -0.09) may explain the poor performance of these techniques.

From the studies in Table 2.1, there appears to be no consistency in which soil
properties were analyzed, which terrain attributes were selected for prediction, the
resolution of sampling of either soil properties or elevation, the scale of analysis, or the
measures of accuracy of prediction, if used at all. Furthermore, the analytical techniques
used in the various studies varied in complexity and in the effort required for analysis
(regression < kriging with external drifi << random field analysis << cokriging.
Increased complexity is only warranted if it leads to significantly improved spatial
prediction. The objective of this study was to evaluate how different analytical
techniques using terrain attributes affect the spatial prediction of soil properties and

whether relative performance of these techniques is affected by the scale of soil sampling.

44

Four analytical techniques were used to generate spatial predictions of soil carbon
obtained on 30 and 100 m regular grids using elevation, slope, and curvature as predictors

and regression, residuals, and prediction efficiency as measures of performance.

MATERIALS AND METHODS

Site description

This study was conducted in a 12.5 ha field (47° 47' 30" N, 83° 52' 30 W) located
6-km south of Durand, Michigan in the Shiawassee River watershed (Figure 2.1). The
field had been in a corn (Zea mays L.)-soybean (Glycine max L. (Merr.)) rotation for
more than lO-yr. The soil color differences in the aerial photograph (Figure 2.1a) were
related primarily to differences in soil organic matter content and drainage but did not
match well with the second soil order survey map unit boundaries (Figure 2.1b). Because
moisture conditions were not optimal in other years when USDA-AF S aerial photographs
were taken (e.g. 1979, 1983, and 1992), the striking visual differences were not captured
as they were in this photo (Figure 2.1a). The soil scientists who created this survey relied
on aerial photography taken prior to 195 8, which may have been of a lower quality or,
taken with a full crop canopy. With better aerial photos, there may have been a better
match between the color differences in Figure 2.1a and the soil survey map unit

boundaries.

45

 

Northing (m)

 

200 300 400 500 600
Easting (m)

Figure 2.1. Study location: (a) areal photograph and (b) soil and hepatic GPS
measurement locations overlain by soil type boundaries.

1‘ The scanned and enlarged aerial photograph (original scale = 1:7,920;
not georectified) taken 6/23/88 was purchased from the USDA-AF S
Aerial Photography Field Office in Salt Lake City, UT. The locations
of sample points for the three soil sampling strategies (G30 = V; Gioo =
E]; VAL = +) and the kinematic survey (0) are overlain by NRCS soil
map unit boundaries (Bt = Breckenridge sandy loam; CtA = Conover
loam with 0 to 2 % slope; MaA = Macomb loam with 0 to 2 % slope;
MsA = Metamora sandy loam with 0 to 2 % slopes; MsB = Metarnora
sandy loam with 2 to 6 % slopes; WeA = Wasepi sandy loam with 0 to
2 % slopes; Threlkeld and Feenstra, 1974).

46

 

Threlkeld and Feenstra (1974) classified and described the soils in Shiawassee
County, MI at a scale of 1:20,000. The soils were mapped (Figure 2. lb) as somewhat
poorly drained Alfisols with the exception of the Breckenridge (Bt; Coarse-loamy,
mixed, nonacid, frigid Mollie Haplaquepts), a poorly drained Inceptisol. The Metamora
sandy loam and Macomb loam (F inc-loamy, mixed, mesic Udollic Ochraqualfs) were
very similar but the Metamora was coarser in texture which means that its surface drained
somewhat faster but they both have slow subsurface drainage. While permeability is
moderately rapid for the Wasepi (Coarse-loamy, mixed, mesic Aquollic Hapludalfs)
series, it has low available water holding capacity. Most of the field had slopes of less

than 2% except for the Metamora map unit with a B slope (MsB), ranging from 2 to 6%.

Soil Sampling Design and Laboratory Analysis

Soil samples were obtained from the field (Figure 2.1b) in May of 1997 using a
30.5-m (G30; 11 = 134; 10.7 samples ha‘l) regular grid, a 100-m(Gloo; n = 12; 1 sample ha'
1) regular grid, and a set of validation points (VAL; n=26; 2.1 samples ha"). The VAL
points were collected using a 200 m unaligned grid with additional random points. The
sample point locations were flagged using a DGPS system with a base station for on-the-
go differential correction. At the each sample locations, 5 sub-samples (l at the grid
point and 4 within a 1.5-m radius) were obtained to a depth of 20-cm using a 2.5-cm
diameter soil core and composited. Soils were dried under forced air at 35° C for 3 days
and pulverized to pass a 2-mm sieve. Sieved soil was finely ground with a roller mill and
then analyzed for total carbon using a Carlo Erba NA 1500 Series 2 N/C/S analyzer (CE
Instruments Milan, Italy). A 61-m grid (G51; 11 = 38; 2.7 samples ha") was extracted

from the G30 data set to be used as a third prediction data set.

47

DEM Creation and Terrain Analysis

A kinematic GPS survey was conducted in January of 1996 using two Z-12
Ashtech GPS sensors. The mobile GPS unit was mounted on an all terrain vehicle (ATV)
traveling at about 17 km hr'l logging GPS location and elevation. Every second, a data
point was logged so that the approximate distance between measurements was 4.7 meters.
The ATV traversed the field in the east-west direction making swaths every 4.6 meters so
the field was sampled at an approximate scale of 463 samples ha". Data that had high
position dilution of precision values (PDOP) and large vertical jumps between
sequentially logged data points were removed. Several swaths were removed from the
northern region of the field because of a systematic error in the GPS data (Figure 2.1b).
Topogrid (Arclnfo ver. 7.1.1, ESRI 1997) was used to create a l x l-m grid without
drainage enforcement. Slope, aspect, and curvature (plan, profile, and tangential) were

calculated with Surfer® (Golden Software, Golden, CO).

Data Analysis

For the G30 data set, normal probability (Q-Q) plots with 95 % confidence
intervals (Friendly, 1991) were used to assess normality. Contour maps of semivariogram
surfaces were created for total carbon and the terrain attributes to determine the direction
of the anisotropic axes if anisotropy existed (Goovaerts, 1997; Isaaks and Srivastava,
1989). Directional (anisotropic) semivariograms were calculated for soil properties and
terrain attributes using angular tolerances of i 22.5° for the soil variables (Goovaerts,
1997) and i 15° for terrain attributes (a smaller angle was used because terrain model

was more densely sampled). For total carbon, nested semivariograms models

48

(combinations of spherical, exponential, and/or Gaussian models) if warranted, were

chosen based on their fit to the empirical variograms. All variogram modeling was
performed with Variowin (Pannatier, 1996). Correlations (or = 0.15) and multiple

regression (or = 0.15) were calculated using SAS (SAS, 1990) for each grid sampling

interval.

The G30 and G61 data sets were interpolated with regression, ordinary kriging,
kriging with a trend model, kriging with an external drift, random field analysis, and
standardized ordinary cokriging. Because there were so few points in the G100 data set,
only regression analysis was performed. All geostatistical interpolation methods were
conducted with GSLib (Deutsch and Joumel, 1998) except for random field analysis

which was performed using SAS (SAS, 1990).

Some theoretical understanding is required to fully appreciate these procedures.
A geostatistical prediction at an unobserved location 3 of an attribute Z is the weighted
average of the observed values of the attribute at spatial locations 5;. The weights wSm
may be restricted to non-zero values in some neighborhood N“) of the target location 3.
If n(s) is the number of observed values in the neighborhood, the predicted value is

calculated as

n(s

ZKoo—m(u) = warms-mm] Eqn. 1

i=1
Goovaerts (1997, 1999) distinguishes between kriging interpolators by the
treatment of the mean term m(u) or m(ui). The mean is constant throughout the study

area for simple kriging (SK), and within each search neighborhood, but varies through

the study area for OK, and varies gradually within each neighborhood N“) for KT, KED,

49

and RFA. The mean component is modeled as a linear combination of the coordinates
for KT, and as a linear function of a smoothly varying secondary variable (e. g. terrain
attribute) for KED, and a linear or nonlinear combination of secondary variables for
RFA. The kriging weights are calculated by solving the kriging system presented in

Chapter 1 of this dissertation.

Measures of map error included MSE, RMSE (root mean squared error), and bias.
Let V, denote the differences between predicted value and observed value at location 5;,
i=1, ..., nv = 62 of the validation data set. Map errors and prediction efficiencies were

calculated as described in the Introduction to the dissertation.

RESULTS AND DISCUSSION

The discussion here focuses on two issues. The first is whether variability in soil
and terrain attributes within the field have the magnitude and structure needed for spatial
prediction. Then interpolation methods that utilize terrain attributes will be evaluate

using measures of prediction error and efficiency.

Nature of the Data

Due to its glacial origin, the elevation and derived terrain attributes within this
field varied considerably (Figure 2.2). While total relief in the field is only 4 m, there is
considerable micro-variability within the field as evidenced by rapid changes in slope and
aspect over short distances. Therefore, micro-variability in terrain attributes may exert
significant influence in the soil and hydrologic properties of this field. Elevation, slope
and aspect (Figure 2.3) were normally distributed while plan curvature, profile curvature,

and tangential curvature were not. Kriging or cokriging with non-Gaussian data is

50

  

20 x vertical
exaggeration

  
 
     

Benefits“ \m\
N N

S\Q§% \Sio\

     
      
  
 
   

   
  

  
      

 
 
 
  

 
 

__-; to izrvrtitrrvtssxxxzfi" <9 i 261
g i Wfi‘ifiifiji‘) «1" ,1"?! LI.
a °°°i "tr-«J» W
1E ‘ 2 iii-3'3“, 3;?
o . 4d:
2 i 3
7007

 

Easting (m)

Figure 2.2. Surface maps (a) of elevation and (b) slope and a contour map (c) of
elevation overlain with arrows indicating aspect.

51

TC
(mg kg")

Elevation

Aspect

Profile Plan Slope
Curvature (angle) (degrees)

Curvature

Tangential

(m)

Curvature

 

404
20-
04

 

268 4
264 a

260 _‘ . o...

 

400 l
200 -

0-...

 

43
2-4

0+0.

 

 

0.3 ~
0.0 .4
-0.3 -

 

0.003 .
0.000 -
-0.003 -

i" -.

 

0.005 -
0.000 -

-0.005 4 '

 

 

 

-3

-2 2 3
Normal Quantiles

-1

Figure 2.3. Normal probability (Q-Q) plots for variables and log transformation with

95% confidence intervals.

52

 

 

Separation distance (meters North

Separation distance (meters East)

Figure 2.4. Semivariogram surfaces for total carbon (G30) and elevation (n = 1000).

permissible but the predictions are not guaranteed to be best linear unbiased estimates

(Cressie, 1993).

Elevation was severely anisotropic (Figure 2.4) with the axis of maximum spatial
continuity 62° East from due North. Intrinsic stationarity was assumed in the orthogonal
direction because the semivariogram did not reach a plateau. The anisotropic axes for
slope were similar to the anisotropic axes for elevation except they were rotated 90° and
were less severe (not shown). Plan curvature, profile curvature, and tangential curvature

were mildly anisotropic (not shown).

The terrain attributes had large RSV values. Elevation had a range of spatial
correlation of 275 m. Slope and aspect had ranges of about 70 m and the curvature
parameters had ranges between 10 and 20 m. Elevation, slope and aspect were suitable to
be used in a geostatistical analysis. The use of curvature in the geostatistical study is
questionable because they were spatially correlated over such a short range and because

they were not normally distributed.

53

Total carbon was normally distributed (Figure 1.3). The average value for the
field was 13 g kg'l, typical for a Michigan landscape. Despite just moderate changes in
relief and slope (Figure 2.2), total carbon content ranged substantially (2 to 29 g kg'l).
The anisotropic axes occurred in the North-South and East—West directions but at
distances of 200-m and greater the axes appear to shift to the same anisotropic axes
system as for elevation (Figure 2.4). The RSV values were not as high as might be
expected for total carbon which tends to be well structured (Table 2.2) indicating a large
nugget effect. The nugget can not be accurately estimated when distances between
sample points are great (e.g. 30.5 and 61-m). Therefore, the RSV could not be
interpreted as a measure of spatial dependence for total carbon. The range parameters for
the anisotropic model also was not interpreted because a technique had been employed to
account for zonal anisotropic accomplished by manipulating the range parameters. The
range of spatial correlation was quite large for the two isotropic models (244 and 249 m)
indicating that the data were spatially well structured. Based on the large range of spatial

correlation, total carbon is suitable for geostatistical analyses.

Ordinary kriging using the geostatistical parameters in Table 2.2 reveals that
carbon values were generally lower in the on hill tops and in the northern region of the
field and greater in the depressions in the southern and southeastern regions of the field
(Figure 2.5). Unfortunately, however, much of the detail apparent in the aerial photo
(Figure 2. l a) was not represented in this interpolation. In short, ordinary kriging did an
adequate job of predicting soil carbon across the landscape but there is room for
improvement. To apply other interpolation techniques, however, additional data

requirements must be satisfied.

54

Table 2.2. Directional (G30) and omnidirectional (G30 and G51) SCmivariogram model
parameters for carbon to be used for ordinary kriging

 

 

 

 

—————-Structure 1—— ———Strueture 2
nugget model sill direction range RSV model sill direction range
(m) (%) EL
Isotropic
(Gm) 6.2 S 31 O 249 75
- - N90° E 146 N80° E 20000
An'sg‘mp'c 10.2 G 15 60 G
(30) N 0°w 118 N10°W 950
Isotropic
(G61) 2.6 S 50 O 244 95

 

i” S = spherical; G = Gaussian; O = omnidirectional

 

600
Easting (m)

15 1 7 20 25 30

Figure 2.5. Total carbon (g kg") contour map created with ordinary kriging using the G30
data set overlain with elevation contours (lines).

 

For kriging with an external drift, the relationship between primary and secondary
variables must be linear (Goovaerts, 1997). For cokriging variables must be both
correlated and have structured cross semivariograms. Ahmed and De Marsily (1987)
state that for cokriging to be of greater predictive value than ordinary kriging, the

absolute value of the correlation coefficients between predicted and measured must

55

exceed 0.70 (Table 2.3). While most of the variables were significantly correlated with
total carbon, only elevation could explain a substantial portion of its variability. In this

study, elevation was the only variable suitable for cokriging and kriging with an external

drift.

Table 2.3. Correlations between total carbon and terrain attributes and explained

 

 

variability (G30 data set) I.
Correlations Variability in total
with carbon explained
Total by terrain
Carbon attributes (%)
Elevation -0.72 * 51
Aspect 0.17 3
Slope -0.40 * 16
Plan curvature 0.19 * 4
Profile curvature 0.20 * 4
Tangential curvature 0.21 * 4

 

l‘n = 134; * indicates significance at or = 0.05

An important requirement for cokriging is that a linear model of coregionalization
be developed that has covariance matrices that are positive semi-definite (Goovaerts,
1997). The parameters for the model at the G30 scale are listed in Table 2.4. Another
important requirement for kriging with a trend and kriging with an external drift is that a

stationary trend exists. The models for the directional trend are presented in Table 2.5.

In summary, total carbon and elevation had sufficiently large correlations and
structured semivariograms and cross semivariograms so ordinary kriging, and cokriging

were appropriate methods. Because total carbon also had a directional trend, kriging with

56

Table 2.4. Linear model of coregionalization (G30) used for cokriging. .

 

 

Structure 1 Structure 2 Structure 3
(Gaussian) (Gaussian) (Gaussian)
svl or NugT Sill Direction Range Sill Direction Range Sill Direction range
Cross-8V (m) (m) (m)
f N 90° E 145 N 90° E 2000 N 62 E 6000
TC 9 15 50 6.3
N 0°W 117 N 0°W 440 N28W 420
TC x N 90° E 145 N 90° E 2000 N 62 E 6000
. -0.23 -1.71 -1.56 -3.61
E'evatlon N 0° w 117 N 0° w 440 N 28 w 420
N 90° E 145 0 09 N 90° E 2000 N 62 E 6000
Elevation 0.02 0.225 ' 2.1
N 0°w 117 3 N 0°w 440 N28W 420

 

1 SV = semivariogram; Nug = nugget; TC = total carbon

Table 2.5. Semivariogram model parameters fo the direction of maximum spatial

continuity (N 90 E) for kriging with a trend model and kriging with an
external drift.

 

scale nugget modelI sill direction range
(m)
G30 9.4 G 16 N 90° E 190

G61 4.4 S 50 N 90° E 213
1' S = spherical; G = Gaussian

a (quadratic) trend and kriging with an external drifi (elevation) were also appropriate
methods. As will be presented in the next section, there were significant regression
relationships between total carbon and the terrain attributes. Because of this, multiple

regression and random field analysis are appropriate prediction methods.

57

Evaluation of Interpolation Methods

Stepwise regression (or = 0.15) was used to predict total carbon each scale. The
regressor variables were Easting (m), Northing (m), elevation (m), slope (%), plan
curvature (m'l), profile curvature (ml), and tangential curvature (m") as independent
variables; however, at each scale only various combinations of these regressors were

selected to be in the model by the stepwise procedure.

Total Carbon (G30) = 1451 - 5.45 x elevation - 0.0145 x Northing - 2.95 x slope +
3.24 x plan curvature + 619 x profile curvature

Total carbon (G61) = 1567 - 5.88 x elevation - 0.0203 x Northing - 4.66 x slope +
18.6 x plan curvature - 993 x profile curvature

Total Carbon (G100) = 59.48 - 0.0589 x Northing

More than half of the variability in total carbon was predicted at the G30 (R2 =
0.66), G61 (R2 = 0.77), and G100 (R2 = 0.74) scales. At the G100 scale, only Northing was
retained in the model and the R2 for this relation was large. Visually, the gradient in TC
is from N to S (Figure 2.1) and the few data points in the G100 grid could only identify
this major trend. Therefore the regression from the G100 data set (n=12) represents
spurious results because the measure of goodness were low. This is evident by a large
RMSE, low prediction efficiency (Figure 2.6), and low 1'2 between predicted and

measured (Figure 2.7) for regression procedure at the 0100 scale.

58

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

c5 60-+'fi'—'—__fl flm
9: ’7
13.9.23
0: -
o ,. , ,
4- —
tub _—
(D'c: j
2.:
“EZWIV—ll—‘Iflflfl fl HHl—lfl
o ,
X D< (9 x 0 <0 (9
I)! x
00 58 §& E 05 S'i‘fific‘ 35

Figure 2.6. Prediction efficiencies and RMSE's.

'l' FA30 = mean value of the G30 data set; OK; = isotropic ordinary kriging; OKa =
anisotropic ordinary kriging; KT = kriging with a trend; COK = cokriging; KED
= kriging with an external drift; RFA = random field analysis.

59

 

 

 

 

 

Total carbon predicted (g kg")

 

 

 

 

 

 

 

 

Total carbon measured (9 kg")

Figure 2.7. Predicted versus measured values for several prediction methods at two
scales of measurement. 7

‘l' OK= isotropic ordinary kriging; KT = kriging with a trend model; COK =
cokriging; KED = kriging with an external drift; RF A = random field analysis;

60

At the G30 intensity, there was little difference between any of the prediction
methods as assessed with measures of prediction efficiency and RMSE (Figure 2.6) or by
deviations from a 1:1 line of predicted versus measured (Figure 2.7). At this scale,
ordinary kriging and kriging with a trend model, methods which rely solely on a
statistical model of the spatial variability, were of similar predictive value as methods that
relied only on the relationship between total carbon and terrain attributes (multiple
regression). Methods that utilized both the spatial variability of carbon and its
relationship with other variables (e. g. kriging with an external drift) only performed

slightly better than these methods.

At the G61 grid intensity, the regression approach had substantially lower RMSEs
and higher prediction efficiencies than any other method at this scale. In fact, the RMSE
for multiple regression was nearly in the same range as the RMSEs for the G30
interpolations (Figure 2.6). The plots of predicted versus measured show that at the 61-m
grid scale, all of the methods that incorporated terrain information out performed those
that used only geostatistical information, despite the large range of spatial correlation. A
great deal of information about total carbon exists in the terrain attributes. The
implication is that by using this information, it is possible to reduce the number of

samples needed to achieve a certain level of accuracy.

The cost of soil sampling is inversely proportional to the square of the grid
increment. For regression, prediction efficiency and grid sampling interval were
inversely related (Figure 2.8) which allow the relationship between prediction efficiency
and cost to be modeled. It should be noted that by considering the prediction efficiency,

a squared loss function is assumed (since prediction efficiency is based on the MSE).

61

This is likely an incorrect assumption but unfortunately is the basis for most statistical
measures of map error. In other words, an accurate understanding of the cost benefit
relationship for grid sampling will require a better understanding of the relationship

between map accuracy and the relative benefit to the farmer.

Recall that in Table 2.1, researchers who used the regression approach used the
coefficients of determination (R2) as an indicator of goodness. Figure 2.8 shows that the
prediction efficiency decreases going from the G30 to G61 scales however the stepwise
regression R2 values increased from 0.66 and 0.77 for the G30 and G61 grids. This

suggests that the R2 may not be a very robust indicator of goodness.

70
65-
60- r2=1.00
55-
502
45..
40-
35 1 —T 1 1

20 40 60 80 100 120

Grid sampling interval (meters)

 

Prediction efficiency (%)

 

 

 

Figure 2.8. The relationship between prediction efficiency and sample point density for
multiple regression interpolation.

Standard measures of prediction accuracy and efficiency are important. Some
researchers use the correlations between predicted and measured while others compare
MSE values (Table 2.1). Unforttmately, these measures are not the same and can give

conflicting results. In Figure 2.7, the r2 for ordinary kriging at G61 is 0.47 and the r2 for

62

regression interpolation at G100 is 0.41 but the RMSE values are lower for ordinary
kriging than for the regression interpolation (Figure 2.6). Clearly some standards are

needed because these approaches can give conflicting results.

Although the predictions at the G30 intensity were not substantially different,
kriging with an external drift performed the best based on higher correlation coefficients
(Figure 2.7), greater prediction efficiencies, and lower RMSEs (Figure 2.6). However,
the best interpolator at the G61 scale was multiple regression. The kriging with an
external drift and multiple regression (G30 and G61) interpolations have been overlain
with elevation contours in Figure 2.9. Of these interpolations, by visual inspection, there
is more correspondence between the kriging with an external drift interpolation and the
aerial photograph (Figure 2.1a). In the southern and northeastern regions of the field, the
kriging with an external drift interpolation matches the darker tones in the aerial photo
better than the two multiple regression interpolations. But regression analysis does a
better job of assigning high carbon values to low areas like the veiny feature in the west
half of the field and low carbon values to ridges the bright red area in the north central
region of the field. The regression approach is successful for these areas because the

ridges and valleys have ether extremely positive or negative curvature values.

63

Northing (m)

 

600
Easting (m)

o 54’10121517202530

 

Figure 2.9. Predicted total carbon (g kg") using kriging with an external drift (G30) and
multiple regression (G30 and G61).

1' KED = Kriging with an external drift; REG = Multiple regression; For display
purposes regressed interpolations were matrix smoothed (8 surrounding cells, central
cell weight was 2.5; other 8 cells weights were I)

CONCLUSIONS

Terrain attributes improved spatial predictions of total carbon particularly with the
coarser sampling intensities. The comparable performance of multiple regression
interpolation procedure suggests that it may not be necessary to use the more
sophisticated geostatistical techniques. At high sampling densities (G30), interpolation
method had little impact on overall map accuracy. There was an interaction between the
scale of measurement and the most appropriate interpolation procedure. At lower
sampling intensity (G61), methods that utilized terrain attributes were more precise than
methods that did not. At this scale, multiple regression analysis yielded the best

predictions, which were nearly as accurate as the methods sampled at the G30 scale.

By using techniques that incorporate terrain attributes, sampling intensity can be
substantially reduced, while maintaining high levels of prediction accuracy and precision.
Since the cost of grid sampling is inversely related to the square of the grid sampling

increment, enhancing spatial estimates with terrain attributes is economically appealing.

SUMMARY AND RECOMENDATIONS

In the first chapter of this dissertation, grid sampling was found to be inadequate
for accurate spatial predictions of soil chemical properties for a field in Clinton County,
MI. The second chapter provides an example of how auxiliary terrain information
(elevation, slope, aspect, and curvature) can successfully be used to enhance spatial
estimates of total carbon. More work is needed to determine if terrain attributes can be
used to enhance the predictions of soil fertility variables (e.g. pH, P, K). Additionally,

other terrain attributes (e. g. wetness index, specific catchment area), high-resolution

65

multi-spectral images, electromagnetic conductivity, ground-penetrating radar could be

used as secondary information.

Soil property sensors are greatly needed for precision agriculture, unfortunately
there are few commercially available sensors that directly measure soil properties of
agronomic importance. Development of on-the-go sensors for soil nutrients has been
limited primarily to nitrogen and has not been very successful. Sensors for soil organic
matter, water content, structure, compaction are still in development but have had greater

SUCCESS.

Until better sensors are developed, a two step approach is recommended. The
first step is directed sampling based on field history, landscape features, remote sensed
imagery, and yield map variability. Directed sampling may not provide accurate soil
property maps. This was the case in the first chapter of this dissertation; however,
management zones were based on NRCS soil types from a 1215,840 survey. However, if
a producer has records of within field manure applications or records of the locations of
old field boundaries, this approach may prove to be quite accurate. The second step is the
composite sampling of areas that have depressed yields or stressed plants as indicated in
yield maps and remote sensed imagery. Plant nitrogen deficiencies can cause chlorosis,
which changes the reflective properties of the plants. Chlorosis can be identified by

aerial images.

66

BIBLIOGRAPHY

BIBLIOGRAPHY

Ahmed, S. and De Marsily, G. 1987. Comparison of geostatistical methods for estimating
transmissivity using data on transmissivity and specific capacity. Water Resources
Research, 23: 1717-1737.

Bell, J .C ., C.A. Butler, and J .A. Thomson. 1995. Soil-terrain modeling for site-specific
agriculture management. p 208-227. In P.C. Robert et al. (ed.) Site —specific
management for agriculture systems. ASA Misc. Publ. ASA, CSSA, and SSSA,
Madison, WI.

Bourennane, H. D. King, P. Chery, and A. Bruand. 1996. Improving the kriging of soil

variables using slope gradient as external drift. European journal of soil science.
47: 473-483.

Brown, J .R. 1998.(ed.) Recommended Chemical Soil Test Procedures for the North
Central Region. Missouri Ag Exp. Sta. Bul 1001.

Cressie, N. 1993. Statistics for spatial data, revised edition. John Wiley & Sons, New
York, 1991.

Deutsch, CV. and AG. JourneI. 1998. GSLIB: geostatistical software library and users
guide, second edition. Oxford University Press, New York.

Franzen, D.W. and TR. Peck. 1995. Field soil sampling density for variable rate
fertilization. J. Prod. Agric. 8(4): 568-574.

Friendly, M. 1991, SAS system for statistical graphics, lst edition, SAS Institute, Inc.',
Cary NC.

Gessler, P.E., I.D. Moore, N.J. McKenzie, and P.J. Ryan. 1995. Soil-landscape modelling
and spatial prediction of soil attributes. Int. J. Geographical information systems.
9(4):421-432.

Goovaerts. P. 1997. Geostatistics for natural resource evaluation. Oxford University
Press, New-York. 483 pp.

68

Goovaerts. P. 1999. Geostatistics in soil science: state—of-the-art prospective. In press.
Geoderrna.

Gotway, C.A. R.B. Ferguson, G.W. Hergert, and TA. Peterson. 1996. Comparison of

kriging and inverse-distance methods for mapping soil parameters. Soil Sci. Soc.
Am. J. 60:1237-1247.

Gotway, CA. and AH. Hartford. 1996. Geostatistical methods for incorporating
auxiliary information in the prediction of spatial variables. Journal of
Agricultural, biological, and Environmental Statistics. 1(1): 17-39.

Isaaks, E. and R. Srivastava. 1989. An introduction to applied geostatistics. Oxford
University Press, New York, NY.

Jaynes, DB. 1996 Improved soil mapping using electromagnetic induction surveys. p.
169-179. In. P.C. Robert et al. (ed.) Proc. 3rd international conference on precision
agriculture. ASA Misc. Publ., ASA, CSSA, and SSSA, Madison, WI.

Juang, K.W. and D.Y. Lee. 1998. A comparison of three kriging methods using auxiliary
variables in heavy-metal contaminated soils. J. Environ. Qual. 27:355-363.

McCann, B.L., D.J. Pennock., C. van Kessel, F.L. Walley. 1996. The development of
management units for site-specific farming. In. P.C. Robert et al. (ed.) Proc. 3rd

international conference on precision agriculture. ASA Misc. Publ., ASA, CSSA,
and SSSA, Madison, WI.

Mohamed, S.B., E.J. Evans, R.S. Shiel. 1996. Mapping techniques and intensity of soil
sampling for precision farming. p. 217 - 226. In. P.C. Robert et al. (ed.) Proc. 3rd
international conference on precision agriculture. ASA Misc. Publ., ASA, CSSA,
and SSSA, Madison, WI.

Moore, I.D., P.E. Gessler, G.A. Nielsen, and GA. Peterson. 1993. Soil attribute
prediction using terrain analysis. Soil Sci.Soc. Am. J. 57:443-452.

Pannatier, Y. 1996. Variowin: Software for spatial data analysis in 2D. Springer, New
York, NY.

Pierce, P.J. and Nowak. 1999. Aspects of precision farming. In. Advances in Agronomy.
In Press.

69

Pocknee, S. B.C. Boydell, H.M. Green, D]. Waters, C.K. Kvien. 1996. Directed Soil
Sampling. In. P.C. Robert et al. (ed.) Proc. 3rd international conference on
precision agriculture. ASA Misc. Publ., ASA, CSSA, and SSSA, Madison, WI.

Pregitzer, K.E. 1978. Soil survey of Clinton county, Michigan. Soil Conservation
Service.

Robertson, G.P., J .R. Crum, B.G. Ellis. 1993. The spatial variability of soil resources
following long-terrn disturbances. Oecologia. 96:451-456.

Robertson, G.P., K.M. Klingensmith, M.J. Klug, E.A. Paul, J.R. Crum, and B.G. Ellis.
1997. Soil Resources, microbial activity, and primary production across an
agricultural ecosystems. Ecological Applications, 7(1): 158-170.

Rosenbaum, MS. and M. Stiderstrtim. 1996. Cokriging of heavy metals as an aid to
biogeochemical mapping. Acta Agric. Scand. Sect. B, Soil and plant Sci. 46: 1-8.

Sadler, E.J., W.J. Busscher, P.J. Bauer, and D.L. Karlen. 1998. Spatial scale requirements
for precision farming: a case study in the southeastern USA. Agron. J. 90: 191-
197.

SAS Institute 1990. SAS/ STAT user's guide. Version 6. SAS Inst., Cary, NC.

Sayer, J.E. 1994. Concepts of variable rate technology with considerations for fertilizer
application. J. Prod. Agric., 7(2): 195-201.

Snyder, C., T. Schroeder, J. Havlin, and G. Kluitenberg. 1996. An economic anaysis of
variable rate nitrogen management. p. 989-998. In P.C. Robert et al. (ed.) Proc. 3rd
international conference on precision agriculture. ASA Misc. Publ., ASA, CSSA,
and SSSA, Madison, WI.

Stroup, W.W., P.S. Baenziger, and D.K. Multize. 1994. Removing spatial variation from
wheat yield trials: a comparison of methods. Crop Science. 34: p. 62-66.

Sudduth, K.A., J .W. Hummel, S.J. Birrell. 1997. Sensors for site-specific management. p.
183-210. In F.J. Pierce and E.J. Sadler (ed.) The state of site-specific management
for agriculture. ASA Misc. Publ., ASA, CSSA, and SSSA, Madison, WI.

70

Thompson, W.H. and RC. Robert. 1995. Evaluation of mapping strategies for variable
rate applications p303-323. In P.C. Robert et al. ed.) Site —specific management
for agriculture systems. ASA Misc. Publ. ASA, CSSA, and SSSA, Madison, WI.

Threlkeld, G.W. and J .E. F eenstra. 1974. Soil survey of Shiawassee county, Michigan.
Soil Conservation service.

Vaughan, P.J., S.M. Lesch, D.L. Corwin, and D.G. Cone. 1995. Water content effect on
soil salinity prediction: a geostatistical study using cokriging. Soil Sci. Soc. Am.
J. 59:1146-1156.

Vitosh, M.L, J .W. Johnson, and DB. Mengel: 1995. Tri-state fertilizer recommendations
for corn, soybeans, wheat and alfalfa. Michigan State University Ext. Bull. E-
2567.

Wibawa, W.D., D.L. Dludlu, L.J. Swenson, D.G. Hopkins, and WC. Dahnke. 1993.
Variable fertilizer application based on yield goal, soil fertility, and soil map unit.
J. Prod. Agric., 6(2): 255-262.

Wollenhaupt, NC, and DD. Buchholz. 1993. Profitability of farming by soils. p. 199-
211. In P.C. Roberts et al. (ed). Soil specific crop management. ASA Misc. Publ.
ASA, CSSA, and SSSA, Madison, WI.

Wollenhaupt, N.C., R.P. Wolkowski, and MK. Clayton. 1994. Mapping Soil Test

phosphorus and potassium for variable-rate fertilizer application. J. Prod. Agric.
7(4): 441 -448.

Zhang, R., D.E. Myers, and A.W. Warrick. 1992. Estimation of the spatial distribution of

soil chemicals using pseudo-cross-variograrns. Soil Sci. Soc. Am. J. 5621444-
1452.

Zhang, R., P. Shouse, and S. Yates. 1997. Use of pseudo-crossvariograms and cokriging
to improve estimates of soil solute concentrations. Soil Sci Soc. Am. J. 61 :1342-
1347.

71

 

MICHIGAN srnre UNIV. LreRnRrEs
lllWilliillWillIllllillillllllllllllllllllllllll
31293017665062