kriging克里金插值以及前端渲染jS代码部分解释
2024-04-23 15:19:26
.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的数组，且每一个元素都被赋值成undefinedreturn 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的矩阵Zvar 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 000006kriging_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; // Singularpivinv = 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//复制未投影的逆矩阵为Kvar 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++) { // Verticesif(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) spacevar xtarget, ytarget;var a = Array(2),b = Array(2);var lxlim = Array(2); // Local dimensionsvar lylim = Array(2); // Local dimensionsvar 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++) { // Verticesif(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 subspacea[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 boundariesvar 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;
}();
kriging克里金插值以及前端渲染jS代码部分解释相关推荐
- 克里金插值详细步骤_openlayers4 入门开发系列之前端动态渲染克里金插值 kriging 篇（附源码下载）...
前言 openlayers4 官网的 api 文档介绍地址 openlayers4 api,里面详细的介绍 openlayers4 各个类的介绍,还有就是在线例子:openlayers4 官网在线例子 ...
- 前端实现克里金插值分析(二)
作者:yangjunlin 在上一篇文章中我们已经使用了像素法实现克里金插值的方式,但是问题也就随之抛了出来.1.第一点,在反距离权重插值的时候,因为处理的数据量大会直接导致主线程卡,导致用户体验 ...
- cesium基于克里金插值实现温度数据渲染
Cesium 是一款面向三维地球和地图的,世界级的JavaScript开源产品.它提供了基于JavaScript语言的开发包,方便用户快速搭建一款零插件的虚拟地球Web应用,并在性能,精度,渲染质量以 ...
- 克里金插值（Kriging）在MATLAB中的实现【优化】
该部分是基于克里金插值(Kriging)在MATLAB中的实现(克里金工具箱),由于在运行过程中有部分问题,基于此做的一些理解+优化. 工具箱的下载见上面的链接,其提供了工具箱. clc clearl ...
- 克里金插值（Kriging）在MATLAB中的实现（克里金工具箱）
一,直接献上克里金插值MATLAB工具箱 链接:https://pan.baidu.com/s/1p0lt2G8KQ-els7Dyw5jPFg 提取码:wcss 下载后将该程序添加到MATLAB安装文 ...
- 基于主动学习和克里金插值的空气质量推测
基于主动学习和克里金插值的空气质量推测 常慧娟, 於志文, 於志勇, 安琦, 郭斌 西北工业大学计算机学院,陕西 西安 710072 福州大学数学与计算机科学学院,福建 福州 350108 摘要 ...
- gstat | 空间插值（三）——克里金插值之泛克里金和简单克里金
本篇接着上篇继续介绍克里金插值.首先加载相关工具包和上篇使用的示例数据: library(gstat) library(sf) library(tidyverse) library(readxl) l ...
- gstat | 空间插值（二）——克里金插值之普通克里金
说明:昨天的推文误把可吸入颗粒物当作PM2.5,实应该为PM10,这里修正后重发. 从本篇开始计划分三篇介绍克里金插值.与反距离权重插值不同,克里金插值是无偏估计,其中也涉及到模型估计.本篇先对普通克 ...
- ArcGIS中使用协同克里金插值（co-kriging interplotation ）对气象数据插值
ArcGIS中如何使用协同克里金插值(co-kriging interplotation )对气象数据插值 ANUSPLIN气象站点数据插值局限性 百度搜索ArcGIS 克里金插值 搭建梯子搜索Arc ...
最新文章
- 配置CentOS的网络环境
- 日常生活 -- 嵌入式再学习前言
- erlang精要(5)-列表推导式
- 【转】什么是ERP、SCM、CRM？
- SublimeText3 初探(工欲善其事,必先利其器)
- Starter Kit for ASP.NET 2.0 家族又添新丁！
- 計算機二級-java10
- python 生成testbench_(Testbench用法总结)1. Testbench中文本数据的存储读取操作对比
- html-菜鸟--书架仿饿了么首页—Html学习（1）
- 解决网页上内容不能复制的几种方法
- IDEA开发工具整合YAPI接口平台
- python中绘制柱形图、饼形图等
- EasyUI上传图片，前台预览，后台读取
- 批处理文件方式 加密windows系统目录
- 大数据、java、python、区块链、人工智能发展前景
- SPSS随机对照研究总结
- 小米电视系统服务器升级,小米电视怎么更新系统 升级步骤图文详解
- 区块链安全和传统安全有什么不同
- 各厂内推整理 | 第三期
- java excel 65535_java导出excel数据超过65535