KRIGING

Description

KRIGING generates a raster image by interpolating points specified in vector segments, using the Kriging method.

Parameters

Name	Type	Caption	Length	Value range
FILV	String	Input vector file name	0 - 192
DBVS *	Integer	Input vector segment(s)	1 - 48
FILE *	String	Output file name	1 - 192
DBOC *	Integer	Output raster channel(s)	1 - 2
DBOW	Integer	Raster output window	0 - 4	Xoffset, Yoffset, Xsize, Ysize
FLDNME	String	Field name for elevation	0 - 64	Default: ATTRIBUTE
BACKVAL	Float	Background gray-level value	0 - 1
VMODE	String	Kriging variogram mode	0 - 11	SPH | EXP | GAU Default: SPH
VPARM	Float	Variogram Nugget, Sill, Range	0 - 3
BLKMD	String	Kriging blocking method	0 - 5
BLKPTS	Integer	Number of points in block	0 - 1	Default: 5
BLKSZ	Float	Block size	0 - 2
MONITOR	String	Monitor mode	0 - 3	ON, OFF Default: ON

* Required parameter

Parameter descriptions

FILV

The name of the file containing the vector segments to rasterize. If no file is specified, it is assumed to be the same as FILE.

DBVS

The vector segments in which the attribute value or the z-coordinate will be used as the basis for interpolation.

If FILV is specified, it should contain this segment. If FILV is not specified, it is assumed to be FILE.

Ranges of channels or segments can be specified with negative values. For example, {1,-4,10} is internally expanded to {1,2,3,4,10}. When you are not specifying a range in this way, only 48 numbers can be specified explicitly.

FILE

The name of the GDB raster file to receive the interpolated image.

DBOC

The image channel to receive the interpolated result.

When you specify two values, the second channel is assigned the variance value produced for each pixel. For more information on variance, see the Algorithm section.

DBOW

The rectangular subset of the image (Xoffset, Yoffset, Xsize, Ysize) to output. If no value is specified for this parameter, the entire channel is output by default.

Xoffset, Yoffset define the upper-left starting pixel coordinates of the window. X-size is the number of pixels that define the window width. Y-size is the number of lines that define the window height.

FLDNME

The field name that contains the elevation values.

When the value of this parameter is specified as ZCOORD, the actual z-coordinates of the vectors are used. Field names are not case-sensitive, and they do not need to be specified in complete form. If more than one match exists, the first name is used.

The Field Name (FLDNME) parameter is required only for vector layers that contain 3-D points or contours. With 3-D lines, the z-coordinates are used automatically. With 2-D layers, the elevation values are not specified.

The same name applies to all contour and 3-D point segments specified in the input vectors (DBVS) parameter. The segments that do not satisfy this requirement can be converted to the required format with ZVALTRNS.

If a field name other than ZCOORD is specified but is not found, VDEMINT attempts to use the z-coordinates of each vertex or point; however, this may cause unexpected results.

BACKVAL

The background, or NoData value, for the output channel.

If no value is specified for this parameter, all pixels in the output window will have interpolated results computed and assigned. If a value is specified, only pixels with the defined value will be replaced by interpolated values. Any values which cannot be computed, due to lack of applicable data, will have a value of BACKVAL assigned.

VMODE

The type of variogram to apply.

For more information on the variogram modes, see the Algorithm section.

VPARM

The Nugget, Sill, and Range parameters of the variogram.

For more information on the variogram modes, see the Algorithm section.

BLKMD

The type of Kriging blocking method to be used.

BLKPTS

If the blocking method (BLKMD) is set to RTREE, this parameter specifies the maximum number of points in a block before splitting it.

If the blocking method is BLOCK, this parameter specifies the initial number of points in the block.

In both cases, the default value is 5.

BLKSZ

The x and y dimensions of the block size.

This value should be specified in the same units as the pixel size; for example, meters. If no value is specified for this parameter, an appropriate value is calculated based on the total number of data points, maximum number of points per block, and the image size.

