How to use Kriging

Kriging is an alternative to many other point interpolation techniques. Unlike straightforward methods, such as Nearest Point, Trend Surface, Moving Average or Moving Surface; Kriging is based on a statistical method. Kriging is the only interpolation method available in ILWIS that gives you an an interpolated map and output error map with the standard errors of the estimates. As Kriging is an advanced technique that heavily relies on statistical theory, the novice user of GIS may need to ask advice from an expert.

Before you are going to use the Kriging method you should have thought about things like:

• Do I really need the Kriging interpolation method?
  When estimates with their errors are required, you should use Kriging instead of another interpolation technique. Examples of situations where Kriging could be very helpful are the mining industry, environmental research where decisions could have major economical and juridical consequences (e.g. is the area under study polluted or not) and so on.
• Is Kriging the most appropriate interpolation method for my sample data?
  Before using an interpolation technique, first the assumptions of the method(s) should be considered carefully. As a user, you should choose the most appropriate method for your job. For example, when you want to do some modelling it is better to choose a straightforward interpolation method. When you are interested in the estimation errors however you should use Kriging. Several articles are available on the comparison of different interpolation techniques and Kriging, which may help the GIS user to decide on the method to use. Examples of such kind of studies can be found in Burrough 1986 and Burrough and McDonnell 1998.

Once you have decided that Kriging is the method you want to use, you should continue with the following steps.

Step 1: Examining the input data

1.1 Visual and statistical data inspection

In the Catalog, select the point map with the sample points. Press the right mouse button and select Properties from the context-sensitive menu. The Point Map - Properties dialog box is opened.

• Check the map's domain and the amount of points and press Cancel.

Double-click the point map in the Catalog. The Display Options - Point Map dialog box is presented. Accept the defaults by clicking the OK button. The point map is opened. Check the distribution of the sample points visually. Another way to investigate whether your points are randomly distributed, or appear clustered, regular, or paired, etc., is by doing a Pattern Analysis.

If you decide that there are enough sample points available and that the distribution of points is good enough to do a Kriging interpolation, measure the shortest and longest distance between two sample points in the point map.

You can measure the length of both the shortest and the longest point-pair vector with the Measure Distance button (the pair of compasses) from the toolbar of the map window. You need this distance information later on for defining the number of lags and the lag spacing in the Variogram Surface operation (step 1.2), the Spatial Correlation or the Cross Variogram operation (step 2.1).

Tip:

If you have sufficient sample points to do a proper Kriging interpolation, split up your data set in two parts. Use one part for interpolation and the other part for verification of the interpolated map.

Next, calculate the variance of the sample data set.

• When the value data that you want to interpolate is stored in a column of the map's attribute table: open the attribute table. When values are stored in the point map itself, open the point map as a table via the context-sensitive menu. From the Columns menu in the table window, choose the Statistics command. In the Column Statistics dialog box:
• select the Variance function, and
• select the input variable for which you want to calculate the variance.
• Press the OK button. You can use this calculated variance as an indication for the sill when modelling the variogram (step 3.1).

1.2 Anisotropy

When the variation of the variable under study is not the same in all directions, then anisotropy is present. In case of suspected anisotropy, calculate a Variogram Surface with the Variogram Surface operation.

• Select the point map in the Catalog. Press the right mouse button and choose Statistics, Variogram Surface from the context-sensitive menu.
• The Variogram Surface dialog box is opened. In case the point map is linked to an attribute table, choose the attribute column with the sample data. Enter the lag spacing, the number of lags and a name for the output map.
• Select the Show check box and press OK. The variogram surface map is calculated.

The output map can best be viewed in a map window using representation Pseudo while a histogram has been calculated. To view the coordinates and the position of the origin in the output raster map, you can add grid lines, where the grid distance equals the specified lag spacing. It is important to recognize the origin of the plot/output map.

• Semi-variogram values close to the origin of the output map are expected to be small (blue in representation Pseudo), as values of points at very short distances to each other are expected to be similar. When there is no anisotropy, semi-variogram values will gradually increase from the origin into all directions. You will thus find circle-like shapes from the origin outwards where the color gradually changes from blue at the origin to green and red further away from the origin.
• Your input data is supposed to be anisotropic when you find an ellipse-like shape of low semi-variogram values (blue in representation Pseudo) in a certain direction going through the origin. In this direction, the semi-variogram values do not increase much. However, in the perpendicular direction, you find a clear increase of semi-variogram values: from blue at the origin to green and red further away from the origin. If anisotropy is present you should use Anisotropic Kriging.

