# Coyote Graphics Routines

## Date: Fri Mar 27 12:14:20 2015

single page | use frames     summary     class     fields     routine details     file attributes

# cgkrig2d.pro

Math, Interpolation, Gridding

The cgKrig2D function interpolates a regularly or irregularly sampled set of points of the form z = f(x, y) to produced a gridded 2D array using a statistical process known as kriging. Kriging is a method of optimal interpolation based on regression against known or observed z values of surrounding data points, weighted according to spatial covariance values by various types of kriging model functions. Each grid location is estimated from observed values at surrounding locations. It is often used with spatial data.

Like all interpolation schemes, kriging can produces spurious results in extreme cases, but has the advantage of being able to compensate for the effects of data clustering and other, similar problems better than other interpolation methods such as inverse distance squared, splines, and triangulation methods. This particular version of Krig2D is orders of magnitude faster than the version of Krig2D that was distributed with IDL through IDL 8.2.3.

An excellent explanation of the kriging process can be found here:

``` http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z00000076000000.htm http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Semivariograms_and_covariance_functions ```
An explanation of the innovation that caused Krig2D to be made faster by several orders of magnitude can be found here:
``` http://www.idlcoyote.com/code_tips/krigspeed.php ```
I've implemented the kriging mathematical models described in the following references:
``` http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z00000076000000.htm http://www.nbb.cornell.edu/neurobio/land/OldStudentProjects/cs490-94to95/clang/kriging.html ```

## Examples

To create a dataset of N random points and determine the surface formed from such points:

``` n = 500 ;# of scattered points seed = -121147L ;For consistency x = RANDOMU(seed, n) y = RANDOMU(seed, n) ; Create a dependent variable in the form a function of (x,y) data = 3 * EXP(-((9*x-2)^2 + (7-9*y)^2)/4) + \$ 3 * EXP(-((9*x+1)^2)/49 - (1-0.9*y)) + \$ 2 * EXP(-((9*x-7)^2 + (6-9*y)^2)/4) - \$ EXP(-(9*x-4)^2 - (2-9*y)^2) params = [0.5, 0.0] interpArray = cgKrig2D(data, x, y, EXPONENTIAL=params, XOUT=xout, YOUT=yout) cgSurf, interpArray, xout, yout, /Save cgPlots, x, y, data, PSYM=2, Color='red', /T3D ```

## Author information

Author

FANNING SOFTWARE CONSULTING:

``` David W. Fanning 1645 Sheely Drive Fort Collins, CO 80526 USA Phone: 970-221-0438 E-mail: david@idlcoyote.com Coyote's Guide to IDL Programming: http://www.idlcoyote.com ```

Copyright (c) 2013, Fanning Software Consulting, Inc.

History

Written, 15 Oct 2013, based on a fast varient of the Krig2D program in the IDL library.

## Routines

`result = cgKrig2d_Exponential(d, t)`

The exponential kriging semivariogram model function.

`result = cgKrig2d_Spherical(d, t)`

The spherical kriging semivariogram model function.

`result = cgKrig2d(z [, x] [, y] [, BOUNDS=array] [, /DOUBLE] [, EXPONENTIAL=array] [, GS=array] [, NX=integer] [, NY=integer] [, /REGULAR] [, SPHERICAL=array] [, XGRID=array], XOUT=XOUT [, XVALUES=array] [, YGRID=array], YOUT=YOUT [, YVALUES=array])`

This function interpolates a regularly or irregularly gridded set of points Z = F(X,Y) using kriging.

## top cgKrig2d_Exponential

`result = cgKrig2d_Exponential(d, t)`

The exponential kriging semivariogram model function. This model should be applied when spatial autocorrelation decreases exponentially with increasing distance.

### Return value

A two-dimensional array containing the covariance model.

### Parameters

d in required type=float

The distance matrix of every point to each other.

t in required type=float

A three-element vector containing, in this order, the values for the range, nugget, and covariance value for the sample population (also called the partial sill), or [A, C0, C]. The sill is properly described as sill = (c0 + c).

## top cgKrig2d_Spherical

`result = cgKrig2d_Spherical(d, t)`

The spherical kriging semivariogram model function. This model show a progressive decrease of spatial autocorrelation until some distance, beyond which the autocorrelation is zero. This is one of the most commonly used models.

### Return value

A two-dimensional array containing the covariance model.

### Parameters

d in required type=float

The distance matrix of every point to each other.

t in required type=float

A three-element vector containing, in this order, the values for the range, nugget, and covariance value for the sample population (also called the partial sill), or [A, C0, C]. The sill is properly described as sill = (c0 + c).

## top cgKrig2d

`result = cgKrig2d(z [, x] [, y] [, BOUNDS=array] [, /DOUBLE] [, EXPONENTIAL=array] [, GS=array] [, NX=integer] [, NY=integer] [, /REGULAR] [, SPHERICAL=array] [, XGRID=array], XOUT=XOUT [, XVALUES=array] [, YGRID=array], YOUT=YOUT [, YVALUES=array])`

This function interpolates a regularly or irregularly gridded set of points Z = F(X,Y) using kriging. One of the kriging models MUST be used in the call to cgKrig2D. These are `CIRCULAR`, `EXPONENTIAL`, `GAUSSIAN`, `LINEAR`, or `SPHERICAL`. The only models currently tested and in use are EXPONENTIAL and SPHERICAL.

