---

一、克里金法

• 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构建气象图层的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。