MATLAB Kriging Toolbox

(version 3.0: février 1998)

Note: During conversion from a Word document to html, the figures were either lost or only partially converted. Please use the original Word document in order to see all equations and other figures.

Caroline Lafleur

INRS-Océanologie, 310 Allée des Ursulines, Rimouski, Qc., Canada, G5L 3A1

The Matlab Kringing Toolbox is free and hence no support or warranty are provided.

Specifications

The Kriging Toolbox version 3.0 is matlab 5.1 and 5.2 compatible under Windows 95. It is an upgrade of version 2.0 which has been compiled under matlab 4.2. Please note that this upgrade only uses 2-D matrices even though the new matlab version supports greater matrix dimensions.

This toolbox can be used "as it is" without other matlab toolbox except the probability functions (chi). However, it may also use the leastsq.m function from Matlab's Opimization toolbox (see fitvario.m). The mex files have been compiled with Microsoft Fortran Power station V4.0 for Matlab 5.1. As far as we have been able to test, they also work with Matlab 5.2.

The Chi Toolbox is necessary for normality study. It is provided with the Kriging Toolbox.

Description

The development of this toolbox is based on the necessity of using objective analysis of scalars in 2 or 3 dimensions in physical oceanography. This type of interpolation usually gives better results than standard interpolation methods. Furthermore, it has the non-negligible advantage of giving estimates of interpolation errors.

This toolbox is almost entirely made up of functions from the book by Deutsch and Journel (1992) and from the paper by Marcotte (1991). The variogram functions are MEX-files compiled from the former while the cokriging functions were published, in Matlab format, in Marcotte's 1991 paper. All the parameters and examples can be found, in English, in the two publications. The book by Journel and Huijbregts (1992) is the best book on semivariograms. A complete example of optimal estimation in physical oceanography can be found in the paper by Denman and Freeland (1985). As well, kridemo shows outlines of a 2-D objective analysis.

Denman, K.L. and H.J. Freeland, 1985. Correlation Scales, Objective Mapping and a Statistical Test of Geostrophy over the Continental Shelf. J. Mar. Res., 43: 517-539.

Deutsch, C. V and A. G. Journel, 1992. GSLIB: Geostatistical Software Library and User's Guide. Oxford University Press, Oxford, 340 p.

Journel, A.G. and C.J. Huijbregts, 1992. Mining Geostatistics. Academic Press, New York, 600 p.

Marcotte, D. 1991. Cokrigeage with MATLAB. Computers & Geosciences. 17(9): 1265-1280.

Comments, suggestions or questions?

Many functions are still not completely tested. Please report any bugs or problems to

Caroline Lafleur, INRS-Océanologie, Tél: (418) 724 1650 poste 1296

caroline_lafleur@uqar.uquebec.ca

Yves Gratton, INRS-Océanologie, Tél: (418) 724 1741

yves_gratton@uqar.uquebec.ca

Kriging Toolbox Contents

Variogram functions

confint Confidence intervals.

fitvario Optimization of variogr.

fitvari2 Optimization of variogr2 (without the Optimization Toolbox).

fun Called from fitvario to estimate variograms.

gam2 MEX-file called from vario2dr.

gam3 MEX-file called from vario3dr.

gamv2 MEX-file called from vario2di.

gamv2uv MEX-file called from var2diuv.

gamv3 MEX-file called from vario3di.

mrqmin Least-square fitting: Levenberg-Marquardt method.

outvario Output of variogram functions.

variogr Models of semivariogram and correlogram.

variogr2 variogr for fitting procedures not using the Optimization Toolbox.

vario2di Variogram of irregularly spaced 2-D data.

vario2dr Variogram of regularly spaced 2-D data.

vario3di Variogram of irregularly spaced 3-D data.

vario3dr Variogram of regularly spaced 3-D data.

var2diuv Variogram of irregularly spaced 2-D vectors.

Kriging functions

barnes 2-D spatial filter for kriged data.

cokri Point or block cokriging in D dimensions.

cokri2 Cokriging function called from cokri.

davis Point kriging using Davis' set of equations.

filresp Barnes' filter response in the wavelength domain.

tintore Application of Barnes' filter with the Tintoré's parameters (1991).

Related functions

covsrt Transform the covariance matrix of mrqmin.

déplie Vector to matrix transformation.

gaussj Linear equation solution by Gauss-Jordan elimination.

kregrid Matrix (m x 2) of 2-D grid coordinates.

kregrid3 Matrix (m x 3) of 3-D grid coordinates.

