// Extend the Array class
|
Array.prototype.max = function () {
|
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.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_default_params = {
|
// krigingModel: "exponential", //model还可选'exponential','gaussian','spherical'
|
// krigingSigma2: 0,
|
// krigingAlpha: 226,
|
// canvasAlpha: 0.9, //canvas图层透明度
|
// // colors: ["#e8f8f5", "#d1f2eb", "#a3e4d7", "#76d7c4", "#48c9b0", "#1abc9c"],
|
// colors: [
|
// "#f1c40f",
|
// "#f4d03f",
|
// "#f7dc6f",
|
// "#f9e79f",
|
// "#fcf3cf",
|
// "#fef9e7",
|
// "#eafaf1",
|
// "#d5f5e3",
|
// "#abebc6",
|
// "#82e0aa",
|
// "#58d68d",
|
// "#2ecc71",
|
// "#28b463",
|
// "#239b56",
|
// "#1d8348",
|
// "#186a3b",
|
// "#145a32",
|
// "#196f3d",
|
// "#1e8449",
|
// ],
|
// };
|
|
var kriging = (function () {
|
var kriging = {};
|
var createArrayWithValues = function (value, n) {
|
var array = [];
|
for (var i = 0; i < n; i++) {
|
array.push(value);
|
}
|
return array;
|
};
|
// Matrix algebra
|
kriging_matrix_diag = function (c, n) {
|
var Z = createArrayWithValues(0, n * n);
|
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;
|
};
|
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) {
|
var i,
|
j,
|
Z = Array(n * m);
|
for (i = 0; i < n; i++)
|
for (j = 0; j < m; j++) Z[i * m + j] = X[i * m + j] + Y[i * m + j];
|
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
|
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
|
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 = 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 = createArrayWithValues(0, lags);
|
var semi = createArrayWithValues(0, 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 = createArrayWithValues(1, 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
|
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
|
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
|
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 = 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
|
kriging.plot = function (canvas, grid, xlim, ylim, colors) {
|
// Clear screen
|
var ctx = canvas.getContext("2d");
|
ctx.clearRect(0, 0, canvas.width, canvas.height);
|
|
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;
|
})();
|
|
if (module && module.exports) {
|
module.exports = kriging;
|
}
|