You can measure the direction of anisotropy with the Measure Distance button from the toolbar of the map window, e.g. by following a 'line' of blue pixels going through the origin of the plot. You need this angle later on in step 2.1 (Spatial Correlation; bidirectional method) and step 4.3 (Anisotropic Kriging).

Tips:

• The semi-variogram values in the output map can be compared to the overall variance of your input data (calculated in step 1.1)
• When no points are encountered in a certain directional distance class, the semi-variogram value of that cell/pixel in the output surface will be undefined.
• When you find many undefined semi-variogram surface values in between few rather large semi-variogram values, you should consider to increase the lag spacing. Mind that results will be more reliable when say more than 30 point pairs are found in individual directional distance classes.
• When you find very many undefined semi-variogram values mainly at the outer parts of the surface, you should consider to reduce the lag spacing.
• When using an input point map with very many points, calculations of a large surface may take quite long. It is advised, to start using the operation with rather few lags and/or a rather small lag spacing.

1.3 Trend

One can perform a Kriging operation (i.e. Universal Kriging) while taking into account a local drift or trend that is supposed to exist within the limiting distance defined around each pixel to be interpolated. Very often you know the trend already:

• when you carried out another interpolation technique like Moving surface or a Trend surface, before you decided that you wanted to use the Kriging interpolation method, or
• when you know the natural behavior of the variable (e.g. soil textures near a river levee are very often more extreme than in a backswamp area).

Tip:

If a global trend is present in the sample set, subtract the trend for the input data with a TabCalc statement. Perform Ordinary Kriging on the de-trended data set and use MapCalc to add both output maps again. This method is an alternative for Universal Kriging. However, a major disadvantage of this alternative is that the error map is incorrect.

1.4 Multiple variables

If you decide that the variable under study is sparsely sampled, find out if there is another variable that is better sampled and has many corresponding sample points (identical XY-coordinates).

Open the attribute table that is linked to the point map. Find out if there are two columns with value domains and corresponding XY- coordinates. If there is a second variable, calculate the variance of both variables individually and the correlation between the two columns.

When the two variables are highly correlated, you can use the better-sampled variable and the relationship between the two variables to help to interpolate the sparsely sampled one with CoKriging.

If the correlation between the two variables is low it is advised to use another interpolation technique or not to interpolate at all.

Tips:

• The correlation between the two variables should make sense. In other words, it should be based on physical relationships/laws (e.g. temperature and relative humidity, temperature and height)
• If the variable is poorly sampled and there is no other variable that can help to interpolate the sparsely sampled one, it is not wise to interpolate because the results do not make much sense. Taking the average of the sample values is as good as interpolating.
• If the variable is sparsely sampled, you may consider to go back into the field and take more measurements.

Step 2: Calculation of the experimental variograms

2.1 Spatial Correlation and Cross Variogram

Kriging assumes a certain degree of spatial correlation between the input point values. To investigate whether your point values are spatially correlated and until which distance from any point this correlation occurs, you can use the Spatial Correlation operation. When you want to use CoKriging, you can investigate spatial correlation with Cross Variogram.

• Select the point map with the sample points in the Catalog. Press the right mouse button and select Statistics, Spatial Correlation (or Cross Variogram when you want to do CoKriging) from the context-sensitive menu. The Spatial Correlation dialog box or the Cross Variogram dialog box is opened.

   

• In the Spatial Correlation dialog box, you can choose to use either the omnidirectional or the bidirectional method:
• The omnidirectional method simply determines all distances between all point pairs, regardless of any direction, i.e. in all directions. Thus, all point pairs that have a certain distance towards each other will be counted in a certain distance class. Then, Moran's I, Geary's c, and the semi-variance are calculated for all point pairs within each distance class.
• The bidirectional method first counts, just like the omnidirectional method, all pairs of points that have a certain distance to each other, and then calculates the Moran's I and Geary's c for these point pairs within each distance class. Furthermore, all point pairs are counted with a certain distance to each other and with a certain direction towards each other. For these point pairs, semi-variance will be calculated. Then, also, for the direction perpendicular to the specified direction, point pairs are counted and semi-variances are calculated.

You must use the bidirectional method when anisotropy is present in your data set and you have to do an Anisotropic Kriging interpolation. Use the same anisotropy angle as measured in step 1.2.

