Calculates kriged estimates using a model fitted to the sample variogram.
Options
PRINT = string token |
Controls printed output (description, search, weights, monitor, data); default desc |
|---|---|
Y = variate |
Y positions (not needed for 2-dimensional regular data i.e. when DATA is a matrix) |
X = variate |
X positions (needed only for 2-dimensional irregular data) |
YOUTER = variate |
Variate containing 2 values to define the Y-bounds of the region to be examined (bottom then top); by default the whole region is used |
XOUTER = variate |
Variate containing 2 values to define the X-bounds of the region to be examined (left then right); by default the whole region is used |
YINNER = variate |
Variate containing 2 values to define the Y-bounds of the interpolated region (bottom then top); no default |
XINNER = variate |
Variate containing 2 values to define the X-bounds of the interpolated region (left then right); no default |
BLOCK = variate |
Dimensions (length and height) of block; default !(0, 0) i.e. punctual kriging |
RADIUS = scalar |
Maximum distance between target point in block and usable data |
SEARCH = string token |
Type of search (isotropic, anisotropic); default isot |
MINPOINTS = scalar |
Minimum number of data points from which to compute elements; default 7 |
MAXPOINTS = scalar |
Maximum number of data points from which to compute elements (2 < MINPOINTS ≤ MAXPOINTS < 41); default 20 |
NSTEP = scalar |
Number of steps for numerical integration; (3 < NSTEP < 11); default 8 |
DRIFT = string token |
Amount of drift (constant, linear, quadratic); default cons |
YXRATIO = scalar |
Ratio of Y interval to X interval; default 1.0 |
INTERVAL = scalar |
Distance between successive interpolations; default 1.0 |
Parameters
DATA = variates or matrices |
Observed measurements as a variate or, for data on a regular grid, as a matrix |
|---|---|
ISOTROPY = string tokens |
Form of variogram (isotropic, Burgess, geometrical); default isot |
MODELTYPE = string tokens |
Model fitted to the variogram (power, boundedlinear, circular, spherical, doublespherical, pentaspherical, exponential, besselk1, gaussian, cubic, stable, cardinalsine, matern); default powe |
NUGGET = scalars |
The nugget variance |
SILLVARIANCES = variates |
Sill variances of the spatially dependent component; default none |
RANGES = variates |
Ranges of the spatially dependent component; default none |
GRADIENT = variates |
Slope of the unbounded component; default none |
EXPONENT = variates |
Power of the unbounded component or power for the stable model; default none |
SMOOTHNESS = scalar |
Value of ν parameter for the Matern model; defalt none |
PHI = variates |
Phi parameters of an anistropic model (ISOTROPY = Burg or geom) |
RMAX = variates |
Maximum gradient or distance parameter of an anistropic model |
RMIN = variates |
Minimum gradient or distance parameter of an anistropic model |
PREDICTIONS = matrices |
Kriged estimates |
VARIANCES = matrices |
Estimation variances |
LAGRANGEMULTIPLIER = matrices or pointers |
Saves the Lagrange multipliers from each kriging solution |
MEASUREMENTERROR = scalar |
Specifies the variance of the measurement error |
SAVE = pointers |
Supplies the model name and estimates, as saved from MVARIOGRAM |
Description
The KRIGE directive computes the ordinary kriging estimates of a variable at positions on a grid from data and a model variogram. The data must be supplied, using the DATA parameter, in one of the two forms as for the FVARIOGRAM procedure: i.e. for data on a regular grid, in a matrix defined with a variate of column labels to provide the x-values and a variate of row labels to provide the y-values or, for irregularly scattered data, as a variate with the X and Y options set to variates to supply the spatial coordinates.
By default all data are considered when forming the kriging system. However, a subset of the data may be selected by limiting the area to a rectangle defined by XOUTER and YOUTER options. Each of these should be set to a variate with two values to define lower and upper limits in the x (East-West) and y (North-South) directions respectively.
The positions at which Z is predicted (estimated) are contained in a rectangle defined by the XINNER, YINNER and INTERVAL options. XINNER and YINNER are set to variates similarly to XOUTER and YOUTER, and their limits should not lie outside those of XOUTER and YOUTER. INTERVAL is set to a scalar to define the distance between the successive positions in the rows and columns of the grid at which kriging is to be done, specified in the same units as the data. However, if the aim is to make a map, INTERVAL should be chosen so that it represents no more than 2 mm on the final printed document. The optimality of the kriging will then not be degraded noticeably by the subsequent contouring.
Kriging may be either punctual, i.e. at “points” which have the same size and shape as the sample support, or on bigger rectangular blocks. The size of the blocks is specified by the BLOCK option, in a variate whose two values define the length of the block first in the x direction (eastings) and then in the y direction (northings). By default the BLOCK variate contains two zero values, to give punctual kriging. The average semivariances between point and block are computed by integrating the variogram numerically over the block. The number of steps in each direction is defined by the NSTEP option. The default of 8 is recommended as a compromise between speed and accuracy. The kriging may be accelerated at the expense of accuracy by reducing NSTEP, or accuracy gained by increasing it. The minimum is 4 and the maximum 10.
The minimum and maximum number of points for the kriging system are set by the MINPOINTS and MAXPOINTS options. There is a minimum limit of 3 for MINPOINTS and a maximum of 40 for MAXPOINTS, and MINPOINTS must be less than or equal to MAXPOINTS. The defaults are 7 and 20 respectively. Data points may be selected around the point or block to be kriged by setting the RADIUS option to the radius within which they must lie. If the variogram is anisotropic, the search may be requested to be anisotropic by setting option SEARCH to anisotropic; by default SEARCH=isotropic.
Universal kriging may be invoked by setting the DRIFT option to linear or to quadratic, i.e. to be of order 1 or 2 respectively. By default is DRIFT=constant, to give ordinary kriging. For data in a regular grid that is not square, the ratio of the spacing in the y direction to that in the x direction is given by the YXRATIO option. The default is 1.0 for square.
The variogram is specified by its type and parameters. The model and estimates can be saved using the SAVE parameter of MVARIOGRAM, and passed on to KRIGE using its SAVE parameter. Alternatively, they can be supplied as follows.
The model can be defined by setting the MODELTYPE option to either power, boundedlinear (one dimension only), circular, spherical, doublespherical, pentaspherical, exponential, besselk1 (Whittle’s function), gaussian, cubic, stable (i.e. powered exponential; see Webster & Oliver 2001), cardinalsine (Chiles & Delfiner 1999) or matern. All models may have a nugget variance, supplied using the NUGGET option; this is the constant estimated by MVARIOGRAM. For punctual kriging, you can specify the variance of any measurement error using the MEASUREMENTERROR parameter. The parameters of the power function (the only unbounded model) are defined by the GRADIENT and EXPONENT parameters. The parameter for the power of the stable model is supplied using the EXPONENT parameter. The parameter ν for the matern function is supplied using the SMOOTHNESS parameter. The simple bounded models, i.e. all other settings of MODELTYPE except doublespherical, require the SILLVARIANCES (the sill of the correlated variance) and RANGES parameters. The latter is strictly the correlation range of the boundedlinear, circular, spherical and pentaspherical models, while for the asymptotic models it is the distance parameter of the model. The doublespherical model requires SILLVARIANCES and RANGES to be set to variates of length two, to correspond to the two components of the model.
The ISOTROPY parameter allows the variation to be defined to be either isotropic or anisotropic in one of two ways: either Burgess anisotropy (Burgess & Webster 1980) or geometric anisotropy (Journel & Huijbregts 1978, Webster & Oliver 1990). The anisotropy is specified by three parameters, namely PHI, the angle in radians of the direction of maximum variation, RMAX, the maximum gradient or distance parameter of the model, and RMIN, the minimum gradient or distance parameter. The power, stable, exponential, Gaussian, pentashperical, spherical, cubic and circular functions may be anisotropic.
KRIGE calculates two matrices, one of predictions (or estimates), which can be saved using the PREDICTIONS parameter, and the other of the prediction (estimation or kriging) variances saved using the VARIANCES parameter. The matrices are arranged with the first row of each matrix at the bottom following geographic rather than mathematical convention. You can save the Lagrange multipliers from the kriging solution using the LAGRANGEMULTIPLIER parameter. For ordinary Kriging the Lagrange multipliers are saved in a matrix (with a multiplier for each point). For universal Kriging a pointer of matrices is saved, where a matrix to save the Lagrange multipliers of each equation term.
The PRINT option can be set to data to print the data (2-dimensional regular data only). It also allows intermediate results to be printed. The setting search lists the results of the search for data around each position to be kriged, weights lists the kriging weights at each position and monitor monitors the formation and inversion of the kriging matrices for each position. These options enable you to check that the kriging is working reasonably. However, they can produce a great deal of output, and should not be requested when kriging large matrices, such as might be wanted for mapping.
Options: PRINT, Y, X, YOUTER, XOUTER, YINNER, XINNER, BLOCK, RADIUS, SEARCH, MINPOINTS, MAXPOINTS, NSTEP, DRIFT, YXRATIO, INTERVAL.
Parameters: DATA, ISOTROPY, MODELTYPE, NUGGET, SILLVARIANCES, RANGES, GRADIENT, EXPONENT, SMOOTHNESS, PHI, RMAX, RMIN, PREDICTIONS, VARIANCES, LAGRANGEMULTIPLIER, MEASUREMENTERROR, SAVE.
Action with RESTRICT
You can restrict any of the DATA variate to do the estimation using only a subset of the units. If more than one of the variates is restricted, they must all be restricted in the same way.
References
Burgess, T.M. & Webster, R. (1980). Optimal interpolation and isarithmic mapping of soil properties. I. The semi-variogram and punctual kriging. Journal of Soil Science, 31, 315-331.
Chiles, J.P. & Delfiner, P. (1999). Geostatistics: Modeling Spatial Uncertainty. Wiley, New York.
Journel, A.G. & Huijbregts, C.J. (1978). Mining Geostatistics. Academic Press, London.
Webster, R. & Oliver, M.A. (1990). Statistical Methods in Soil and Land Resource Survey. Oxford University Press, Oxford.
Webster, R. & Oliver, M.A. (2001). Geostatistics for Environmental Scientists. Wiley, Chichester.
See also
Directives: FVARIOGRAM, FCOVARIOGRAM, MCOVARIOGRAM, COKRIGE.
Procedures: MVARIOGRAM, DVARIOGRAM, DCOVARIOGRAM, DHSCATTERGRAM, KCROSSVALIDATION.
Commands for: Spatial statistics.
Example
" Example KRIG-1: Kriging
Form a variogram for levels of potassium at Brooms Barn Experimental Station
(see Webster & Oliver, 1990, Statistical Methods in Soil and
Land Resource Survey, Oxford University Press, pages 267-269)."
FILEREAD [NAME='%gendir%/examples/KRIG-1.DAT'] East,North,K
" Analyse on the log scale because of skewness of distribution"
CALCULATE LogK = LOG10(K)
" Form variograms in four directions, at 45 degree intervals, each
summarizing the semivariance across a 45-degree segment"
VARIATE [VALUES=0,45,90,135] Angles
& [VALUES=45,45,45,45] Segments
FVARIOGRAM [PRINT=statistics; Y=North; X=East; STEP=1; XMAX=13;\
DIRECTIONS=Angles; SEGMENTS=Segments]\
LogK; VARIOGRAM=LogKvar; COUNTS=Kcounts; DISTANCES=Midpoints
" Display the calculated variograms"
VARIATE Vgram[#Angles],Lag[#Angles],Count[#Angles]
CALCULATE Vgram[] = LogKvar$[*; 1...4]
& Lag[] = Midpoints$[*; 1...4]
& Count[] = Kcounts$[*; 1...4]
PRINT Lag[0],Vgram[0],Count[0],Lag[45],Vgram[45],Count[45]
& Lag[90],Vgram[90],Count[90],Lag[135],Vgram[135],Count[135]
AXES 1; YLOWER=0; XLOWER=0 :
PEN 1...4; COLOUR=1; SYMBOL=1...4
DGRAPH Vgram[]; Lag[]; PEN=1...4
" Model the variogram, trying three different models "
CALCULATE Kcounts=Kcounts*(Midpoints<11.75)
FOR Mod='LINEAR','SPHERICAL','EXPONENTIAL'
MVARIOGRAM [MODELTYPE=#Mod; PRINT=model,summary,estimates; \
WEIGHTING=counts; WINDOW=3; TITLE=Mod; XUPPER=15]\
LogKvar; COUNTS=Kcounts; DISTANCES=Midpoints
ENDFOR
" Produce matrices of predictions Kest and prediction variances Kvar
on a coarse grid, with interval 2 units (on scale of input coordinates) "
KRIGE [PRINT=d; X=East; Y=North; YOUTER=!(1,30); XOUTER=!(1,18);\
YINNER=!(1,30); XINNER=!(1,18); BLOCK=!(1.0,1.0); RADIUS=4.75;\
MINPOINTS=7; MAXPOINTS=20; INTERVAL=2]\
LogK; ISOTROPY=isotropic; MODELTYPE=spherical; NUGGET=0.0046;\
SILL=0.01528; RANGE=10.81; PREDICTIONS=Kest; VARIANCES=Kvar
PRINT Kest,Kvar; FIELD=7; DECIMALS=4
" Produce a finer grid with interval 0.5 units. This takes considerably
longer to calculate"
KRIGE [PRINT=d; X=East; Y=North; YOUTER=!(1,30); XOUTER=!(1,18);\
YINNER=!(1,30); XINNER=!(1,18); BLOCK=!(1.0,1.0); RADIUS=4.75;\
MINPOINTS=7; MAXPOINTS=20; INTERVAL=0.5]\
LogK; ISOTROPY=isotropic; MODELTYPE=spherical; NUGGET=0.0046;\
SILL=0.01528; RANGE=10.81; PREDICTIONS=Egrid; VARIANCES=Vgrid
" Reflect the calaculate grid so thath it is suitable for plotting
rather than printing as an array"
GETATTRIBUTE [ATTRIBUTE=rows,columns] Egrid; SAVE=Dim
CALCULATE Dim['rows'] = REVERSE(Dim['rows'])
& Nrow = NVALUES(Dim['rows'])
MATRIX [ROWS=Dim['rows']; COLUMNS=Dim['columns']] Ergrid,Vrgrid
CALCULATE (Ergrid,Vrgrid)$[Nrow...1;*] = (Egrid,Vgrid)$[1...Nrow;*]
" Produce a contour map of the predictions"
FRAME WINDOW=1,2; YLOWER=0; YUPPER=0.97,0.9; XLOWER=0,0.65; XUPPER=0.65,0.99
AXES 1; YLOWER=0.5; YUPPER=30.5; XLOWER=0.5; XUPPER=18.5
PEN [RESET=y] 2,3; CFILL=2,3
DCONTOUR [WINDOW=1; TITLE='Brooms Barn LogK'] Ergrid;INTERVAL=0.05; PENFILL=2
DSURFACE [WINDOW=1; TITLE='Brooms Barn LogK'] Ergrid;INTERVAL=0.05; PENFILL=2
" Map the variance of the predictions"
DCONTOUR [WINDOW=1; TITLE='LogK estimation variance'] Vrgrid;INTERVAL=0.0005