OpenLayers实现克里金插值渲染图-Kriging

0.OpenLayers

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

1.克里金插值

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

2.算法实现

克里金插值算法有开源的实现，项目的GitHub地址:kriging.js

如果你想通过es6的方式使用这个工具，这里是我修改过的源码:http://c.trojx.me/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.使用示例

上面说到，使用kriging.plot可以将格网数据渲染至canvas上。但是要将canvas叠加到地图上，还是需要借助OpenLayers的ImageCanvas组件。

效果预览

#map{
  height:100%;
  width:100%;
}

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 {
    let extent = dragBox.getGeometry().getExtent();
    WFSVectorSource.forEachFeatureIntersectingExtent(extent, (feature)=> {
        selectedFeatures.push(feature);
    });
    drawKriging(extent);
});
dragBox.on('boxstart', ()=>{
    selectedFeatures.clear();
});

//绘制kriging插值图
let canvasLayer=null;
const drawKriging=(extent)=>{
    let values=[],lngs=[],lats=[];
    selectedFeatures.forEach(feature=>{
        values.push(feature.values_.value);
        lngs.push(feature.values_.geometry.flatCoordinates[0]);
        lats.push(feature.values_.geometry.flatCoordinates[1]);
    });
    if (values.length>3){
        let variogram=kriging.train(values,lngs,lats,
            params.krigingModel,params.krigingSigma2,params.krigingAlpha);

        let polygons=[];
        polygons.push([[extent[0],extent[1]],[extent[0],extent[3]],
            [extent[2],extent[3]],[extent[2],extent[1]]]);
        let grid=kriging.grid(polygons,variogram,(extent[2]-extent[0])/200);

        let dragboxExtent=extent;
        //移除已有图层
        if (canvasLayer !== null){
            map.removeLayer(canvasLayer);
        }
        //创建新图层
        canvasLayer=new ol.layer.Image({
            source: new ol.source.ImageCanvas({
                canvasFunction:(extent, resolution, pixelRatio, size, projection) =>{
                    let canvas = document.createElement('canvas');
                    canvas.width = size[0];
                    canvas.height = size[1];
                    canvas.style.display='block';
                    //设置canvas透明度
                    canvas.getContext('2d').globalAlpha=params.canvasAlpha;                          

                    //使用分层设色渲染
                    kriging.plot(canvas,grid,
                        [extent[0],extent[2]],[extent[1],extent[3]],params.colors);

                    return canvas;
                },
                projection: 'EPSG:4326'
            })
        })
        //向map添加图层
        map.addLayer(canvasLayer);
    }else {
        alert("有效样点个数不足，无法插值");
    }
}
//首次加载，自动渲染一次差值图
let extent = [params.mapCenter[0]-.005,params.mapCenter[1]-.005,params.mapCenter[0]+.005,params.mapCenter[1]+.005];
    WFSVectorSource.forEachFeatureIntersectingExtent(extent, (feature)=> {
        selectedFeatures.push(feature);
    });
drawKriging(extent);

4.在线调试

在线编辑示例代码：https://codepen.io/jianxunrao/pen/oadBPq