% MAPINFO "Info" for mapping roitines in SaGA. % This is "INFO" file for the following % programs contained in SaGA toolbox: % OBJMAP, KRIGING, MINCURVI % It can be called from MATLAB by typing: % >> type mapinfo or % >> objmap info % (or kriging info or mincurvi info) % Kirill K. Pankratov, kirill@plume.mit.edu OBJMAP, KRIGING and MINCURVI programs stand for 2-D interpolation from irregular points by objective mapping, kriging, and minimum curvature methods respectively. These programs are quite similar in structure and syntax (especially OBJMAP and KRIGING). All these methods are based on the inversion of the "Green's function" matrix G where G_ij element depends on the distance r_ij between i-th and j-th points of a set. Yet the Green's function itself is different in these functions: OBJMAP uses a 2-point correlation function KRIGING - 2-point semi-variogram: 1/2*<(z1-z2)^2> MINCURVI is most closely related to GRIDDATA and uses r^2*(log(r)-1) as a Green's function. ADAPTIVE QUADTREE DIVISION: All these routines use QUADTREE adaptive "tree" division for memory efficiency. The domain with limits [min(x) max(x)] and [min(y) max(y)] encompassing data is divided succesively into subdomains (blocks) until no subdomain contains more than a specified number of points N. Then the interpolation procedure is performed as follows: For each subdomain we find basis and interpolation points inside it and its immediate neighbours (sharing an edge or at least a vertex). We use the known (basis) points in these regions for gridding to interpolation points within these regions also. The resulting interpolation is taken with weight equal to 1 for "center" subdomain and weights rapidly decreasing away from it for neighbouring subdomains. Therefore final estimate at each interpolation points is obtained as a combination of weighted estimates from several sets of surrounding known points. This ensures smooth interpolated values at the boundaries of regions. If a neighbouring region contains fewer than a specified "depopulation threshold" number of points, then its own "secondary" neighbours are found and basis points inside them are also used. This implies that interpolation points will most likely be "completely surrounded" by known points and at the same time very distant points are not used, thus easing memory and computational load. Because subdomains has variable sizes, each subdomain has a different number of neighbours (typically about 10) and different number of known points are used in the interpolation. The parameters of the Quadtree division can be specified as the last input argument into MINCURVI, OBJMAP or KRIGING. They must compose a vector of length from 1 to 3 with following integer parameters: OPT = [N ND NMAX], where N is max. number of points in a block (subdomain), ND - "depopulation threshold" below which new "secondary" neighbours are sought, NMAX - maximum number of points in a set below which the QUADTREE procedure is not used. Defaults are OPT=[32 8 500]. See TREEDEMO for demonstration of QUADTREE division. SYNTAX, OPTIONS AND PARAMETERS: All these 3 routines must have at least 5 arguments similar to the "standard" MATLAB interpolation routines such as GRIDDATA or INTERP2. For example: MINCURVI(X,Y,Z,XI,YI) where X,Y,Z (vectors of the same lengths) are "basis" (known) coordinates and values, XI, YI (vectors or matrices of the same size) - coordinates of interpolation points. If XI is a row vector and YI - column vector of (possibly) different sizes, they specify matrices of constant rows and columns respectively. The rest are optional parameters, which are different for these three routines. MINCURVI(...,N) also specifies upper threshold number for Quatree division - the maximum number of points allowed inside one block. MINCURVI(...,OPT) where OPT = [N ND NMAX] also allows to input depopulation threshhold and "whole set" threshold for Quadtree division (see above section about QUADTREE). For example OPT = [20 5] specifies max. block population 2 and minimal (depopulation threshold) 5. In this case the whole set threshold has default value (500). OBJMAP(...,[LX LY],E) specifies horizontal scales LX, LY and relative measurement error E for the classical gaussian correlation function in the form C(x,y) = E*D(x,y)+(1-E)*exp(-(x/LX)^2-(y/LY)^2), where D - Dirac delta function. OBJMAP(...,L,E) uses the same scale l for LX and LY. OBJMAP(...,[LX LY],E,OPT) also include OPTions: the complete syntax is OPT = [N ND NMAX PB VERBOSE] where the first 3 numbers specify parameters of Quadtree division (see above). PB - lengthscales as parts of the average block size (when they are not supplied directly, in this case this agument can be empty). Default is 1/3 of the block size. VERBOSE - 1 or 0 - whether to display some information (scales, number of blocks and progress in processing them) during computations. default value is 1 - "verbose" mode. OBJMAP(...,FUN,P,E,OPT) uses another approach to correlation function: it is specified as a function name or expression (string) FUN. This function or expression can depend on arguments x, y, r=sqrt(x^2+y^2) and also parameters (such as lengthscales) stored in the vector P (next input argument). For example: objmap(...,'myfun(x,y,z,p)',p) objmap(...,'1/(1+(x/p(1))^2+(y/p(2))^2)',p) E and OPT are the same arguments as was discussed above. !!! Note that the results can be sensitive to the choice of the correlation function. Generally speaking, it must be "possible" correlation function, the one which can be realized for some random field. More rigorously, matrix G must be positive-definite, otherwise the results can be meaningless. For example, the function such as '1/(1+(r/p)^2)' will likely behave reasonably well, while similar-looking '1/(1+(r/p)^4)' will give extremely bad results. One should be careful in choosing the functional form of a correlation function. OBJMAP can have from 1 to 4 output arguments. If number is 1 or 3 then the form ZI or [XI,YI,ZI] is assumed. If number is 2 or 4 then it also returns error map EM: [ZI,EM] = ... or [XI,YI,ZI,EM] = ... Error map is a matrix of the same size as ZI. It shows a relative error of interpolation at the same points as ZI. If the correlaion function is "feasible", then EM has values between 0 and 1, if not - its values can be very unrealistic - <0 or >>1. When input argument Z is empty Z==[], then only error map is calculated and returned. KRIGING is similar to the OBJMAP in syntax. The only difference is that it uses a semi-variogram as opposed to the correlation matrix. KRIGING(...,[LX LY],E,OPT) assumes gaussian form of a semi-variogram: V(x,y) = V0*(1-(1-E)*exp(-(x/LX)^2-(y/LY)^2)-E*D(x,y)), where D - Dirac delta function, V0 - variance of Z (). E is a value of a semi-variogram at point (0,0) and is equivalent to a relative error for objective mapping. OPT is again the same vector of options as in OBJMAP. KRIGING(...,FUN,P,E,OPT) uses functional form of semi- variogram (expression of function name), similar to the OBJMAP. INTERPOLATION OF VECTOR FIELDS. In case when several variables measured at the same points X, Y need to be interpolated this can be done in one call to the OBJMAP of KRIGING functions. The advantage of this is that only the most computationally intensive procedure of calculating the Green's function matrix is performed only once and the weights of the basis points are the same for each variable. This is reasonable only for the case when all variables have the same coorelation function (and correspondingly, semi-variogram), otherwise each variable should be calculated separately. In case of several variables (or vector fields) all variables must be combined in one matrix with number of columns or number of rows equal to the number of known points (length of X, Y). For example, if velocities U, V and temperature T are measured in the ocean at points X, Y, the call should be constructed as [ZI,EM] = objmap(X,Y,[U(:) V(:) T(:)],XI,YI,...) In this case the output matrix ZI has the form ZI = [UI VI TI], where each column has the length of prod(XI). Errormap in this case is the same for all variables.