克里金插值算法文件

news/2024/9/29 15:40:52/

由于作者时间发布的比较久,导致一些cesium api和当前版本api不同,故改之:
文件名:kriging.js

/*!* author: [object Object]* @sakitam-gis/kriging v0.1.0* build-time: 2019-7-6 20:41* LICENSE: MIT* (c) 2019-2019 https://github.com/sakitam-gis/kriging.js*/
function max(source) {return Math.max.apply(null, source);
}
function min(source) {return Math.min.apply(null, source);
}
function rep(source, n) {let array = [];for (let i = 0; i < n; i++) {array.push(source);}return array;
}
function pip(source, x, y) {let i = 0;let j = source.length - 1;let c = false;let length = source.length;for (; i < length; j = i++) {if (((source[i][1] > y) !== (source[j][1] > y))&& (x < (source[j][0] - source[i][0]) * (y - source[i][1]) / (source[j][1] - source[i][1]) + source[i][0])) {c = !c;}}return c;
}
function matrixDiag(c, n) {let i = 0;let Z = rep(0, n * n);for (; i < n; i++) {Z[i * n + i] = c;}return Z;
}
function matrixTranspose(X, n, m) {let i = 0;let j;let Z = Array(m * n);for (; i < n; i++) {j = 0;for (; j < m; j++) {Z[j * n + i] = X[i * m + j];}}return Z;
}
function matrixAdd(X, Y, n, m) {let i = 0;let j;let Z = Array(n * m);for (; i < n; i++) {j = 0;for (; j < m; j++) {Z[i * m + j] = X[i * m + j] + Y[i * m + j];}}return Z;
}
function matrixMultiply(X, Y, n, m, p) {let i = 0;let j;let k;let Z = Array(n * p);for (; i < n; i++) {j = 0;for (; j < p; j++) {Z[i * p + j] = 0;k = 0;for (; k < m; k++) {Z[i * p + j] += X[i * m + k] * Y[k * p + j];}}}return Z;
}
function matrixChol(X, n) {let i;let j;let k;let 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;
}
function matrixChol2inv(X, n) {let i;let j;let k;let 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];
}
function matrixSolve(X, n) {let m = n;let b = Array(n * n);let indxc = Array(n);let indxr = Array(n);let ipiv = Array(n);let i;let icol = 0;let irow = 0;let j;let k;let l;let ll;let big;let dum;let pivinv;let temp;for (i = 0; i < n; i++) {for (j = 0; j < n; j++) {if (i === j)b[i * n + j] = 1;elseb[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;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;
}
function variogramGaussian(h, nugget, range, sill, A) {return nugget + ((sill - nugget) / range)* (1.0 - Math.exp(-(1.0 / A) * Math.pow(h / range, 2)));
}
function variogramExponential(h, nugget, range, sill, A) {return nugget + ((sill - nugget) / range)* (1.0 - Math.exp(-(1.0 / A) * (h / range)));
}
function variogramSpherical(h, nugget, range, sill) {if (h > range)return nugget + (sill - nugget) / range;return nugget + ((sill - nugget) / range)* (1.5 * (h / range) - 0.5 * Math.pow(h / range, 3));
}function train(t, x, y, model, sigma2, alpha) {let variogram = {t: t,x: x,y: y,nugget: 0.0,range: 0.0,sill: 0.0,A: 1 / 3,n: 0,model: variogramExponential,K: [],M: [],};switch (model) {case 'gaussian':variogram.model = variogramGaussian;break;case 'exponential':variogram.model = variogramExponential;break;case 'spherical':variogram.model = variogramSpherical;break;default:variogram.model = variogramExponential;}let i;let j;let k;let l;let n = t.length;let 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];let lags = ((n * n - n) / 2) > 30 ? 30 : (n * n - n) / 2;let tolerance = variogram.range / lags;let lag = rep(0, lags);let semi = rep(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;}n = l;variogram.range = lag[n - 1] - lag[0];let X = rep(1, 2 * n);let Y = Array(n);let 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];}let Xt = matrixTranspose(X, n, 2);let Z = matrixMultiply(Xt, X, 2, n, 2);Z = matrixAdd(Z, matrixDiag(1 / alpha, 2), 2, 2);let cloneZ = Z.slice(0);if (matrixChol(Z, 2)) {matrixChol2inv(Z, 2);}else {matrixSolve(cloneZ, 2);Z = cloneZ;}let W = matrixMultiply(matrixMultiply(Z, Xt, 2, 2, n), Y, 2, n, 1);variogram.nugget = W[0];variogram.sill = W[1] * variogram.range + variogram.nugget;variogram.n = x.length;n = x.length;let 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);}let C = matrixAdd(K, matrixDiag(sigma2, n), n, n);let cloneC = C.slice(0);if (matrixChol(C, n)) {matrixChol2inv(C, n);}else {matrixSolve(cloneC, n);C = cloneC;}let K1 = C.slice(0);let M = matrixMultiply(C, t, n, n, 1);variogram.K = K1;variogram.M = M;return variogram;
}
function predict(x, y, variogram) {let i;let 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 matrixMultiply(k, variogram.M, 1, variogram.n, 1)[0];
}
function variance(x, y, variogram) {let i;let 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);}let val = matrixMultiply(matrixMultiply(k, variogram.K, 1, variogram.n, variogram.n), k, 1, variogram.n, 1)[0];return variogram.model(0, variogram.nugget, variogram.range, variogram.sill, variogram.A) + val;
}
function grid(polygons, variogram, width) {let i;let j;let k;let n = polygons.length;if (n === 0)return;let xlim = [polygons[0][0][0], polygons[0][0][0]];let ylim = [polygons[0][0][1], polygons[0][0][1]];for (i = 0; i < n; i++) {for (j = 0; j < polygons[i].length; j++) {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];}}let xtarget;let ytarget;let a = Array(2);let b = Array(2);let lxlim = Array(2);let lylim = Array(2);let x = Math.ceil((xlim[1] - xlim[0]) / width);let y = Math.ceil((ylim[1] - ylim[0]) / width);let A = Array(x + 1);for (i = 0; i <= x; i++)A[i] = Array(y + 1);for (i = 0; i < n; 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++) {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];}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 (pip(polygons[i], xtarget, ytarget)) {A[j][k] = predict(xtarget, ytarget, variogram);}}}}return {xlim: xlim,ylim: ylim,width: width,data: A,zlim: [min(variogram.t), max(variogram.t)],};
}
function plot(canvas, grid, xlim, ylim, colors) {let ctx = canvas.getContext('2d');let data = grid.data, zlim = grid.zlim, width = grid.width;if (ctx) {ctx.clearRect(0, 0, canvas.width, canvas.height);let range = [xlim[1] - xlim[0], ylim[1] - ylim[0], zlim[1] - zlim[0]];let i = void 0;let j = void 0;let x = void 0;let y = void 0;let z = void 0;let n = data.length;let m = data[0].length;let wx = Math.ceil(width * canvas.width / (xlim[1] - xlim[0]));let wy = Math.ceil(width * canvas.height / (ylim[1] - ylim[0]));for (i = 0; i < n; i++) {for (j = 0; j < m; j++) {if (data[i][j] === undefined)continue;x = canvas.width * (i * width + grid.xlim[0] - xlim[0]) / range[0];y = canvas.height * (1 - (j * width + grid.ylim[0] - ylim[0]) / range[1]);z = (data[i][j] - 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.fillStyle = getColor(colors,data[i][j]);ctx.fillRect(Math.round(x - wx / 2), Math.round(y - wy / 2), wx, wy);}}}
}//自定义色彩
function getColor(colors, z) {let l = colors.length;for (let i = 0; i < l; i++) {if (z >= colors[i].min && z < colors[i].max) {    return colors[i].color;}}if (z < 0) {return colors[0].color;}
}let index = {train: train,predict: predict,variance: variance,grid: grid,plot: plot,max: max,min: min,pip: pip,rep: rep,matrixDiag: matrixDiag,matrixTranspose: matrixTranspose,matrixAdd: matrixAdd,matrixMultiply: matrixMultiply,matrixChol: matrixChol,matrixChol2inv: matrixChol2inv,matrixSolve: matrixSolve,variogramGaussian: variogramGaussian,variogramExponential: variogramExponential,variogramSpherical: variogramSpherical,
};export default index;
export { grid, matrixAdd, matrixChol, matrixChol2inv, matrixDiag, matrixMultiply, matrixSolve, matrixTranspose, max, min, pip, plot, predict, rep, train, variance, variogramExponential, variogramGaussian, variogramSpherical };
//# sourceMappingURL=kriging.esm.js.map

http://www.ppmy.cn/news/1531983.html

相关文章

grafana加载缓慢解决方案

背景 目前随着数据和图表的逐渐增多&#xff0c;Grafana 页面加载速度明显变慢&#xff0c;严重影响了用户体验&#xff0c;几次都有骂娘的冲动.&#xff0c;因此我们需要对 Grafana 进行优化&#xff0c;以提升加载性能。 对于速度优化&#xff0c;我们可以从以下方面进行入…

标准化和归一化的定义、公式、作用、示例、区别

标准化&#xff08;Standardization&#xff09; 和 **归一化&#xff08;Normalization&#xff09;**是数据预处理中常用的两种技术&#xff0c;目的是调整数据的尺度&#xff0c;使得不同特征的数据可以在同一水平上进行比较或处理。这两种方法在形式和用途上有所不同&#…

海滨体育馆管理系统:SpringBoot实现技巧与案例

2系统关键技术 2.1JAVA技术 Java是一种非常常用的编程语言&#xff0c;在全球编程语言排行版上总是前三。在方兴未艾的计算机技术发展历程中&#xff0c;Java的身影无处不在&#xff0c;并且拥有旺盛的生命力。Java的跨平台能力十分强大&#xff0c;只需一次编译&#xff0c;任…

mac怎么设置ip地址映射

最近开发的项目分为了两种版本&#xff0c;一个自己用的&#xff0c;一个是卖出去的。 卖出的域名是和自己的不一样的&#xff0c;系统中有一些功能是只有卖出去的版本有的&#xff0c;但我们开发完之后还得测试&#xff0c;那就需要给自己的电脑配置一个IP地址映射了&#xf…

【Sentinel-2简介】

Sentinel-2简介 Sentinel-2是欧洲空间局&#xff08;European Space Agency, ESA&#xff09;全球环境和安全监视&#xff08;即哥白尼计划&#xff09;系列卫星的重要组成部分&#xff0c;由Sentinel-2A和Sentinel-2B两颗卫星组成。以下是关于Sentinel-2的详细介绍&#xff1…

留学生如何适应海外生活以及应对文化差异

对于即将出国学习和生活的留学生来说&#xff0c;文化差异和生活方式的变化常常是一个紧迫的问题。那么&#xff0c;如何应对这些文化差异&#xff0c;以及如何适应新的学习环境和社交生活呢&#xff1f;本文将分享一些具体可行的建议和方法&#xff0c;助您顺利跨越这道难关&a…

ENV | 5步安装 npm node(homebrew 简洁版)

1. 操作步骤 1.1 安装 homebrew /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"1.2 安装 node # 安装最新版 brew install node # 安装指定版本&#xff0c;如18 brew install node181.3 安装 nvm&#xff08…

Craft:年度 Mac 应用,卡片式笔记新星

今年的年度 Mac 应用大奖颁给了Craft&#xff0c;这是一款集笔记、文档和个人管理于一体的独特工具。Craft 最大的亮点在于其卡片式的交互设计&#xff0c;这种设计让信息组织变得更加直观且高效。 尽管它仅上线了一年时间&#xff0c;但已经展现出了不输于许多老牌笔记应用的…