克里金插值c程序_【OpenLayer 实战】实现克里金插值渲染图Kriging

0.OpenLayers

OpenLayers是一个开源的Javascript库，用来在Web浏览器显示地图。它提供API接口创建类似于Google Maps和Bing Maps的基于web的地理应用。

1.克里金插值

克里金插值也称作空间局部插值法，或空间自协方差最佳插值法，它以变异函数理论和结构分析为基础， 在特定区域内对区域化变量进行无偏最优估计，被广泛应用于土壤制图领域，是一种非常重要的地统计学方法。
具体算法解释见xg1990的文章《克里金(Kriging)插值的原理与公式推导》(https://zhuanlan.zhihu.com/p/25377842)

2.算法实现

克里金插值算法有开源的实现，项目的GitHub地址:kriging.js(https://github.com/oeo4b/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上。

3.使用ImageCanvas

上面说到，使用kriging.plot可以将格网数据渲染至canvas上。但是要将canvas叠加到地图上，还是需要借助OpenLayers的ImageCanvas组件(https://openlayers.org/en/latest/apidoc/module-ol_source_ImageCanvas-ImageCanvasSource.html)。

效果预览

let params={
        mapCenter:[114.360456, 30.538622],    maxValue:100,    krigingModel:'exponential',//model还可选'gaussian','spherical'    krigingSigma2:0,    krigingAlpha:100,    canvasAlpha:0.9,//canvas图层透明度    colors:["#006837", "#1a9850", "#66bd63", "#a6d96a", "#d9ef8b", "#ffffbf",        "#fee08b", "#fdae61", "#f46d43", "#d73027", "#a50026"],};let baseLayer = new ol.layer.Tile({
        title: "base",    source: new ol.source.OSM()});let map = new ol.Map({
        target: 'map',    layers: [baseLayer],    view: new ol.View({
            center: params.mapCenter,        projection: 'EPSG:4326',        zoom: 16    })});let WFSVectorSource=new ol.source.Vector();let WFSVectorLayer=new ol.layer.Vector(    {
            source:WFSVectorSource,    });map.addLayer(WFSVectorLayer);//添加选择和框选控件，按住Ctrl/⌘键，使用鼠标框选采样点let select = new ol.interaction.Select();map.addInteraction(select);let dragBox = new ol.interaction.DragBox({
        condition: ol.events.condition.platformModifierKeyOnly});map.addInteraction(dragBox);//创建10个位置随机、属性值随机的特征点for (let i = 0; i < 10; i++) {
        let feature = new ol.Feature({
            geometry: new ol.geom.Point([params.mapCenter[0]+