cesium可视域分析实现方案

通过度娘资料，拼接来的，使用起来性能不是很流畅，有点卡顿，为了应付一下项目交互，望请各位大神指点迷津！

 <!-- //可视域分析 -->
      <canvas id="canvasMap" style="display: none"> </canvas>
      <div class="vs-tools">
      <el-button @click="start">创建可视域分析</el-button>
      <el-button @click="clear">清空</el-button>


<script setup>

import ViewShedAnalysis from '../public/labs/createViewshed.js'
 const start = ()=> {

      vaObj.createViewshed(10);
    };
   const clear=()=> {
    
      vaObj.clearAll();
    };
   const initModel=()=> {
      let tilesetModel = new Cesium.Cesium3DTileset({
        url: `${window.location.origin}/data/model/osgb/tileset.json`,
      });
      viewer.value.scene.primitives.add(tilesetModel);
      viewer.value.flyTo(tilesetModel);
    };
  const getLocation=()=> {
      let handler = new Cesium.ScreenSpaceEventHandler(viewer.canvas);
      handler.setInputAction(function (event) {
        let earthPosition = viewer.value.scene.pickPosition(event.position);
        if (Cesium.defined(earthPosition)) {
          let cartographic = Cesium.Cartographic.fromCartesian(earthPosition);
          let lon = Cesium.Math.toDegrees(cartographic.longitude).toFixed(5);
          let lat = Cesium.Math.toDegrees(cartographic.latitude).toFixed(5);
          let height = cartographic.height.toFixed(2);
          console.log(earthPosition, {
            lon: lon,
            lat: lat,
            height: height,
          });
        }
      }, Cesium.ScreenSpaceEventType.LEFT_CLICK);
    };
  
let vaObj

onMount(()=>{
/所有的API都是从这个实例开始,3D视图
  viewer.value = new Cesium.Viewer("cesiumContainer", {
    imageryProvider: esri,
    infoBox: false,
    selectionIndicator: false,
    shadows: true,
    shouldAnimate: true,
    //  清楚所有控件
    timeline: false, //时间控件
    animation: false, //动画控件
    geocoder: false, //搜索控件
    homeButton: false, //主页控件
    sceneModePicker: false, //投影方式按钮
    baseLayerPicker: false, //图层选择
    navigationHelpButton: false, //帮助手势
    fullscreenButton: false, //全屏按钮
    //地形图层
    terrainProvider: Cesium.createWorldTerrain({
      // 水面特效,消耗电脑gpu
      // requestWaterMask:true,
    }),
  });


//可视域
viewer.value._cesiumWidget._creditContainer.style.display = "none"; // 隐藏版权

vaObj = new ViewShedAnalysis(viewer.value, "canvasMap");

initModel()
getLocation()
})
</script>

找一个工程文件夹，把下面两个js代码分别保存起来，我的命名是createViewshed.js和keiging.js,直接在组件引入使用
注意：需要安装turf=>yarn add  @turf/turf

// Extend the Array class
Array.prototype.max = function () {
    return Math.max.apply(null, this);
};
Array.prototype.min = function () {
    return Math.min.apply(null, this);
};
Array.prototype.mean = function () {
    var i, sum;
    for (i = 0, sum = 0; i < this.length; i++)
        sum += this[i];
    return sum / this.length;
};
Array.prototype.pip = function (x, y) {
    var i, j, c = false;
    for (i = 0, j = this.length - 1; i < this.length; j = i++) {
        if (((this[i][1] > y) != (this[j][1] > y)) &&
            (x < (this[j][0] - this[i][0]) * (y - this[i][1]) / (this[j][1] - this[i][1]) + this[i][0])) {
            c = !c;
        }
    }
    return c;
};

// var kriging = function () {




//     return train, predict, variance, grid, contour, plot;
// }();



var createArrayWithValues = function (value, n) {
    var array = [];
    for (var i = 0; i < n; i++) {
        array.push(value);
    }
    return array;
};