### Return value

A two-dimensional array containing the interpolated surface, sampled at input locations.

### Parameters

z in required type=float

An array containing the values of the data points as a function of X and Y. If X and Y are provided, this vector should be the same length. If X and Y are not provided, this array must be a 2D array. In this case the output grid is determined by `XGRID` (or `XVALUES`) and `YGRID` (or `YVALUES`) keywords, and default values for `NX` and `NY` are determined by the 2D dimensions of the input data array. If X and Y are not provided, regular gridding is assumed. Otherwise, the input data is assumed to be irregularly gridded, unless the 'REGULAR` keyword is set.`

x in optional type=float

An array containing the x locations of the surface to be gridded. If use, the `Y` data parameter must also be used and all three positional parameters must be the same length.

y in optional type=float

An array containing the y locations of the surface to be gridded. If use, the `X` data parameter must also be used and all three positional parameters must be the same length.

### Keywords

BOUNDS in optional type=array

Set this keyword to a four-element array [xmin, ymin, xmax, ymax] containing the grid limits of the output grid. If not provided, the grid limits are set to the extent of the X and Y vectors.

DOUBLE in optional type=boolean default=0

Set this keyword to perform all calculations in double precision floating point math. Otherwise, the calculations are done in since precision floating point math.

EXPONENTIAL in optional type=array

Set this keyword to a two- or three-element vector containing the kriging model parameters [A, C0, C] for the kriging exponential model. The parameter A is the range. At distances beyond A, the semivariogram or covariance remains essentially constant. The parameter C0 is the nugget. Theoretically, a zero separation distance, the semivariogram model should be zero. But, sometimes the semivariogram model displays a "lag" where the model function intercepts that Y axis at a location other than zero. This is called the "nugget". The parameter C, if it is present, is the value at which autocorrelation ceases to exist. If it is not present, it is calculated as the sample variance. The value C0+C is called the sill, which is the variogram value for very large distances. One of the kriging model keywords, `EXPONENTIAL` or `SPHERICAL`, must be used in the call to cgKrig2D.

For exponential models, the semivariagram at distance d is given as: C(d) = C1 * Exp(-3 * (d/A) if d not equal 0. C(d) = C1 + C0 if d equal 0.

GS in optional type=array

A two-element array [xs, ys] giving the grid spacing of the output grid, where xs is the spacing in the horizontal spacing between grid points, and ys is the vertical spacing. The default is based on the extents of x and y. If the grid starts at x value xmin and ends at xmax, then the default horizontal spacing is (xmax - xmin)/(`NX`-1). The ys parameter is computed in the same way. The default grid size, if neither `NX` or `NY` are specified, is 51 by 51.

NX in optional type=integer default=51

The output grid size in the X direction. If not specified, it can be be inferred from the `GS` and `BOUNDS` keywords. If not specified, and required by the code, a value of 51 is used.

NY in optional type=integer default=51

The output grid size in the Y direction. If not specified, it can be be inferred from the `GS` and `BOUNDS` keywords. If not specified, and required by the code, a value of 51 is used.

REGULAR in optional type=boolean default=0

Set this keyword to indicate the `Data` parameter is a 2D array containing measurements on a regular grid. It is rare to set this keyword, as it is set automatically under many circumstances.

SPHERICAL in optional type=array

Set this keyword to a two- or three-element vector containing the exponential model parameters [A, C0, C] for the kriging spherical model. The parameter A is the range. At distances beyond A, the semivariogram or covariance remains essentially constant. The parameter C0 is the nugget. Theoretically, a zero separation distance, the semivariogram model should be zero. But, sometimes the semivariogram model displays a "lag" where the model function intercepts that Y axis at a location other than zero. This is called the "nugget". The parameter C, if it is present, is the value at which autocorrelation ceases to exist. If it is not present, it is calculated as the sample variance. The value C0+C is called the sill, which is the variogram value for very large distances. One of the kriging model keywords, `EXPONENTIAL` or `SPHERICAL`, must be used in the call to cgKrig2D.

For spherical models, the semivariagram at distance d is given as: C(d) = c0 + C * ( ( 1.5 * (d/a) ) - ( 0.5 * (d/a)^3) ) if d less than a. C(d) = C + C0 if d greater than a. C(d) = 0 if d equals 0.

XGRID in optional type=array

Set this keyword to a two-element array, [xstart, xspacing] to indicate where the output grid starts and what the horizontal spacing will be. Do not specify the `XVALUES` keyword if this keyword is used.

XOUT
XVALUES in optional type=array

Set this keyword to a vector of X location values corresponding to the equivalent Z values in the `Data` parameter. Do not use this keyword if using the `XGRID` keyword.

YGRID in optional type=array

Set this keyword to a two-element array, [ystart, yspacing] to indicate where the output grid starts and what the vertical spacing will be. Do not specify the `YVALUES` keyword if this keyword is used.

YOUT
YVALUES in optional type=array

Set this keyword to a vector of Y location values corresponding to the equivalent Z values in the `Data` parameter. Do not use this keyword if using the `YGRID` keyword.

## File attributes

 Modification date: Fri Mar 27 11:07:39 2015 Lines: 455 Docformat: rst rst