cesium 实现克里金生成矢量等值面，使用worker浏览器线程

需求背景

需要实现一个上百个栅格存在数值的等值面，为了提高用户体验
1.需要插值，考虑使用 kriging 算法，原有的kriging算法，未优化前，计算300点位就会造成栈溢出
2.kriging前端计算量大，使用worker浏览器线程，避免阻塞主线程

解决效果

index.vue

<!--/**
* @author: liuk
* @date: 2024-07-29
* @describe:克里金生成矢量等值面
*/-->
<template>
  <div class="ybPanel-wrap">
    <div class="operation-area">
      <div class="ctrl-row">
        <el-row style="width: 100%;height: 100%">
          <el-col :span="24">
            <el-form-item label="" style="margin-top: 5px">
              <el-checkbox v-model="gridShow" style="position: relative;top:-5px;margin-right: 20px"
                           @change="gridShowChange">
                <span style="font-size: 14px;color: #fff">--</span>
              </el-checkbox>
            </el-form-item>
          </el-col>
          <el-col :span="24">
            <el-form-item label="">
              <el-radio-group v-model="formData.type">
                <el-radio-button value="h">--</el-radio-button>
                <el-radio-button value="total">--</el-radio-button>
              </el-radio-group>
            </el-form-item>
          </el-col>
        </el-row>
      </div>
      <div class="time-player">
        <span class="iconfont icon-icon_prev" @click="doPrev"></span>
        <span class="play-icon iconfont" @click="togglePlay" :class="playFlag ? 'icon-pause' : 'icon-play'"/>
        <span class="iconfont icon-icon_next2" @click="doNext"></span>
        <el-slider v-model="curTime" :show-tooltip="false" :min="timeList[0]" :max="timeList[timeList.length - 1]"
                   tooltip-class="toolTip"
                   :step="60*60*1e3" style="width: 700px;margin-left: 10px"/>
        <div class="tip-area">
          <div class="tip" ref="tipRef">
            <span>{{ moment(curTime).format(`YYYY-MM-DD HH:mm:ss`) }}</span>
          </div>
          <div class="min-val">{{ moment(timeList[0]).format(`YYYY-MM-DD HH:mm:ss`) }}</div>
          <div class="max-val">{{ moment(timeList[timeList.length - 1]).format(`YYYY-MM-DD HH:mm:ss`) }}</div>
        </div>
        <el-select class="tip-area-select" v-model="speedType" size="small" placeholder="" style="width: 40px">
          <el-option v-for="item in speedOptions" :key="item.value" :label="item.label" :value="item.value"/>
        </el-select>
      </div>
    </div>
  </div>
</template>

<script lang="ts" setup>
import {reactive, toRefs, onMounted, defineProps, defineEmits, watch, onUnmounted, ref, computed} from "vue"
import {theHourTime, formatToFixed, cartesianToWgs84} from "@/utils/dictionary";
import moment from "moment";
import {ElMessage, ElLoading} from "element-plus";
import {getDrainageBasinGridByCY, getLon2LatByCY, getTimeByCY} from "@/api/forecast";
// Props
const props = defineProps(['modelValue'])
// Emits
const emit = defineEmits(['update:modelValue'])
// Refs
const tipRef = ref(null)


const model = reactive({
  formData: {
    type: 'h',
    maxEndTime: moment().valueOf(),
  },
  gridShow: false,// 是否展示栅格
  playFlag: false,// 是否播放
  curTime: 0, // 当前时间
  speedType: 1,// 播放数据
  timeList: [],// 所有数据时间
  curLoadingIndex: 0,
})
const {
  formData,
  gridShow,
  playFlag,
  curTime,
  timeList,
  speedType,
  curLoadingIndex
} = toRefs(model)

watch(() => model.curTime, (curTime) => {
  if (!tipRef.value) return
  tipRef.value.style.left = (curTime - model.timeList[0]) / (model.timeList[model.timeList.length - 1] - model.timeList[0]) * 100 + "%"
  krigingModel.loading || updateGridColor()
})

onMounted(() => {
  getlist()
})