// Matrix algebra
var kriging_matrix_diag = function (c, n) {
    var Z = createArrayWithValues(0, n * n);
    for (let i = 0; i < n; i++) Z[i * n + i] = c;
    return Z;
};
var kriging_matrix_transpose = function (X, n, m) {
    var i, j, Z = Array(m * n);
    for (i = 0; i < n; i++)
        for (j = 0; j < m; j++)
            Z[j * n + i] = X[i * m + j];
    return Z;
};
var kriging_matrix_scale = function (X, c, n, m) {
    var i, j;
    for (i = 0; i < n; i++)
        for (j = 0; j < m; j++)
            X[i * m + j] *= c;
};
var kriging_matrix_add = function (X, Y, n, m) {
    var i, j, Z = Array(n * m);
    for (i = 0; i < n; i++)
        for (j = 0; j < m; j++)
            Z[i * m + j] = X[i * m + j] + Y[i * m + j];
    return Z;
};
// Naive matrix multiplication
var kriging_matrix_multiply = function (X, Y, n, m, p) {
    var i, j, k, Z = Array(n * p);
    for (i = 0; i < n; i++) {
        for (j = 0; j < p; j++) {
            Z[i * p + j] = 0;
            for (k = 0; k < m; k++)
                Z[i * p + j] += X[i * m + k] * Y[k * p + j];
        }
    }
    return Z;
};
// Cholesky decomposition
var kriging_matrix_chol = function (X, n) {
    var i, j, k, sum, p = Array(n);
    for (i = 0; i < n; i++) p[i] = X[i * n + i];
    for (i = 0; i < n; i++) {
        for (j = 0; j < i; j++)
            p[i] -= X[i * n + j] * X[i * n + j];
        if (p[i] <= 0) return false;
        p[i] = Math.sqrt(p[i]);
        for (j = i + 1; j < n; j++) {
            for (k = 0; k < i; k++)
                X[j * n + i] -= X[j * n + k] * X[i * n + k];
            X[j * n + i] /= p[i];
        }
    }
    for (i = 0; i < n; i++) X[i * n + i] = p[i];
    return true;
};
// Inversion of cholesky decomposition
var kriging_matrix_chol2inv = function (X, n) {
    var i, j, k, sum;
    for (i = 0; i < n; i++) {
        X[i * n + i] = 1 / X[i * n + i];
        for (j = i + 1; j < n; j++) {
            sum = 0;
            for (k = i; k < j; k++)
                sum -= X[j * n + k] * X[k * n + i];
            X[j * n + i] = sum / X[j * n + j];
        }
    }
    for (i = 0; i < n; i++)
        for (j = i + 1; j < n; j++)
            X[i * n + j] = 0;
    for (i = 0; i < n; i++) {
        X[i * n + i] *= X[i * n + i];
        for (k = i + 1; k < n; k++)
            X[i * n + i] += X[k * n + i] * X[k * n + i];
        for (j = i + 1; j < n; j++)
            for (k = j; k < n; k++)
                X[i * n + j] += X[k * n + i] * X[k * n + j];
    }
    for (i = 0; i < n; i++)
        for (j = 0; j < i; j++)
            X[i * n + j] = X[j * n + i];

};
// Inversion via gauss-jordan elimination
var kriging_matrix_solve = function (X, n) {
    var m = n;
    var b = Array(n * n);
    var indxc = Array(n);
    var indxr = Array(n);
    var ipiv = Array(n);
    var i, icol, irow, j, k, l, ll;
    var big, dum, pivinv, temp;

    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++) {
            if (i == j) b[i * n + j] = 1;
            else b[i * n + j] = 0;
        }
    for (j = 0; j < n; j++) ipiv[j] = 0;
    for (i = 0; i < n; i++) {
        big = 0;
        for (j = 0; j < n; j++) {
            if (ipiv[j] != 1) {
                for (k = 0; k < n; k++) {
                    if (ipiv[k] == 0) {
                        if (Math.abs(X[j * n + k]) >= big) {
                            big = Math.abs(X[j * n + k]);
                            irow = j;
                            icol = k;
                        }
                    }
                }
            }
        }
        ++(ipiv[icol]);

        if (irow != icol) {
            for (l = 0; l < n; l++) {
                temp = X[irow * n + l];
                X[irow * n + l] = X[icol * n + l];
                X[icol * n + l] = temp;
            }
            for (l = 0; l < m; l++) {
                temp = b[irow * n + l];
                b[irow * n + l] = b[icol * n + l];
                b[icol * n + l] = temp;
            }
        }
        indxr[i] = irow;
        indxc[i] = icol;

        if (X[icol * n + icol] == 0) return false; // Singular

        pivinv = 1 / X[icol * n + icol];
        X[icol * n + icol] = 1;
        for (l = 0; l < n; l++) X[icol * n + l] *= pivinv;
        for (l = 0; l < m; l++) b[icol * n + l] *= pivinv;

        for (ll = 0; ll < n; ll++) {
            if (ll != icol) {
                dum = X[ll * n + icol];
                X[ll * n + icol] = 0;
                for (l = 0; l < n; l++) X[ll * n + l] -= X[icol * n + l] * dum;
                for (l = 0; l < m; l++) b[ll * n + l] -= b[icol * n + l] * dum;
            }
        }
    }
    for (l = (n - 1); l >= 0; l--)
        if (indxr[l] != indxc[l]) {
            for (k = 0; k < n; k++) {
                temp = X[k * n + indxr[l]];
                X[k * n + indxr[l]] = X[k * n + indxc[l]];
                X[k * n + indxc[l]] = temp;
            }
        }

    return true;
};

// Variogram models
var kriging_variogram_gaussian = function (h, nugget, range, sill, A) {
    return nugget + ((sill - nugget) / range) *
        (1.0 - Math.exp(-(1.0 / A) * Math.pow(h / range, 2)));
};
var kriging_variogram_exponential = function (h, nugget, range, sill, A) {
    return nugget + ((sill - nugget) / range) *
        (1.0 - Math.exp(-(1.0 / A) * (h / range)));
};
var kriging_variogram_spherical = function (h, nugget, range, sill, A) {
    if (h > range) return nugget + (sill - nugget) / range;
    return nugget + ((sill - nugget) / range) *
        (1.5 * (h / range) - 0.5 * Math.pow(h / range, 3));
};

