codepen demo

https://codepen.io/jianxunrao/pen/oadBPq

<div id="map">
</div>
#map{
  height:100%;
  100%;
}
let params={
    mapCenter:[114.360456, 30.538622],
    maxValue:100,
    krigingModel:'exponential',//model还可选'gaussian','spherical'
    krigingSigma2:0,
    krigingAlpha:100,
    canvasAlpha:0.75,//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]+Math.random()*0.01-.005,params.mapCenter[1]+Math.random()*0.01-.005]), value: Math.round(Math.random()*params.maxValue)
    });
    feature.setStyle(new ol.style.Style({
        image: new ol.style.Circle({
            radius: 6,
            fill: new ol.style.Fill({color: "#00F"})
        })
    }));
    WFSVectorSource.addFeature(feature);
}

//设置框选事件
let selectedFeatures = select.getFeatures();
dragBox.on('boxend', ()=>{
    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);
 /**
     * 算法来源:https://www.cnblogs.com/imgss/p/10707150.html 经纬度转笛卡尔坐标
     * lat,lng为弧度表示的经纬度,r为地球半径,由于是算夹角,r是多少不重要
     *  *
     * @param lat
     * @param lng
     * @param r
     * @returns {{x: number, y: number, z: number}}
     */
    ball2xyz(lat, lng, r = 6400) {
        try {
            return {
                x: r * Math.cos(lat) * Math.cos(lng),
                y: r * Math.cos(lat) * Math.sin(lng),
                z: r * Math.sin(lat)
            };
        } catch (e) {

        }
    }

    /**
     *  将地理经纬度转换成笛卡尔坐标系
     *  https://blog.csdn.net/reborn_lee/article/details/82497577
     *
     * @param lng
     * @param lat
     * @returns {*}
     */
    geo2xyz({lng, lat}) {
        try {
            let thera = (Math.PI * lat) / 180;
            let fie = (Math.PI * lng) / 180;
            return this.ball2xyz(thera, fie);
        } catch (e) {

        }
    }

    /**
     * 计算3个地理坐标点之间的夹角
     *
     * @param l1
     * @param l2
     * @param l3
     * @returns {number}
     */
    angleOflocation(l1, l2, l3) {
        try {
            let p1 = this.geo2xyz(l1);
            let p2 = this.geo2xyz(l2);
            let p3 = this.geo2xyz(l3);
            let {x: x1, y: y1, z: z1} = p1;
            let {x: x2, y: y2, z: z2} = p2;
            let {x: x3, y: y3, z: z3} = p3;
            // 计算向量 P2P1 和 P2P3 的夹角 https://www.zybang.com/question/3379a30c0dd3041b3ef966803f0bf758.html
            let _P1P2 = Math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2);
            let _P2P3 = Math.sqrt((x3 - x2) ** 2 + (y3 - y2) ** 2 + (z3 - z2) ** 2);
            let P = (x1 - x2) * (x3 - x2) + (y1 - y2) * (y3 - y2) + (z1 - z2) * (z3 - z2); //P2P1*P2P3
            return (Math.acos(P / (_P1P2 * _P2P3)) / Math.PI) * 180;
        } catch (e) {

        }
    }
原文地址:https://www.cnblogs.com/hustshu/p/14769286.html