Trauth M.H., MATLAB® Recipes for Earth Sciences, Third edition

Подождите немного. Документ загружается.

228 7 SPATIAL DATA

111

Elevation Mode

Flow Direction Flow Accumulation

Fig. 7.14 Schematic of calculation of  ow accumulation by the D8 method.

7.10 ANALYSIS OF DIGITAL ELEVATION MODELS (BY R. GEBBERS) 229

7 SPATIAL DATA

tation of hydrological characteristics from the DEM.

Flow accumulation, also called specific catchment area or upslope

contributing area, is de ned as the number of cells, or area, that contrib-

ute runo to a particular cell (Fig. 7.14). In contrast to the local parameters

slope and aspect,  ow accumulation can only be determined from the glob-

al neighborhood.  e principal operation is to add cell in ows from higher

neighboring cells, starting from the speci ed cell and working up to the

watersheds. Before adding together the out ows from each cell, we need

to determine the gradient of each individual cell towards each neighbor-

ing cell, indexed by

N.  e array N contains indices for the eight adjacent

cells according to the MATLAB convention as shown in Figure 7.13. We

make use of the

circshift function to access the neighboring cells. For

a two-dimensional matrix

Z, the function circshift(Z,[r c]) circu-

larly shi s the values in the matrix

Z by an amount of rows and columns

given by

r and c, respectively. For example, circshift(Z,[1 1]) will

circularly shi

Z one row down and one column to the right.  e individual

gradients are calculated by

for the eastern, southern, western, and northern neighbors (the so-called

rook’s case) and by

for the diagonal neighbors (bishop’s case). In these formulae, h is the

cell size, z

r,c

is the elevation of the central cell and z

r+y,c+x

the eleva-

tion of a neighbor.  e cell indices x and y are obtained from the matrix

N.  e gradients are stored in a three-dimensional matrix grads, where

grads(:,:,1) contains the gradients towards the neighboring cells to

the east,

grads(:,:,2) contains the gradients towards the neighboring

cells to the south-east, and so on. Negative gradients indicate out ow from

the central cell towards the relevant neighboring cell. To obtain the sur-

face  ow between cells, gradients are transformed by the inverse tangent of

grads divided by 0.5π.

N = [0 -1;-1 -1;-1 0;+1 -1;0 +1;+1 +1;+1 0;-1 +1];

[a b] = size(SRTM);

grads = zeros(a,b,8);

for c = 2 : 2 : 8

grads(:,:,c) = (circshift(SRTM,[N(c,1) N(c,2)]) ...

230 7 SPATIAL DATA

-SRTM)/sqrt(2*90);

end

for c = 1 : 2 : 7

grads(:,:,c) = (circshift(SRTM,[N(c,1) N(c,2)]) ...

-SRTM)/90;

end

grads = atan(grads)/pi*2;

Since a central cell can have several downslope neighbors, water can  ow

in several directions.  is phenomenon is called divergent flow. Early  ow

accumulation algorithms were based on the single- ow-direction method

(know as the D8 method, Fig. 7.14), which allows  ow to only one of the

cell's eight neighboring cells.  is method cannot model divergence in

ridge areas and tends to produce parallel  ow lines in some situations. Here,

we illustrate the use of a multiple- ow-direction method, which allows  ow

from a cell to multiple neighboring cells.  e proportion of the total  ow

that is assigned to a neighboring cell is dependent on the gradient between

the central cell and that particular neighboring cell, and is therefore a frac-

tion of the total out ow. Even though multiple- ow methods produce more

realistic results in most situations, they tend to result in dispersion in val-

leys, where the  ow should be more concentrated. A weighting factor w is

therefore introduced, which controls the relationship between the out ows.

A recommended value for w is 1.1; higher values would concentrate the  ow

in the direction of the steepest slope, while w=0 would result in an extreme

dispersion. In the following sequence of commands, we  rst select those

gradients that are less than zero and then multiply the gradients by the

weighting.

w = 1.1;

flow = (grads.*(-1*grads<0)).^w;

We then sum up the upslope gradients along the third dimension of the

flow matrix. Replacing all upslope gradient values of 0 by a value of 1

avoids the problems created by division by zero.

upssum = sum(flow,3);

upssum(upssum==0) = 1;

We divide the  ows by upssum to obtain fractional weights that add up to a

7.10 ANALYSIS OF DIGITAL ELEVATION MODELS (BY R. GEBBERS) 231

7 SPATIAL DATA

3200

3400

3600

3800

100

150

200

250

300

350

4000

Slope (°)Elevation (m)

Watershed

Flow Accumulation Wetness Index

Aspect (°)

Fig. 7.15 Display of a subset of the SRTM data set used in Section 7.5 and primary and

secondary attributes of the digital elevation model; a elevation, b slope, c aspect, d watershed,

e  ow accumulation and f wetness index.