onUnmounted(() => {
  window.clearInterval(timer)
  viewer.entities.remove(heatmapEntity)
  viewer.dataSources.remove(pointGridDatasouce);
})

const getlist = async () => {
  // const {data: data1} = await getLon2LatByCY()
  // krigingModel.lonLats = data1
  createHeatMap()
  // const {data: data2} = await getTimeByCY()
  // model.timeList = data2.map(item => item * 1e3)
  // const longitudes = krigingModel.lonLats.map(item => item.split("_").map(Number)[0])
  // const latitude = krigingModel.lonLats.map(item => item.split("_").map(Number)[1])
  const timeList = [1724641200000, 1724644800000, 1724648400000, 1724652000000, 1724655600000, 1724659200000, 1724662800000, 1724666400000, 1724670000000, 1724673600000, 1724677200000, 1724680800000, 1724684400000, 1724688000000, 1724691600000, 1724695200000, 1724698800000, 1724702400000, 1724706000000, 1724709600000, 1724713200000, 1724716800000, 1724720400000, 1724724000000]
  const loadingInstance = ElLoading.service({
    fullscreen: true,
    text: '模型生成中...',
    maskColor: 'rgba(41, 12, 12, .5)'
  })
  model.timeList = timeList
  krigingModel.dataMap.clear()
  const worker = new Worker('/worker/worker.js');
  worker.postMessage({timeList});

  const colors = [
    {min: 0, max: 1, color: "rgba(233,233,233,.03)"},
    {min: 1, max: 1.5, color: "rgba(163,243,146,.03)"},
    {min: 1.6, max: 6.9, color: "rgba(57,165,1,.03)"},
    {min: 7, max: 14.9, color: "rgba(99,183,255,.03)"},
    {min: 15, max: 39.9, color: "rgba(244,171,24,.03)"},
    {min: 40, max: 49.9, color: "rgba(234,27,231,.03)"},
    {min: 50, max: 100, color: "rgba(157,0,79,.03)"},
  ]
  worker.onmessage = function (e) {
    const {time, grid, data, index} = e.data;
    const canvasDom = document.createElement('canvas');
    canvasDom.className = `heatMap-box-${time}`
    kriging.plot(canvasDom, grid, [112.82999504183772, 113.19000495776265], [25.50999504223745, 25.950004957462852], colors)
    krigingModel.dataMap.set(time, {grid, data, canvasDom, time: new Date(time).toLocaleString()})
    model.curLoadingIndex = index
    if (krigingModel.loading) {
      model.curTime = model.timeList[0]
      updateGridColor()
      krigingModel.loading = false
      loadingInstance.close()
    }
  };
  worker.onerror = function (error) {
    console.error('Worker error:', error);
  };
}

let timer
const togglePlay = () => {
  model.playFlag = !model.playFlag
  if (model.playFlag) {
    if (model.curTime >= model.timeList[model.timeList.length - 1]) model.curTime = model.timeList[0]
    timer = window.setInterval(() => {
      if (model.curTime >= model.timeList[model.timeList.length - 1]) {
        window.clearInterval(timer)
        return
      }
      model.curTime += 60 * 60 * 1e3 / 60
    }, 1000 / model.speedType / 60)
  } else {
    window.clearInterval(timer)
  }
}

const doPrev = () => {
  model.curTime -= 60 * 60 * 1e3
}

const doNext = () => {
  model.curTime += 60 * 60 * 1e3
}

const gridShowChange = (bool) => {
  pointGridDatasouce.show = bool
}

const speedOptions = [
  {label: "0.5x", value: 0.5},
  {label: "1x", value: 1},
  {label: "2x", value: 2},
  {label: "3x", value: 3},
]
// 地图模块
import {usemapStore} from "@/store/modules/cesiumMap";
import boundaryDataOld from "@/assets/data/boundaryOld20240807.json";
import pointGrid from "@/assets/data/pointGrid20240815.json"
import * as turf from '@turf/turf'
import {kriging} from "@/utils/kriging";

const krigingModel = reactive({
  lonLats: [],// 所有点位坐标
  dataMap: new Map(),// 所有点位坐标相对于的降雨量数值
  loading: true
})
const {dataMap} = toRefs(krigingModel)