kridemo Kriging Toolbox demo.

ksone MEX-file called from kstest.

kstest Kolmogorov-Smirnov normality test.

mat3dp Get a value out of a pseudo 3-D matrix.

mat4dp Get a value out of a pseudo 4-D matrix.

means Mean function called from cokri2.

mrqcof M-file called from mrqmin.

trans Translation function called from cokri2.

Variogram Options

Available variograms are:

1. Semivariogram:

2. Cross semivariogram:

3. Covariance:

4. Correlogram:

5. General relative semivariogram:

6. Pairwise relative semivariogram:

7. Semivariogram of logarithms:

8. Semirodogram:

9. Semimadogram:

10. Semivariogramme indicator:

Variable descriptions:

C(h) = covariance as a function of distance h; C(0) gives the variance.

g(h) = semivariance = [ C(0) - C(h) ].

h = separation vector.

N(h) = number of sample pairs

xi,yi = value of the sample pair separated by vector h: xi is the value at the start (or tail) and yi is the value at the end (or head) of interval h.

zi,zi' et yi,yi' = same as (xi,yi) in the cross semivariogram: yi and zi are the values at the start and yi' and zi' are the values at the end of interval h.

For more information, please consult Deutsch and Journel (1992).

Kriging Options

Available kriging options are:

simple cokriging

ordinary cokriging with one nonbias condition (Isaaks and Srivastava)

ordinary cokriging with p nonbias condition

universal cokriging with drift of order 1

universal cokriging with drift of order 2

99. cokriging is not performed, only sample variance sv is computed

Cokriging means kriging with more than one variable. When the cokriging program is called with only one variable at a time, the results will be those of simple kriging, ordinary kriging, universal kriging, point kriging or block kriging. More details can be found in the paper of Marcotte (1991).

Chi Toolbox Contents

This toolbox is used by the normality test. The included functions compute the Chi-squared probability function and the percentage points of the probability function. The m-files were downloaded from the matlab public site:

http://www.mathworks.com/ftp/statv4/shtml

Probability function c2

chiprob Probability of observing a given c2value.

chitable c2 value for a given probability.

chiaux Function called from chitable.

Toolbox author:

Peter R. Shaw

Woods Hole Oceanographic Institution, Woods Hole, MA 02543

(508) 457-2000 ext. 2473

pshaw@whoi.edu

Note: Optimization Toolbox needed

barnes

Purpose

2-D spatial filter for kriged data.

Synopsis

F = barnes (xi, yi, zi, c, g)

Description

The Barnes' filter is a low-pass 2-D filter whose mathematical description is:

where , ,

and .

Input variables:

xi: column grid coordinates

yi: row grid coordinates

zi: grid data

c, g: filter parameters

Output variable:

F: filtered data

Example

The tintore function provides a good example of spatial filters in physical oceanography.

References

Maddox, R.A., 1980. An Objective Technique for Separating Macroscale and Mesoscale Features in Meteorological Data. Monthly Weather Rev., 108: 1108-1121.

Tintoré et al., 1991. Mesoscale Dynamics and Vertical Motion in the Alboran Sea. J. Phys. Oceanogr., 21: 811-823.

See also

tintore, filresp

cokri / cokri2

Purpose

cokri: Point or block cokriging in D dimensions.

cokri2: Cokriging function called from cokri.

Synopsis

[x0s, s, sv, id, l] = cokri (x, x0, model, c, itype, avg, block, nd, ival, nk, rad, ntok, d)

[x0s, s, id, l, k0] = cokri2 (x, x0, id, model, c, sv, itype, avg, ng, d)

Description

Cokriging means kriging with more than one variable. When cokri is called with only one variable at a time, the results will be those of simple kriging, ordinary kriging, universal kriging, point kriging or block kriging. More details can be found in the paper of Marcotte (1991).

Available kriging options are:

simple cokriging

ordinary cokriging with one nonbias condition (Isaaks and Srivastava)

ordinary cokriging with p nonbias condition

universal cokriging with drift of order 1

universal cokriging with drift of order 2

99.cokriging is not performed, only sv is computed

Available semivariogram models are:

nugget effect

exponential model

gaussian model

spherical model

linear model

quadratic model

power (hd) model

logarithmic

sinc(h)

Bessel [ Jo (h) ]

exp(-h) * cos(dh)

exp(-h) * Jo(dh)

exp(-h2) * cos(dh)

exp(-h2) * Jo(dh)

