前言
程序生成的插值分析图需要人工去制作，结合后端接口后可使用网页自动绘制插值分析图，网上搜索发现大多用 kriging.js 绘制，本文也使用 kriging.js 绘制。
一、kriging.js
kriging.js 是一个用于地统计插值的 JavaScript 库，它提供了多种插值方法，并且可以在 Web 浏览器中运行。
二、程序简介
1. 定义测试数据
代码如下（示例数据，随机生成的）：
let testDataJson = [
{"city":"江苏","latitude":"118.8062","longitude":"31.9208","value":"22"},
{"city":"黑龙江","latitude":"127.9688","longitude":"45.368","value":"25"},
{"city":"内蒙古","latitude":"110.3467","longitude":"41.4899","value":"26"},
{"city":"吉林","latitude":"125.8154","longitude":"44.2584","value":"17"},
{"city":"北京市","latitude":"116.4551","longitude":"40.2539","value":"18"},
{"city":"辽宁","latitude":"123.1238","longitude":"42.1216","value":"23"},
{"city":"河北","latitude":"114.4995","longitude":"38.1006","value":"24"},
{"city":"天津","latitude":"117.4219","longitude":"39.4189","value":"28"},
{"city":"山西","latitude":"112.3352","longitude":"37.9413","value":"32"},
{"city":"陕西","latitude":"109.1162","longitude":"34.2004","value":"29"},
{"city":"甘肃","latitude":"103.5901","longitude":"36.3043","value":"33"},
{"city":"宁夏","latitude":"106.3586","longitude":"38.1775","value":"31"},
{"city":"青海","latitude":"101.4038","longitude":"36.8207","value":"20"},
{"city":"新疆","latitude":"87.9236","longitude":"43.5883","value":"22"},
{"city":"四川","latitude":"103.9526","longitude":"30.7617","value":"22"},
{"city":"重庆","latitude":"108.384366","longitude":"30.439702","value":"25"},
{"city":"山东","latitude":"117.1582","longitude":"36.8701","value":"26"},
{"city":"河南","latitude":"113.4668","longitude":"34.6234","value":"17"},
{"city":"安徽","latitude":"117.29","longitude":"32.0581","value":"18"},
{"city":"湖北","latitude":"114.3896","longitude":"30.6628","value":"23"},
{"city":"浙江","latitude":"119.5313","longitude":"29.8773","value":"24"},
{"city":"福建","latitude":"119.4543","longitude":"25.9222","value":"28"},
{"city":"江西","latitude":"116.0046","longitude":"28.6633","value":"32"},
{"city":"湖南","latitude":"113.0823","longitude":"28.2568","value":"29"},
{"city":"贵州","latitude":"106.6992","longitude":"26.7682","value":"33"},
{"city":"云南","latitude":"102.9199","longitude":"25.4663","value":"31"},
{"city":"广东","latitude":"113.12244","longitude":"23.009505","value":"28"},
{"city":"广西","latitude":"108.479","longitude":"23.1152","value":"32"},
{"city":"海南","latitude":"110.3893","longitude":"19.8516","value":"29"},
{"city":"上海","latitude":"121.4648","longitude":"31.2891","value":"33"},
{"city":"西藏","latitude":"91.11","longitude":"29.97","value":"31"},
];
2. 定义 kriging 默认值
代码如下（示例数据，随机生成的）：
let params = {
mapCenter: [103.5901, 36.3043], // 地图中心点
// maxValue: 100,
krigingModel: 'exponential', // 'exponential':指数模型 'gaussian':高斯模型,'spherical':球形模型
krigingSigma2: 0, // 方差参数，通常可设为0
krigingAlpha: 100, // 变差函数模型的先验，通常可设为100
canvasAlpha: 1, // canvas插值图层透明度
colors:[ // 颜色区域集合（插值分析会按照当前颜色集进行渲染）
"rgb(0,0,255)",
"rgb(54,97,255)",
"rgb(56,172,255)",
"rgb(0,255,255)",
"rgb(145,255,180)",
"rgb(210,255,105)",
"rgb(255,255,0)",
"rgb(255,183,0)",
"rgb(255,111,0)",
"rgb(255,0,0)"
]
};
3. 绘制地图
代码如下（示例），当前使用天地图底图，需要单独申请key，或者可参考前一篇文章 openlayers显示指定区域，其它区域隐藏 使用高德底图影像：
let vec_c = getTdtLayer(ol, 'vec_w') // 默认
let img_c = getTdtLayer(ol, 'img_w') // 卫星影像
let cva_c = getTdtLayer(ol, 'cva_w') // 路线
let map = new ol.Map({
target: 'map',
layers: [img_c, getChina()],
view: new ol.View({
center: params.mapCenter,
projection: 'EPSG:4326',
zoom: 5
})
});
// 添加各省市点位信息
let obj = { "type": "FeatureCollection", "features": [] };
for (let i = 0; i < testDataJson.length; i++) {
// 地图站点
obj.features.push({
"type": "Feature",
"properties": {
"name": testDataJson[i].city,
"value": testDataJson[i].value
},
"geometry": { "type": "Point", "coordinates": [testDataJson[i].longitude, testDataJson[i].latitude] }
})
}
let myfeatures = (new ol.format.GeoJSON({
featureProjection: 'EPSG:4326',
})).readFeatures(obj);
let mySource = new ol.layer.Vector({ //(矢量图层)
source: new ol.source.Vector({
features: myfeatures
}),
style: function (feature, resolution) {
return new ol.style.Style({
text: new ol.style.Text({
font: 'bold 12px "Open Sans", "Arial Unicode MS", "sans-serif"',
offsetX: 20,
offsetY: 20,
text: feature.get('name') + '\n' + feature.get('value') + '℃',
textAlign: 'left',
fill: new ol.style.Fill({
color: '
#fff
'
})
})
})
}
})
mySource.setZIndex(1);
map.addLayer(mySource);
function getTdtLayer(ol, lyr) {
let key = 'key';
let url = 'http://t{0-7}.tianditu.com/DataServer?T=' + lyr + '&X={x}&Y={y}&L={z}&tk=' + key
let layer = new ol.layer.Tile({
source: new ol.source.XYZ({
url: url
})
})
return layer
}
function getChina() {
return new ol.layer.Vector({
source: new ol.source.Vector({
url: "https://geo.datav.aliyun.com/areas_v3/bound/100000_full.json",
format: new ol.format.GeoJSON(),
}),
style: function (feature) {
return new ol.style.Style({
zIndex: 0,
stroke: new ol.style.Stroke({
color: "
#6FE9FF
",
width: 2
}),
text: new ol.style.Text({
font: 'normal 12px "Open Sans", "Arial Unicode MS", "sans-serif"',
text: feature.get('name'),
textAlign: 'left',
fill: new ol.style.Fill({
color: '
#6FE9FF
'
})
}),
})
},
})
}
4. 加载 kriging 绘制区域及剪切底图
$.getJSON('https://geo.datav.aliyun.com/areas_v3/bound/100000.json', function(data){
let coord = data.features[0].geometry.coordinates[0];
drawKriging(coord);
// 屏蔽除了加载的中国地图外的其他区域（默认裁切大的地图，其它区域可以使用concat拼接，或者修改上述底图链接，可生成省级、市级数据示例）
var coords = data['features'][0]['geometry']['coordinates'];
var f = new ol.Feature(new ol.geom.MultiPolygon(coords));
var crop = new ol.filter.Crop({
feature: f,
wrapX: true,
inner: false
});
img_c.addFilter(crop);
var mask = new ol.filter.Mask({
feature: f,
wrapX: true,
inner: false,
fill: new ol.style.Fill({
color: [8, 20, 36, 0.8]
})
});
img_c.addFilter(mask);
})
5. 绘制插值分析
const drawKriging = (coord) => {
let kriging = test.kriging;
let values = [],
lngs = [],
lats = [];
for (let i = 0; i < testDataJson.length; i++) {
values.push(testDataJson[i].value)
lngs.push(testDataJson[i].longitude)
lats.push(testDataJson[i].latitude)
}
// 使用kriging绘制插值分析起码得三个点位数据
if (values.length >= 3){
let variogram = kriging.train(
values, // 样本值数组（如每个站点的降雨量）
lngs, // 样本点的经度坐标数组
lats, // 样本点的纬度坐标数组
params.krigingModel, // 变差函数模型，可选 'gaussian'（高斯模型）、'exponential'（指数模型）或 'spherical'（球形模型）
params.krigingSigma2, // 高斯过程的方差参数，通常可设为 0
params.krigingAlpha //变差函数模型的先验，通常可设为 100
);
let grid = kriging.grid(
coord, // 定义插值区域的边界多边形坐标数组
variogram, // 训练好的变差函数模型
0.1// 每个栅格单元的宽度（度），值越小，生成的分辨率越高；数值越小,轮廓越细腻，但是占浏览器内存越大（很容易崩溃）
);
//移除已有图层
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("有效样点个数不足，无法插值");
}
}
三、完整程序
完整网页代码如下，可直接运行：
<!DOCTYPE html>
<html>
<head>
<title>结合openlayers、ol-ext、kriging.js绘制全国气温插值分析图</title>
<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/ol@v10.6.0/ol.css">
<script src="https://cdn.jsdelivr.net/npm/jquery@3.6.4/dist/jquery.min.js"></script>
<script src="https://cdn.jsdelivr.net/npm/ol@v10.6.0/dist/ol.js"></script>
<script src="https://cdn.jsdelivr.net/npm/ol-ext@4.0.35/dist/ol-ext.min.js"></script>
<script src="https://cdn.jsdelivr.net/npm/kriging@0.1.12/dist/kriging.min.js"></script>
<style>
* {
margin: 0;
padding: 0;
}
#map
{
width: 100vw;
height: 100vh;
}
</style>
</head>
<body>
<div id="map" class="map"></div>
<script>
$(function() {
let params = {
mapCenter: [103.5901, 36.3043],
// maxValue: 100,
krigingModel: 'exponential', // 'exponential':指数 'gaussian':高斯,'spherical':球体
krigingSigma2: 0,
krigingAlpha: 100,
canvasAlpha: 1, //canvas图层透明度
colors:[
"rgb(0,0,255)",
"rgb(54,97,255)",
"rgb(56,172,255)",
"rgb(0,255,255)",
"rgb(145,255,180)",
"rgb(210,255,105)",
"rgb(255,255,0)",
"rgb(255,183,0)",
"rgb(255,111,0)",
"rgb(255,0,0)"
]
};
let testDataJson = [
{"city":"江苏","longitude":"118.8062","latitude":"31.9208","value":"22"},
{"city":"黑龙江","longitude":"127.9688","latitude":"45.368","value":"25"},
{"city":"内蒙古","longitude":"110.3467","latitude":"41.4899","value":"26"},
{"city":"吉林","longitude":"125.8154","latitude":"44.2584","value":"17"},
{"city":"北京市","longitude":"116.4551","latitude":"40.2539","value":"18"},
{"city":"辽宁","longitude":"123.1238","latitude":"42.1216","value":"23"},
{"city":"河北","longitude":"114.4995","latitude":"38.1006","value":"24"},
{"city":"天津","longitude":"117.4219","latitude":"39.4189","value":"28"},
{"city":"山西","longitude":"112.3352","latitude":"37.9413","value":"32"},
{"city":"陕西","longitude":"109.1162","latitude":"34.2004","value":"29"},
{"city":"甘肃","longitude":"103.5901","latitude":"36.3043","value":"33"},
{"city":"宁夏","longitude":"106.3586","latitude":"38.1775","value":"31"},
{"city":"青海","longitude":"101.4038","latitude":"36.8207","value":"20"},
{"city":"新疆","longitude":"87.9236","latitude":"43.5883","value":"22"},
{"city":"四川","longitude":"103.9526","latitude":"30.7617","value":"22"},
{"city":"重庆","longitude":"108.384366","latitude":"30.439702","value":"25"},
{"city":"山东","longitude":"117.1582","latitude":"36.8701","value":"26"},
{"city":"河南","longitude":"113.4668","latitude":"34.6234","value":"17"},
{"city":"安徽","longitude":"117.29","latitude":"32.0581","value":"18"},
{"city":"湖北","longitude":"114.3896","latitude":"30.6628","value":"23"},
{"city":"浙江","longitude":"119.5313","latitude":"29.8773","value":"24"},
{"city":"福建","longitude":"119.4543","latitude":"25.9222","value":"28"},
{"city":"江西","longitude":"116.0046","latitude":"28.6633","value":"32"},
{"city":"湖南","longitude":"113.0823","latitude":"28.2568","value":"29"},
{"city":"贵州","longitude":"106.6992","latitude":"26.7682","value":"33"},
{"city":"云南","longitude":"102.9199","latitude":"25.4663","value":"31"},
{"city":"广东","longitude":"113.12244","latitude":"23.009505","value":"28"},
{"city":"广西","longitude":"108.479","latitude":"23.1152","value":"32"},
{"city":"海南","longitude":"110.3893","latitude":"19.8516","value":"29"},
{"city":"上海","longitude":"121.4648","latitude":"31.2891","value":"33"},
{"city":"西藏","longitude":"91.11","latitude":"29.97","value":"31"},
];
let vec_c = getTdtLayer(ol, 'vec_w') // 默认
let img_c = getTdtLayer(ol, 'img_w') // 卫星影像
let cva_c = getTdtLayer(ol, 'cva_w') // 路线
// 绘制地图
let map = new ol.Map({
target: 'map',
layers: [img_c, getChina()],
view: new ol.View({
center: params.mapCenter,
projection: 'EPSG:4326',
zoom: 5
})
});
// 添加各省市点位信息
let obj = { "type": "FeatureCollection", "features": [] };
for (let i = 0; i < testDataJson.length; i++) {
// 地图站点
obj.features.push({
"type": "Feature",
"properties": {
"name": testDataJson[i].city,
"value": testDataJson[i].value
},
"geometry": { "type": "Point", "coordinates": [testDataJson[i].longitude, testDataJson[i].latitude] }
})
}
let myfeatures = (new ol.format.GeoJSON({
featureProjection: 'EPSG:4326',
})).readFeatures(obj);
let mySource = new ol.layer.Vector({ //(矢量图层)
source: new ol.source.Vector({
features: myfeatures
}),
style: function (feature, resolution) {
return new ol.style.Style({
text: new ol.style.Text({
font: 'bold 12px "Open Sans", "Arial Unicode MS", "sans-serif"',
offsetX: 20,
offsetY: 20,
text: feature.get('name') + '\n' + feature.get('value') + '℃',
textAlign: 'left',
fill: new ol.style.Fill({
color: '
#fff
'
})
}),
// image: new ol.style.Icon(({
// anchor: [0.5, 0],
// anchorXUnits: 'fraction',
// anchorYUnits: 'pixels',
// src: './yun.png',
// scale: .5
// }))
})
}
})
mySource.setZIndex(1);
map.addLayer(mySource);
// 绘制kriging插值图
let canvasLayer = null;
const drawKriging = (coord) => {
let kriging = test.kriging;
// testDataJson
let values = [],
lngs = [],
lats = [];
for (let i = 0; i < testDataJson.length; i++) {
values.push(testDataJson[i].value)
lngs.push(testDataJson[i].longitude)
lats.push(testDataJson[i].latitude)
}
// 使用kriging绘制插值分析起码得三个点位数据
if (values.length >= 3){
let variogram = kriging.train(
values, // 样本值数组（如每个站点的降雨量）
lngs, // 样本点的经度坐标数组
lats, // 样本点的纬度坐标数组
params.krigingModel, // 变差函数模型，可选 'gaussian'（高斯模型）、'exponential'（指数模型）或 'spherical'（球形模型）
params.krigingSigma2, // 高斯过程的方差参数，通常可设为 0
params.krigingAlpha //变差函数模型的先验，通常可设为 100
);
let grid = kriging.grid(
coord, // 定义插值区域的边界多边形坐标数组
variogram, // 训练好的变差函数模型
0.1 // 每个栅格单元的宽度（度），值越小，生成的分辨率越高；数值越小,轮廓越细腻，但是占浏览器内存越大（很容易崩溃）
);
//移除已有图层
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("有效样点个数不足，无法插值");
}
}
// 绘制插值图、添加遮罩
$.getJSON('https://geo.datav.aliyun.com/areas_v3/bound/100000.json', function(data){
let coord = data.features[0].geometry.coordinates[0];
drawKriging(coord);
// 屏蔽除了加载的中国地图外的其他区域
var coords = data['features'][0]['geometry']['coordinates'];
var f = new ol.Feature(new ol.geom.MultiPolygon(coords));
var crop = new ol.filter.Crop({
feature: f,
wrapX: true,
inner: false
});
img_c.addFilter(crop);
var mask = new ol.filter.Mask({
feature: f,
wrapX: true,
inner: false,
fill: new ol.style.Fill({
color: [8, 20, 36, 0.8]
})
});
img_c.addFilter(mask);
})
})
function getTdtLayer(ol, lyr) {
let key = 'key';
let url = 'http://t{0-7}.tianditu.com/DataServer?T=' + lyr + '&X={x}&Y={y}&L={z}&tk=' + key
let layer = new ol.layer.Tile({
source: new ol.source.XYZ({
url: url
})
})
return layer
}
function getChina() {
return new ol.layer.Vector({
source: new ol.source.Vector({
url: "https://geo.datav.aliyun.com/areas_v3/bound/100000_full.json",
format: new ol.format.GeoJSON(),
}),
style: function (feature) {
return new ol.style.Style({
zIndex: 0,
stroke: new ol.style.Stroke({
color: "
#6FE9FF
",
width: 2
}),
text: new ol.style.Text({
font: 'normal 12px "Open Sans", "Arial Unicode MS", "sans-serif"',
text: feature.get('name'),
textAlign: 'left',
fill: new ol.style.Fill({
color: '
#6FE9FF
'
})
}),
})
},
})
}
</script>
</body>
</html>
四、绘制效果
？问题
跪求大神指导🙋
- 使用 kriging.js 绘制时，点位数据不能多余 500 个，数据量太多会报错且绘制不出来效果，怎么实现类似 arcmap 程序生成的多数据插值分析？
- 使用 kriging.js 绘制时，如何实现类似 arcmap 程序按照图例数据级颜色级别进行渲染？
- 使用 kriging.js 绘制时，如何保障绘制准确性，用上述相同数据及 10 层颜色使用 arcmap 进行克里金及反距离权重法生成的图结果完全不对等？
