kriging克里金插值以及前端渲染jS代码部分解释
韩靖琪
.html中的方法调用
//训练使用高斯过程与贝叶斯先验
let variogram=kriging.train(positions.map(pos=>pos[2]),positions.map(pos=>pos[0]),positions.map(pos=>pos[1]),params.krigingModel,params.krigingSigma2,params.krigingAlpha);
//网格矩阵或轮廓路径
let grid=kriging.grid(polygons,variogram,(extent[2]-extent[0])/1000);
//在DOM上绘图.
//Canvas是HTML5提供的一个标签,我们可以在这个盒子区域绘画
kriging.plot(canvas,grid,[extent[0],extent[2]],[extent[1],extent[3]],colors);
kriging-original.js文件中的部分代码解释
// Extend the Array class
// Array.prototype.max重写数组原型链
//表示取得最大值
Array.prototype.max = function() {
//apply()方法接受的是一个参数数组
//返回一个最大值的数组
return Math.max.apply(null, this);
};
//这里表示取得最小值
Array.prototype.min = function() {
//返回一个最小值
return Math.min.apply(null, this);
};
//这里表示算平均数
Array.prototype.mean = function() {
var i, sum;
for(i = 0, sum = 0; i < this.length; i++)
sum += this[i];
return sum / this.length;
};
Array.prototype.rep = function(n) {
//返回一个长度为n的数组，且每一个元素都被赋值成undefined
return Array.apply(null, new Array(n))
.map(Number.prototype.valueOf, this[0]);
//Number.prototype.valueOf()方法返回数值对象的原始值
};
Array.prototype.pip = function(x, y) {
var i, j, c = false;
for(i = 0, j = this.length - 1; i < this.length; j = i++) {
if(((this[i][1] > y) != (this[j][1] > y)) &&
(x < (this[j][0] - this[i][0]) * (y - this[i][1]) / (this[j][1] - this[i][1]) + this[i][0])) {
c = !c;
}
}
return c;
}
var kriging = function() {
var kriging = {};
// Matrix algebra矩阵代数
kriging_matrix_diag = function(c, n) {
//新建一个n*n的矩阵
var i, Z = [0].rep(n * n);
//循环赋值c给Z矩阵的每一元素
for(i = 0; i < n; i++) Z[i * n + i] = c;
return Z;
};
//将这个矩阵变为转置阵，也就是将元素颠倒顺序
kriging_matrix_transpose = function(X, n, m) {
var i, j, Z = Array(m * n);
for(i = 0; i < n; i++)
for(j = 0; j < m; j++)
Z[j * n + i] = X[i * m + j];
return Z;
};
//再次改变数值，把c给每一个二维元素赋值
kriging_matrix_scale = function(X, c, n, m) {
var i, j;
for(i = 0; i < n; i++)
for(j = 0; j < m; j++)
X[i * m + j] *= c;
};
//添加的方法
kriging_matrix_add = function(X, Y, n, m) {
//新建一个m*n的矩阵Z
var i, j, Z = Array(n * m);
for(i = 0; i < n; i++)
for(j = 0; j < m; j++)
//将X和Y矩阵相加合并成一个矩阵
Z[i * m + j] = X[i * m + j] + Y[i * m + j];
//返回一个Z矩阵
return Z;
};
// Naive matrix multiplication
//简单的矩阵乘法，矩阵和矩阵的乘法
//也就是前一个矩阵中的行乘以后一个矩阵中的列
kriging_matrix_multiply = function(X, Y, n, m, p) {
var i, j, k, Z = Array(n * p);
for(i = 0; i < n; i++) {
for(j = 0; j < p; j++) {
Z[i * p + j] = 0;
for(k = 0; k < m; k++)
Z[i * p + j] += X[i * m + k] * Y[k * p + j];
}
}
return Z;
};
// Cholesky decomposition
//柯列斯基分解，这是一种将正定矩阵分解为上三角矩阵和下三角矩阵的方法，
//在优化矩阵计算的时候会用到的一种技巧
//也就是，下面左边为下三角，右边为上三角
//100000 123456
//120000 023456
//123000 003456
//123400 000456
//123450 000056
//123456 000006
kriging_matrix_chol = function(X, n) {
var i, j, k, sum, p = Array(n);
for(i = 0; i < n; i++) p[i] = X[i * n + i];
for(i = 0; i < n; i++) {
for(j = 0; j < i; j++)
p[i] -= X[i * n + j] * X[i * n + j];
if(p[i] <= 0) return false;
p[i] = Math.sqrt(p[i]);
for(j = i + 1; j < n; j++) {
for(k = 0; k < i; k++)
X[j * n + i] -= X[j * n + k] * X[i * n + k];
X[j * n + i] /= p[i];
}
}
for(i = 0; i < n; i++) X[i * n + i] = p[i];
return true;
};
// Inversion of cholesky decomposition
//用斯基分解求矩阵的逆
kriging_matrix_chol2inv = function(X, n) {
var i, j, k, sum;
for(i = 0; i < n; i++) {
X[i * n + i] = 1 / X[i * n + i];
for(j = i + 1; j < n; j++) {
sum = 0;
for(k = i; k < j; k++)
sum -= X[j * n + k] * X[k * n + i];
X[j * n + i] = sum / X[j * n + j];
}
}
for(i = 0; i < n; i++)
for(j = i + 1; j < n; j++)
X[i * n + j] = 0;
for(i = 0; i < n; i++) {
X[i * n + i] *= X[i * n + i];
for(k = i + 1; k < n; k++)
X[i * n + i] += X[k * n + i] * X[k * n + i];
for(j = i + 1; j < n; j++)
for(k = j; k < n; k++)
X[i * n + j] += X[k * n + i] * X[k * n + j];
}
for(i = 0; i < n; i++)
for(j = 0; j < i; j++)
X[i * n + j] = X[j * n + i];
};
// Inversion via gauss-jordan elimination
//用高斯-约当消去法求逆，它的速度不是最快的，但是它非常稳定
//如果A是求解矩阵，那么求A的逆矩阵则为
//用A矩阵右边乘以单位矩阵I（与A同行同列值为1的单位矩阵）
//公式为A*I=I*B，(等号右边要同时变化)，也就是一个矩阵右边乘以单位矩阵化为，
//左边单位矩阵乘以B，则B就是A矩阵的逆
kriging_matrix_solve = function(X, n) {
var m = n;
var b = Array(n * n);
var indxc = Array(n);
var indxr = Array(n);
var ipiv = Array(n);
var i, icol, irow, j, k, l, ll;
var big, dum, pivinv, temp;
for(i = 0; i < n; i++)
for(j = 0; j < n; j++) {
if(i == j) b[i * n + j] = 1;
else b[i * n + j] = 0;
}
for(j = 0; j < n; j++) ipiv[j] = 0;
for(i = 0; i < n; i++) {
big = 0;
for(j = 0; j < n; j++) {
if(ipiv[j] != 1) {
for(k = 0; k < n; k++) {
if(ipiv[k] == 0) {
if(Math.abs(X[j * n + k]) >= big) {
big = Math.abs(X[j * n + k]);
irow = j;
icol = k;
}
}
}
}
}
++(ipiv[icol]);
if(irow != icol) {
for(l = 0; l < n; l++) {
temp = X[irow * n + l];
X[irow * n + l] = X[icol * n + l];
X[icol * n + l] = temp;
}
for(l = 0; l < m; l++) {
temp = b[irow * n + l];
b[irow * n + l] = b[icol * n + l];
b[icol * n + l] = temp;
}
}
indxr[i] = irow;
indxc[i] = icol;
if(X[icol * n + icol] == 0) return false; // Singular
pivinv = 1 / X[icol * n + icol];
X[icol * n + icol] = 1;
for(l = 0; l < n; l++) X[icol * n + l] *= pivinv;
for(l = 0; l < m; l++) b[icol * n + l] *= pivinv;
for(ll = 0; ll < n; ll++) {
if(ll != icol) {
dum = X[ll * n + icol];
X[ll * n + icol] = 0;
for(l = 0; l < n; l++) X[ll * n + l] -= X[icol * n + l] * dum;
for(l = 0; l < m; l++) b[ll * n + l] -= b[icol * n + l] * dum;
}
}
}
for(l = (n - 1); l >= 0; l--)
if(indxr[l] != indxc[l]) {
for(k = 0; k < n; k++) {
temp = X[k * n + indxr[l]];
X[k * n + indxr[l]] = X[k * n + indxc[l]];
X[k * n + indxc[l]] = temp;
}
}
return true;
}
// Variogram models
//变差函数模型
//变差函数高斯
kriging_variogram_gaussian = function(h, nugget, range, sill, A) {
return nugget + ((sill - nugget) / range) *
(1.0 - Math.exp(-(1.0 / A) * Math.pow(h / range, 2)));
};
//变差函数指数
kriging_variogram_exponential = function(h, nugget, range, sill, A) {
return nugget + ((sill - nugget) / range) *
(1.0 - Math.exp(-(1.0 / A) * (h / range)));
};
//变差函数的球形
kriging_variogram_spherical = function(h, nugget, range, sill, A) {
if(h > range) return nugget + (sill - nugget) / range;
return nugget + ((sill - nugget) / range) *
(1.5 * (h / range) - 0.5 * Math.pow(h / range, 3));
};
// Train using gaussian processes with bayesian priors
//训练使用高斯过程与贝叶斯先验
//kriging.train(t, x, y, model, sigma2, alpha):
//使用gaussian、exponential或spherical模型对数据集进行训练，返回的是一个variogram对象；
kriging.train = function(t, x, y, model, sigma2, alpha) {
var variogram = {
t: t,
x: x,
y: y,
nugget: 0.0,
range: 0.0,
sill: 0.0,
A: 1 / 3,
n: 0
};
switch(model) {
case "gaussian":
variogram.model = kriging_variogram_gaussian;
break;
case "exponential":
variogram.model = kriging_variogram_exponential;
break;
case "spherical":
variogram.model = kriging_variogram_spherical;
break;
};
// Lag distance/semivariance
// 滞后距离/半方差
var i, j, k, l, n = t.length;
var distance = Array((n * n - n) / 2);
for(i = 0, k = 0; i < n; i++)
for(j = 0; j < i; j++, k++) {
distance[k] = Array(2);
distance[k][0] = Math.pow(
Math.pow(x[i] - x[j], 2) +
Math.pow(y[i] - y[j], 2), 0.5);
distance[k][1] = Math.abs(t[i] - t[j]);
}
distance.sort(function(a, b) { return a[0] - b[0]; });
variogram.range = distance[(n * n - n) / 2 - 1][0];
// Bin lag distance
//本滞后距离
var lags = ((n * n - n) / 2) > 30 ? 30 : (n * n - n) / 2;
var tolerance = variogram.range / lags;
var lag = [0].rep(lags);
var semi = [0].rep(lags);
if(lags < 30) {
for(l = 0; l < lags; l++) {
lag[l] = distance[l][0];
semi[l] = distance[l][1];
}
} else {
for(i = 0, j = 0, k = 0, l = 0; i < lags && j < ((n * n - n) / 2); i++, k = 0) {
while(distance[j][0] <= ((i + 1) * tolerance)) {
lag[l] += distance[j][0];
semi[l] += distance[j][1];
j++;
k++;
if(j >= ((n * n - n) / 2)) break;
}
if(k > 0) {
lag[l] /= k;
semi[l] /= k;
l++;
}
}
if(l < 2) return variogram; // Error: Not enough points错误:分数不够
}
// Feature transformation功能转换
n = l;
variogram.range = lag[n - 1] - lag[0];
var X = [1].rep(2 * n);
var Y = Array(n);
var A = variogram.A;
for(i = 0; i < n; i++) {
switch(model) {
case "gaussian":
X[i * 2 + 1] = 1.0 - Math.exp(-(1.0 / A) * Math.pow(lag[i] / variogram.range, 2));
break;
case "exponential":
X[i * 2 + 1] = 1.0 - Math.exp(-(1.0 / A) * lag[i] / variogram.range);
break;
case "spherical":
X[i * 2 + 1] = 1.5 * (lag[i] / variogram.range) -
0.5 * Math.pow(lag[i] / variogram.range, 3);
break;
};
Y[i] = semi[i];
}
// Least squares最小平方
var Xt = kriging_matrix_transpose(X, n, 2);
var Z = kriging_matrix_multiply(Xt, X, 2, n, 2);
Z = kriging_matrix_add(Z, kriging_matrix_diag(1 / alpha, 2), 2, 2);
var cloneZ = Z.slice(0);
if(kriging_matrix_chol(Z, 2))
kriging_matrix_chol2inv(Z, 2);
else {
kriging_matrix_solve(cloneZ, 2);
Z = cloneZ;
}
var W = kriging_matrix_multiply(kriging_matrix_multiply(Z, Xt, 2, 2, n), Y, 2, n, 1);
// Variogram parameters变差函数参数
variogram.nugget = W[0];
variogram.sill = W[1] * variogram.range + variogram.nugget;
variogram.n = x.length;
// Gram matrix with prior有先验Gram矩阵
n = x.length;
var K = Array(n * n);
for(i = 0; i < n; i++) {
for(j = 0; j < i; j++) {
K[i * n + j] = variogram.model(Math.pow(Math.pow(x[i] - x[j], 2) +
Math.pow(y[i] - y[j], 2), 0.5),
variogram.nugget,
variogram.range,
variogram.sill,
variogram.A);
K[j * n + i] = K[i * n + j];
}
K[i * n + i] = variogram.model(0, variogram.nugget,
variogram.range,
variogram.sill,
variogram.A);
}
// Inverse penalized Gram matrix projected to target vector
//反向，，克矩阵投影到目标向量
var C = kriging_matrix_add(K, kriging_matrix_diag(sigma2, n), n, n);
var cloneC = C.slice(0);
if(kriging_matrix_chol(C, n))
kriging_matrix_chol2inv(C, n);
else {
kriging_matrix_solve(cloneC, n);
C = cloneC;
}
// Copy unprojected inverted matrix as K
//复制未投影的逆矩阵为K
var K = C.slice(0);
var M = kriging_matrix_multiply(C, t, n, n, 1);
variogram.K = K;
variogram.M = M;
return variogram;
};
// Model prediction
//模型预测，预测新的坐标点p=(xnew,ynew)的新的值（如高度，温度等）
kriging.predict = function(x, y, variogram) {
var i, k = Array(variogram.n);
for(i = 0; i < variogram.n; i++)
k[i] = variogram.model(Math.pow(Math.pow(x - variogram.x[i], 2) +
Math.pow(y - variogram.y[i], 2), 0.5),
variogram.nugget, variogram.range,
variogram.sill, variogram.A);
return kriging_matrix_multiply(k, variogram.M, 1, variogram.n, 1)[0];
};
//模型方差
kriging.variance = function(x, y, variogram) {
var i, k = Array(variogram.n);
for(i = 0; i < variogram.n; i++)
k[i] = variogram.model(Math.pow(Math.pow(x - variogram.x[i], 2) +
Math.pow(y - variogram.y[i], 2), 0.5),
variogram.nugget, variogram.range,
variogram.sill, variogram.A);
return variogram.model(0, variogram.nugget, variogram.range,
variogram.sill, variogram.A) +
kriging_matrix_multiply(kriging_matrix_multiply(k, variogram.K,
1, variogram.n, variogram.n),
k, 1, variogram.n, 1)[0];
};
// Gridded matrices or contour paths
//网格矩阵或轮廓路径
//kriging.grid(polygons,variogram,width);
//使用刚才的variogram对象使polygons描述的地理位置内的格网元素具备不一样的预测值；
//使用一个边界区域按间距生成grid格网数据数组
//polygons为区域的坐标数组,可以为多个polygon，variogram为第一步train产生的结果，width为生成grid格网的间距
kriging.grid = function(polygons, variogram, width) {
var i, j, k, n = polygons.length;
if(n == 0) return;
// Boundaries of polygons space
//多边形空间的边界
var xlim = [polygons[0][0][0], polygons[0][0][0]];
var ylim = [polygons[0][0][1], polygons[0][0][1]];
for(i = 0; i < n; i++) // Polygons多边形
for(j = 0; j < polygons[i].length; j++) { // Vertices
if(polygons[i][j][0] < xlim[0])
xlim[0] = polygons[i][j][0];
if(polygons[i][j][0] > xlim[1])
xlim[1] = polygons[i][j][0];
if(polygons[i][j][1] < ylim[0])
ylim[0] = polygons[i][j][1];
if(polygons[i][j][1] > ylim[1])
ylim[1] = polygons[i][j][1];
}
// Alloc for O(n^2) space
var xtarget, ytarget;
var a = Array(2),
b = Array(2);
var lxlim = Array(2); // Local dimensions
var lylim = Array(2); // Local dimensions
var x = Math.ceil((xlim[1] - xlim[0]) / width);
var y = Math.ceil((ylim[1] - ylim[0]) / width);
var A = Array(x + 1);
for(i = 0; i <= x; i++) A[i] = Array(y + 1);
for(i = 0; i < n; i++) {
// Range for polygons[i]
lxlim[0] = polygons[i][0][0];
lxlim[1] = lxlim[0];
lylim[0] = polygons[i][0][1];
lylim[1] = lylim[0];
for(j = 1; j < polygons[i].length; j++) { // Vertices
if(polygons[i][j][0] < lxlim[0])
lxlim[0] = polygons[i][j][0];
if(polygons[i][j][0] > lxlim[1])
lxlim[1] = polygons[i][j][0];
if(polygons[i][j][1] < lylim[0])
lylim[0] = polygons[i][j][1];
if(polygons[i][j][1] > lylim[1])
lylim[1] = polygons[i][j][1];
}
// Loop through polygon subspace
a[0] = Math.floor(((lxlim[0] - ((lxlim[0] - xlim[0]) % width)) - xlim[0]) / width);
a[1] = Math.ceil(((lxlim[1] - ((lxlim[1] - xlim[1]) % width)) - xlim[0]) / width);
b[0] = Math.floor(((lylim[0] - ((lylim[0] - ylim[0]) % width)) - ylim[0]) / width);
b[1] = Math.ceil(((lylim[1] - ((lylim[1] - ylim[1]) % width)) - ylim[0]) / width);
for(j = a[0]; j <= a[1]; j++)
for(k = b[0]; k <= b[1]; k++) {
xtarget = xlim[0] + j * width;
ytarget = ylim[0] + k * width;
if(polygons[i].pip(xtarget, ytarget))
A[j][k] = kriging.predict(xtarget,
ytarget,
variogram);
}
}
A.xlim = xlim;
A.ylim = ylim;
A.zlim = [variogram.t.min(), variogram.t.max()];
A.width = width;
return A;
};
kriging.contour = function(value, polygons, variogram) {
};
// Plotting on the DOM
//在DOM上绘图
//kriging.plot(canvas,grid,xlim,ylim,colors);将得到的格网grid渲染至canvas上
kriging.plot = function(canvas, grid, xlim, ylim, colors) {
// Clear screen
var ctx = canvas.getContext("2d");
ctx.clearRect(0, 0, canvas.width, canvas.height);
// Starting boundaries
var range = [xlim[1] - xlim[0], ylim[1] - ylim[0], grid.zlim[1] - grid.zlim[0]];
var i, j, x, y, z;
var n = grid.length;
var m = grid[0].length;
var wx = Math.ceil(grid.width * canvas.width / (xlim[1] - xlim[0]));
var wy = Math.ceil(grid.width * canvas.height / (ylim[1] - ylim[0]));
for(i = 0; i < n; i++)
for(j = 0; j < m; j++) {
if(grid[i][j] == undefined) continue;
x = canvas.width * (i * grid.width + grid.xlim[0] - xlim[0]) / range[0];
y = canvas.height * (1 - (j * grid.width + grid.ylim[0] - ylim[0]) / range[1]);
z = (grid[i][j] - grid.zlim[0]) / range[2];
if(z < 0.0) z = 0.0;
if(z > 1.0) z = 1.0;
ctx.fillStyle = colors[Math.floor((colors.length - 1) * z)];
ctx.fillRect(Math.round(x - wx / 2), Math.round(y - wy / 2), wx, wy);
}
};
return kriging;
}();
类似资料：