R - "Kriging" (Ordinary Kriging)

Olivia爱数学

没有

R语言里面已经有Kriging的包，可以直接下载使用，library("kriging")。

关于kriging的基础部分，我会在后面一篇文章单独写，具体的推导一下ordinary kriging的来源。

图片来自于ArcGIS help，为了表述方便，统一变量：块金（Nugget）、变程(range)、基台(sill)、偏基台值(partial sill)。

在CRAN中下载kriging, CRAN - Package kriging 为下载地址。直接下载源码解压之后的目录显示为如下所示：

先了解一下文件夹里面的内容。其中man/文件里面包含的文件是类似于解释文档之类的，文件中的内容以.Rd为后缀，这个文件夹是用来编译的时候转为pdf或者http所示；在R/文件夹中，存储的是以.r为结尾R文件；src文件夹中存储的是其他语言的代码，通常为C或者C++，这个的意义所在是因为用C语言或者C++所运行的代码速度比R快很多，因此可以通过这个代码来获取更好的效率。

接下来来描述kriging包里面的具体函数来阐述一下kriging插值（ordinary kriging）

首先是描述两个函数，一个是 plot.kriging ，一个是 image.kriging 。plot.kriging是用来得到plot semivariogram。以plot图的形式来得到距离和值的差的关系。这个部分的函数仅仅是绘图的内容，具体的可以参考一下别的文档里面对于如何定义绘图坐标轴，最大值，最小值之类的。image.kriging是用来得到map，在经过插值之后。这两个函数存在于methods.R中。

除开这两个函数，最重要的是kriging函数，主函数如下所示：

kriging(x, y, response, model = 'spherical', lags = 10, pixels = 100, polygons = NULL)

x 是空间点的x坐标，以vector的格式（ 有书上说过这个R语言里面对数据的处理是以vector的格式处理 ）

y是空间点的y坐标，以vector的格式

response是空间点的属性值，以vector的格式存储，这个属性值通常与x,y的长度相同

model是选择克里金插值的variogram的内容，可以选择的方式包含："spherical", "exponential", "gaussian"，默认的方式是spherical

lags，（为了将插值区域以距离的方式分割成为块状）

pixel，是沿着任何一个轴的（x,y坐标）所允许的最大的点数，默认值为100

polygons是可以用于来插值的多边形

Kriging主代码块（R）

