ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

Cesium+克里金插值:从离散站点到连续温度场的WebGIS可视化实战

Cesium+克里金插值:从离散站点到连续温度场的WebGIS可视化实战 简介这是一份面向 WebGIS 开发者与三维可视化学习者的 Cesium Kriging 克里金插值实现资料针对离散温度测点生成连续温度场并在 Cesium 三维地球上进行渲染。资料基于原生 Cesium 与克里金插值算法涵盖变异函数分析、空间插值到三维表面生成的核心流程并附有完整的源码、示例数据与样式配置。压缩包内共 585 个文件以 JavaScript 文件为主搭配 JSON 配置数据、CSS 样式及大量 PNG/JPG 图片资源覆盖核心逻辑、场景构建、界面美化与 Cesium 运行依赖等模块整体仅 7.27MB轻量易部署。内置的可运行示例页面配合清晰的目录结构便于对照源码理解插值参数设置与渲染细节也能快速替换为其他气象要素进行迁移实践。已有 165 人学习下载适合希望快速掌握温度场可视化实现思路并做功能扩展的中高级 GIS 开发者。 接到一个温度分布图需求的时候我一开始真没把它当回事。气象站点的观测数据就摆在那里几百个点散落在整个省域范围内需求方想要的是类似天气App里那种连续平滑的温度色斑图而不是一个一个点标注的温度值。我以为套个热力图插件就能交差结果放大到站点稀疏的区域画面全是虚影数据方看完直接摇头。后来转到cesium kriging这条路线把克里金插值和三维地球渲染串起来才真正把离散站点变成了能看、能分析、能切换时次的连续温度场。这篇文章就围绕“在Cesium里用克里金插值渲染温度数据”这一整条链路展开从插值原理、参数选择到网格转GeoJSON多边形、Entity/Primitive两种渲染方式再到动态更新和性能调优的实战坑位。适合正在做气象、海洋、环境、土壤等WebGIS可视化手里有离散站点数据、想在三维地球上呈现连续场的开发者。不管你是刚接触Cesium还是已经写过一些实体渲染这份流程都可以直接拿过去改。1. 离散站点到连续温度场为什么绕不开插值这一步1.1 需求到底长什么样温度这类气象要素的本质是一个在空间上连续分布的场。任意一个经纬度点理论上都有对应的气温值。但气象站的观测只能覆盖极少数位置全省几百个站已经是比较密的站网了放在Cesium那个全球尺度上看依然稀疏得像撒了一把芝麻。直接把这些点用Cesium.Entity渲染成圆点或者标签能表达的信息只有“这个站此刻多少度”无法回答“任意一点今天大概多少度”“哪片区域在降温”“等温线大致走向”这类问题。所以必须走空间插值用已知站点的温度值按照某种空间关系去推算未知位置的温度得到一个规则网格面再把这个网格面渲染到地球上。这一步不做后续一切可视化都无从谈起。1.2 三种常规插值方案怎么选提到空间插值WebGIS圈子里最先冒出来的一般是三个名字IDW反距离加权、样条插值、克里金插值。做之前我把它们放在一张表里对比了一遍插值方法核心原理优点缺点典型场景IDW反距离加权距离越近权重越大直接用距离倒数加权平均实现简单、计算快极易出现“牛眼效应”站点附近形成一圈圈异常闭合圈离群点影响很大快速预览、数据量极小时临时用样条插值构造通过所有采样点的最小曲率曲面表面光滑连续容易过拟合出现超过实际可能范围的极值区地形趋势面、无需严格误差控制克里金插值通过变异函数刻画空间自相关性再按统计最优原则加权预测结果更平稳、符合空间规律、还能给出估计误差参数选择有门槛变异函数模型不合适时同样会出怪图气象、环境、土壤、地下水等专业领域纯粹从统计角度讲克里金被称为空间插值的“最佳线性无偏估计”因为它不止看距离还考虑了数据在空间上的结构关系。比如一个站点和它东边1km的点与它和北边5km的点温度相关性是完全不同的IDW不管这些只管距离。克里金通过变异函数把这种方向性、尺度性的自相关捕捉出来再参与加权计算。但在真实项目里我选克里金还有一个更朴素的理由这个算法在浏览器里有个成熟的kriging.js实现一套流程走下来生成的色斑图边界平滑、颜色过渡自然视觉上最接近气象业务里常用的温度产品图。领导不看你用啥算法只看图顺不顺眼。1.3 克里金是怎么做到的克里金的核心动作可以拆成三步先根据所有站点对的温度差和空间距离拟合一个变异函数模型描述“距离越远温度差异越大”这个关系具体有多强然后对任意一个待预测点以它和周围站点的距离关系为约束求解一组最优权重让预测误差的方差最小最后用这组权重对站点值做加权求和得到该点的预测温度。整个过程在数学上保证了对每个未知点的估计都是无偏且方差最小的。实际使用中我们感受不到这么细但理解这两句话对后面调参很重要。因为变异函数模型选得不对或者参数给得不合适温度场会出现“站点周围一个色块、离开站点就乱跳”的怪现象这不是Cesium的锅是插值本身就没算对。2. 克里金插值的核心逻辑与参数选择2.1 训练变异函数kriging.js的入口前端做克里金插值目前最顺手的库是kriging.jsnpm上叫sakitam-gis/kriging保持了经典kriging.js的APItrain负责训练变异函数grid负责按范围生成网格predict做单点预测。如果只是做批量网格插值用前两个就够。站点数据一般长这样一个数组里带经纬度和温度值const stations [ { lng: 116.35, lat: 39.86, temp: 24.3 }, { lng: 116.52, lat: 39.91, temp: 25.1 }, { lng: 116.28, lat: 39.73, temp: 23.8 }, // ... 几百个站点 ];然后调用trainimport { train, grid } from sakitam-gis/kriging; const xs stations.map(s s.lng); const ys stations.map(s s.lat); const zs stations.map(s s.temp); const variogram train(zs, xs, ys, exponential, 0, 100);train的参数里exponential是变异函数模型后面两个数值分别代表测量误差sigma2和模型的形状参数alpha。sigma2一般给0就行表示认为站点观测是准的alpha给100是经典案例里的默认值它控制变异函数达到平稳时的整体尺度实际项目里可以先保持默认出图后再根据效果微调。2.2 三种变异函数模型的取舍kriging.js支持exponential、gaussian、spherical三种模型我三种都用同一批数据跑过差别非常直观模型空间连续性特征温度场效果我的使用建议exponential指数衰减短距离内快速变化表面偏“锐”站点间过渡明显能看出局部冷暖中心默认首选大多数气象数据够用gaussian高斯型光滑度高整体非常平滑但容易把冷/热中心抹掉一部分极值站点很密、想得到极致光滑图时试spherical球状模型带有明显的变程概念过渡介于两者之间边缘可能出现轻微波动专业地统计人员按区域特征选择温度场的可视化个人经验是exponential最稳妥。gaussian虽然好看但对站点稀疏区域的预测会过度平滑把山区和平原的温差糊在一起数据方一旦对比真实站点值就会发现偏差偏大。2.3 网格生成的范围与步长train跑完拿到variogram接下来定义插值范围和分辨率。第一步是确定边界别用站点的最小最大经纬度直接圈要外扩一点不然边缘区域会因缺少邻域站点而出现明显的插值跳变。我习惯外扩0.1到0.2度视站点分布密度而定。const bounds { minLng: Math.min(...xs) - 0.2, maxLng: Math.max(...xs) 0.2, minLat: Math.min(...ys) - 0.2, maxLat: Math.max(...ys) 0.2 }; // 网格步长单位是度 const step 0.02; const kGrid grid( [ [bounds.minLng, bounds.minLat], [bounds.maxLng, bounds.minLat], [bounds.maxLng, bounds.maxLat], [bounds.minLng, bounds.maxLat] ], variogram, step );这里最关键的是step它直接决定最终生成的格子数量也就是Cesium要渲染的多边形数量。grid会在边界范围内按step切出n * g个网格点对每个点调用predict算温度值。step越小格网越细腻但多边形数量按平方级上升。我实测过在覆盖约100km乘100km范围的省级局部区域内不同step对应的情况如下step值估算格子数Cesium渲染方式建议视觉效果0.01度100 x 100 10000需要PrimitiveEntity会卡细腻接近气象产品0.02度50 x 50 2500Entity勉强能跑频繁交互略掉帧较细日常够用0.05度20 x 20 400Entity流畅能看出块状边缘3. 网格结果如何转成Cesium能认的图形3.1 网格数据结构拆解kGrid返回的对象里面的data是一个Float64Array存了所有网格点的温度预测值排列顺序是按列优先的。它同时给了几个关键的字段t是最小经度d是网格步长n是经度方向格子数g是纬度方向格子数xlim和ylim是边界范围。这意味着我可以从kGrid里拿到每一个网格中心的经纬度和温度值但Cesium不认识这种网格数据。要让它在三维地球上显示最直接的方式就是把每个格网转成一个规则矩形多边形根据温度值给它算一个颜色然后塞给Cesium渲染。温度值的范围可以从kGrid.zlim拿也可以自己在站点数据里算。3.2 色带映射函数颜色映射是出图好不好看的关键。我常用的是一套蓝-青-绿-黄-红的连续渐变色带低温深蓝、高温橙红和人感知冷暖的直觉一致。核心思路是把温度值标准化到0到1之间再在两个色标之间做线性插值function temperatureColor(value, minTemp, maxTemp) { const t Cesium.Math.clamp((value - minTemp) / (maxTemp - minTemp), 0, 1); const stops [0, 0.25, 0.5, 0.75, 1]; const colors [ Cesium.Color.fromBytes(30, 70, 190), Cesium.Color.fromBytes(70, 160, 220), Cesium.Color.fromBytes(110, 200, 120), Cesium.Color.fromBytes(250, 210, 80), Cesium.Color.fromBytes(220, 50, 50) ]; let i 0; while (i stops.length - 2 t stops[i 1]) { i; } const local (t - stops[i]) / (stops[i 1] - stops[i]); return Cesium.Color.lerp(colors[i], colors[i 1], local, new Cesium.Color()); }注意minTemp和maxTemp一定要用整个数据集里的合理全局范围别用单个渲染批次里的动态最小最大值。否则做时间序列时每个时次都在重新拉伸色带上一帧的蓝色可能对应20度下一帧的蓝色变成10度体感就是画面整体在忽明忽暗地闪烁。3.3 核心转换循环有了色带函数从kGrid到Cesium实体就只剩一个双循环const entities []; const minTemp kGrid.zlim[0]; const maxTemp kGrid.zlim[1]; for (let ix 0; ix kGrid.n; ix) { for (let iy 0; iy kGrid.g; iy) { const value kGrid.data[iy * kGrid.n ix]; if (Number.isNaN(value)) continue; const west kGrid.t ix * kGrid.d; const south kGrid.ylim[0] iy * kGrid.d; const east west kGrid.d; const north south kGrid.d; const positions Cesium.Cartesian3.fromDegreesArray([ west, south, east, south, east, north, west, north ]); entities.push({ polygon: { hierarchy: positions, material: temperatureColor(value, minTemp, maxTemp), perPositionHeight: false } }); } } viewer.entities.add(entities);perPositionHeight这里必须设成false让多边形贴地跟随地形走势。如果设成true所有格子都会停在0海拔在一个山区项目里会直接看到一层网格悬浮在空中或者钻进山体里面。3.4 Entity、GeoJsonDataSource还是Primitive上面用的是Entity优点是代码短、好调试鼠标点击和样式覆盖都很方便。但格网一多就露馅。grid生成的格子大约2500个时Entity整体加载还勉强一旦浏览器窗口缩放或者相机频繁旋转帧率会明显下滑。超过5000个格子交互就开始卡了。GeoJsonDataSource是我最早试的方案把格子拼成一个GeoJSON FeatureCollection再一次性加载。这个写法数据结构上很干净但GeoJsonDataSource内部要自己对每个Feature做样式解析、颜色生成大量Feature时的性能和逐Entity添加没有本质差别反而多了字符串序列化和坐标转换的开销我后来基本不这么干。真正的性能进阶是用Primitive PerInstanceColorAppearance把几千个格子的几何和颜色合成一次绘制调用const instances []; for (let ix 0; ix kGrid.n; ix) { for (let iy 0; iy kGrid.g; iy) { const value kGrid.data[iy * kGrid.n ix]; if (Number.isNaN(value)) continue; const west kGrid.t ix * kGrid.d; const south kGrid.ylim[0] iy * kGrid.d; const positions Cesium.Cartesian3.fromDegreesArray([ west, south, west kGrid.d, south, west kGrid.d, south kGrid.d, west, south kGrid.d ]); instances.push( new Cesium.GeometryInstance({ geometry: new Cesium.PolygonGeometry({ polygonHierarchy: new Cesium.PolygonHierarchy(positions) }), attributes: { color: Cesium.ColorGeometryAttribute.fromColor( temperatureColor(value, minTemp, maxTemp) ) } }) ); } } const primitive new Cesium.Primitive({ geometryInstances: instances, appearance: new Cesium.PerInstanceColorAppearance({ flat: true, translucent: false }) }); viewer.scene.primitives.add(primitive);Entity总数超过两三千是我切换Primitive的临界点。Primitive对实例颜色和几何的管理更贴近底层几百上千个实例丢给一次绘制比建几千个Entity对象省太多开销。4. 性能调优、动态更新与实战避坑4.1 分辨率与渲染性能的平衡写到这里你已经能跑出一张静态的温度色斑图了。但真实项目里问题往往不是“能不能画出来”而是“在交互状态下能不能保持流畅”。step的选择要结合站点密度来。站点本身只有几十个那step再小也没有意义空地方全靠插值猜网格再细腻也掩盖不了数据稀疏站点有几百个step0.02是个比较甜点的值保留细节的同时2500个格子对多数机器压力不大。如果要全省范围、几千个站点、还要3km左右的高分辨率格网前端算和前端渲染都别硬撑直接走服务端方案。我踩过一次比较深的坑是单纯追求细腻度把step设成0.005生成了两三万个格子Cesium页面基本处于“拖一下卡三秒”的状态。后来把格网分辨率降到0.02肉眼几乎看不出差别流畅度却回来了。4.2 离群值和色带归一化气象观测里偶尔会有站点的数据异常比如传感器故障导致的40度高温毛刺。如果用包含离群值的全局最大最小值去做色带拉伸整个温度场会变成一片暗色真正的温度差异全部被压缩在色带底部图面信息全丢了。处理办法是分级定色带范围。我通常取所有站点温度值的某个分位数比如2%到98%超出范围的统一截断到边界色画出来之后用GroudTr真值诡异点再单独排查。这比让一个异常值毁掉整张图的色带要好得多。4.3 时间序列温度场的动态切换气象可视化大概率不只做一个时次24小时、72小时的温度数据一批批过来需要在Cesium里做时间轴的切换动画。我最早用CallbackProperty实现动态换色想法是每个格子的颜色都动态计算结果2500个格子配上频繁更新CPU直接爆了色彩切换一顿一顿的。后来改成预计算方案在数据到达时把所有时次的kGrid先算好每个时次生成好对应的Primitive切换时销毁旧的、添加新的。虽然内存换时间但切换过程干净利落不卡顿。// 切换时次 function showTimeLayer(index) { if (currentPrimitive) { viewer.scene.primitives.remove(currentPrimitive); currentPrimitive.destroy(); } currentPrimitive timeLayers[index].primitive; viewer.scene.primitives.add(currentPrimitive); }销毁时不要只从场景里remove要记得调用destroy()释放显存和几何资源。Vue3项目里尤其要注意组件卸载时如果只销毁viewer而不清理这些Primitive刷新页面或者重新进入模块时浏览器显存会被旧资源越堆越多。4.4 高纬度区域和投影问题克里金插值本身是基于经纬度距离直接计算的这在低纬度、中小范围区域问题不大。但项目覆盖范围到了东北、西北这种高纬度地区直接用经纬度距离算权重会有偏差经度1度在赤道地区约111km在北纬60度地区只剩约55km同样跨1度经度实际物理距离差了一倍。站点温度的空间相关性本质是按真实物理距离衰减的所以计算结果会失真。处理思路有两种。一种是插值前把站点坐标投影到Web Mercator或者Albers等距圆锥投影在投影坐标系里算距离、做插值得到网格后再把网格坐标反算回经纬度渲染。另一种是干脆把这一步挪到服务端用pykrige这类库在真实的投影坐标系里算好把格网结果以GeoJSON或者GeoTIFF形式发给前端前端只负责渲染。4.5 服务端计算与前端渲染的架构取舍说到服务端这其实是我现在做温度场项目更推荐的分工方式。白天当项目站点数量超过1000、网格精度要求又高时浏览器端做克里金插值的耗时和内存占用都很吓人。kriging网格生成的过程是在前端主线程里跑的算一次可能就要一两秒再叠加一个时间序列里几十个时次预计算页面直接白屏好几秒体验很差。服务端用pykrige或者gstatR语言算好每个时次的规则格网输出一个压缩过的JSON或者GeoTIFF前端只负责解析和渲染姿态就从容很多。拿JSON举例服务端返回的可以是一个{ n, g, minLng, minLat, maxLng, maxLat, mean, max, data }结构前端拿到后同样走上面第3章的转换循环。这样前端计算压力降到最低插值参数、投影方式都可以由数据团队统一控制出问题也好追溯。4.6 几个容易翻车的细节写渲染代码时几个小细节值得单独记一下。网格循环里一定要跳过NaN值。插值范围外扩后边缘区域可能有无法预测的点kriging会返回NaN不跳过的话多边形层级会出错渲染直接异常。温度色斑图叠在有地形和3D Tiles的球面上时多边形贴地可能会出现闪烁和建筑物纹理互相“打架”。简单的处理是把温度层放在独立的一个数据源里并且给多边形材质一个极小的透明值或者使用classificationType控制只压平到地形不参与3D Tiles的贴合。项目只在大屏上看整体趋势时直接用0高度平面反而干净。大规模格网一次渲染时不要用viewer.entities.add(entities)这种全文覆盖的写法重复添加。每次切换时次先判断是否是同一个Entity数组避免对同一个数组反复执行add导致重复实体累积把浏览器拖垮。5. 项目落地时几个能少走弯路的判断最后说一些项目里沉淀下来的个人判断不一定所有项目通用但多数情况能帮你少走一步弯路。如果你的范围控制在市级、站点数量在几十到一百个前端用kriging.js加Entity方案足够不用上服务端也别让架构复杂化。省市级几百个站、多方数据源、时间序列多时次的场景建议前端只做显示层插值放服务端算格网分辨率统一由数据侧控制。色带边界要固定时间序列切换时画面才能稳定这是体感差异最大的一个细节。网格解析和渲染的核心流程前端可以做成一个独立模块输入统一的数据结构输出独立的Entity集合或Primitive实例这样换数据源、换插值算法都不至于动到渲染核心。温度数据这类空间插值可视化本质上是“离散点-连续场-图形渲染”三层结构把每一层的边界切清楚后期加湿度、气压、风速渲染时会轻松非常多。我就是这样一次次迭代过来的这套结构每一步都不是多余的。本文还有配套的精品资源点击获取
返回列表