232 7 SPATIAL DATA

total of one.  is is achieved separately for each layer of the 3D flow matrix

by a

for loop:

for i = 1:8

flow(:,:,i) = flow(:,:,i).*(flow(:,:,i)>0)./upssum;

end

 e 2D matrix inflowsum will store the intermediate in ow totals for each

iteration.  ese intermediate totals are then summed to reach a  gure for

the total accumulated  ow

flowac at the end of each iteration.  e initial

values for

inflowsum and flowac are obtained through upssum.

inflowsum = upssum;

flowac = upssum;

Another 3D matrix inflow is now needed, in which to store the total inter-

mediate in ow from all neighbors:

inflow = grads*0;

Flow accumulation is terminated when there is no in ow, or translated

into MATLAB code, we use a conditional

while loop that terminates if

sum(inflowsum(:))==0.  e number of non-zero entries in inflow-

sum will decrease during each loop.  is is achieved by alternately updating

inflow and inflowsum. Here, inflowsum is updated with the interme-

diate

inflow of the neighboring cells weighted by flow under the condi-

tion that the neighboring cells are contributing cells, i.e., where

grads are

positive. Where not all neighboring cells are contributing cells, the inter-

mediate

inflowsum is reduced, as also is inflow.  e  ow accumula-

tion

flowac increases through consecutive summation of the intermediate

inflowsum.

while sum(inflowsum(:))>0

for i = 1:8

inflow(:,:,i) = circshift(inflowsum,[N(i,1) N(i,2)]);

end

inflowsum = sum(inflow.*flow.*grads>0,3);

flowac = flowac + inflowsum;

end

We display the result as a pseudocolor map with log-scaled values

(Fig. 7.15e).

h = pcolor(log(1+flowac));

colormap(flipud(jet)), colorbar

set(h,'LineStyle','none')

axis equal

7.10 ANALYSIS OF DIGITAL ELEVATION MODELS (BY R. GEBBERS) 233

7 SPATIAL DATA

title('Flow accumulation')

[r c] = size(flowac);

axis([1 c 1 r])

set(gca,'TickDir','out');

 e plot displays areas with high  ow accumulation in shades of blue, and

areas with low  ow accumulation, usually corresponding to ridges, in

shades of red. We used a logarithmic scale for mapping the  ow accumula-

tion in order to obtain a better representation of the results.  e simpli ed

algorithm introduced here for calculating  ow accumulation can be used

to analyze sloping terrains in DEMs. In  at terrains, where the slope ap-

proaches zero, no  ow direction can be generated by our algorithm and thus

 ow accumulation stops. Such situations require more sophisticated algo-

rithms to perform analyses on completely  at terrain.  ese more advanced

algorithms also include sink- lling routines to avoid spurious sinks that

interrupt  ow accumulation. Small depressions can be  lled by smoothing,

as we did at the beginning of this section.

 e  rst part of this section was about primary relief attributes.

Secondary attributes of a DEM are functions of two or more primary at-

tributes. Examples of secondary attributes are the wetness index and the

stream power index.  e wetness index for a cell is the log of the ratio be-

tween the area of the catchment for that particular cell and the tangent of

its slope:

 e term 1+flowac avoids the problems associated with calculating the loga-

rithm of zero when

flowac=0.  e wetness index is used to predict the soil

water content ( saturation) resulting from lateral water movement.  e po-

tential for waterlogging is usually highest in the lower parts of catchments,

where the slopes are more gentle. Flat areas with a large upslope area have

a high wetness index compared to steep areas with small catchments.  e

wetness index

weti is computed and displayed by

weti = log((1+flowac)./tand(slp));

h = pcolor(weti);

colormap(flipud(jet)), colorbar

set(h,'LineStyle','none')

axis equal

title('Wetness index')

[r c] = size(weti);

axis([1 c 1 r])

234 7 SPATIAL DATA

set(gca,'TickDir','out');

In this plot, blue colors indicate high values for the wetness index, while red

colors represent low values (Fig. 7.15f). In our example, soils in the south-

east are the most likely to have a high water content due to the runo from

the large central valley and the  atness of the terrain.

 e stream power index is another important secondary relief attribute

which is frequently used in hillslope hydrology, geomorphology, soil sci-

ence, and related disciplines. As a measure of stream power it provides an

indication of the potential for sediment transport and erosion by water. It

is de ned as the product of the area of catchment for a speci c cell and the

tangent of the slope of that cell:

 e potential for erosion is high when large quantities of water (calculated

by  ow accumulation) are fast  owing due to an extreme slope.  e follow-

ing series of commands compute and display the stream power index:

spi = flowac.*tand(slp);

h = pcolor(log(1+spi));

colormap(jet), colorbar

set(h,'LineStyle','none')

axis equal

title('Stream power index')

[r c] = size(spi);

axis([1 c 1 r])

