coyote: CGKRIG2D

Description
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 
Categories
Math, Interpolation, Gridding
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
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
History
Written, 15 Oct 2013, based on a fast varient of the Krig2D program in the IDL library.
Copyright
Copyright (c) 2013, Fanning Software Consulting, Inc.
exponential kriging semivariogram model function. This model should be applied when spatial
correlation decreases exponentially with increasing distance.
Returns
A two-dimensional array containing the interpolated surface, sampled at input locations.
Params
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.
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.
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.