exp(-h2) * (1 - dh2)

1 - 3*min(h,1) + 2*min(h,1)

h * log(max(h,eps))

New models can be added quite easily since models are calculated using the eval function.

Input variables:

x: data matrix [x y z var1 var2 ...]

x0: grid coordinates [xi yi zi]

model: [models, a (h=r/a), rotation angles].

No rotation angle is required for an isotropic distribution.

Example: model = [1 10; 4 30] means that the distribution is isotropic and that it is represented by a nugget effect of range 10 plus a spherical model of range 30.

c: amplitudes of the models

itype: kriging option

block: vector (1 x D), giving the size of the block to estimate; for point cokriging: block = ones(1,D)

nd: Vector (1 x D), giving the discretization grid for block cokriging; for point cokriging: nd = ones(1,D)

ival: Code for cross-validation.

0: no cross-validation

1: cross-validation is performed by removing one variable at a

time at a given location.

2: cross-validation is performed by removing all variables at a

given location.

nk: number of nearest neighbors in x matrix to use in the cokriging.

rad: search radius.

ntok: points in x0 will be kriged by groups of ntok grid points.

d: model coefficients. This coefficient has been added to the original Marcotte's function. Warning: In cokri, models are defined in terms of h = r/a. In variogr, the dependent variable is r and hence d = b*a (see variogr).

Output variables

x0s: kriged data matrix at x0 locations.

s: kriged data variance matrix at x0 locations.

sv: variance of each variable.

id, l see Marcotte (1991)

Reference

Marcotte, D. 1991. Cokrigeage with matlab. Computers & Geosciences. 17(9): 1265-1280.

See also

cokri2, variogr, trans, means

confint

Purpose

Confidence intervals.

Synopsis

[k2, k1] = confint (g, m, S2)

Description

Confidence intervals for the structure function

CONF {k2 variance k1} (1)

The structure function is a measure of the variance of a given variable as a function of distance. The estimation of the confidence intervals in such a case is given by (1).

k1 = (n-1) * S2 / c1

k2 = (n-1) * S2 / c2 (2)

where n = sample size = m+1

m = number of degrees of freedom

S2 = variance of the sample

c1 and c2 are determined by the solution to the equations

F(c1) = (1-g) /2

F(c2) = (1+g) /2 (3)

where g = confidence level (95%, 99% or the like)

F = distribution

Solutions are obtained by function chitable (Chi Toolbox).

References

Denman, K. L. and H. J. Freeland (1985). Correlation Scales, Objective Mapping and a Statistical Test of Geostrophy over the Continental Shelf. J. of Mar. Res., 43: 517-539.

Kreyszig, E., 1988. Advanced Engineering Mathematics, sixth ed., John Wiley & Sons, New York, p.1252

See also

chitable

davis

Purpose

Point kriging using Davis' set of equations A W = B

where , and .

The g(hik) is the semivariance of sample pairs separated by distance hik. Non bias conditions require the sum of Wi to be equal to 1. In that case, one more degree of freedom must be introduced with the use of a Lagrange multiplier in order to minimize the estimation error.

Synopsis

[Zp, Sp] = davis (data, x0, model, a, d, c, A)

Description

Available models are:

nugget effect

exponential model

gaussian model

spherical model

linear model

quadratic model

power model (hd)

logarithmic model

sinc(h)

Bessel [ Jo (h) ]

exp(-h) * cos(dh)

exp(-h) * Jo(dh)

exp(-h2) * cos(dh)

exp(-h2) * Jo(dh)

exp(-h2) * (1 - dh 2)



Input variables:

data: data [x y variable]

x0: grid coordinates [xi yi]

model: semivariogram model

a: semivariogram range

d: model coefficient (different from coefficient b of variogr, same as d in cokri)

c: model amplitude

A: g(hik) matrix if already calculated; if not, ignore that input.

Output variables:

Zp: kriged data matrix at x0 positions

Sp: variance of kriged data at x0 positions

Reference

Davis, J.C. 1986. Statistics and Data Analysis in Geology, 2nd ed., John Wiley & Sons, New York, 289 p.

deplie

Purpose

Vector to matrix transformation.

Synopsis

mat = deplie (y, nx, ny)

Description

Transformation of a vector y into a matrix mat of size ny x nx.

Input variables:

y: row or column vector

nx: number of columns in matrix mat

ny: number of rows in matrix mat

Output variable:

mat: matrix

Example

y = [(1:10) (1:10) (1:10)];