MONITOR

The program progress can be monitored by printing the percentage of processing completed. A system parameter, MONITOR, controls this activity.

Available options are:

• ON: Turns the monitor ON (default)
• OFF: Turns the monitor OFF (recommended if running in batch or background mode)

Details

KRIGING generates a raster image by interpolating points specified in vector segments, using the Kriging method.

The Kriging method for estimating a surface over an area of interest is a statistical technique that minimizes the error variance for the estimated values. It calculates the value at a specific point as a weighted sum of the known points. The weights are determined by considering the covariance between each pair of the known points, and the covariance between the specific point and each of the known points.

KRIGING implements point Kriging with user control over a variety of input parameters. Points stored in vector segments (DBVS) with a scalar value at each point (FLDNME) have their value interpolated into an output raster (DBOC) over the indicated subwindow (DBOW).

Example

Use KRIGING to generate a digital elevation model (DEM) by interpolating a set of sample points.

EASI>filv	=	''	! read vector from FILE
EASI>dbvs	=	2	! vector segment with sample points
EASI>file	=	'demfile.pix'
EASI>dboc	=	1	! output channel
EASI>dbow	=		! process entire image
EASI>fldnme	=	'ZCOORD'	! elevation stored in z-coordinates
EASI>backval	=		! process all pixels
EASI>vmod	=	''	! default variogram mode (SPHERICAL)
EASI>vparm	=		! default Nugget, Sill, and Range 
EASI>blkmd	=	''	! default block mode (RTREE)
EASI>blkpts	=		! default points in block (5)
EASI>blksz	=		! default block size

EASI>RUN KRIGING
        

Algorithm

The Kriging technique for estimating a surface over an area of interest is a statistical technique that minimizes the error variance for the estimated values. It calculates the value at a specific point as a weighted sum of the known points. The weights are determined by considering the covariance between each pair of the known points, and the covariance between the specific point and each of the known points.

The covariance between any two points is determined by evaluating a variogram for the distance between the points. There are three variograms available through this function. Each of these requires three parameters: the nugget, the sill and the range. The preferred variogram and parameters can be determined by building an experimental variogram. This is achieved by considering all pairs of points.

For various distance increments, the variance–that is, the average of the squares of the differences in z–are computed and plotted. The nugget is the value at zero, which will usually be zero. As the distance gets larger, the variances are likely to level out. The distance at which they level out is the range. The value at which they level out is the sill.

The choice of variogram depends on the shape of the curve. The spherical and exponential variograms are both linear at the origin. The exponential variogram reaches the sill asymptotically; that is, it approaches the sill gradually without ever reaching it. The Gaussian variogram is parabolic at the origin and reaches the sill asymptotically.

The default variogram is spherical. This was chosen because it is linear at the origin, and linear interpolation is preferable when building elevation models to avoid a terraced effect in the surface. It was chosen in preference to the exponential variogram because the influence of closer points is significantly greater than further points and the spherical variogram reached the sill at the range whereas the exponential variogram does not.

The default nugget value is zero. At zero, the Kriging technique will produce the exact input z-value when evaluated at the input point. In general, when evaluating at each raster center, the precise input position is not evaluated; therefore, a sample input point may fall within a pixel that has a slightly different associated value, even if the nugget is zero. If the nugget is non-zero, the surface value of all data points will be closer to the mean. A non-zero value for the nugget should only be used when the confidence in the data is low.

The value used for distances greater than the range is the nugget plus the sill. Any distance less than the range will produce a value between the nugget and the sill plus the nugget. Note that if the sill is zero, the covariance for all distances will be the same. This will cause the matrix to be singular. If the nugget is zero, the value of the sill has virtually no influence over the result of the interpolation. If it becomes too small, round-off error can occur. A value of one is the default.

