前言

最近项目上需要前端实现一个降雨量图，一开始一头蒙，这不是gis干的事吗？挠了一天头，后面经过几番曲折，找到了解决办法，就是kriging.js。先给大家看下降雨量图大致的样子，我找了深圳市气象局的一张图：

什么是克里金插值

克里金法（Kriging）是依据协方差函数对随机过程/随机场进行空间建模和预测（插值）的回归算法。在特定的随机过程，例如固有平稳过程中，克里金法能够给出最优线性无偏估计（Best Linear Unbiased Prediction, BLUP），因此在地统计学中也被称为空间最优无偏估计器（spatial BLUP）。

对克里金法的研究可以追溯至二十世纪60年代，其算法原型被称为普通克里金（Ordinary Kriging, OK），常见的改进算法包括泛克里金（Universal Kriging, UK）、协同克里金（Co-Kriging, CK）和析取克里金（Disjunctive Kriging, DK）；克里金法能够与其它模型组成混合算法。

若协方差函数的形式等价，且建模对象是平稳高斯过程，普通克里金的输出与高斯过程回归（Gaussian Process Regression, GPR）在正态似然下输出的均值和置信区间相同，有稳定的预测效果。克里金法是典型的地统计学算法，被应用于地理科学、环境科学、大气科学等领域。一大堆官方解释我也看的头晕，这些都是官方说法，我们只需要知道克里金算法在我们实现如降雨量图时起到了什么作用,我总结下来就是克里金算法在特定区域内对区域化变量进行无偏最优估计,能够顺滑衔接。

kriging.js源码主要模块分析

克里金插值算法有开源的js，GitHub地址kriging.js

kriging.train(t, x, y, model, sigma2, alpha)：使用gaussian、exponential或spherical模型对数据集进行训练，返回的是一个variogram对象；

kriging.grid(polygons,variogram,width)：使用刚才的variogram对象使polygons描述的地理位置内的格网元素具备不一样的预测值；

kriging.plot(canvas,grid,xlim,ylim,colors): 将得到的格网grid渲染至canvas上。

数据准备

1、站点数据

站点可以理解成监测点的位置或范围，多站点的数据形成一个特定范围,样例如下：

const sitesData = [
  { name: "广坤街道", lon: 113.7878622, lat: 22.68899985 },
  { name: "刘能街道", lon: 113.8817965, lat: 22.76980065 },
  { name: "老四街道", lon: 114.0154123999, lat: 22.7211864899 },
  { name: "晓峰街道", lon: 113.9847207999, lat: 22.5880693099 },
  { name: "大拿街道", lon: 113.961574, lat: 22.7522136099 },
  { name: "老徐街道", lon: 113.7898898, lat: 22.73634492 }
]

1、站点值

每个站点值和上面的站点一一对应,样例如下：

const siteValues = [52.4, 0.5, 5.5, 11.7, 0, 13]

2、定义色带

根据自己的具体需求定义每个范围的最大最小值，并设置一个颜色

const colors = [
  { min: 0, max: 0.1, color: "#008FFF" },
  { min: 0.1, max: 10, color: "#72D66B" },
  { min: 10, max: 20, color: "#3DB83D" },
  { min: 20, max: 30, color: "#3DB83D" },
  { min: 30, max: 40, color: "#3DB83D" },
  { min: 70, max: 300, color: "#3DB83D" }
]

3、修改源码，以获取自定义色值

//自定义色彩
kriging.getColor = function (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;
  }
};

4、添加站点到Cesium中

通过遍历站点数据添加Cesium的实体entity中，并用数组保存下来，并且调用绘制方法

import { drawKriging } from "./drawKriging";
let sitesEntity = []
sitesData.forEach((item) => {
    let entity = viewer.entities.add({
      id: item.id,
      position: Cesium.Cartesian3.fromDegrees(item.lon, item.lat, height),
      label: {},
      monitoItems: {
        ...item
      }
    });
    sitesEntity.push(entity);
});
drawKriging(sitesEntity, viewer, values, colors, "KrigingRain");