deplie (y, 10, 3) = 1 2 3 4 5 6 7 8 9 10

1 2 3 4 5 6 7 8 9 10

1 2 3 4 5 6 7 8 9 10

See also

kregrid

filresp

Purpose

Barnes' filter response in the wavelength domain.

Synopsis

R = filresp (c, g)

Description

Barnes' filter response in the wavelength domain is given by:

where , c and g are the filter parameters and l is the wavelength.

Example

f2 = filresp(200,0.6);

References

Barnes, S.L.,1973. Mesoscale Objective Map Analysis Using Weighted Time Series Observations. NOAA Tech. Memo. ERL NSSL-62, 60 p.

Maddox, R.A., 1980. An Objective Technique for Separating Macroscale and Mesoscale Features in Meteorological Data. Mon. Wea. Rev., 108, 1108:1121.

See also

barnes, tintore

fitvario / fitvari2 / fun

Purpose

fitvario: Optimization of variogr.

fitvari2: Optimization of variogr2 (without the Optimization Toolbox).

fun: Called from fitvario to estimate variograms.

Synopsis

fitvario (model, data, a)

fitvario (model, data, a, b)

fitvari2 (model, data, a, b, C)

f = fun (lam, data)

Description

Least-square fitting of semivariogram model coefficients a, b and C.

Input variables:

model: model type (see variogr)

data(:,1): x-axis (distance) (gam(:,1))

data(:,2): y-axis (variance) (gam(:,2))

a, b et C: starting values of coefficients a, b and C of variogr

Output:

The graphical output shows the plot of the semivariance as a function of the distance. Experimental results appear as symbols. The plain line gives the best fit for the model chosen. The optimal values of a, b and C are also displayed on the graph.

Example

r = (0:10)';

x = 1 - exp(-r);

fitvario(2,[r x],2)

fitvari2(2,[r x],2,0,1.5)



See also

leastsq into the Optimization Toolbox, variogr, variogr2

gam2, gamv2, gamv2uv gam3, gamv3

Purpose

MEX-file called from vario2dr, vario2di, var2diuv, vario3dr and vario3di.

Synopsis

[np, gam, hm, tm, hv, tv] = gam2 (nlag, nx, ny, ndir, ixd, iyd, vr, tmin, tmax, nvarg, ivtail, ivhead, ivtype, nvar)

[np, dis, gam, hm, tm, hv, tv] = gamv2 (nd, x, y, vr, tmin, tmax, nlag, xlag, xltol, ndir, azm, atol, bandw, nvarg, ivtail, ivhead, ivtype, nvar)

[np, dis, gam, hm, tm, hv, tv] = gamv2uv (nd, x, y, u, v, tmin, tmax, nlag, xlag, xltol, ndir, azm, atol, bandw, nvarg, ivtail, ivhead, ivtype, nvar, option)

[np, gam, hm, tm, hv, tv] = gam3 (nlag, nx, ny, nz, ndir, ixd, iyd, izd, vr, tmin, tmax, nvarg, ivtail, ivhead, ivtype, nvar)

[np, dis, gam, hm, tm, hv, tv] = gamv3 (nd, x, y, z, vr, tmin, tmax, nlag, xlag, xltol, ndir, azm, atol, bandwh, dip, dtol, bandwd, nvarg, ivtail, ivhead, ivtype, nvar)

Description

The description of inputs and outputs is given in vario2dr, vario2di, var2diuv, vario3dr and vario3di. These functions are MEX-files compiled from source Fortran codes gam2.for, gamv2.for, gam3.for and gamv3.for of GSLIB (Deutsch et Journel, 1992).

Reference

Deutsch, C. V and A. G. Journel,1992. GSLIB: Geostatistical Software Library and User's Guide. Oxford University Press, Oxford, 340 p.

See also

vario2dr, vario2di, vario3dr, vario3di, var2diuv

kregrid / kregrid3

Purpose

Matrix (m x 2) of 2-D grid coordinates.

Matrix (m x 3) of 3-D grid coordinates.

Synopsis

y = kregrid (x0, dx, xf, y0, dy, yf)

y = kregrid3 (x0, dx, xf, y0, dy, yf, z0, dz, zf)

Description

Grid (x,y) coordinates are reshape in an m x 2 matrix.

Grid (x,y,z) coordinates are reshape in an m x 3 matrix.

Input variables:

(x0,y0,z0): (x,y,z) position of the lower left corner of the grid

(xf, yf, zf): (x,y,z) position of the upper right corner of the grid

