Forms experimental variograms.
Options
PRINT = string token |
Controls printed output (statistics); default stat |
|---|---|
Y = variate |
Y positions (needed only for 2-dimensional irregular data) |
X = variate |
X positions or interval (not needed for 2-dimensional regular data i.e. when DATA is a matrix) |
YMAX = scalar |
Maximum lag in the y direction (2-dimensional regular data only) |
XMAX = scalar |
Maximum lag in the x direction |
STEPLENGTH = scalar or variate |
Length(s) of the steps in which lag is incremented |
METHOD = string token |
How to estimate the variogram (moments, cressiehawkins, dowd, genton); default mome |
DIRECTIONS = scalar or variate |
Directions (degrees) along which to form the variogram (relevant only for 2-dimensional irregular data) |
SEGMENTS = scalar or variate |
Angles subtended by the segments (degrees) over which averaging is to be done (relevant only for 2-dimensional irregular data) |
Parameters
DATA = variates or matrices |
Measurements as a variate or, for data on a regular grid, as a matrix |
|---|---|
VARIOGRAMS = variates or matrices |
Structure to store the sample variogram |
COUNTS = variates or matrices |
Numbers of comparisons involved in the calculation of each variogram |
DISTANCES = variates or matrices |
Mean lag distances at each step |
LAGPOINTS = pointer |
Saves lag classes, indexes to observations and directions to plot in an h-scattergram |
Description
The FVARIOGRAM directive forms an experimental variogram from a set of values of a variable, Z, distributed in one or two dimensions. By default the variogram is calculated by Matheron’s method of moments, as
γ(h) = (1 / (2 × m(h))) × ∑ i = 1 … m(h) { z(xi) – z(xi + h) }2 ,
where z(xi) and z(xi + h) are the values at positions xi + h, and m(h) is the number of paired comparisons contributing to the estimate. For data on a regular grid or transect h is an integer multiple of the sampling interval. For irregularly scattered data h is discretized so that for each nominal lag there is a range of distance equal to the increment and an angular range set by the user. The nominal lag is at the centre of both ranges. However, you can set the METHOD option to calculate robust estimates instead. The cressiehawkins setting uses the estimator of Cressie & Hawkins (1980), which essentially damps the effect of outliers from the secondary process. Dowd’s (1984) and Genton (1978) methods, which estimate the variogram for a dominant intrinsic process in the presence of outliers, can be requested by the dowd and genton settings respectively. For further details see Webster & Oliver (2007) pages 67-68 and 115-116.
The data are specified using the DATA parameter. If they are on a regular grid, they should be supplied 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. Alternatively, if they are irregularly scattered, then they should be supplied in a variate, and the X and Y options should be set to variates to supply their spatial coordinates.
The experimental variogram is controlled by five options. For irregular data the maximum distance to which the variogram is calculated is set by the XMAX option for all directions. For regular data XMAX defines the maximum lag distance in the X direction, and YMAX must also be given to limit the distance in the Y direction. The increments in distance are set by the STEPLENGTH option, where you can supply a scalar to define equally-spaced steps or a variate to specify the steps themselves. The variogram may be computed in one or more directions. These are given by the DIRECTIONS option in degrees counterclockwise from east in the usual convention. Each direction is at the centre of an angular range, which is defined by the SEGMENTS option. DIRECTIONS and SEGMENTS should be set to scalars if the variogram is to be calculated for only one direction, or to variates if there are to be several.
A variogram can be computed without regard to direction by setting DIRECTIONS to 0 and SEGMENTS to 180. This is advisable if variation seems to be isotropic, i.e. the same in all directions, or if there are too few data to compute (h) for two or more directions separately. The lag then becomes a scalar |h| = h in distance only. Experience suggests that some 300 data are needed to distinguish anisotropy.
By default some statistics are printed concerning the variogram, but these can be supressed by setting option PRINT=*. Other information can be saved using the various parameters, in variates if there is a single direction, or in matrices with one column for each direction if there are several: VARIOGRAMS stores the ordered set of semivariances; DISTANCES stores the mean lag distances at which the semivariances have been computed; and COUNTS stores the numbers of paired comparisons from which the semivariances have been computed.
The LAGPOINTS parameter allows you to save a pointer containing lag classes, indexes to observations and directions that can be used to plot an h-scattergram.
Options: PRINT, Y, X, YMAX, XMAX, STEPLENGTH, METHOD, DIRECTIONS, SEGMENTS.
Parameters: DATA, VARIOGRAMS, COUNTS, DISTANCES, LAGPOINTS.
Action with RESTRICT
You can restrict a DATA variate to form the variogram from only a subset of its units.
References
Cressie, N. & Hawkins, D.M. (1980). Robust estimation of the variogram. Journal of the International Association of Mathematical Geology, 12, 115-125.
Dowd, P.A. (1984). The variogram and kriging: robust and resistant estimators. In: Geostatistics for Natural Resources Characterization (ed. G. Verly, M. David, A.G. Journel & A. Marechal), 91-106. D. Reidel, Dordrecht.
Genton, M.G. (1998). Highly robust variogram estimation. Mathematical Geology, 30, 213-221.
Webster, R. & Oliver, M.A. (2007). Geostatistics for Environmental Scientists, 2nd Edition. Wiley, Chichester.
See also
Directives: KRIGE, 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