5、drawKriging方法

function getCanvas(siteValue, lngs, lats, ex, minx, maxy, maxx, miny, colors) {
  //1.用克里金训练一个variogram对象
  let variogram = kriging.train(siteValue, lngs, lats, "exponential", 0, 100);
  //2.使用刚才的variogram对象使polygons描述的地理位置内的格网元素具备不一样的预测值；
  let grid = kriging.grid(ex, variogram, (maxy - miny) / 500);
  let canvas = document.createElement("canvas");
  canvas.width = 1000;
  canvas.height = 1000;
  canvas.style.display = "block";
  canvas.getContext("2d").globalAlpha = 1; //设置透明度
  //3.将得到的格网预测值渲染到canvas画布上
  kriging.plot(canvas, grid, [minx, maxx], [miny, maxy], colors);
  return canvas;
}
export function drawKriging(sites, viewer, values, colors, layerName) {
  //如果存在雨量图则删除雨量图
  let KrigingRain = viewer.entities.getById(layerName);
  viewer.entities.remove(KrigingRain);
  let lngs = []; //经度集合
  let lats = []; //纬度集合
  let siteValue = []; //站点数值集合
  let coords = []; //绘制面的所有点
  let ex = []; //绘制面的jeojson
  ex = json.features[0].geometry.coordinates[0];
  for (let i = 0; i < sites.length; i++) {
    // sites[i].label.text = values[i]; // 展示每个站点名称
    let ellipsoid = viewer.scene.globe.ellipsoid;
    let cartesian3 = new Cesium.Cartesian3(
      sites[i]._position._value.x,
      sites[i]._position._value.y,
      sites[i]._position._value.z
    );
    let cartographic = ellipsoid.cartesianToCartographic(cartesian3);
    let lat = Cesium.Math.toDegrees(cartographic.latitude);
    let lng = Cesium.Math.toDegrees(cartographic.longitude);
    siteValue.push(values[i]);
    lngs.push(lng);
    lats.push(lat);
  }
  let coordinates = json.features[0].geometry.coordinates[0][0];
  for (let deg = 0; deg < coordinates.length; deg++) {
    coords.push(coordinates[deg][0]);
    coords.push(coordinates[deg][1]);
  }
  if (siteValue.length > 2) {
    new Cesium.PolygonGeometry({
      polygonHierarchy: new Cesium.PolygonHierarchy(Cesium.Cartesian3.fromDegreesArray(coords))
    }); //构造面，方便计算范围
    let extent = Cesium.PolygonGeometry.computeRectangle({
      polygonHierarchy: new Cesium.PolygonHierarchy(Cesium.Cartesian3.fromDegreesArray(coords))
    }); //范围（弧度）
    let minx = Cesium.Math.toDegrees(extent.west); //转换为经纬度
    let miny = Cesium.Math.toDegrees(extent.south);
    let maxx = Cesium.Math.toDegrees(extent.east);
    let maxy = Cesium.Math.toDegrees(extent.north);
    let canvas = null; //画布
    canvas = getCanvas(siteValue, lngs, lats, ex, minx, maxy, maxx, miny, colors);
    if (canvas !== null) {
      viewer.entities.add({
        id: layerName,
        polygon: {
          clampToGround: true,
          height: 0,
          hierarchy: {
            positions: Cesium.Cartesian3.fromDegreesArray(coords)
          },
          material: new Cesium.ImageMaterialProperty({
            image: canvas //使用贴图的方式将结果贴到面上
          })
        }
      });
    }
  }
}

实现效果

总结

通过克里金我们实现了降雨量图，我们还能实现气温图、湿度图、气压图等等，也了解了克里金函数的作用，以及和Cesium如何相结合使用。实际开发中，大家可以结合具体的业务需求开发自己的功能。

Cesium通过kriging.js实现降雨量图

前言