set(gca,'TickDir','out');

 e wetness and stream power indices are particularly useful in high res-

olution terrain analysis, i.e., digital elevation models sampled at intervals

of less than 30 meters. In our terrain analysis example we have calculated

weti and spi from a medium resolution DEM, and must expect a degree

of scale dependency in these attributes.

 is section has illustrated the use of basic tools for analyzing digital

elevation models. A more detailed introduction to digital terrain modeling

is given in the book by Wilson & Galant (2002). Furthermore, the article

by Freeman (1991) provides a comprehensive summary of digital terrain

analysis including an introduction to the use of advanced algorithms for

 ow accumulation.

7.11 GEOSTATISTICS AND KRIGING (BY R. GEBBERS) 235

7 SPATIAL DATA

7.11 Geostatistics and Kriging (by R. Gebbers)

Geostatistics describes the autocorrelation of one or more variables in 1D,

2D, and 3D space, or even in 4D space-time, to make predictions for un-

observed locations, to give information about the accuracy of the predic-

tions and to reproduce spatial variability and uncertainty.  e shape, range,

and direction of the spatial autocorrelation are described by a variogram,

which is the main tool in linear geostatistics.  e origins of geostatistics can

be traced back to the early 1950s when the South African mining engineer

Daniel G. Krige  rst published an interpolation method based on the spa-

tial dependency of samples. In the 60s and 70s, the French mathematician

George Matheron developed the theory of regionalized variables, which

provides the theoretical foundations for Krige’s more practical methods.

 is theory forms the basis of several procedures for the analysis and esti-

mation of spatially dependent variables, which Matheron called geostatis-

tics. Matheron as well coined the term kriging for spatial interpolation by

geostatistical methods.

Theorical Background

A basic assumption in geostatistics is that a spatiotemporal process is com-

posed of both deterministic and stochastic components (Fig. 7.16).  e de-

terministic components can be global and local trends (sometimes called

drifts).  e stochastic component is composed of a purely random part and

an autocorrelated part.  e autocorrelated component implies that on aver-

age, closer observations are more similar to each other than more widely

separated observations.  is behavior is described by the variogram where

squared di erences between observations are plotted against their separa-

tion distances. Krige's fundamental idea was to use the variogram for inter-

polation, as means of determining the amount of in uence that neighboring

observations have when predicting values for unobserved locations. Basic

linear geostatistics includes two main procedures: variography for model-

ing the variogram, and kriging for interpolation.

Preceding Analysis

Because linear geostatistics as presented here is a parametric method, the

underlying assumptions need to be checked by a preceding analysis. As

with other parametric methods, linear geostatistics is sensitive to outli-

ers and deviations from a normal distribution. We  rst open the data  le

236 7 SPATIAL DATA

0 100 200 300 400 500 600 0 100 200 300 400 500 600

0 100 200 300 400 500 600

0 102030405060

3.0

1.5

0.0

-1.5

-3.0

3.0

1.5

0.0

-1.5

-3.0

3.0

1.5

0.0

-1.5

-3.0

3.0

1.5

0.0

-1.5

-3.0

3.0

1.5

0.0

-1.5

-3.0

0.8

0.6

0.4

0.2

1.0

Lag Distance

Spatiotemporal Process

Global Trend Component

Local Trend Component

Random Component

Autocorrelation Component

Variogram

Fig. 7.16 Components of a spatiotemporal process and the variogram.  e variogram (f)

should only be derived from the autocorrelated component.

7.11 GEOSTATISTICS AND KRIGING (BY R. GEBBERS) 237

7 SPATIAL DATA

geost_dat.mat containing xyz data triplets, and then plot the sampling lo-

cations. By doing this, we can check the point distribution and detect any

major errors in the data coordinates, x and y.

clear

load geost_dat.mat

plot(x,y,'.')

 e range of the observations z can be checked by

min(z)

ans =

3.7199

max(z)

ans =

7.8460

For linear geostatistics, the observations z should be Gaussian distributed.

 is is usually only tested by visual inspection of the histogram because sta-

tistical tests are o en too sensitive if the number of samples exceeds ca. 100.

One can also calculate the skewness and kurtosis of the data.

hist(z)

skewness(z)

ans =

0.2568

kurtosis(z)

ans =

2.5220

A  at-topped or multiple-peaked distribution suggests that there is more

than one population present in the data set. If these populations can be

related to particular areas they should be treated separately. Another rea-

son for multiple peaks can be preferential sampling of areas with high and/

or low values.  is usually happens as a result of some a priori knowledge

and is known as a cluster e ect. Dealing with a cluster e ect is described in

Deutsch and Journel (1998) and Isaaks and Srivastava (1998).

Most problems arise from positive skewness, i.e., if the distribution has

a long tail to the right. According to Webster and Oliver (2001), one should

consider root transformation if the skewness is between 0.5 and 1, and loga-