克里金插值
组件库简介
kriging.js是一个开源的克里金插值算法组件库,核心方法如下:
// 使用gaussian、exponential或spherical模型对数据集进行训练,返回的是一个variogram对象;
kriging.train(t, x, y, model, sigma2, alpha)
// 使用variogram对象使polygons范围描述的地理位置内的格网元素具备不一样的预测值,width是插值格点精度大小
kriging.grid(polygons,variogram,width)
// 将得到的格网grid渲染至canvas上。
kriging.plot(canvas,grid,xlim,ylim,colors):
https://github.com/sakitam-gis/kriging.js
kriging-contour组件库是一个基于克里金插值算法,根据离散点位置及其权重,生成等值面矢量数据(GeoJSON格式)和栅格数据(Canvas绘制图片),这些数据在任何WebGIS客户端上都可通用展示。
https://github.com/FreeGIS/kriging-contour
安装依赖库
# 安装克里金插值所需的插件
npm i @turf/turf
# 方法一
npm i @sakitam-gis/kriging
# 方法二
npm i kriging-contour
# 引入kriging.js
import kriging from '@sakitam-gis/kriging';
# 引入kriging-contour
import { getVectorContour, drawCanvasContour } from 'kriging-contour/dist/kriging-contour.js'
实现思路
首先划定一个区域利用turfJS随机一部分数据(也可使用真实点状数据集),然后插件渲染空间插值至canvas或生成面状geojson,然后添加到地图容器中进行展示。
turf生成随机数据
// 随机点的边界(折线的最大包围盒坐标)
const boundaries = turf.lineString([[110, 32], [118, 40], [120, 35]]);
// 随机50个点状要素数据
let positionData = turf.randomPoint(50, { bbox: turf.bbox(boundaries) });
// 再生成些随机数做属性
turf.featureEach(positionData, function (currentFeature, featureIndex) {
currentFeature.properties = { value: (Math.random() * 100).toFixed(2) };
});
方法一:kriging库插值
使用体验:在大量数据进行插值时,处理速度较慢;色带颜色如果较少,格网间颜色差异较大;插值的范围要按照特定的格式填入,否则无法加载。
<template>
<div class="app-contain">
<!-- leaflet 地图容器 -->
<canvas id="canvasMap" style="display: none;"></canvas>
<div id="myMap"></div>
<div class="controls">
<el-button color="#626aef" @click="startKriging('kriging')">普通克里金</el-button>
<el-button color="#626aef" @click="startKriging('Vector')">克里金矢量</el-button>
<el-button color="#626aef" @click="startKriging('Image')">克里金图像</el-button>
<el-button color="#626aef" @click="clearKriging()">清空</el-button>
</div>
</div>
</template>
<script setup>
// 引入样式
import L from 'leaflet';
import 'leaflet/dist/leaflet.css'
// 引入turf
import * as turf from '@turf/turf'
// 引入kriging.js
import kriging from '@sakitam-gis/kriging';
// 引入kriging-contour
import { getVectorContour, drawCanvasContour } from 'kriging-contour/dist/kriging-contour.js'
// VUE组件
import { onMounted, reactive, ref } from "vue"
// 天地图TK
let tdtKey = 'YOURS_TK'
let map = null;
let featureLayerGroup = null;
let imageLayerGroup = null;
const initMap = () => {
// 矢量地图
const tiandituMap = new L.TileLayer(`http://t0.tianditu.gov.cn/cva_c/wmts?layer=cva&style=default&tilematrixset=c&Service=WMTS&Request=GetTile&Version=1.0.0&Format=tiles&TileMatrix={z}&TileCol={x}&TileRow={y}&tk=${tdtKey}`,
{
tileSize: 512,
noWrap: true,
bounds: [[-90, -180], [90, 180]]
})
// 文字注记
const tiandituText = new L.TileLayer(`http://t0.tianditu.com/vec_c/wmts?layer=vec&style=default&tilematrixset=c&Service=WMTS&Request=GetTile&Version=1.0.0&Format=tiles&TileMatrix={z}&TileCol={x}&TileRow={y}&tk=${tdtKey}`,
{
tileSize: 512,
noWrap: true,
bounds: [[-90, -180], [90, 180]]
})
const layers = L.layerGroup([tiandituText, tiandituMap])
map = L.map('myMap', { //需绑定地图容器div的id
center: [39.56, 116.20], //初始地图中心
crs: L.CRS.EPSG4326,
zoom: 5, //初始缩放等级
maxZoom: 18, //最大缩放等级
minZoom: 0, //最小缩放等级
zoomControl: true, //缩放组件
attributionControl: false, //去掉右下角logol
scrollWheelZoom: true, //默认开启鼠标滚轮缩放
// 限制显示地理范围
maxBounds: L.latLngBounds(L.latLng(-90, -180), L.latLng(90, 180)),
layers: [layers] // 图层
})
// 矢量图层组
featureLayerGroup = new L.FeatureGroup().addTo(map).bringToFront()
// 图像图层组
imageLayerGroup = new L.FeatureGroup().addTo(map).bringToFront()
}
const startKriging = (krigingType) => {
// 随机点的边界(折线的最大包围盒坐标)
const boundaries = turf.lineString([[110, 32], [118, 40], [120, 35]]);
// 随机50个点状要素数据
let positionData = turf.randomPoint(50, { bbox: turf.bbox(boundaries) });
// 再生成些随机数做属性
turf.featureEach(positionData, function (currentFeature, featureIndex) {
currentFeature.properties = { value: (Math.random() * 100).toFixed(2) };
});
if ('Vector' == krigingType) {
showKrigingVector(boundaries, positionData);
} else if ('Image' == krigingType) {
showKrigingImage(boundaries, positionData)
} else if ('kriging' == krigingType) {
showKriging(boundaries, positionData)
}
}
const showKriging = (boundaries, positionData) => {
// 清空图层
clearKriging();
// 完全透明
let scope = L.geoJSON(boundaries, {
style: function () {
return {
fillColor: '6666ff',
color: 'red',
weight: 2,
opacity: 0,
fillOpacity: 0,
};
}
}).addTo(imageLayerGroup);
map.fitBounds(scope.getBounds());
//根据scope边界线,生成范围信息
let xlim = [scope.getBounds()._southWest.lng, scope.getBounds()._northEast.lng];
let ylim = [scope.getBounds()._southWest.lat, scope.getBounds()._northEast.lat];
function loadkriging(points) {
let canvas = document.getElementById("canvasMap");
canvas.width = 2000;
canvas.height = 1000;
// 数量
let pointLength = points.features.length;
let t = [];// 数值
let x = [];// 经度
let y = [];// 纬度
// 加载点数过多的话,会出现卡顿
for (let i = 0; i < pointLength; i++) {
x.push(points.features[i].geometry.coordinates[0]);
y.push(points.features[i].geometry.coordinates[1]);
t.push(points.features[i].properties.value);
// 将插值点展示到地图中,并添加提示文字(可删除,非核心代码)
L.circle([y[i], x[i]], { radius: 1 })
.addTo(imageLayerGroup)
.bindTooltip(points.features[i].properties.value, {
permanent: true, //是永久打开还是悬停打开
direction: 'top' //方向
}).openTooltip();
}
// 克里金插值参数
const params = {
krigingModel: 'exponential',//model还可选'gaussian','spherical'
krigingSigma2: 0,
krigingAlpha: 100,
canvasAlpha: 0.8,//canvas图层透明度-0.75
colors: ["#00A600", "#01A600", "#03A700", "#04A700", "#05A800", "#07A800", "#08A900", "#09A900", "#0BAA00", "#0CAA00", "#0DAB00", "#0FAB00", "#10AC00", "#12AC00", "#13AD00", "#14AD00", "#16AE00", "#17AE00", "#19AF00", "#1AAF00", "#1CB000", "#1DB000", "#1FB100", "#20B100", "#22B200", "#23B200", "#25B300", "#26B300", "#28B400", "#29B400", "#2BB500", "#2CB500", "#2EB600", "#2FB600", "#31B700", "#33B700", "#34B800", "#36B800", "#37B900", "#39B900", "#3BBA00", "#3CBA00", "#3EBB00", "#3FBB00", "#41BC00", "#43BC00", "#44BD00", "#46BD00", "#48BE00", "#49BE00", "#4BBF00", "#4DBF00", "#4FC000", "#50C000", "#52C100", "#54C100", "#55C200", "#57C200", "#59C300", "#5BC300", "#5DC400", "#5EC400", "#60C500", "#62C500", "#64C600", "#66C600", "#67C700", "#69C700", "#6BC800", "#6DC800", "#6FC900", "#71C900", "#72CA00", "#74CA00", "#76CB00", "#78CB00", "#7ACC00", "#7CCC00", "#7ECD00", "#80CD00", "#82CE00", "#84CE00", "#86CF00", "#88CF00", "#8AD000", "#8BD000", "#8DD100", "#8FD100", "#91D200", "#93D200", "#95D300", "#97D300", "#9AD400", "#9CD400", "#9ED500", "#A0D500", "#A2D600", "#A4D600", "#A6D700", "#A8D700", "#AAD800", "#ACD800", "#AED900", "#B0D900", "#B2DA00", "#B5DA00", "#B7DB00", "#B9DB00", "#BBDC00", "#BDDC00", "#BFDD00", "#C2DD00", "#C4DE00", "#C6DE00", "#C8DF00", "#CADF00", "#CDE000", "#CFE000", "#D1E100", "#D3E100", "#D6E200", "#D8E200", "#DAE300", "#DCE300", "#DFE400", "#E1E400", "#E3E500", "#E6E600", "#E6E402", "#E6E204", "#E6E105", "#E6DF07", "#E6DD09", "#E6DC0B", "#E6DA0D", "#E6D90E", "#E6D710", "#E6D612", "#E7D414", "#E7D316", "#E7D217", "#E7D019", "#E7CF1B", "#E7CE1D", "#E7CD1F", "#E7CB21", "#E7CA22", "#E7C924", "#E8C826", "#E8C728", "#E8C62A", "#E8C52B", "#E8C42D", "#E8C32F", "#E8C231", "#E8C133", "#E8C035", "#E8BF36", "#E9BE38", "#E9BD3A", "#E9BC3C", "#E9BB3E", "#E9BB40", "#E9BA42", "#E9B943", "#E9B945", "#E9B847", "#E9B749", "#EAB74B", "#EAB64D", "#EAB64F", "#EAB550", "#EAB552", "#EAB454", "#EAB456", "#EAB358", "#EAB35A", "#EAB35C", "#EBB25D", "#EBB25F", "#EBB261", "#EBB263", "#EBB165", "#EBB167", "#EBB169", "#EBB16B", "#EBB16C", "#EBB16E", "#ECB170", "#ECB172", "#ECB174", "#ECB176", "#ECB178", "#ECB17A", "#ECB17C", "#ECB17E", "#ECB27F", "#ECB281", "#EDB283", "#EDB285", "#EDB387", "#EDB389", "#EDB38B", "#EDB48D", "#EDB48F", "#EDB591", "#EDB593", "#EDB694", "#EEB696", "#EEB798", "#EEB89A", "#EEB89C", "#EEB99E", "#EEBAA0", "#EEBAA2", "#EEBBA4", "#EEBCA6", "#EEBDA8", "#EFBEAA", "#EFBEAC", "#EFBFAD", "#EFC0AF", "#EFC1B1", "#EFC2B3", "#EFC3B5", "#EFC4B7", "#EFC5B9", "#EFC7BB", "#F0C8BD", "#F0C9BF", "#F0CAC1", "#F0CBC3", "#F0CDC5", "#F0CEC7", "#F0CFC9", "#F0D1CB", "#F0D2CD", "#F0D3CF", "#F1D5D1", "#F1D6D3", "#F1D8D5", "#F1D9D7", "#F1DBD8", "#F1DDDA", "#F1DEDC", "#F1E0DE", "#F1E2E0", "#F1E3E2", "#F2E5E4", "#F2E7E6", "#F2E9E8", "#F2EBEA", "#F2ECEC", "#F2EEEE", "#F2F0F0", "#F2F2F2"
]
}
// 对数据集进行训练
let variogram = kriging.train(t, x, y, params.krigingModel, params.krigingSigma2, params.krigingAlpha);
// 将插值范围封装成特定格式
let bbox = turf.bbox(boundaries); // 外包矩形范围
// 根据外包矩形范围生成外包矩形面Polygon
let bboxPolygon = turf.bboxPolygon(bbox);
let positions = [];
bboxPolygon.geometry.coordinates[0].forEach((v) => {
positions.push([v[0], v[1]])
})
// 将边界封装成特定的格式
let range = [positions]
// 使用variogram对象使polygons描述的地理位置内的格网元素具备不一样的预测值,最后一个参数,是插值格点精度大小
let grid = kriging.grid(range, variogram, 0.05);
// 将得到的格网grid渲染至canvas上
kriging.plot(canvas, grid, [xlim[0], xlim[1]], [ylim[0], ylim[1]], params.colors);
}
//将canvas对象转换成image的URL
function returnImgae() {
let mycanvas = document.getElementById("canvasMap");
return mycanvas.toDataURL("image/png");
}
// 执行克里金插值函数
loadkriging(positionData);
let imageBounds = [[ylim[0], xlim[0]], [ylim[1], xlim[1]]];
L.imageOverlay(returnImgae(), imageBounds, { opacity: 0.8 }).addTo(imageLayerGroup);
}
// 生成矢量等值面并渲染
const showKrigingVector = (boundaries, positionData) => {
// 清空图层
clearKriging();
// 展点(可删除)
L.geoJSON(positionData, {
pointToLayer: function (feature, latlng) {
return L.circleMarker(latlng, {
radius: 5,
fillColor: '#6666ff',
fillOpacity: 1,
color: "#fff",
weight: 2,
});
}, onEachFeature(feature, layer) {
// 显示文字
let content = feature.properties.value
// marker的icon文字
let myIcon = L.divIcon({
html: `<div style="white-space: nowrap;color:#6666ff;">${content}</div>`,
iconAnchor: [0, 0],
className: 'my-div-icon',
iconSize: 120
});
let featureCenter = L.latLng(feature.geometry.coordinates[1], feature.geometry.coordinates[0]);
featureLayerGroup.addLayer(L.marker(featureCenter, { icon: myIcon }));
}
}).addTo(featureLayerGroup)
// 颜色色带
let colors = [{ fill: "#ffdc84" }, { fill: "#ffd782" },
{ fill: "#ffd281" }, { fill: "#ffcd7f" }, { fill: "#ffc87e" }, { fill: "#ffc37c" }, { fill: "#ffbe7a" }, {fill: "#ffb979"}, {fill: "#feb477"},{fill: "#feaf76"},
{fill: "#feaa74"}, {fill: "#fea573"}, {fill: "#fea071"}, {fill: "#fe9b6f"},
{fill: "#fe966e"}, {fill: "#fe906c"}, {fill: "#fe8b6b"}, {fill: "#fe8669"},
{fill: "#fe8167"}, {fill: "#fe7c66"}, {fill: "#fe7764"}, {fill: "#fe7263"},
{fill: "#fd6d61"}, {fill: "#fd6860"}, {fill: "#fd635e"}, {fill: "#fd5e5c"},
{fill: "#fd595b"}, {fill: "#fd5459"}, {fill: "#fd4f58"}, {fill: "#fd4a56"}]
// 等级分级
let levelV = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 250, 260, 270, 280, 290, 300, 400];
let kriging_contours = getVectorContour(positionData, 'value', {
model: 'exponential',
sigma2: 0,
alpha: 100
}, levelV, boundaries);
// 展示生成的矢量等值面
L.geoJSON(kriging_contours, {
style: function (feature) {
return {
fillColor: hotColor(feature.properties.value),
weight: 0,
fillOpacity: 0.3,
};
}
}).addTo(featureLayerGroup);
// 根据值来配色
function hotColor(d) {
let index = levelV.findIndex((item) => item >= d);
if (index > -1) {
return colors[index].fill
} else {
return colors[colors.length - 1].fill
}
}
}
// 生成图像等值面并渲染
const showKrigingImage = (boundaries, positionData) => {
// 清空图层
clearKriging();
// 完全透明
let scope = L.geoJSON(boundaries, {
style: function () {
return {
fillColor: '6666ff',
color: 'red',
weight: 2,
opacity: 0,
fillOpacity: 0,
};
}
}).addTo(imageLayerGroup);
map.fitBounds(scope.getBounds());
//根据scope边界线,生成范围信息
let xlim = [scope.getBounds()._southWest.lng, scope.getBounds()._northEast.lng];
let ylim = [scope.getBounds()._southWest.lat, scope.getBounds()._northEast.lat];
// 色带
let colors = ["#006837", "#1a9850", "#66bd63", "#a6d96a", "#d9ef8b", "#ffffbf", "#fee08b", "#fdae61", "#f46d43", "#d73027", "#a50026"]
// 画布
let canvas = document.getElementById("canvasMap");
canvas.width = 1000;
canvas.height = 1000;
let kriging_contours = drawCanvasContour(positionData, 'value', {
model: 'exponential',
sigma2: 0,
alpha: 100
}, canvas, [xlim[0], xlim[1]], [ylim[0], ylim[1]], colors);
//将canvas对象转换成image的URL
function returnImgae() {
let mycanvas = document.getElementById("canvasMap");
return mycanvas.toDataURL("image/png");
}
let imageBounds = [[ylim[0], xlim[0]], [ylim[1], xlim[1]]];
L.imageOverlay(returnImgae(), imageBounds, { opacity: 0.9 }).addTo(imageLayerGroup);
}
// 清空图层
const clearKriging = () => {
imageLayerGroup.clearLayers();
featureLayerGroup.clearLayers();
}
onMounted(() => {
initMap();
})
</script>
<style scoped>
#myMap {
width: 96vw;
height: 96vh;
}
.controls {
position: absolute;
top: 0px;
left: 200px;
padding: 15px;
z-index: 1000;
}
</style>
leaflet-canvas-marker插件,已解决了缩放时飘的问题,增加了添加标签的功能
<script src="../lib/leaflet/plugins/leaflet.canvas-markers.js"></script>
<script>
var map = L.map('map', {
center: [39.905963, 116.390813],
zoom: 15,
preferCanvas: true //使用canvas模式渲染矢量图形
});
//添加底图
var tiles = L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png').addTo(map);
//使用canvas模式渲染marker
var ciLayer = L.canvasIconLayer({}).addTo(map);
for (var i = 0; i < 100000; i++) {
var lat = 39.905963 + (Math.random() - Math.random()) * 3;
var lng = 116.390813 + (Math.random() - Math.random()) * 3;
//在这里设置图片和文字,调整位置,让文字显示在图片上
var icon = L.icon({
iconUrl: './img/bg.png',
iconSize: [84, 34],
iconAnchor: [40, 20],
//添加文字
text:i.toString(),
textAnchor: [5, 0],
textFont:'14px bold', //设置字体大小和样式
textFillStyle:'#FFFFFF' //设置字体颜色
});
var marker = L.marker([lat, lng], {
icon: icon
})
.bindPopup("I Am " + i); //绑定气泡窗口
ciLayer.addLayer(marker);
}
</script>