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方法的色彩赋值,更改后就可以根据数据范围赋予自定义的色彩了。