Trauth M.H., MATLAB® Recipes for Earth Sciences, Third edition
Подождите немного. Документ загружается.
248 7 SPATIAL DATA
where λ
i
are the weights that have to be estimated. e sum of the weights
should be equal to one in order to guarantee that the estimates are unbi-
ased:
e expected (average) error for the estimation must be zero. at is:
where z
x0
is the true, but unknown value. We can use the above equations to
compute algebraically the mean-squared error in terms of the variogram:
where E is the estimation or kriging variance, which must be minimized,
γ(x
i
, x
0
) is the variogram (semivariance) between the data points and the
unobserved points, γ(x
i
, x
j
) is the variogram between the data points x
i
and
x
j
, and λ
i
and λ
j
are the weights of the ith and jth data point.
For kriging we must minimize this equation (a quadratic objective func-
tion), satisfying the condition that the sum of the weights should be equal
to one (linear constraint). is optimization problem can be solved using a
Lagrange multiplier ν resulting in a linear kriging system of N+1 equations
and N+1 unknowns:
A er obtaining the weights λ
i
, the kriging variance is given by
e kriging system can be presented in a matrix notation:
7.11 GEOSTATISTICS AND KRIGING (BY R. GEBBERS) 249
7 SPATIAL DATA
where
is the matrix of the coe cients: these are the modeled variogram values for
the pairs of observations. Note that on the diagonal of the matrix, where
separation distance is zero, the value of γ disappears.
is the vector of the unknown weights and the Lagrange multiplier.
is the right-hand-side vector. To obtain the weights and the Lagrange multi-
plier the matrix G_mod is inverted:
e kriging variance is given by
250 7 SPATIAL DATA
For our calculations with MATLAB we need the matrix of coe cients de-
rived from the distance matrix D and a variogram model. D was calculated
in the variography section above and we use the exponential variogram
model with a nugget, sill and range from the previous section:
G_mod = (nugget + sill*(1 - exp(-3*D/range))).*(D>0);
We then take the number of observations and add a column and row vector
of all ones to the
G_mod matrix and a zero in the lower le corner:
n = length(x);
G_mod(:,n+1) = 1;
G_mod(n+1,:) = 1;
G_mod(n+1,n+1) = 0;
Now the G_mod matrix has to be inverted:
G_inv = inv(G_mod);
A grid with the locations of the unknown values is needed. Here we use a
grid cell size of ve within a quadratic area ranging from 0 to 200 in x and y
direction. e coordinates are created in matrix form by:
R = 0 : 5 : 200;
[Xg1,Xg2] = meshgrid(R,R);
and converted to vectors by:
Xg = reshape(Xg1,[],1);
Yg = reshape(Xg2,[],1);
We then allocate memory for the kriging estimates Zg and the kriging vari-
ance
s2_k by:
Zg = Xg * NaN;
s2_k = Xg * NaN;
We now krige the unknown values at each grid point:
for k = 1 : length(Xg)
DOR = ((x - Xg(k)).^2 + (y - Yg(k)).^2).^0.5;
G_R = (nugget + sill*(1 - exp(-3*DOR/range))).*(DOR>0);
G_R(n+1) = 1;
E = G_inv * G_R;
Zg(k) = sum(E(1:n,1).*z);
s2_k(k) = sum(E(1:n,1).*G_R(1:n,1))+E(n+1,1);
end
e rst command computes the distance between the grid points (Xg,Yg)
7.11 GEOSTATISTICS AND KRIGING (BY R. GEBBERS) 251
7 SPATIAL DATA
and the observation points (x,y). en we construct the right-hand-side
vector of the kriging system by using the variogram model
G_R and adding
one to the last row. We next obtain the matrix
E with the weights and the
Lagrange multiplier. e estimate
Zg at each point k is the weighted sum of
the observations
z. Finally, the kriging variance s2_k of the grid point is
computed and we can plot the results. First, we create a grid of the kriging
estimate and the kriging variance:
r = length(R);
Z = reshape(Zg,r,r);
SK = reshape(s2_k,r,r);
A subplot on the le presents the kriged values:
subplot(1,2,1)
h = pcolor(Xg1,Xg2,Z);
set(h,'LineStyle','none')
axis equal
ylim([0 200])
title('Kriging Estimate')
xlabel('x-Coordinates')
ylabel('y-Coordinates')
colorbar
e le hand subplot presents the kriging variance:
subplot(1,2,2)
h = pcolor(Xg1,Xg2,SK);
set(h,'LineStyle','none')
axis equal
ylim([0 200])
title('Kriging Variance')
xlabel('x-Coordinates')
ylabel('y-Coordinates')
colorbar
hold on
and we overlay overlaying the sampling positions:
plot(x,y,'ok')
hold off
e kriged values are shown in Figure 7.22a. e kriging variance depends
only on the distance from the observations and not on the observed values
(Fig. 7.22b). Kriging reproduces the population mean when observations
are beyond the range of the variogram; at the same time, the kriging vari-
ance increases (lower right corner of the maps in Figure 7.22). e kriging
variance can be used as a criterion to improve sampling design and it is
needed for back-transformation in lognormal kriging. Back-transformation
252 7 SPATIAL DATA
4.555.566.57
0 50 100 150 200
0
50
100
150
200
0 50 100 150 200
0
50
100
150
200
0.3 0.4 0.5 0.6 0.7 0.8 0.9
x−Coordinates x−Coordinates
y−Coordinates
y−Coordinates
Kriging Estimate
Kriging Variance
ab
Fig. 7.22 Interpolated values on a regular grid by ordinary point kriging using a an ex-
ponential variogram model; b kriging variance as a function of the distance from the
observations (empty circles).
for lognormal kriging is achieved by:
Discussion of Kriging
Point kriging as presented here is an exact interpolator. It reproduces ex-
actly the values at an observation point, even though a variogram with a
nugget e ect is used. Smoothing can be achieved by including the variance
of the measurement errors (see Kitanidis 1997), and by block kriging which
averages the observations within a certain neighborhood (or block). While
kriging variance depends only on the distance between the observed and
the unobserved locations it is primarily a measure of the density of infor-
mation (Wackernagel 2003). e accuracy of kriging is better evaluated
by cross-validation using a resampling method or surrogate test (Sections
4.6 and 4.7). e in uence of the neighboring observations on the estima-
tion depends on their con guration, as summarized by Webster and Oliver
RECOMMENDED READING 253
7 SPATIAL DATA
(2001): "Near points carry more weight than more distant ones; the relative
weight of a point decreases when the number of points in the neighborhood
increases; clustered points carry less weight individually than isolated ones
at the same distance; data points can be screened by ones lying between
them and the target." Sampling design for kriging is di erent from the de-
sign that might be optimal for variography. A regular grid, triangular or
quadratic, can be regarded as optimal.
e MATLAB code presented here is a straightforward implementa-
tion of the above formulae. In professional programs the number of data
points entering the G_mod matrix is restricted and the inversion of G_mod
is avoided by working with the covariances instead of the variograms
(Webster and Oliver 2001, Kitanidis 1997). For those who are interested in
programming and in a deeper understanding of algorithms, Deutsch and
Journel (1992) is essential reading. e best internet source is the homepage
for AI-GEOSTATISTICS:
http://www.ai-geostats.org
Recommended Reading
Cressie N (1993) Statistics for Spatial Data, Revised Edition. John Wiley & Sons, New York
Davis JC (2002) Statistics and Data Analysis in Geology, third edition. John Wiley & Sons,
New York
Deutsch CV, Journel AG (1998) GSLIB – Geostatistical So ware Library and User’s Guide,
Second edition. Oxford University Press, Oxford
Freeman TG (1991) Calculating Catchment Area with Divergent Flow Based on a Regular
Grid. Computers and Geosciences 17:413–422
Gringarten E, Deutsch CV (2001) Teacher’s Aide Variogram Interpretation and Modeling.
Mathematical Geology 33:507–534
Isaaks E, Srivastava M (1989) An Introduction to Applied Geostatistics. Oxford University
Press, Oxford
Gringarten E, Deutsch CV (2001) Teacher’s Aide Variogram Interpretation and Modeling.
Mathematical Geology 33:507–534
Kitanidis P (1997) Introduction to Geostatistics – Applications in Hydrogeology. Cambridge
University Press, Cambridge
Olea RA (1984) Systematic Sampling of Spatial Functions. Kansas Series on Spatial Analysis
7, Kansas Geological Survey, Lawrence, KS
Pannatier Y (1996) VarioWin – So ware for Spatial Data Analysis in 2D, Springer, Berlin
Heidelberg New York
Pardo-Igúzquiza E, Dowd PA (1997) AMLE3D: A Computer Program for the Interference
of Spatial Covariance Parameters by Approximate Maximum Likelihood Estimation.
Computers and Geosciences 23:793–805
Rendu JM, Readdy L (1982) Geology and Semivariogram – A Critical Relationship. In:
Johnson TB, Barns RJ (eds) Application of Computer & Operation Research in the
254 7 SPATIAL DATA
Mineral Industry. 17th Intern. Symp. American Institute of Mining. Metallurgical and
Petroleum Engineers, New York, pp. 771–783
Sandwell DT (1987) Biharmonic Spline Interpolation of GEOS-3 and SEASAT Altimeter
data. Geophysical Research Letters 2:139–142
Swan ARH, Sandilands M (1995) Introduction to Geological Data Analysis. Blackwell
Sciences, Oxford
e Mathworks (2010) Mapping Toolbox – User’s Guide. e MathWorks, Natick, MA
Golden So ware, Inc. (2002) Surfer 8 (Surface Mapping System). Golden, Colorado
Wackernagel H. (2003) Multivariate Geostatistics: An Introduction with Applications.
ird, completely revised edition. Springer, Berlin Heidelberg New York
Webster R, Oliver MA (2001) Geostatistics for Environmental Scientists. John Wiley &
Sons, New York
Wessel P, Bercovici D (1998) Gridding with Splines in Tension: A Green Function Approach.
Mathematical Geology 30:77–93
Wilson JP, Gallant JC (2000) Terrain Analysis, Principles and Applications. John Wiley &
Sons, New York
8 IMAGE PROCESSING
8 Image Processing
8.1 Introduction
Computer graphics are stored and processed as either vector or raster data.
Most of the data types that were encountered in the previous chapter were
vector data, i. e., points, lines and polygons. Drainage networks, the outlines
of geologic units, sampling locations and topographic contours are all ex-
amples of vector data. In Chapter 7, coastlines are stored in a vector format
while bathymetric and topographic data are saved in a raster format. Vector
and raster data are o en combined in a single data set, for instance, to dis-
play the course of a river on a satellite image. Raster data are o en converted
to vector data by digitizing points, lines or polygons. Conversely, vector data
are sometimes transformed to raster data.
Images are generally represented as raster data, i. e., as a 2D array of
color intensities. Images are everywhere in geosciences. Field geologists use
aerial photos and satellite images to identify lithologic units, tectonic struc-
tures, landslides and other features within a study area. Geomorphologists
use such images for the analysis of drainage networks, river catchments,
and vegetation or soil types. e analysis of images from thin sections, the
automated identi cation of objects, and the measurement of varve thick-
nesses all make use of a great variety of image processing methods.
is chapter is concerned with the analysis and display of image data.
Firstly, the various ways that raster data can be stored on the computer are
explained (Section 8.2). e main tools for importing, manipulating and
exporting image data are presented in Section 8.3. is knowledge is then
used to process and to georeference satellite images (Sections 8.4 and 8.5).
On-screen digitization techniques are discussed in Section 8.6. Finally, sec-
tions 8.7 to 8.9 are concerned with color-intensity transects of laminated
sediments, and with automated grain size analysis and charcoal quanti-
cation in microscope images. e Image Processing Toolbox is used for
the speci c examples throughout this chapter. While the MATLAB User’s
Guide to the Image Processing Toolbox provides an excellent general in-
M.H. Trauth, MATLAB
®
Recipes for Earth Sciences, 3rd ed.,
DOI 10.1007/978-3-642-12762-5_8, © Springer-Verlag Berlin Heidelberg 2010
256 8 IMAGE PROCESSING
troduction to the analysis of images, this chapter provides an overview of
typical applications in earth sciences.
8.2 Data Storage
Vector and raster graphics are the two fundamental methods for storing
pictures. e typical format for storing vector data has already been intro-
duced in the previous chapter. In the following example, the two columns
in the le coastline.txt represent the longitudes and latitudes of the points
of a polygon.
NaN NaN
42.892067 0.000000
42.893692 0.001760
NaN NaN
42.891052 0.001467
42.898093 0.007921
42.904546 0.013201
42.907480 0.016721
42.910414 0.020828
42.913054 0.024642
(cont'd)
e NaNs help to identify break points in the data (Section 7.2).
e raster data are stored as 2D arrays. e elements of these arrays
represent variables such as the altitude of a grid point above sea level, the
annual rainfall or, in the case of an image, the color intensity values.
174 177 180 182 182 182
165 169 170 168 168 170
171 174 173 168 167 170
184 186 183 177 174 176
191 192 190 185 181 181
189 190 190 188 186 183
Raster data can be visualized as 3D plot. e x and y are the indices of
the 2D array or any other reference frame, and z is the numerical value of
the elements of the array (see also Chapter 7). Alternatively, the numerical
values contained in the 2D array can be displayed as a pseudocolor plot,
which is a rectangular array of cells with colors determined by a colormap.
A colormap is an m-by-3 array of real numbers between 0.0 and 1.0. Each
row de nes a red, green, blue (RGB) color. An example is the above array
that could be interpreted as grayscale intensities ranging from 0 (black) to
255 (white). More complex examples include satellite images that are stored
in 3D arrays.
8.2 DATA STORAGE 257
8 IMAGE PROCESSING
As previously discussed, a computer stores data as bits that have one of
two states, one and zero (Chapter 2). If the elements of the 2D array repre-
sent the color intensity values of the pixels (short for picture elements) of an
image, 1-bit arrays contain only ones and zeros.
0 0 1 1 1 1
1 1 0 0 1 1
1 1 1 1 0 0
1 1 1 1 0 1
0 0 0 0 0 0
0 0 0 0 0 0
is 2D array of ones and zeros can be simply interpreted as a black-and-
white image, where the value of one represents white and zero corresponds
to black. Alternatively, the 1-bit array could be used to store an image con-
sisting of any two di erent colors, such as red and blue.
In order to store more complex types of data, the bits are joined together
to form larger groups, such as bytes consisting of eight bits. Since the earli-
est computers could only process eight bits at a time, early computer code
was written in sets of eight bits, which came to be called bytes. Hence, each
element of the 2D array or pixel contains a vector of eight ones or zeros.
1 0 1 0 0 0 0 1
ese 8 bits (or 1 byte) allow 2
8
=256 possible combinations of the eight ones
or zeros. erefore, 8 bits are enough to represent 256 di erent intensities
such as grayscales. e 8 bits can be read in the following way reading from
right to le : a single bit represents two numbers, two bits give four numbers,
three bits show eight numbers, and so forth up to a byte, or eight bits, which
represents 256 numbers. Each added bit doubles the count of numbers. Here
is a comparison of binary and decimal representations of the number 161:
128 64 32 16 8 4 2 1 (value of the bit)
1 0 1 0 0 0 0 1 (binary)
128 + 0 + 32 + 0 + 0 + 0 + 0 + 1 = 161 (decimal)
e end members of the binary representation of grayscales are
0 0 0 0 0 0 0 0
which is black, and
1 1 1 1 1 1 1 1
which is pure white. In contrast to the above 1-bit array, the one-byte array