const mapStore = usemapStore()
const viewer = mapStore.getCesiumViewer();

let heatmapEntity, ctx, heatDoc, pointGridDatasouce;
const polyTurf = turf.polygon([boundaryDataOld.geometry.coordinates[0][0].map(item => ([item[0], item[1]]))]);

const createHeatMap = async () => {
  heatDoc = document.createElement("canvas");
  heatDoc.className = "heatMap-box"
  ctx = heatDoc.getContext("2d");
  heatDoc.width = 180;
  heatDoc.height = 220;
  heatmapEntity = viewer.entities.add({
    rectangle: {
      coordinates: Cesium.Rectangle.fromDegrees(112.82999504183772, 25.50999504223745, 113.19000495776265, 25.950004957462852),
      material: heatDoc
    },
  });
  pointGridDatasouce = await Cesium.GeoJsonDataSource.load(pointGrid, {
    fill: Cesium.Color.fromCssColorString("transparent"),
    clampToGround: true,
  });
  viewer.dataSources.add(pointGridDatasouce);
  pointGridDatasouce.show = false
}

const updateGridColor = async () => {
  if (moment(model.curTime).minutes() !== 0) return
  const temp = krigingModel.dataMap.get(model.curTime)
  if (!temp) {
    ElMessage.warning(`${moment(model.curTime).format("MM月DD日 HH时")}计算中，请稍候`)
    timer && window.clearInterval(timer)
    model.playFlag = false
    model.curTime -= 60 * 60 * 1e3
    return
  }
  // const params = {
  //   forecastTime: moment(model.curTime).unix(),
  //   forecastStartTime: model.formData.type === 'total' ? moment(model.formData.startTime).unix() : "",
  // }
  pointGridDatasouce.entities.values.forEach((entity, i) => {
    const entityPolygon = turf.polygon([entity.polygon.hierarchy.getValue().positions.map(pos => cartesianToWgs84(pos).slice(0, 2))])
    if (!turf.intersect(entityPolygon, polyTurf)) return
    entity.position = Cesium.Cartesian3.fromDegrees(+entity.properties.longitude.getValue(), +entity.properties.latitude.getValue())
    entity.label = new Cesium.LabelGraphics({
      text: formatToFixed(temp.data?.[i] || 0, 1), // 可以根据需要调整标签的文本
      // text: formatToFixed(i, 1),
      font: '14px PingFangSC-Regular, PingFang SC',
      pixelOffset: new Cesium.Cartesian2(0, 5),
      outlineColor: new Cesium.Color.fromCssColorString('black'),
      horizontalOrigin: Cesium.HorizontalOrigin.CENTER,
      verticalOrigin: Cesium.VerticalOrigin.BOTTOM,
      fillColor: Cesium.Color.fromCssColorString("#fff").withAlpha(1),
      disableDepthTestDistance: Number.POSITIVE_INFINITY,
      heightReference: Cesium.HeightReference.CLAMP_TO_GROUND,
      distanceDisplayCondition: new Cesium.DistanceDisplayCondition(0.0, 1e5)
    });
    entity.polyline = new Cesium.PolylineGraphics({
      positions: entity.polygon.hierarchy.getValue().positions,
      material: Cesium.Color.fromCssColorString("#fff").withAlpha(1),
      width: 1,
      clampToGround: true
    })
  });
  heatmapEntity.rectangle.material = temp.canvasDom
}
</script>

woker.js

/**
 * @author: liuk
 * @date: 2024-08-26
 * @describe: 克里金插值
 */
const longitudes = [...]
const latitudes = [...]

importScripts('/worker/kriging.ts')
self.onmessage = async (e) => {
    const {timeList} = e.data;
    const res = await fetch("/worker/---.json")
    const boundaryDataOld = await res.json();
    const ex = boundaryDataOld.geometry.coordinates[0]
    timeList.forEach(time => {
        const data = new Array(396).fill().map(() => Math.random() * 50);
        const variogram = kriging.train(data, longitudes, latitudes, 'exponential', 0, 100);
        const grid = kriging.grid(ex, variogram, 0.0005)
        self.postMessage({time, grid, data});
    });
};