dx: x grid interval

dy: y grid interval

dz: z grid interval

Output variable:

y: m x 2 or m x 3 matrix of grid coordinates

Example

Let us say that we want to generate the following grid:

y = kregrid (1, 1, 3, 0, 5, 5)

y =

1 0

2 0

3 0

1 5

2 5

3 5

ksone

Purpose

MEX-file called from kstest.

Synopsis

[d, prob] = ksone (sample, n, normal)

Description

MEX-file compiled from the source Fortran code ksone.for, a Numerical Recipes subroutine. A comparison is made between the sample cumulative distribution and a normal cumulative distribution.

Input variables:

sample: standardized sample (mean(sample) = 0 et std(sample) = 1)

n: number of data in the sample

normal: normal cumulative distribution ranging from 0 to 1

Output variable:

d: K-S statistic

prob: significance level. Small values show that the sample cumulative distribution is significantly different from normal.

Reference

Press, W et al. 1992. Numerical Recipes in Fortran, The Art of Scientific Computing, second ed., Cambridge University Press, Cambridge, p 619.

See also

kstest

kstest

Purpose

Kolmogorov-Smirnov normality test.

Synopsis

[d, prob] = kstest (sample)

Description

Normality test with the ksone MEX-file.

Input variables:

sample: sample

Output variable:

d: K-S statistic

prob: significance level. Small values show that the sample cumulative distribution is significantly different from normal.

Example

Normality test of a vector of normally distributed random numbers:

data = randn(100,1);

[d,prob] = kstest(data)

d =

0.0456

prob =

0.9854

References

Kreyszig, E., 1988. Advanced Engineering Mathematics, sixth ed., John Wiley & Sons, New York, p.1211.

Legendre, L. and Legendre, P. (1983) Numerical Ecology. Developments in Environment Modeling 3. Elsevier, New York, 419p.

See also

ksone

mat3dp

Purpose

Get a value out of a pseudo 3-D matrix.

Synopsis

r = mat3d ( mat, d3, i, j, k )

Description

A pseudo 3-D matrix is a 2-D matrix made up of successive 2-D slices.

Input variables:

mat: pseudo 3-D matrix

d3: number of levels in the third dimension

i: position in the first dimension

j: position in the second dimension

k: position in the third dimension

Output variable:

r: value corresponding to mat(i,j,k)

Example

x =

1.0 1.2

1.5 1.7

2.0 2.2

2.5 2.7

3.0 3.2

3.5 3.7

mat3dp (x, 3, 1, 2, 2) = 2.2

mat4dp

Purpose

Get a value out of a pseudo 4-D matrix.

Synopsis

r = mat4d (mat, d3, d4, ii, jj, kk, ll)

Description

A pseudo 4-D matrix is a 2-D matrix made up of successive pseudo 3-D matrices.

Input variables:

mat: pseudo 4-D matrix

d3: number of levels in the third dimension

d4: number of levels in the fourth dimension

ii: position in the first dimension

jj: position in the second dimension

kk: position in the third dimension

ll: position in the fourth dimension

Output variable:

r: value corresponding to mat(ii,jj,kk,ll)

Example

x =

1.0 1.2

1.5 1.7

2.0 2.2

2.5 2.7

3.0 3.2

3.5 3.7

4.0 4.2

4.5 4.7

mat4dp (x, 2, 2, 1, 2, 2, 1) = 2.2

means

Purpose

Mean function called from cokri2.

Synopsis

y = means (x)

Description

Average or mean value. The only difference with matlab function mean is for a row vector. Means returns the same row vector instead of the mean value of the elements of the row.

Input variable:

x: a vector or a matrix

Output variable:

y: row vector of the mean of each column

Example

Let us consider the following matrix:

x = 1.0 1.5 2.0 2.5 3.0 3.5

1.2 1.7 2.2 2.7 3.2 3.7

means (x) = 1.1 1.6 2.1 2.6 3.1 3.6

mean (x) = 1.1 1.6 2.1 2.6 3.1 3.6

whereas

means (x(1,:)) = 1.0 1.5 2.0 2.5 3.0 3.5

mean (x(1,:)) = 2.25

Reference

Marcotte, D., 1991. Cokrigeage with MATLAB. Computers & Geosciences. 17(9): 1265-1280.

mrqmin / mrqcof / covsrt / gaussj

Purpose

mrqmin: Least-square fitting: Levenberg-Marquardt method.

mrqcof: M-file called from mrqmin.

