Cesium基础-表面面积量算（依地形量算、依模型表面量算）

上一篇Cesium基础-表面距离量算简单介绍了表面距离的计算原理。这一篇想介绍一下表面面积量算的实现过程。

和表面距离一样，本质还是“微积分”的原理。将整块面无限细化为很小的小面，以提升计算精度。

面细化方法，我用了两种方法。如果大家有更好的方法，欢迎评论区留言。

方法一：将多边形面对应屏幕区域，构建泰森多边形（我比较懒，直接用turf.js的API计算的）

  private calculateSurfaceArea(positions:Cartesian2[]):number{
    let result = 0;
    const bounds = getBounds(positions);
    const points = turf.randomPoint(50,{bbox:[bounds[0],bounds[1],bounds[2],bounds[3]]});
    const mainPoly = Cartesian2turfPolygon(positions);
    const voronoiPolygons = turf.voronoi(points,{bbox:[bounds[0],bounds[1],bounds[2],bounds[3]]});
    voronoiPolygons.features.forEach(element => {
        const intersectPoints = intersect(mainPoly,element.geometry)
        result += this.calculateDetailSurfaceArea(intersectPoints);
    });
    return result;
  }


  private calculateDetailSurfaceArea(positions:Cartesian2[]):number{
    let worldPositions:Cesium.Cartesian3[] = [];
    positions.forEach(element => {
      worldPositions.push(pickCartesian(this.viewer,element).cartesian);
    });
    return this.getArea(worldPositions);
  }

private getArea(positions: Cartesian3[]) {
    const x: number[] = [0];
    const y: number[] = [0];
    const geodesic = new Cesium.EllipsoidGeodesic();
    const radiansPerDegree = Math.PI / 180.0; //角度转化为弧度(rad)
    //数组x,y分别按顺序存储各点的横、纵坐标值
    for (let i = 0; i < positions.length - 1; i++) {
      const p1 = positions[i];
      const p2 = positions[i + 1];
      const point1cartographic = Cesium.Cartographic.fromCartesian(p1);
      const point2cartographic = Cesium.Cartographic.fromCartesian(p2);
      geodesic.setEndPoints(point1cartographic, point2cartographic);
      const s = Math.sqrt(Math.pow(geodesic.surfaceDistance, 2) +
        Math.pow(point2cartographic.height - point1cartographic.height, 2));
      // console.log(s, p2.y - p1.y, p2.x - p1.x)
      const lat1 = point2cartographic.latitude * radiansPerDegree;
      const lon1 = point2cartographic.longitude * radiansPerDegree;
      const lat2 = point1cartographic.latitude * radiansPerDegree;
      const lon2 = point1cartographic.longitude * radiansPerDegree;
      let angle = -Math.atan2(
        Math.sin(lon1 - lon2) * Math.cos(lat2),
        Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(lon1 - lon2)
      );
      if (angle < 0) {
        angle += Math.PI * 2.0;
      }
      console.log('角度：' + (angle * 180) / Math.PI);

      y.push(Math.sin(angle) * s + y[i]);
      x.push(Math.cos(angle) * s + x[i]);
    }

    let sum = 0;
    for (let i = 0; i < x.length - 1; i++) {
      sum += x[i] * y[i + 1] - x[i + 1] * y[i];
    }
    // console.log(x, y)

    return Math.abs(sum + x[x.length - 1] * y[0] - x[0] * y[y.length - 1]) / 2;
  }

export function pickCartesian(viewer:Cesium.Viewer,windowPosition:Cesium.Cartesian2):PickResult{
    //根据窗口坐标，从场景的深度缓冲区中拾取相应的位置，返回笛卡尔坐标。
    const cartesianModel = viewer.scene.pickPosition(windowPosition); 
    
    //场景相机向指定的鼠标位置（屏幕坐标）发射射线
    const ray = viewer.camera.getPickRay(windowPosition);
    //获取射线与三维球相交的点（即该鼠标位置对应的三维球坐标点，因为模型不属于球面的物体，所以无法捕捉模型表面）
    const cartesianTerrain = viewer.scene.globe.pick(ray,viewer.scene);

    const result = new PickResult();
    if(typeof(cartesianModel) !== 'undefined' && typeof(cartesianTerrain) !== 'undefined'){
      result.cartesian = cartesianModel || cartesianTerrain;
      result.CartesianModel = cartesianModel;
      result.cartesianTerrain = cartesianTerrain as Cesium.Cartesian3;
      result.windowCoordinates = windowPosition.clone();
      //坐标不一致，证明是模型，采用绝对高度。否则是地形，用贴地模式。
      result.altitudeMode = cartesianModel.z.toFixed(0) !== cartesianTerrain!.z.toFixed(0) ? Cesium.HeightReference.NONE:Cesium.HeightReference.CLAMP_TO_GROUND;
    }
    return result;
  }

export function getBounds(points:Cesium.Cartesian2[]):number[]{
    let bounds: number[] = [];
    let left = Number.MAX_VALUE;
    let right = Number.MIN_VALUE;
    let top = Number.MAX_VALUE;
    let bottom = Number.MIN_VALUE;
    points.forEach(element => {
        left = Math.min(left,element.x);
        right = Math.max(right,element.x);
        top = Math.min(top,element.y);
        bottom = Math.max(bottom,element.y);
    });
    bounds.push(left);
    bounds.push(top);
    bounds.push(right);
    bounds.push(bottom);
    return bounds;
}


export function Cartesian2turfPolygon(positions:Cesium.Cartesian2[]):turf.helpers.Polygon{
    var coordinates:number[][][] = [[]];
    positions.forEach(element => {
        coordinates[0].push([element.x,element.y]);
    });
    coordinates[0].push([positions[0].x,positions[0].y]);
    const polygon =  turf.polygon(coordinates);
    return polygon.geometry;
}

export function intersect(poly1:turf.helpers.Polygon,poly2:turf.helpers.Polygon):Cesium.Cartesian2[]{
    var intersection = turf.intersect(poly1, poly2);
    if(intersection?.geometry !== undefined){
        return turfPloygon2CartesianArr(intersection?.geometry as turf.Polygon);
    }else{
        return [];
    }
}

 方法二：在填挖方计算文章中，我借鉴了大神的思路，发现PolygonGeometry可以设置granularity用以对多边形进行颗粒度细分。默认是三角网。

可以通过计算三角网中每个三角形的面积，进而累加得到表面面积。

该思路我还没有去写代码验证。如果小伙伴们有兴趣，可以验证一下呀！