// Train using gaussian processes with bayesian priors
var train = function (t, x, y, model, sigma2, alpha) {
    var variogram = {
        t: t,
        x: x,
        y: y,
        nugget: 0.0,
        range: 0.0,
        sill: 0.0,
        A: 1 / 3,
        n: 0
    };
    switch (model) {
        case "gaussian":
            variogram.model = kriging_variogram_gaussian;
            break;
        case "exponential":
            variogram.model = kriging_variogram_exponential;
            break;
        case "spherical":
            variogram.model = kriging_variogram_spherical;
            break;
    };

    // Lag distance/semivariance
    var i, j, k, l, n = t.length;
    var distance = Array((n * n - n) / 2);
    for (i = 0, k = 0; i < n; i++)
        for (j = 0; j < i; j++, k++) {
            distance[k] = Array(2);
            distance[k][0] = Math.pow(
                Math.pow(x[i] - x[j], 2) +
                Math.pow(y[i] - y[j], 2), 0.5);
            distance[k][1] = Math.abs(t[i] - t[j]);
        }
    distance.sort(function (a, b) { return a[0] - b[0]; });
    variogram.range = distance[(n * n - n) / 2 - 1][0];

    // Bin lag distance
    var lags = ((n * n - n) / 2) > 30 ? 30 : (n * n - n) / 2;
    var tolerance = variogram.range / lags;
    var lag = createArrayWithValues(0, lags);
    var semi = createArrayWithValues(0, lags);
    if (lags < 30) {
        for (l = 0; l < lags; l++) {
            lag[l] = distance[l][0];
            semi[l] = distance[l][1];
        }
    }
    else {
        for (i = 0, j = 0, k = 0, l = 0; i < lags && j < ((n * n - n) / 2); i++, k = 0) {
            while (distance[j][0] <= ((i + 1) * tolerance)) {
                lag[l] += distance[j][0];
                semi[l] += distance[j][1];
                j++; k++;
                if (j >= ((n * n - n) / 2)) break;
            }
            if (k > 0) {
                lag[l] /= k;
                semi[l] /= k;
                l++;
            }
        }
        if (l < 2) return variogram; // Error: Not enough points
    }

    // Feature transformation
    n = l;
    variogram.range = lag[n - 1] - lag[0];
    var X = createArrayWithValues(1, 2 * n);
    var Y = Array(n);
    var A = variogram.A;
    for (i = 0; i < n; i++) {
        switch (model) {
            case "gaussian":
                X[i * 2 + 1] = 1.0 - Math.exp(-(1.0 / A) * Math.pow(lag[i] / variogram.range, 2));
                break;
            case "exponential":
                X[i * 2 + 1] = 1.0 - Math.exp(-(1.0 / A) * lag[i] / variogram.range);
                break;
            case "spherical":
                X[i * 2 + 1] = 1.5 * (lag[i] / variogram.range) -
                    0.5 * Math.pow(lag[i] / variogram.range, 3);
                break;
        };
        Y[i] = semi[i];
    }

    // Least squares
    var Xt = kriging_matrix_transpose(X, n, 2);
    var Z = kriging_matrix_multiply(Xt, X, 2, n, 2);
    Z = kriging_matrix_add(Z, kriging_matrix_diag(1 / alpha, 2), 2, 2);
    var cloneZ = Z.slice(0);
    if (kriging_matrix_chol(Z, 2))
        kriging_matrix_chol2inv(Z, 2);
    else {
        kriging_matrix_solve(cloneZ, 2);
        Z = cloneZ;
    }
    var W = kriging_matrix_multiply(kriging_matrix_multiply(Z, Xt, 2, 2, n), Y, 2, n, 1);

    // Variogram parameters
    variogram.nugget = W[0];
    variogram.sill = W[1] * variogram.range + variogram.nugget;
    variogram.n = x.length;

    // Gram matrix with prior
    n = x.length;
    var K = Array(n * n);
    for (i = 0; i < n; i++) {
        for (j = 0; j < i; j++) {
            K[i * n + j] = variogram.model(Math.pow(Math.pow(x[i] - x[j], 2) +
                Math.pow(y[i] - y[j], 2), 0.5),
                variogram.nugget,
                variogram.range,
                variogram.sill,
                variogram.A);
            K[j * n + i] = K[i * n + j];
        }
        K[i * n + i] = variogram.model(0, variogram.nugget,
            variogram.range,
            variogram.sill,
            variogram.A);
    }

    // Inverse penalized Gram matrix projected to target vector
    var C = kriging_matrix_add(K, kriging_matrix_diag(sigma2, n), n, n);
    var cloneC = C.slice(0);
    if (kriging_matrix_chol(C, n))
        kriging_matrix_chol2inv(C, n);
    else {
        kriging_matrix_solve(cloneC, n);
        C = cloneC;
    }

    // Copy unprojected inverted matrix as K
    var K = C.slice(0);
    var M = kriging_matrix_multiply(C, t, n, n, 1);
    variogram.K = K;
    variogram.M = M;

    return variogram;
};

