Coyote Graphics Routines

Date: Fri Mar 27 12:14:20 2015

Overview Directory File Etc Categories Search Index Warnings Help Developer documentation

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

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.

Routine details

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
Produced by IDLdoc 3.5.0 on Fri Mar 27 12:15:01 2015 IDLdoc project information