Cesium结合kriging.js制作降雨等值面

前因

因工作需要使用cesium制作降雨等值面，所以在参考众多博客后。终于是成功实现了降雨等值面。
参考博客：
https://blog.csdn.net/weixin_44265800/article/details/103106321
https://www.jianshu.com/p/c8473ac0b08a
https://zhuanlan.zhihu.com/p/73769138
https://blog.csdn.net/jessezappy/article/details/108891162

实现效果图

使用克里金插值

1.kriging.js
克里金法有开源的克里金插值算法实现，可以将kriging.js直接导入项目进行使用。
克里金项目的GitHub地址：https://github.com/oeo4b/kriging.js
2.使用kriging.js
相关代码如下：

import kriging from "@/utils/kriging.js";
colors: [
        { min: 0, max: 5, color: "#A9F08E" }, 
        { min: 5, max: 10, color: "#72D66B" }, 
        { min: 10, max: 25, color: "#3DB83D" }, 
        { min: 25, max: 50, color: "#61B7FC" }, 
        { min: 50, max: 100, color: "#0001FE" }, 
        { min: 100, max: 250, color: "#FD00FA" }, 
        { min: 250, max: 1000, color: "#7F013E" }, 
],
map.drawKriging = function(viewer,values,colors){
    //如果存在雨量图则删除雨量图
    var KrigingRain = viewer.entities.getById('KrigingRain');
    viewer.entities.remove(KrigingRain);
    var lngs = [];//经度集合
    var lats = [];//纬度集合
    var siteValue = [];//站点数值集合
    var coords = [];//绘制面的所有点
    var ex = [];//绘制面的jeojson
    ex = [LYBJJSON.features[0].geometry.coordinates[0]];
    for(var i=0; i < sites.length; i++){
      if(sites[i].type == "1"){
        for(var j=0; j < values.length; j++){
          if(sites[i].id == values[j].code && sites[i].state == "normal"){
            sites[i].label.text = values[j].value+"";
            var ellipsoid=viewer.scene.globe.ellipsoid;
            var cartesian3=new Cesium.Cartesian3(sites[i]._position._value.x,sites[i]._position._value.y,sites[i]._position._value.z);
            var cartographic=ellipsoid.cartesianToCartographic(cartesian3);
            var lat=Cesium.Math.toDegrees(cartographic.latitude);
            var lng=Cesium.Math.toDegrees(cartographic.longitude);
            // var alt=cartographic.height;
            siteValue.push(values[j].value);
            lngs.push(lng);
            lats.push(lat);
          }
        }
      }
    }
    for(let item of LYBJJSON.features[0].geometry.coordinates[0]){
      coords.push(...item)
    }
    if (siteValue.length > 2) {
        const polygon = 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;//画布
        function getCanvas() {
            //1.用克里金训练一个variogram对象
            let variogram = kriging.train ( siteValue, lngs, lats, 'exponential', 0, 100 );
            //2.使用刚才的variogram对象使polygons描述的地理位置内的格网元素具备不一样的预测值；
            let grid = kriging.grid ( ex, variogram, (maxy - miny) / 500 );
            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 );
        }
        getCanvas ();
        if (canvas != null) {
          KrigingObj = viewer.entities.add ({
            id: "KrigingRain",
            polygon: {
              show: ShowName == "drawKriging"?true:false,
              clampToGround: true,
              hierarchy: {
                  positions: Cesium.Cartesian3.fromDegreesArray ( coords )
              },
              material: new Cesium.ImageMaterialProperty ({
                  image: canvas//使用贴图的方式将结果贴到面上
              })
            }
        });
      }
    }
};

3.计算使用数据

var lngs = [];//经度集合 结构：[102.12186, 101.97394, 102.0567]
var lats = [];//纬度集合 结构：[22.3089599, 22.34198, 22.1919399]
var siteValue = [];//站点数值集合 结构：[12, 15.5]
var coords = [];//绘制面的所有点 结构：[102.2158, 20.05916, 102.2175, 20.05916, 102.22083, 20.05583]
var ex = [];//绘制面的jeojson 结构：[[[102.2158, 20.05916], [102.2175, 20.05916], [102.22083, 20.05583]]]

4.更改kriging.js的plot方法

kriging.plot = function(canvas, grid, xlim, ylim, colors) {
    // Clear screen
    let ctx = canvas.getContext("2d");
    ctx.clearRect(0, 0, canvas.width, canvas.height);
    // Starting boundaries
    let range = [xlim[1]-xlim[0], ylim[1]-ylim[0], grid.zlim[1]-grid.zlim[0]];
    let i, j, x, y, z;
    let n = grid.length;
    let m = grid[0].length;
    let wx = Math.ceil(grid.width*canvas.width/(xlim[1]-xlim[0]));
    let 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 = kriging.getColor(colors,grid[i][j]);
            ctx.fillRect(Math.round(x-wx/2), Math.round(y-wy/2), wx, wy);
        }
    }
};
//自定义色彩
kriging.getColor = function (colors, z) {
    var l = colors.length;
    for (var 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;
    }
};

kriging.js使用方法解析

kriging.train ( siteValue, lngs, lats, ‘exponential’, 0, 100 );
使用gaussian、exponential或spherical模型对数据集进行训练，返回的是一个variogram对象。
0对应高斯过程的方差参数，一般设置为 0。
100对应方差函数的先验值，默认设置为100。

kriging.grid ( ex, variogram, (maxy - miny) / 500 );
使用train返回的variogram对象使ex描述的地理位置内的格网元素具备不一样的预测值。
(maxy - miny) / 500的500是返回的网格数量,越大越细处理越慢。

kriging.plot ( canvas, grid, [minx, maxx], [miny, maxy], colors );
将得到的格网grid渲染至canvas上。
这里我们更改了plot方法的色彩赋值，更改后就可以根据数据范围赋予自定义的色彩了。