// Model prediction
var predict = function (x, y, variogram) {
    var i, k = Array(variogram.n);
    for (i = 0; i < variogram.n; i++)
        k[i] = variogram.model(Math.pow(Math.pow(x - variogram.x[i], 2) +
            Math.pow(y - variogram.y[i], 2), 0.5),
            variogram.nugget, variogram.range,
            variogram.sill, variogram.A);
    return kriging_matrix_multiply(k, variogram.M, 1, variogram.n, 1)[0];
};
var variance = function (x, y, variogram) {
    var i, k = Array(variogram.n);
    for (i = 0; i < variogram.n; i++)
        k[i] = variogram.model(Math.pow(Math.pow(x - variogram.x[i], 2) +
            Math.pow(y - variogram.y[i], 2), 0.5),
            variogram.nugget, variogram.range,
            variogram.sill, variogram.A);
    return variogram.model(0, variogram.nugget, variogram.range,
        variogram.sill, variogram.A) +
        kriging_matrix_multiply(kriging_matrix_multiply(k, variogram.K,
            1, variogram.n, variogram.n),
            k, 1, variogram.n, 1)[0];
};

// Gridded matrices or contour paths
var grid = function (polygons, variogram, width) {
    var i, j, k, n = polygons.length;
    if (n == 0) return;

    // Boundaries of polygons space
    var xlim = [polygons[0][0][0], polygons[0][0][0]];
    var ylim = [polygons[0][0][1], polygons[0][0][1]];
    for (i = 0; i < n; i++) // Polygons
        for (j = 0; j < polygons[i].length; j++) { // Vertices
            if (polygons[i][j][0] < xlim[0])
                xlim[0] = polygons[i][j][0];
            if (polygons[i][j][0] > xlim[1])
                xlim[1] = polygons[i][j][0];
            if (polygons[i][j][1] < ylim[0])
                ylim[0] = polygons[i][j][1];
            if (polygons[i][j][1] > ylim[1])
                ylim[1] = polygons[i][j][1];
        }

    // Alloc for O(n^2) space
    var xtarget, ytarget;
    var a = Array(2), b = Array(2);
    var lxlim = Array(2); // Local dimensions
    var lylim = Array(2); // Local dimensions
    var x = Math.ceil((xlim[1] - xlim[0]) / width);
    var y = Math.ceil((ylim[1] - ylim[0]) / width);

    var A = Array(x + 1);
    for (i = 0; i <= x; i++) A[i] = Array(y + 1);
    for (i = 0; i < n; i++) {
        // Range for polygons[i]
        lxlim[0] = polygons[i][0][0];
        lxlim[1] = lxlim[0];
        lylim[0] = polygons[i][0][1];
        lylim[1] = lylim[0];
        for (j = 1; j < polygons[i].length; j++) { // Vertices
            if (polygons[i][j][0] < lxlim[0])
                lxlim[0] = polygons[i][j][0];
            if (polygons[i][j][0] > lxlim[1])
                lxlim[1] = polygons[i][j][0];
            if (polygons[i][j][1] < lylim[0])
                lylim[0] = polygons[i][j][1];
            if (polygons[i][j][1] > lylim[1])
                lylim[1] = polygons[i][j][1];
        }

        // Loop through polygon subspace
        a[0] = Math.floor(((lxlim[0] - ((lxlim[0] - xlim[0]) % width)) - xlim[0]) / width);
        a[1] = Math.ceil(((lxlim[1] - ((lxlim[1] - xlim[1]) % width)) - xlim[0]) / width);
        b[0] = Math.floor(((lylim[0] - ((lylim[0] - ylim[0]) % width)) - ylim[0]) / width);
        b[1] = Math.ceil(((lylim[1] - ((lylim[1] - ylim[1]) % width)) - ylim[0]) / width);
        for (j = a[0]; j <= a[1]; j++)
            for (k = b[0]; k <= b[1]; k++) {
                xtarget = xlim[0] + j * width;
                ytarget = ylim[0] + k * width;
                if (polygons[i].pip(xtarget, ytarget))
                    A[j][k] = predict(xtarget,
                        ytarget,
                        variogram);
            }
    }
    A.xlim = xlim;
    A.ylim = ylim;
    A.zlim = [variogram.t.min(), variogram.t.max()];
    A.width = width;
    return A;
};
var contour = function (value, polygons, variogram) {

};

// Plotting on the DOM
var plot = function (canvas, grid, xlim, ylim, colors) {
    // Clear screen
    var ctx = canvas.getContext("2d");
    ctx.clearRect(0, 0, canvas.width, canvas.height);

    // Starting boundaries
    var range = [xlim[1] - xlim[0], ylim[1] - ylim[0], grid.zlim[1] - grid.zlim[0]];
    var i, j, x, y, z;
    var n = grid.length;
    var m = grid[0].length;
    var wx = Math.ceil(grid.width * canvas.width / (xlim[1] - xlim[0]));
    var wy = Math.ceil(grid.width * canvas.height / (ylim[1] - ylim[0]));
    for (i = 0; i < n; i++)
        for (j = 0; j < m; j++) {
            if (grid[i][j] == undefined) continue;
            x = canvas.width * (i * grid.width + grid.xlim[0] - xlim[0]) / range[0];
            y = canvas.height * (1 - (j * grid.width + grid.ylim[0] - ylim[0]) / range[1]);
            z = (grid[i][j] - grid.zlim[0]) / range[2];
            if (z < 0.0) z = 0.0;
            if (z > 1.0) z = 1.0;

            ctx.fillStyle = colors[Math.floor((colors.length - 1) * z)];
            ctx.fillRect(Math.round(x - wx / 2), Math.round(y - wy / 2), wx, wy);
        }

};
// console.log(kriging);
// if (module && module.exports) {
//     // module.exports = kriging;
// }
export { train, grid ,plot};