• Enter the lag spacing and a name for the output table. Select the Show check box and press OK.
• The Spatial Correlation table is calculated.

   

• In a Cross Variogram dialog box:
• Enter the correct attribute columns, the lag spacing and a name for the output table. Select the Show check box and press OK.
• The Cross Variogram table is calculated.

Then, display the output table of the Spatial Correlation operation or the output table of the Cross variogram operation in a table window. Inspect in the output table the following columns:

• when you used Spatial Correlation omnidirectional method,

columns Distance, AvgLag, NrPairs;

• when you used Spatial Correlation bidirectional method,

columns Distance, AvgLag1/AvgLag2, NrPairs1/NrPairs2.

• when you used Cross Variogram,

columns: Distance, AvgLag, NrPairsA/NrPairsB/NrPairsAB and verify the Cauchy-Schwarz condition: g²_AB � g_A * g_B , for all values of the distance h.

Tips:

• For columns Distance/AvgLags: usually not more than half of the total sample distance should be taken into account; the larger the distance between the point pairs, the less point pairs, and the less reliable the outcome;
• For column NrPairs: for reliable semi- or cross-variogram values, distance classes should at least contain 30 point pairs. When the amount of point pairs in each distance class becomes too much you should choose a smaller lag spacing to get more information.

2.2 Displaying experimental variograms

The next step is to create point graphs, i.e. experimental semi- and cross-variogram(s), from the columns of the Spatial Correlation or Cross Variogram output table.

• From the File menu in the table window, choose the Create Graph command.
• In the Create Graph dialog box:
• choose for the X-axis, the Distance column or the column with the average distance between point pairs within lags (AvgLag);
• choose for the Y-axis, the column with the experimental semi- or cross-variogram values.

When you used the omnidirectional Spatial Correlation method, you can draw one graph:

• Distance or AvgLag , against SemiVar.

The discrete experimental semi- or cross-variogram values will by default be displayed as points in a new graph window.

When you used the bidirectional Spatial Correlation method, you can draw two graphs:

• Distance or AvgLag1 against SemiVar1 (semi-variogram values in the specified direction), and
• Distance or AvgLag2, against SemiVar2 (semi-variogram values in the perpendicular direction).

When you used the Cross Variogram operation, you can draw three graphs:

• Distance or AvgLag, against SemiVarA (semi-variogram values of the sparsely sampled variable),
• Distance or AvgLag, against SemiVarB (semi-variogram values of the other variable), and
• Distance or AvgLag, against CrossVarAB (cross-variogram values).

To add a second, third etc. graph layer to the graph window, choose Add Graph from Columns from the Edit menu in the graph window, or click the Add Graph button in the toolbar of the graph window.

In the graph window, you may wish to adapt the boundaries of the X-axis from 0 to more or less half of the total distance between the samples. You may also wish to adapt the Y-axis from 0 to more or less the expected variance of your input sample values (s², calculated in step 1.1).

Step 3: Modelling variograms

3.1 Adding models

Then, you need to model the discrete values of your experimental semi- or cross-variogram by a continuous function, which will give an expected value for any desired distance h.

• From the Edit menu in the graph window, choose the Add Semi-variogram Model command.
• In the Add Graph Semi-variogram Model dialog box, you can choose a type of semi- or cross-variogram model, and you can fill out values for the sill, range and nugget (or when using the Power model values for the sill, slope and power). A line will then be drawn according to the model you selected and the values you selected for sill, range and nugget.

Remark:

When you used the bidirectional Spatial Correlation method keep in mind that the model in the specified direction should be the same as the one in the perpendicular direction. ILWIS supports only geometric anisotropy: only the ranges of the models can be different.

3.2 Editing models

You are advised to visually experiment with models and with the values for sill, range, and nugget to find the best line through your discrete experimental semi- or cross variogram values.

• To edit an existing model, double-click the semi-variogram model you wish to edit in the Graph Management pane.
• The Graph Options - Semivariogram Model dialog box will appear in which you can adapt the model and the values for the nugget, sill and range until you are satisfied with the result.

If you like, you can also add other layers, i.e. other semi-variaogram models, for instance by choosing Add Graph, Semi-variogram Model from the Edit menu in the graph window, or by clicking the arrow next to the Add Graph button and selecting Semi-variogram Model.

3.3 Checking models

