Calculates kriged estimates using a model fitted to the sample variogram.
|Controls printed output (
||Y positions (not needed for 2-dimensional regular data i.e. when
||X positions (needed only for 2-dimensional irregular data)|
||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|
||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|
||Variate containing 2 values to define the Y-bounds of the interpolated region (bottom then top); no default|
||Variate containing 2 values to define the X-bounds of the interpolated region (left then right); no default|
||Dimensions (length and height) of block; default !(0, 0) i.e. punctual kriging|
||Maximum distance between target point in block and usable data|
||Type of search (
||Minimum number of data points from which to compute elements; default 7|
||Maximum number of data points from which to compute elements (2 <
||Number of steps for numerical integration; (3 <
||Amount of drift (
||Ratio of Y interval to X interval; default 1.0|
||Distance between successive interpolations; default 1.0|
||Observed measurements as a variate or, for data on a regular grid, as a matrix|
||Form of variogram (
||Model fitted to the variogram (
||The nugget variance|
||Sill variances of the spatially dependent component; default none|
||Ranges of the spatially dependent component; default none|
||Slope of the unbounded component; default none|
||Power of the unbounded component or power for the
||Value of ν parameter for the Matern model; defalt none|
||Phi parameters of an anistropic model (
||Maximum gradient or distance parameter of an anistropic model|
||Minimum gradient or distance parameter of an anistropic model|
||Saves the Lagrange multipliers from each kriging solution|
||Specifies the variance of the measurement error|
||Supplies the model name and estimates, as saved from
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
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
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
YINNER are set to variates similarly to
YOUTER, and their limits should not lie outside those of
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
MAXPOINTS options. There is a minimum limit of 3 for
MINPOINTS and a maximum of 40 for
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
anisotropic; by default
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
boundedlinear (one dimension only),
besselk1 (Whittle’s function),
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
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
doublespherical, require the
SILLVARIANCES (the sill of the correlated variance) and
RANGES parameters. The latter is strictly the correlation range of the
pentaspherical models, while for the asymptotic models it is the distance parameter of the model. The
doublespherical model requires
RANGES to be set to variates of length two, to correspond to the two components of the model.
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
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.
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.
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.
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.
Commands for: Spatial statistics.
" 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,Vgram,Count,Lag,Vgram,Count & Lag,Vgram,Count,Lag,Vgram,Count 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