Vue+leaflet笔记之克里金插值

文章目录

Vue+leaflet笔记之克里金插值
- 开发环境
- 代码简介
- - 组件库简介
  - 安装依赖库
  - 实现思路
  - - turf生成随机数据
    - 方法一：kriging库插值
    - 方法二：kriging-contour插值(矢量)
    - 方法三：kriging-contour插值(栅格)

本文介绍了Web端使用 Leaflet开发库进行克里金插值的三种方法 (底图来源:天地图)，分别结合了kriging、kriging-contour组件库实现克里金插值功能，效果如下图所示。

开发环境

Vue开发库：3.2.37 & Leaflet开发库：1.9.3

Leaflet主要插件：turf、kriging.js、kriging-contour

代码简介

组件库简介

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库插值

使用体验：在大量数据进行插值时，处理速度较慢；色带颜色如果较少，格网间颜色差异较大；插值的范围要按照特定的格式填入，否则无法加载。

核心代码

效果预览

方法二：kriging-contour插值(矢量)

使用体验：插值效果较为平滑；可以对插值效果进行裁剪，仅保留研究区的范围；颜色样式定制更方便一些。

核心代码

效果预览

方法三：kriging-contour插值(栅格)

使用体验：相比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>