leflet使用kriging.js构建气象图层

一、克里金法

kriging.js
kriging.js是一个 Javascript 库，通过普通克里金算法提供空间预测和映射功能。
克里金法是一种高斯过程，其中使用核回归将二维坐标映射到某个目标变量。该算法经过专门设计，可通过为变异函数参数分配先验参数来准确地对较小的数据集进行建模。

二、构图函数调用

canvas在地图实例标签同列，定义一个canvas标签

-POINT数据

world数据


import kriging from "./kirging.js";// 定义一个WEATHER_LAYERS图层
// 气象图层
// let WEATHER_LAYERS = L.layerGroup([], {//   zIndex: 9,
// }).addTo(L_MAP);/*** 气温图层* @param {*} canvas：canvas = window.document.querySelector("#canvasMap");* @param {*} POINT* @param {*} world：边界*/
export const ADD_WEATHER_LAYER = (canvas, POINT, world) => {//WEATHER_LAYERS图层清空WEATHER_LAYERS.clearLayers();//遍历world边界数据，生成scope边界线const positions = [];world[0].forEach(function (point) {positions.push([point[1], point[0]]);});const scope = L.polyline(positions, { color: "rgba(0,0,0,0)" });WEATHER_LAYERS.addLayer(scope);//根据scope边界线的范围，计算范围变量const xlim = [scope.getBounds()._southWest.lng,scope.getBounds()._northEast.lng,];const ylim = [scope.getBounds()._southWest.lat,scope.getBounds()._northEast.lat,];//进行克里金插值function loadkriging() {canvas.width = 1000;canvas.height = 1000;const n = POINT.length;const t = []; //数值const x = []; //经度const y = []; //纬度for (let i = 0; i < n; i++) {t.push(POINT[i].properties.air_temp);x.push(POINT[i].geometry.coordinates[0]);y.push(POINT[i].geometry.coordinates[1]);// WEATHER_LAYERS.addLayer(//   L.circle([y[i], x[i]], { color: "rgba(0,0,0,0)" })// );//定义文本内容const divText = L.divIcon({html: `${POINT[i].properties.air_temp}℃`,className: "weather_number",bgPos: [0, 40],});const text = L.marker([POINT[i].geometry.coordinates[1], POINT[i].geometry.coordinates[0]],{icon: divText,});WEATHER_LAYERS.addLayer(text);}//对数据集进行训练const variogram = kriging.train(t, x, y, "exponential", 0, 100);//使用variogram对象使polygons描述的地理位置内的格网元素具备不一样的预测值,最后一个参数，是插值格点精度大小const grid = kriging.grid(world, variogram, (ylim[1] - ylim[0]) / 150);const colors = ["#002f86","#1a5ca6","#2075d2","#3ca0f0","#abe7f6","#abe7f6","#d2fbff","#f3ffee","#d0ffd0","#bffe8b","#fcfe9c","#fff3c4","#fecfa5","#fa9589","#ff5500","#e60000",];//将得到的格网grid渲染至canvas上kriging.plot(canvas, grid, [xlim[0], xlim[1]], [ylim[0], ylim[1]], colors);}//将canvas对象转换成image的URLfunction returnImage() {return canvas.toDataURL("image/png");}loadkriging();const imageBounds = [[ylim[0], xlim[0]],[ylim[1], xlim[1]],];const weather_layer = L.imageOverlay(returnImage(), imageBounds, {opacity: 0.8,});WEATHER_LAYERS.addLayer(weather_layer);
};

三、实例展示

四、kirging.js

// 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 = (function () {var kriging = {};var createArrayWithValues = function (value, n) {var array = [];for (var i = 0; i < n; i++) {array.push(value);}return array;};// Matrix algebraconst kriging_matrix_diag = function (c, n) {var Z = createArrayWithValues(0, n * n);for (let i = 0; i < n; i++) Z[i * n + i] = c;return Z;};const 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;};const 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;};const 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 multiplicationconst 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 decompositionconst 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 decompositionconst 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 eliminationconst 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 modelsconst 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))));};const kriging_variogram_exponential = function (h, nugget, range, sill, A) {return (nugget +((sill - nugget) / range) * (1.0 - Math.exp(-(1.0 / A) * (h / range))));};const 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 priorskriging.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/semivariancevar 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 distancevar 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 transformationn = 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 squaresvar 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 parametersvariogram.nugget = W[0];variogram.sill = W[1] * variogram.range + variogram.nugget;variogram.n = x.length;// Gram matrix with priorn = 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 vectorvar 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 Kvar K = C.slice(0);var M = kriging_matrix_multiply(C, t, n, n, 1);variogram.K = K;variogram.M = M;return variogram;};// Model predictionkriging.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 pathskriging.grid = function (polygons, variogram, width) {var i,j,k,n = polygons.length;if (n == 0) return;// Boundaries of polygons spacevar 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 DOMkriging.plot = function (canvas, grid, xlim, ylim, colors) {// Clear screenvar 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;
})();if (module && module.exports) {module.exports = kriging;
}

leflet使用kriging.js构建气象图层相关推荐

control层alert弹出框乱码_【ArcGIS for JS】动态图层的属性查询（11）
在真实需求中,我们不仅仅是将shp在地图中显示那么简单,我们往往要查询该图层的属性信息,我们在前面代码的基础上添加上属性查询. 1.1方法1(通过click直接获取) 1.1.1代码实现给要素图层添 ...
Openlayers利用kriging.js实现纯前端插值
关于空间插值这一块,也是GIS一个重要的分析.前端可以通过turf.js后端可以采用contour,如果仅有少量数据需要进行插值分析并进行展示,那应当如何做?可以采用turf.js但是turf.js插 ...
OpenLayers6（5）：基于kriging.js实现插值渲染图
1 版本 OpenLayers:6.4.3 2 相关配置 import kriging from '@sakitam-gis/kriging'; 3 使用kriging.js 3.1 数据准备待准备 ...
meetup_如何使用标准库和Node.js构建Meetup Slack机器人
meetup by Janeth Ledezma 简妮丝·莱德兹玛(Janeth Ledezma) 如何使用标准库和Node.js构建Meetup Slack机器人 (How to build a M ...
构建node.js基础镜像_我如何使用Node.js构建工作抓取网络应用
构建node.js基础镜像 by Oyetoke Tobi Emmanuel 由Oyetoke Tobi Emmanuel 我如何使用Node.js构建工作抓取网络应用 (How I built a ...
python django vue_Django+Vue.js构建项目
本文主要讲述如何从0开始,用Django和Vue.js构建一个项目.文章提要:Django与vue.js整合开发原理从头新建一个Django项目新建一个前端页面,使用vue应用在Django中设 ...
初识区块链——用JS构建你自己的区块链
初识区块链--用JS构建你自己的区块链区块链太复杂,那我们就讲点简单的.用JS来构建你自己的区块链系统,寥寥几行代码就可以说明区块链的底层数据结构.POW挖矿思想和交易过程等.当然了,真实的场景远远 ...
NativeScript - JS 构建跨平台的原生 APP
使用 NativeScript,你可以用现有的 JavaScript 和 CSS 技术来编写 iOS.Android 和 Windows Phone 原生移动应用程序.由原生平台的呈现引擎呈现界面而不 ...
react apollo_Apollo GraphQL：如何使用React和Node Js构建全栈应用
react apollo Apollo Client is a complete state management library for JavaScript apps. It's a powerf ...

leflet使用kriging.js构建气象图层

一、克里金法

二、构图函数调用

三、实例展示

四、kirging.js

leflet使用kriging.js构建气象图层相关推荐

最新文章

热门文章