# Ordinary Kriging kriging <- function(x, y, response, model = "spherical", lags = 10, pixels = 100, polygons = NULL) { #先定义基础的参数 n <- length(response) #n是这个属性值的个数，不是多少种属性 p <- 2 D <- twodim(x, y, n) #twodim函数会在后面单独说到 V <- onedim(response, n) #onedim函数具体内容在kriging.c文件中 cutoff <- sqrt((max(x)-min(x))^2+(max(y)-min(y))^2)/3 #从后面来看，只是给个一个范围的阈值 #由krig.vg函数(variogram)得到一个矩阵W，输入的参数分别是D <- as.matrix(D,n,n)是讲D以矩阵的方式作为parameter输入，目的是计算方便，同样的是V，这一句是得到一个对应的W W <- krig.vg(as.matrix(D, n, n), as.matrix(V, n, n), cutoff, lags) # fit.vg是一组值(系数们),这个函数的目的是为了得到一堆数据，那这些数据是为了得到拟合的值 # 这个系数由矩阵W的第二列和第三列得到，其实就是拟合已知点和已知点之间的相关关系。 lm()是一个线性拟合的函数。 fit.vg <- lm(W[,3]~W[,2])$coefficients #nugget是块金值，可以看作是lm函数拟合之后的intercept nugget <- as.numeric(fit.vg[1]) #range是变程，通俗来说就是对于某一个点其附近对他存在有影响的点的最大的距离即为range，所以做插值的时候，会选择在这个range内的点去插值也是因为这个距离内的点对插值点有影响。这里的range取矩阵W的第二列(mean/distance)的最大值 range <- max(W[,2]) #sill是基台 = partial sill + nugget，类比于线性方程，可以理解下面的内容，就是fit.vg[2]给出的系数(beta)，这个系数乘上range就可以得到在当前方程下，对应的range下的partial sill的值。 sill <- nugget+as.numeric(fit.vg[2])*range a <- 1/3 #krig.fit函数见后面部分描述，A即为对应于不同的模型下的拟合，返回的应该是一组值 A <- krig.fit(D, nugget, range, sill, model, n) # G返回的是krig.grid的值，见后面描述，这个函数返回的是空值（是没有值看，看作是一个啥也没有的矩阵），所以在后面的G$x,G$y都是对空的矩阵取值，取出填充新的值进去。 G <- krig.grid(min(x), min(y), max(x), max(y), pixels) #取G的pixel导出为pixel，格式为vector，返回的是不同值 pixel <- unique(G$pixel) # 见后面的krig.polygon函数 # 这个部分是将之前得到的G这个预先设定好的krig网格（可以看作是没有值）以初始值的形式发送到krig.polygon函数中，G$x,G$y以新的vector形式输入，得到新的kriging格网（注意，这个部分得到的G格网还是一个没有预测值的格网）。 if(!is.null(polygons)) { G <- krig.polygons(G$x, G$y, polygons) # 得到新的值，里面有预测值的已经。 G.pred <- krig.pred(x, y, response, G$x, G$y, c(solve(matrix(A, n+1, n+1))), nugget, range, sill, model, n) # o以list的格式输出。 o <- list( model = model, nugget = nugget, range = range, sill = sill, pixel = pixel, # 将输出的地图以数据框的格式存入。 map = data.frame( x = G$x, y = G$y, pred = G.pred semivariogram = data.frame( distance = W[,2], semivariance = W[,3] class(o) <- "kriging" #定义krig.vg函数 krig.vg <- function(D, V, cutoff, lags) { W <- matrix(0, lags, 3) # 新建一个matrix，命名为W，这个matrix有lags个行，3个列，被0填充为原始矩阵 for(i in 1:lags) { # 这个矩阵的构造：lags个行，我们首先假设lags=10，因此矩阵为10行3列，i是以行计数，因此[i,1],[i,2],[i,3]分别表示矩阵第一列，第二列，第三列。这里的cutoff/lags可以看作常数，D<(i*cutoff/lags)) & (D>((i-1)*cutoff/lags))是为了找到满足这个距离条件的点。存于第一列。之所以用lags的目的，是将插值区域划分为lags个块。由此可以推断，W[i,1]是表示满足这个条件的点的个数，W[i,2]是对应的均值，W[i,3]是对应的半变异函数(V...可以参见kriging.c文件中的onedimdist的函数定义)。Variogram是希望根据距离分成一段一段的，并且计算每一段里面的内容，因此需要有lags。 W[i,1] <- length(D[(D<(i*cutoff/lags)) & (D>((i-1)*cutoff/lags))])/2 W[i,2] <- mean(D[(D<(i*cutoff/lags)) & (D>((i-1)*cutoff/lags))]) # Average W[i,3] <- sum(V[(D<(i*cutoff/lags)) & (D>((i-1)* cutoff / lags ))]) / ( 4 * W[i,1]) # Semivariance return(W)

twodim代码块（onedim类似）

代码块是存在于src/文件夹下的kriging.c里面。用来得到n维下的欧式距离矩阵

