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