The range represents the distance greater than that at which the influence of one point on the other is negligible. This value has a marked effect on the result. For example, if the value is very small, the surface produced will all have the mean value except for small circles surrounding each of the data points.

If all data points being considered are further from the point of interest than the range, the point of interest will be assigned the average of the points. If there are data points on one side within range, but none on the other, the effect can result in extrapolated values beyond the range of the input data points, providing unexpected anomalies. The range should be at least twice the maximum critical distance, where the critical distance is the distance between two adjacent points.

For example, in scan-line data, it should be the maximum distance between scan lines. The range is data-dependent; hence, a default value for all data sets cannot be given. If the range is not specified, the software determines it from the average spacing. The difference between adjacent points in the line and adjacent points in adjacent lines skews the result. From experimentation with the data sets used in production, a value of eight times the average distance was found to be preferred. If this presents anomalies, the value should be increased.

When there are very few points, all points can be used in the calculation of each raster position. In most conditions, however, this is not practical. From the perspective of computation time, a matrix must be inverted that is one greater in each dimension than the number of points being considered.

This process is n-cubed and would be prohibitive for large numbers of points. In addition, distant points should have no influence over the value of the current point in many cases. Note that an input value more distant from the point being evaluated than the range, but still included in the calculations, will have an indirect influence over the value, as it will influence the mean of the data used to calculate the covariance.

This issue is managed by blocking the data into smaller subsets, and including in the calculation for each block some number of surrounding points. The classic means of defining the blocks is to set a reasonable block size and include points that are within half the block distance in either direction. This technique is adequate when the data is evenly distributed, but may cause problems when it is not.

The number of surrounding points to include is critical to avoid the blocking being visible in the output surface. The preceding technique, were the data points are evenly distributed, suggests adding three times the number of points in the block to ensure a smooth surface. Experimentation indicates that when less than twice the number of points is used, the blocking produces unusual anomalies in the resulting surface. Therefore, three times is recommended.

This algorithm provides two techniques. The first, the block technique, groups the data by the specified block size. It then looks for surrounding points to minimize a blocked result according to the search algorithm described later in this section. The issue in this case is the choice of block size. If the data distribution is widely divergent, a choice of a small block size can cause some blocks to have no points, while a large block size will cause some blocks to have many points.

The second technique, based on r-tree indexing, produces blocks with a maximum number of points in each block. Whenever a block becomes too large, it is split along its major axis into two, with the points distributed among them. Classically, the two sub-blocks are chosen such that each has the same number of points.

While this appears to be a reasonable choice, it can cause difficulties when points are brought in along scan lines. Because the majority of the points are clustered together, forcing one half into one block and the other half into a second block tends to produce long, narrow blocks. To minimize this, the block is divided along the major axis at its center.

When a block is being processed, a search is required to find the nearest points surrounding the block. The maximum number of points added is given by a predefined parameter, nMaxPointsPerMatrix. Within that limit, a minimum of nMinPointsPerOctant points are searched for in each octant. This avoids pulling in points from one direction only, thus producing a step effect near, for example, the center of two contour lines.

This process, however, could pull in points that were too far away. An example of this is when the data comes from flight lines almost parallel to the x-axis. A slightly lower point could have no data in the octant just south of east, or have the closest point in that octant very far away. To avoid this, the octant requirement is made "fuzzy" by only requiring that each pair of adjacent octants have twice the nMinPointsPerOctant points.

In addition, once the limit of nMaxPointsPerMatrix and the "fuzzy" minimum of the quadrant in any direction have been met, the search in any direction is limited to the current furthest distance from the block being processed.

The time of execution depends heavily on the number of points per block and the number of points added, as a square matrix one greater than the sum of these must be inverted for each block. This is an n-cubed process. The default value of five points per block and 15 points added runs quite efficiently while maintaining the parameter requirements. Doubling the number of points per block, and hence, the number of points added, more than doubles the execution time. Using 20 points per block, and hence 60 points added, takes over six times as long.