import * as Cesium from 'cesium';
import { train, grid ,plot} from "./keiging";
console.log(train, grid);
import * as turf from "@turf/turf";
let pos;
let earthPosition;
let ViewShedAnalysis = function (viewer, canvasEleId) {
  // console.log(viewer, canvasEleId);
  if (!viewer) throw new Error("no viewer object!");
  let canvasEle = document.getElementById(canvasEleId);
  if (!canvasEle) throw new Error("the canvas element is not exist");
  this.canvasEle = canvasEle;
  this.viewer = viewer;
  this.handler = undefined;
  this.lightCamera;
  this.pyramid;
  this.frustumPrimitive;
  this.viewershedPolygon;
};
ViewShedAnalysis.prototype = {
  /**
   * 初始化handler
   */
  initHandler() {
    if (this.handler) {
      this.handler.destroy();
      this.handler = undefined;
    }
  },
  /**
   * 开始执行视域分析
   * @param {number} precision 精度，值越大创建耗时越长，建议在10~20之间
   */
  createViewshed: function (precision) {
    let $this = this;
    let scene = $this.viewer.scene;
    $this.initHandler();
    $this.clearAll();
    $this.handler = new Cesium.ScreenSpaceEventHandler($this.viewer.canvas);
    $this.handler.setInputAction((event) => {
      // 禁止地球旋转和缩放，地球的旋转会对鼠标移动监听有影响，所以需要禁止
      scene.screenSpaceCameraController.enableRotate = false;
      scene.screenSpaceCameraController.enableZoom = false;
      scene.globe.depthTestAgainstTerrain = false;
      earthPosition = scene.pickPosition(event.position);
      pos = $this.cartesian3ToDegree(earthPosition);
      $this.handler.setInputAction(function (event) {
        let newPosition = scene.pickPosition(event.endPosition);
        if (Cesium.defined(newPosition)) {
          let pos1 = $this.cartesian3ToDegree(newPosition);
          let distance = Cesium.Cartesian3.distance(newPosition, earthPosition);
          let angle = $this.getAngle(pos[0], pos[1], pos1[0], pos1[1]);
          let pitch = $this.getPitch(earthPosition, newPosition);
          $this.ViewShedOptions = {
            viewPosition: earthPosition, //观测点 笛卡尔坐标
            endPosition: newPosition, //目标点 笛卡尔坐标
            direction: angle, //观测方位角 默认为`0`,范围`0~360`
            pitch: pitch, //俯仰角,radius,默认为`0`
            horizontalViewAngle: 90, //可视域水平夹角,默认为 `90`,范围`0~360`
            verticalViewAngle: 60, //可视域垂直夹角,默认为`60`,范围`0~180`
            visibleAreaColor: Cesium.Color.GREEN, //可视区域颜色,默认为`green`
            invisibleAreaColor: Cesium.Color.RED, //不可见区域颜色,默认为`red`
            visualRange: distance, //距离,单位`米`
          };
          $this.updateViewShed();
        }
      }, Cesium.ScreenSpaceEventType.MOUSE_MOVE);
    }, Cesium.ScreenSpaceEventType.LEFT_DOWN);

    $this.handler.setInputAction(() => {
      $this.initHandler();
      // 开启地球旋转和缩放
      scene.screenSpaceCameraController.enableRotate = true;
      scene.screenSpaceCameraController.enableZoom = true;
      $this.drawViewershed(precision);
    }, Cesium.ScreenSpaceEventType.LEFT_UP);
  },

  ReturnDistance(pos0, pos1) {
    let distance = 0;
    let point1cartographic = Cesium.Cartographic.fromCartesian(pos0);
    let point2cartographic = Cesium.Cartographic.fromCartesian(pos1);
    /**根据经纬度计算出距离**/
    let geodesic = new Cesium.EllipsoidGeodesic();
    geodesic.setEndPoints(point1cartographic, point2cartographic);
    let s = geodesic.surfaceDistance;
    return s;
  },
  getHeight(x, y, objectsToExclude) {
    let endCartographic = Cesium.Cartographic.fromDegrees(x, y);
    let endHeight = this.viewer.scene.sampleHeight(
      endCartographic,
      objectsToExclude
    );
    return endHeight;
  },

  cartesian3ToDegree: function (Cartesian3) {
    // console.log(Cartesian3);
    let _ellipsoid = this.viewer.scene.globe.ellipsoid;
    let _cartographic = _ellipsoid.cartesianToCartographic(Cartesian3);
    let _lat = Cesium.Math.toDegrees(_cartographic.latitude);
    let _lng = Cesium.Math.toDegrees(_cartographic.longitude);
    let _alt = _cartographic.height;
    return [_lng, _lat, _alt];
  },
  getAngle: function (lng1, lat1, lng2, lat2) {
    let dRotateAngle = Math.atan2(Math.abs(lng1 - lng2), Math.abs(lat1 - lat2));
    if (lng2 >= lng1) {
      dRotateAngle = lat2 < lat1 ? Math.PI - dRotateAngle : dRotateAngle;
    } else {
      dRotateAngle =
        lat2 >= lat1 ? 2 * Math.PI - dRotateAngle : Math.PI + dRotateAngle;
    }
    dRotateAngle = (dRotateAngle * 180) / Math.PI;
    return dRotateAngle;
  },

  getPitch(pointA, pointB) {
    let transfrom = Cesium.Transforms.eastNorthUpToFixedFrame(pointA);
    const vector = Cesium.Cartesian3.subtract(
      pointB,
      pointA,
      new Cesium.Cartesian3()
    );
    let direction = Cesium.Matrix4.multiplyByPointAsVector(
      Cesium.Matrix4.inverse(transfrom, transfrom),
      vector,
      vector
    );
    Cesium.Cartesian3.normalize(direction, direction);
    return Cesium.Math.PI_OVER_TWO - Cesium.Math.acosClamped(direction.z);
  },

  updateViewShed: function () {
    this.clear();
    this.setLightCamera();
    this.addVisualPyramid();
    this.createFrustum();
  },
  clear: function () {
    if (this.pyramid) {
      this.viewer.entities.removeById(this.pyramid.id);
      this.pyramid = undefined;
    }
    if (this.frustumPrimitive) {
      this.viewer.scene.primitives.remove(this.frustumPrimitive);
      this.frustumPrimitive = undefined;
    }
    if (this.debugModelMatrixPrimitive) {
      this.viewer.scene.primitives.remove(this.debugModelMatrixPrimitive);
      this.debugModelMatrixPrimitive = undefined;
    }
  },
  clearAll: function () {
    this.clear();
    if (this.viewershedPolygon) {
      this.viewer.scene.primitives.remove(this.viewershedPolygon);
      this.viewershedPolygon = undefined;
    }
  },
  addVisualPyramid: function () {
    let options = this.ViewShedOptions;
    let position = options.viewPosition;
    let visualRange = Number(options.visualRange);
    let transform = Cesium.Transforms.eastNorthUpToFixedFrame(position);
    this.debugModelMatrixPrimitive = this.viewer.scene.primitives.add(
      new Cesium.DebugModelMatrixPrimitive({
        modelMatrix: transform,
        length: 5.0,
      })
    );
    const halfClock = options.horizontalViewAngle / 2;
    const halfCone = options.verticalViewAngle / 2;
    const pitch = Cesium.Math.toDegrees(options.pitch);
    const ellipsoid = new Cesium.EllipsoidGraphics({
      radii: new Cesium.Cartesian3(visualRange, visualRange, visualRange),
      minimumClock: Cesium.Math.toRadians(90 - options.direction - halfClock),
      maximumClock: Cesium.Math.toRadians(90 - options.direction + halfClock),
      minimumCone: Cesium.Math.toRadians(90 - pitch - halfCone),
      maximumCone: Cesium.Math.toRadians(90 - pitch + halfCone),
      fill: false,
      outline: true,
      subdivisions: 256,
      stackPartitions: 64,
      slicePartitions: 64,
      outlineColor: Cesium.Color.YELLOWGREEN.withAlpha(0.5),
    });
    const pyramidEntity = new Cesium.Entity({
      position: position,
      ellipsoid,
    });
    this.pyramid = this.viewer.entities.add(pyramidEntity);
  },
  setLightCamera: function () {
    if (!this.lightCamera) {
      this.lightCamera = new Cesium.Camera(this.viewer.scene);
    }
    let options = this.ViewShedOptions;
    let visualRange = Number(options.visualRange);
    this.lightCamera.position = options.viewPosition;
    this.lightCamera.frustum.near = 0.1;
    this.lightCamera.frustum.far = visualRange;
    const hr = Cesium.Math.toRadians(options.horizontalViewAngle);
    const vr = Cesium.Math.toRadians(options.verticalViewAngle);
    this.lightCamera.frustum.aspectRatio =
      (visualRange * Math.tan(hr / 2) * 2) /
      (visualRange * Math.tan(vr / 2) * 2);
    this.lightCamera.frustum.fov = hr > vr ? hr : vr;
    this.lightCamera.setView({
      destination: options.viewPosition,
      orientation: {
        heading: Cesium.Math.toRadians(options.direction || 0),
        pitch: options.pitch || 0,
        roll: 0,
      },
    });
  },
  createFrustum: function () {
    const scratchRight = new Cesium.Cartesian3();
    const scratchRotation = new Cesium.Matrix3();
    const scratchOrientation = new Cesium.Quaternion();
    const direction = this.lightCamera.directionWC;
    const up = this.lightCamera.upWC;
    let right = this.lightCamera.rightWC;
    right = Cesium.Cartesian3.negate(right, scratchRight);
    let rotation = scratchRotation;
    Cesium.Matrix3.setColumn(rotation, 0, right, rotation);
    Cesium.Matrix3.setColumn(rotation, 1, up, rotation);
    Cesium.Matrix3.setColumn(rotation, 2, direction, rotation);
    let orientation = Cesium.Quaternion.fromRotationMatrix(
      rotation,
      scratchOrientation
    );
    let instanceOutline = new Cesium.GeometryInstance({
      geometry: new Cesium.FrustumOutlineGeometry({
        frustum: this.lightCamera.frustum,
        origin: this.ViewShedOptions.viewPosition,
        orientation: orientation,
      }),
      id: "视椎体轮廓线" + Math.random().toString(36).substr(2),
      attributes: {
        color: Cesium.ColorGeometryInstanceAttribute.fromColor(
          new Cesium.Color(0.0, 1.0, 0.0, 1.0)
        ),
        show: new Cesium.ShowGeometryInstanceAttribute(true),
      },
    });
    this.frustumPrimitive = this.viewer.scene.primitives.add(
      new Cesium.Primitive({
        geometryInstances: instanceOutline,
        appearance: new Cesium.PerInstanceColorAppearance({
          flat: true,
          translucent: false,
          closed: true,
        }),
      })
    );
  },
  createPoint: function (firstPos, secondPos) {
    let entity4FirstPos = new Cesium.Entity({
      name: "firstPos",
      show: true,
      position: firstPos,
      point: {
        show: true,
        pixelSize: 20,
        color: Cesium.Color.RED,
        outlineColor: Cesium.Color.YELLOW,
        outlineWidth: 5,
      },
      description: `
              <p>这是绘制的视椎体起点</p>`,
    });
    this.viewer.entities.add(entity4FirstPos);
    let entity4SecondPos = new Cesium.Entity({
      name: "secondPos",
      show: true,
      position: secondPos,
      point: {
        show: true,
        pixelSize: 30,
        color: Cesium.Color.YELLOW,
        outlineColor: Cesium.Color.RED,
        outlineWidth: 8,
      },
      description: `
              <p>这是绘制的视椎体视角终点</p>`,
    });
    this.viewer.entities.add(entity4SecondPos);
  },

  //绘制可视域
  add(positionArr) {
    let polygon = new Cesium.PolygonGeometry({
      polygonHierarchy: new Cesium.PolygonHierarchy(
        Cesium.Cartesian3.fromDegreesArray(positionArr)
      ),
      height: 0.0,
      extrudedHeight: 0.0,
      vertexFormat: Cesium.PerInstanceColorAppearance.VERTEX_FORMAT,
      stRotation: 0.0, // 纹理的旋转坐标（以弧度为单位）,正旋转是逆时针方向
      ellipsoid: Cesium.Ellipsoid.WGS84,
      granularity: Cesium.Math.RADIANS_PER_DEGREE, // 每个纬度和经度之间的距离（以弧度为单位）,确定缓冲区中的位置数
      perPositionHeight: false, // 每个位置点使用的高度
      closeTop: true,
      closeBottom: true,
      // NONE 与椭圆表面不符的直线;GEODESIC 遵循测地路径;RHUMB	遵循大黄蜂或恶魔般的道路。
      arcType: Cesium.ArcType.GEODESIC, // 多边形边缘线型
    });

    let polygonInstance = new Cesium.GeometryInstance({
      geometry: polygon,
      name: "ViewershedPolygon",
      attributes: {
        color: Cesium.ColorGeometryInstanceAttribute.fromColor(
          Cesium.Color.BLUE.withAlpha(0.6)
        ),
        show: new Cesium.ShowGeometryInstanceAttribute(true), //显示或者隐藏
      },
    });
    this.viewershedPolygon = this.viewer.scene.primitives.add(
      new Cesium.GroundPrimitive({
        geometryInstances: polygonInstance,
        appearance: new Cesium.EllipsoidSurfaceAppearance({
          aboveGround: true,
          material: new Cesium.Material({
            fabric: {
              type: "Image",
              uniforms: {
                image: this.returnImgae(),
              },
            },
          }),
        }),
      })
    );
  },
  drawViewershed(precision) {
    const pos = this.cartesian3ToDegree(this.ViewShedOptions.viewPosition);
    const radius = this.ViewShedOptions.visualRange;
    const direction = this.ViewShedOptions.direction;
    let boundary = this.computeBoundaryOptions(pos, radius, direction);
    const bbox = boundary.bbox;
    let mask = turf.polygon([boundary.boundaryPoints]);
    const dis = this.ViewShedOptions.visualRange / (precision * 1000);
    let gridPoints = turf.pointGrid(bbox, dis, { mask: mask });
    // console.log(this.ViewShedOptions.visualRange, precision, dis);

    let pointsResult = this.createTargetPoints(gridPoints, dis, pos);
    let variogram = train(
      pointsResult.values,
      pointsResult.lngs,
      pointsResult.lats,
      "exponential",
      0,
      100
    );
    console.log(variogram);
    var _grid = grid([boundary.boundaryPoints], variogram, dis / 1000);
    const colors = [
      "#ff000080",
      "#ff000080",
      "#ff000080",
      "#ff000080",
      "#ff000080",
      "#ff000080",
      "#00ff0080",
      "#00ff0080",
      "#00ff0080",
      "#00ff0080",
      "#00ff0080",
      "#00ff0080",
    ];

    this.canvasEle.width = 3840;
    this.canvasEle.height = 2160;
    plot(
      this.canvasEle,
      _grid ,
      [bbox[0], bbox[2]],
      [bbox[1], bbox[3]],
      colors
    );
    this.add(boundary.positionArr);
  },
  computeBoundaryOptions(pos, radius, angle) {
    let Ea = 6378137; //   赤道半径
    let Eb = 6356725; // 极半径
    const lng = pos[0],
      lat = pos[1];
    const bbox = [lng, lat, lng, lat]; //[minX, minY, maxX, maxY]
    let positionArr = [];
    let boundaryPoints = [];
    positionArr.push(lng, lat);
    boundaryPoints.push([lng, lat]);
    //正北是0°
    let start = angle + 45 > 360 ? angle - 45 - 360 : angle - 45;
    let end = start + 90;
    for (let i = start; i <= end; i++) {
      let dx = radius * Math.sin((i * Math.PI) / 180.0);
      let dy = radius * Math.cos((i * Math.PI) / 180.0);
      let ec = Eb + ((Ea - Eb) * (90.0 - lat)) / 90.0;
      let ed = ec * Math.cos((lat * Math.PI) / 180);
      let BJD = lng + ((dx / ed) * 180.0) / Math.PI;
      let BWD = lat + ((dy / ec) * 180.0) / Math.PI;
      positionArr.push(BJD, BWD);
      boundaryPoints.push([BJD, BWD]);
      this.refreshBBox(bbox, BJD, BWD);
    }
    boundaryPoints.push([lng, lat]);
    return {
      positionArr,
      boundaryPoints,
      bbox,
    };
  },
  /**
   * 更新外围矩形 Bbox
   * @param {Array} result 外围矩形Bbox-[minX, minY, maxX, maxY]
   * @param {Number} x 经度
   * @param {Number} y 纬度
   */
  refreshBBox(result, x, y) {
    result[0] = x < result[0] ? x : result[0];
    result[1] = y < result[1] ? y : result[1];
    result[2] = x > result[2] ? x : result[2];
    result[3] = y > result[3] ? y : result[3];
  },
  /**
   * 插值点用射线判断通视性
   * @param {*} gridPoints 网格点
   * @param {*} step 步长，可以理解成是精度
   * @param {*} sourcePos 视域分析起点
   * @returns kriging插值所需的参数对象{ values:[], lngs:[], lats:[]}
   */
  createTargetPoints(gridPoints, step, sourcePos) {
    let positionArr = [];
    let objectsToExclude = [
      this.frustumPrimitive,
      this.pyramid,
      this.debugModelMatrixPrimitive,
    ];
    let values = [],
      lngs = [],
      lats = [];
    let height = this.getHeight(sourcePos[0], sourcePos[1], objectsToExclude);
    positionArr.push({
      x: sourcePos[0],
      y: sourcePos[1],
      z: height,
    });
    let viewPoint = this.ViewShedOptions.viewPosition;
    for (let index = 0; index < gridPoints.features.length; index++) {
      const feature = gridPoints.features[index];
      const coords = feature.geometry.coordinates;
      const x = coords[0],
        y = coords[1];
      let h = this.getHeight(x, y, objectsToExclude);
      let endPoint = Cesium.Cartesian3.fromDegrees(x, y, h);
      let direction = Cesium.Cartesian3.normalize(
        Cesium.Cartesian3.subtract(
          endPoint,
          viewPoint,
          new Cesium.Cartesian3()
        ),
        new Cesium.Cartesian3()
      );
      // 建立射线
      let ray = new Cesium.Ray(viewPoint, direction);
      let result = this.viewer.scene.pickFromRay(ray, objectsToExclude); // 计算交互点，返回第一个
      if (result) {
        let buffer = this.ReturnDistance(endPoint, result.position);
        // let M_color = Cesium.Color.GREEN;
        if (buffer > step) {
          // M_color = Cesium.Color.RED;
          values.push(0);
        } else {
          values.push(1);
        }
        lngs.push(x);
        lats.push(y);
        // this.viewer.entities.add(
        //   new Cesium.Entity({
        //     name: "插值点哦",
        //     show: true,
        //     position: endPoint,
        //     point: {
        //       show: true,
        //       pixelSize: 10,
        //       color: M_color,
        //       outlineWidth: 2,
        //       outlineColor: Cesium.Color.YELLOW,
        //     },
        //   })
        // );
      }
    }
    return {
      values,
      lngs,
      lats,
    };
  },
  /**
   * canvas转image图片
   * @returns base64图片
   */
  returnImgae() {
    return this.canvasEle.toDataURL("image/png");
  },
};

export default ViewShedAnalysis;