You can use the Column SemiVariogram operation to find which model fits your experimental semi- or cross-variogram values best. This operation calculates semi- or cross-variogram values according to a user-specified model and stores calculated semi- or cross-variogram values in an output column.

• Open the output table of the Spatial Correlation or Cross Variogram operation and choose the SemiVariogram command from the Columns menu in the table window.
• In the Column SemiVariogram dialog box, select the same distance column and use the same model and values for the nugget, sill and range as you specified when creating the semi-variogram model (step 3.1). Name the output column for instance SemiCol.
• You can test which semi- or cross-variogram model fits best, by calculating the Goodness of Fit (R²).

where:

The numerator of the fraction gives the sum of the squared differences between the experimental semi- or cross-variogram values and the semi- or cross-variogram values calculated by a user-specified semi- or cross-variogram model.

The denominator of the fraction gives the sum of the squared differences between the experimental values and the average experimental semi- or cross-variogram value of all distance classes/intervals.

The Goodness of Fit indicator can be calculated using some TabCalc statements:

• The terms in the numerator are computed by first making a new column in the Spatial Correlation or Cross Variogram table with the command:

g_gam_sqr = SQ(SemiVar - SemiCol)

where:

	SemiVar	is the name of the column with the experimental semi- or cross-variogram values
	SemiCol	is the name of the column with the semi- or cross-variogram values calculated with the Column SemiVariogram operation

   

• The terms in the denominator first require the calculation of S (h_i)/N (i.e. the average experimental semi- or cross-variogram values) in a new column:

mean_g = AVG(SemiVar)

   

• This average must be used in a new column with squared terms:

g_avg_sqr = SQ(SemiVar - mean_g)

   

• The 'Goodness of Fit' indicator R² is now obtained by:

R2 = 1 - SUM(g_gam_sqr) / SUM(g_avg_sqr)

The maximum value of R² is 1, meaning an exact match of semi-or cross-variogram values calculated using a certain semi- or cross-variogram model and parameters, and experimental semi- or cross-variogram values. Once, you have decided which model fits your data best, you can continue with the Kriging operation.

Step 4: Kriging interpolation

Perform one of the 5 Kriging operations available in ILWIS:

4.1 Simple Kriging

Simple Kriging assumes that the randomized spatial function is stationary and that the mean is constant over the area of interest. This assumption is often unrealistic; the more robust Ordinary Kriging is thus frequently used instead. In Simple Kriging, all input points are used to calculate each output pixel value. For more information, see Kriging : functionality.

4.2 Ordinary Kriging

In Ordinary Kriging the randomized spatial function is non-stationary and the mean varies over the area of interest. Ordinary Kriging amounts to re-estimating the mean at each new location. In Ordinary Kriging, you can influence the number of points that should be taken into account in the calculation of an output pixel value by specifying a limiting distance and a minimum and maximum number of points. Only the points that fall within the limiting distance to an output pixel will be used in the calculation for that output pixel value. For more information, see Kriging : functionality.

4.3 Anisotropic Kriging

Anisotropic Kriging is a variant of Ordinary Kriging. Anisotropic Kriging incorporates the influence of direction dependency , i.e. Anisotropy in the Kriging operation. For more information, see Anisotropic Kriging : functionality.

4.4 Universal Kriging

Universal Kriging is another variant of the Ordinary Kriging operation: Universal Kriging is Kriging with a local trend. The local trend or drift is a continuous and slowly varying trend surface on top of which the variation to be interpolated is superimposed. The assumption of a 'constant' mean in Ordinary Kriging is replaced by a prior trend model that is incorporated in the predictions. For more information, see Universal Kriging : functionality.

4.5 CoKriging

CoKriging is a multivariate extension of the Ordinary Kriging method. One parameter, which is too sparsely sampled to carry out an accurate interpolation, is highly correlated with another parameter of which enough samples are available to do Kriging. Based on the well sampled data set Kriging is carried out. For more information, see CoKriging : functionality.

Step 5: Output

5.1 Confidence interval maps

From the combination of a Kriged output map containing the estimates and its output error map, you can create confidence interval maps by using some MapCalc statements. For more information, see How to calculate confidence interval maps.

References:

• Burrough, P.A. 1986. Principles of Geographical Information Systems for Land Resources Assessment. Oxford University Press. 194 pp.
• Burrough, P.A. and R.A. McDonnell 1998. Principles of Geographical Information Systems. Oxford University Press. 330 pp.