covsrt: Transform the covariance matrix of mrqmin.

gaussj: Linear equation solution by Gauss-Jordan elimination.

Synopsis

[covari, alpha, chisq, alamda, a] = mrqmin (x ,y ,sig ,ndata, a, ia, ma, nca, funcs, alamda, model)

[alpha, beta, chisq] = mrqcof (x, y, sig, ndata, a, ia, ma, nalp, funcs, model)

covari = covsrt (npc, ma, ia, mfit, covari)

[a, b] = gaussj (a, n, np, b, m, mp)

Description

Levenberg-Marquardt method, attempting to reduce the value c of a fit between a set of data points x(ndata), y(ndata) with standard deviations sig(ndata), and a nonlinear function dependent on ma coefficients a(1:ma).

Reference

Press, W et al. 1992. Numerical Recipes in Fortran, The Art of Scientific Computing, second ed., Cambridge University Press, Cambridge, p 680.

See also

fitvario2

outvario

Purpose

Output of variogram functions.

Synopsis

[np, gam, hm, tm, hv, tv] = outvario (nlg, in7, ndir, nvarg, in1, in2, in3, in4, in5,... in6, ivtype)

Description

This function is only used to reorganize the outputs of gam2, gamv2, gam3 and gamv3.

See also

vario2di, vario2dr, var2diuv, vario3di, vario3dr

tintore

Purpose

Application of Barnes' filter with the Tintoré's parameters (1991).

Synopsis

[F2, Fb] = tintore (xi, yi, zi)

Description

This function is an example of Barnes' filter to separate mesoscale from macroscale features in physical oceanography. The filter parameters are those established by Tintoré et al. (1991).

Input variables:

xi: vector of the x grid coordinates (positions of the columns of zi)

yi: vector of the y grid coordinates (positions of the rows of zi)

zi: grid data (size(zi) = [length(yi), length(xi)])

Output variables:

F2: macroscale structure

Fb: mesoscale structure

Reference

Tintoré et al., 1991. Mesoscale Dynamics and Vertical Motion in the Alboran Sea, J. Phys. Oceanogr., 21:811-823.

See also

barnes, filresp

trans

Purpose

Translation function called from cokri2.

Synopsis

cx = trans (cx, model, im)

Description

This function takes as input original coordinates and returns the rotated and reduced coordinates following specifications described in the semivariogram model.

Reference

Marcotte, D., 1991. Cokrigeage with MATLAB. Computers & Geosciences. 17 (9): 1265-1280.

See also

cokri2

variogr / variogr2

Purpose

variogr: Models of semivariogram and correlogram.

variogr2: variogr for fitting procedures not using the Optimization Toolbox.

Synopsis

y = variogr (type, r, a, C, b) for semivariograms

y = variogr (-type, r, a, C, b) for correlograms

y = variogr2 (type, r, a, C, b) for semivariograms

y = variogr2 (-type, r, a, C, b) for correlograms

Description

This function is used to obtain the theoretical curve of the variance of a sample as a function of the sample pair distance.

Available model options are:

With a sill

1. spherical

2. exponential

3. gaussian

4. quadratic (in preparation)

Without a sill

5. linear

10. logarithmic (in preparation)

11. power of r (in preparation)

Hole effects

20. C * ( 1 - sin(b*r) / r )

21. C * ( 1 - exp(-r/a) * cos(b*r) )

22. C * ( 1 + exp(-r/a) * cos(b*r) )

23. C * ( 1 - exp(-r/a) * cos(r*b) )

24. C * (1 - exp(-(r/a)2) * cos(b*r) )

25. C * ( 1 - J 0 (b*r) )

26. C * ( 1 - J0 (b*r) * exp(-r/a) )

27. C * ( 1 - exp(-(r/a)2) * J0 (b*r) )

