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组件。
<div id="map"></div>
#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 < 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);