/* n-dimensional euclidian distance matrices */ void twodimdist(double *x, double *y, int *n, double *d) { # *x,*y是指的点的x,y坐标，以vector形式存储，*n是指的有多少个点， int i, j; for(i=0;i<*n;i++) { for(j=0;j<*n;j++) { d[j+(*n)*i] = pow(pow(x[i]-x[j], 2) + pow(y[i]-y[j], 2), 0.5); # d[j+(*n)*i]是一个标签，d为D(在R文件中，最开始定义的部分)的小写，类似的见onedimdist部分的v对应的部分。 } } } # 这段代码就是计算任意两点之间的欧氏距离，onedimension类似于twodimen

kriging.fit代码块

void krigfit(double *d, double *nugget, double *range, double *sill, int *model, int *n, double *a) { /* Models: 0: Spherical (default) 1: Exponential 2: Gaussian int i, j; a[(*n+1)*(*n+1) - 1] = 0; # 也是对于区域内任意两点之间进行拟合。这个spehrical/exponential/gaussian的模型方式见后面具体阐述。 if(*model==0) { for(i=0;i<(*n);i++) { for(j=0;j<(*n);j++) { a[j+(*n+1)*i] = sphermodel(d[j+(*n)*i], *nugget, *range, *sill); } } }.......

variogram模型方式（只描述了默认的spherical，其余的类似，只是具体的方程不一样）

/* Variogram model functions */ # 其他的是以不同的模型构建距离和值之间的关系 double sphermodel(double h, double nugget, double range, double sill) { if(h>=range) return(sill); # h是任意的距离，当距离大于range变程的时候，没有影响，因此返回对应的y值即为基台 return( (sill-nugget) * (1.5*(h/range) - 0.5*pow(h/range, 3)) + nugget );

krig.grid 函数代码块

# 以格网的方式，这个部分是给出一个空的格网 void kriggrid(double *blx, double *bly, double *pixel, int *xpixels, int *ypixels, double *gx, double *gy) { int i, j; int k = 0; # xpixels和ypixels分别是指的x和y轴对应的最大的点数 # *blx和*bly是min(x),min(y),*xpixel,*ypixels是x,y轴的最大值 for(i=0;i<*xpixels;i++) { for(j=0;j<*ypixels;j++) { # i,j即共同构成点，一个是点的x坐标，一个是y坐标，gx[k]，gy[k]是对点进行初始值的赋值，给定初始值。 gx[k] = *blx + i*(*pixel); gy[k] = *bly + j*(*pixel); } } }

krigpolygon()函数

# 给出多边形，对上一步的空的格网给定初始值。 void krigpolygons(int *n, double *X, double *Y, int *nlist, int *npoly, double *polygonsx, double *polygonsy, int *Gn, double *Gx, double *Gy) { int i, j, k; for(i=0, k=0; i<*nlist;i++) { for(j=0; j<*n;j++) { # 这里的pointinpoly函数就是给出0，1来判断是否继续后面的赋值的代码 if(pointinpoly(npoly[i+1]-npoly[i], polygonsx+npoly[i], polygonsy+npoly[i], X[j], Y[j])) { Gx[k] = X[j]; Gy[k] = Y[j]; (*Gn)++; } } }}

mmult函数代码块

/* Mathematical matrix functions */ /* n x p, p x m */ void mmult(double *x, double *y, int *n, int *p, int *m, double *z) { /* x: n x p y: p x m 0: n x m int i, j, k; for(i=0;i<*n;i++) { for(j=0;j<*m;j++) { z[i+j*(*n)] = 0; for(k=0;k<*p;k++) { z[i+j*(*n)] += x[i+k*(*n)] * y[k+j*(*p)]; } } }}

krig.pred函数，给定插值的值

void krigpred(double *x, double *y, double *response, double *gx, double *gy, double *inva, double *nugget, double *range, double *sill, int *model, int *n, int *ng, int *nplus, int *one, double *gpred) { int i, j, k; double xdist; double R[*nplus]; double pred[1]; double invaXR[*nplus]; if(*model==0) { for(i=0;i<(*ng);i++) { # *ng是length(Gx)，其次gx[]是要插值的grid的x位置，是之前得到的grid的内容。xdist是要插值的点和已知的点之间的这个距离。（通过这个以及通过已知的模型方程，可以反推得到这个要插值点的值） for(j=0;j<(*n);j++) { xdist = pow(pow(x[j]-gx[i], 2) + pow(y[j]-gy[i], 2), 0.5); R[j] = sphermodel(xdist, *nugget, *range, *sill); R[*n] = 1; /* Y' * (A^-1 * R) */ mmult(inva, R, nplus, nplus, one, invaXR); mmult(response, invaXR, one, n, one, pred); gpred[i] = pred[0]; } }

上面的部分就是CRAN的"kriging"包的主要代码（忽略了一个cbind.r文件）这个文件其实很重要，是为了将c语言的代码和r语言的代码绑定在一起，所以如果要看源码，需要通过这个文件来查看在变量在不同语言下的命名。

kriging示例运行代码块

Examples # Krige random data for a specified area using a list of polygons library(maps) usa <- map("usa", "main", plot = FALSE) p <- list(data.frame(usa$x, usa$y)) # 运行p之后会发现成为[[1]]，然后接着一堆数据，以数据框的形式存在 # Create some random data x <- runif(50, min(p[[1]][,1]), max(p[[1]][,1])) # 随机产生50个点，这个点在给定的USA的面积的范围内，p[[1]][,1]即为原先数据框的第一列，为x，min,max为这个x的值的最小值和最大值，目的就是使得新产生的点是出于USA内，同时不同于原先的给出的x值，接下来的y同x所示。 y <- runif(50, min(p[[1]][,2]), max(p[[1]][,2])) z <- rnorm(50)