28. C * ( 1 - exp(-(r/a)2) * (1 - br2)



Any new types can be easily added to the list.

Input variables:

r: distance vector

a: range

C: amplitude

b: coefficient used in the models 20

type: variogram type

Output variable:

y: variance as a function of r

Reference

Journel, A.G. and C.J. Huijbregts, 1992. Mining Geostatistics. Academic Press,

New York, 600 p.

See also

fitvario

vario2di

Purpose

Variogram of irregularly spaced 2-D data.

Synopsis

[np, gam, hm, tm, hv, tv] = vario2di (nd, x, y, vr, nlag, xlag, xltol, ndir, azm, atol, bandw, nvarg, ivtail, ivhead, ivtype, nvar)

Description

This function is used to calculate the variograms of a sample which is irregularly distributed on a plane.

Input variables:

nd: number of data (no missing values)

x(nd): x coordinates of the data set

y(nd): y coordinates of the data set

vr(nd,nvar): data values

nlag: number of lags to calculate

xlag: length of the unit lag

xltol: distance tolerance (if <0 then set to xlag/2)

ndir: number of directions to consider

azm(ndir): azimuth angle of direction (measured positive degrees clockwise from NS).

atol(ndir): azimuth (half window) tolerances

bandw(ndir): maximum bandwidth

nvarg: number of variograms to compute

ivtail(nvarg): variable for the tail of each variogram

ivhead(nvarg): variable for the head of each variogram

ivtype(nvarg): type of variogram to compute

nvar: number of variables

Output variables:

np: number of pairs

gam: semivariogram, covariance, correlogram,... values

hm: mean of the tail data

tm: mean of the head data

hv: variance of the tail data

tv: variance of the head data

The first column of each output variable is the vector of the corresponding distances. The first row of each output variable is the vector of the corresponding variogram numbers. If more than one direction is calculated, the values of the output variables for the other directions are added at the end of the outputs of the first one. More details can be found in Deutsch et Journel (1992).

Reference

Deutsch, C. V and A. G. Journel (1992) GSLIB: Geostatistical Software Library and User's Guide. Oxford University Press, Oxford, 340 p.

vario2dr

Purpose

Variogram of regularly spaced 2-D data.

Synopsis

[np, gam, hm, tm, hv, tv] = vario2dr (nlag, nx, ny, ndir, ixd, iyd, vr, nvarg, ivtail, ivhead, ivtype, nvar)

Description

This function is used to calculate the variograms of a sample which is regularly distributed on a plane.

Input variables:

nlag: number of lags to calculate

nx: number of units in x (number of columns)

ny: number of units in y (number of lines)

ndir: number of directions to consider

ixd(ndir): x indicator of direction - number of x grid columns that must be shifted to move from a node on the grid to the next nearest node on the grid which lies on the directional vector.

iyd(ndir): y indicator of direction - similar to ixd, number of grid lines that must be shifted to nearest node which lies on the directional vector

vr(nx,ny*nvar): 3-D array of data.

nvarg: number of variograms to compute

ivtail(nvarg): variable for the tail of each variogram

ivhead(nvarg): variable for the head of each variogram

ivtype(nvarg): type of variogram to compute

nvar: number of variables

Output variables:

np: number of pairs

gam: semivariogram, covariance, correlogram,... values

hm: mean of the tail data

tm: mean of the head data

hv: variance of the tail data

tv: variance of the head data

The first column of each output variable is the vector of the corresponding distances. The first row of each output variable is the vector of the corresponding variogram numbers. If more than one direction is calculated, the values of the output variables for the other directions are added at the end of the outputs of the first one. More details can be found in Deutsch et Journel (1992).

Reference

Deutsch, C. V and A. G. Journel (1992) GSLIB: Geostatistical Software Library and User's Guide. Oxford University Press, Oxford, 340 p.

vario3di

Purpose

Variogram of irregularly spaced 3-D data.

Synopsis

[np, gam, hm, tm, hv, tv] = vario3di (nd, x, y, z, vr, nlag, xlag, xltol, ndir, azm, atol, bandwh, dip, dtol, bandwd, nvarg, ivtail, ivhead, ivtype, nvar)

Description

This function is used to calculate the variograms of a sample which is irregularly distributed in space.

Input variables:

nd: number of data

x(nd): x coordinates of the data set

y(nd): y coordinates of the data set

z(nd): z coordinates of the data set

vr(nd,nvar): data values

nlag: number of lags to calculate

xlag: length of the unit lag

xltol: distance tolerance (if <0 then set to xlag/2)

ndir: number of directions to consider

azm(ndir): azimuth angle of direction (measured positive degrees clockwise from NS)

atol(ndir): azimuth (half window) tolerances

bandwh: maximum horizontal bandwidth

dip(ndir): dip angle of direction (measured in negative degrees down from horizontal)

dtol(ndir): dip (half window) tolerances

bandwd: maximum vertical bandwidth

nvarg: number of variograms to compute

ivtail(nvarg): variable for the tail of each variogram

ivhead(nvarg): variable for the head of each variogram

ivtype(nvarg): type of variogram to compute

nvar: number of variables

Output variables:

np: number of pairs

gam: semivariogram, covariance, correlogram,... values

hm: mean of the tail data

tm: mean of the head data

hv: variance of the tail data

tv: variance of the head data

The first column of each output variable is the vector of the corresponding distances. The first row of each output variable is the vector of the corresponding variogram numbers. If more than one direction is calculated, the values of the output variables for the other directions are added at the end of the outputs of the first one. More details can be found in Deutsch et Journel (1992).

Reference

Deutsch, C. V and A. G. Journel (1992) GSLIB: Geostatistical Software Library and User's Guide. Oxford University Press, Oxford, 340 p.

vario3dr

Purpose

Variogram of regularly spaced 3-D data.

Synopsis

[np, gam, hm, tm, hv, tv] = var2diuv (nd, x, y, u, v, nlag, xlag, xltol, ndir, azm, atol, bandw, nvarg, ivtail, ivhead, ivtype, nvar,option)

Description

This function is used to calculate the variograms of a sample which is regularly distributed in space.

Input variables:

nlag: maximum number of lags to be calculated

nx: number of units in x (number of columns)

ny: number of units in y (number of lines)

nz: number of units in z (number of levels)

ndir: number of directions to consider

ixd(ndir): x (column) indicator of direction - number of grid columns that must be shifted to move from a node on the grid to the next nearest node on the grid which lies on the directional vector

iyd(ndir): y (line) indicator of direction - similar to ixd, number of grid lines that must be shifted to nearest node which lies on the directional vector

izd(ndir): z (level) indicator of direction - similar to ixd, number of grid levels that must shifted to nearest node of directional vector

nv: the number of variables

vr(nx,ny*nz*nvar): 4-D array of data

nvarg: number of variograms to compute

ivtail(nvarg): variable for the tail of the variogram

ivhead(nvarg): variable for the head of the variogram

ivtype(nvarg): type of variogram to compute

nvar: number of variables

Output variables:

np: number of pairs

gam: semivariogram, covariance, correlogram,... values

hm: mean of the tail data

tm: mean of the head data

hv: variance of the tail data

tv: variance of the head data

The first column of each output variable is the vector of the corresponding distances. The first row of each output variable is the vector of the corresponding variogram numbers. If more than one direction is calculated, the values of the output variables for the other directions are added at the end of the outputs of the first one. More details can be found in Deutsch et Journel (1992).

Reference

Deutsch, C. V and A. G. Journel (1992) GSLIB: Geostatistical Software Library and User's Guide. Oxford University Press, Oxford, 340 p.

var2diuv

Purpose

Variogram of irregularly spaced 2-D vectors.

Synopsis

[np, gam, hm, tm, hv, tv] = var2diuv (nd, x, y, u, v, nlag, xlag, xltol, ndir, azm, atol, bandw, nvarg, ivtail, ivhead, ivtype, nvar, option)

Description

This function is used to calculate the variograms of a vector sample which is irregularly distributed on a plane. The option 'parallel' or 'perpendi' allows you to chose the vector component. The 'parallel' component is the one parallel to the plane whereas the 'perpendi' component is perpendicular to the plane. This function is used to krige vertical sections of horizontal currents in physical oceanography.

Input variables:

nd: number of data (no missing values)

x(nd): x coordinates of the data set

y(nd): y coordinates of the data set

u(nd,nvar): east-west current

v(nd,nvar): north-south current

nlag: number of lags to calculate

xlag: length of the unit lag

xltol: distance tolerance (if <0 then set to xlag/2)

ndir: number of directions to consider

azm(ndir): azimuth angle of direction (measured positive degrees clockwise from NS).

atol(ndir): azimuth (half window) tolerances

bandw(ndir): maximum bandwidth

nvarg: number of variograms to compute

ivtail(nvarg): variable for the tail of each variogram

ivhead(nvarg): variable for the head of each variogram

ivtype(nvarg): type of variogram to compute:

nvar: number of variables

option: 'parallel' or 'perpendi'

Output variables:

np: number of pairs

gam: semivariogram, covariance, correlogram,... values

hm: mean of the tail data

tm: mean of the head data

hv: variance of the tail data

tv: variance of the head data

References

Denman, K.L. and H.J. Freeland, 1985. Correlation Scales, Objective Mapping and a Statistical Test of Geostrophy over the Continental Shelf. J. Mar. Res., 43: 517-539.

Deutsch, C. V and A. G. Journel (1992) GSLIB: Geostatistical Software Library and User's Guide. Oxford University Press, Oxford, 340 p.