Cesium空间分析-填挖方计算

Cesium提供的接口很基础，很专业。因此扩展起来非常灵活，但同时也有一定的门槛。而且目前关于Cesium相关的参考资料很少，因此实现功能后准备自己记录一下过程。也希望为后续的开发者们提供一些参考资料。

填挖方计算

填挖方：填方+挖方

过程：绘制一个多边形区域，设置一个填挖基准面，基于这个区域的地形及基准面高度，计算需要填多少立方的土，挖掉多少立方的土。

原理：将多边形划分为无穷多个很小很小的三角体，这一特性可利用Cesium.PolygonGeometry.granularity特性来实现。然后遍历每一个小三角体，判断每一个三角体是挖掉，还是填平。（注意：因为地形表面起伏的原因，体积计算会有误差。所以整体计算结果仅供参考）

效果图中，绿色的是Cesium.Entity中的围墙-wall,实现时参考Cesium沙盒示例。蓝色的是带有extrudedHeight属性的多边形，高度为基准面高度。

注：整体思路来源于另一篇文章Cesium-地形挖方分析，感谢大神提供思路。

关键代码（typeScript）：（该方法不适用与模型上的分析，后续我会介绍模型上填挖方的实现）

    private computeCutAndFillVolume1(positions:Cesium.Cartesian3[]):CutAndFillResult{
        //result是一个结果类，用于存储填挖方结果
        const result = new CutAndFillResult();
        if(!this.viewer.terrainProvider.availability){
            return result ;
        }

        //用来记录整块区域的最小高程
        let minHeight = 15000;

        for(let i = 0; i< positions.length;i ++){
            const cartographic = Cesium.Cartographic.fromCartesian(positions[i]);

            const height = this.viewer.scene.globe.getHeight(cartographic);
            minHeight  = Math.min(minHeight,height as number);
        }
        //由于不好给出基准面高度，所以默认取绘制节点的最低点作为基准。
        this.baseHeight = this.baseHeight ==0 ? minHeight:this.baseHeight;

        //设置颗粒度（来源于大神的思路）
        let granularity = Math.PI/Math.pow(2,11);
        granularity = granularity/64;

        const polygonGeometry = Cesium.PolygonGeometry.fromPositions({
            positions:positions,
            vertexFormat:Cesium.PerInstanceColorAppearance.FLAT_VERTEX_FORMAT,
            granularity:granularity,
        });
        const geom = Cesium.PolygonGeometry.createGeometry(polygonGeometry);
        let totalCutVolume = 0;
        let totalFillVolume = 0;
        let maxHeight = 0;

        let i0,i1,i2;
        let height1,height2,height3;
        let bottomP1,bottomP2,bottomP3;
        const scratchCartesian = new Cesium.Cartesian3();
        let cartographic;
        let bottomArea = 0.0;
        let totalBottomArea = 0.0;
        let subTrianglePositions;

        for(let i = 0; i< (geom as Cesium.Geometry).indices.length;i +=3){
            i0 = geom?.indices[i];
            i1 = geom?.indices[i + 1];
            i2 = geom?.indices[i + 2];

            subTrianglePositions = geom?.attributes.position.values;
            if(subTrianglePositions){
                scratchCartesian.x = subTrianglePositions[i0 * 3]; 
                scratchCartesian.y = subTrianglePositions[i0 * 3 + 1];
                scratchCartesian.z = subTrianglePositions[i0 * 3 + 2];
            }

            cartographic = Cesium.Cartographic.fromCartesian(scratchCartesian);
            height1 = this.viewer.scene.globe.getHeight(cartographic);

            bottomP1 = Cesium.Cartesian3.fromRadians(cartographic.longitude,cartographic.latitude,0);

            maxHeight = Math.max(maxHeight,height1 as number);
            minHeight = Math.min(minHeight,height1 as number);

            if(subTrianglePositions){
                scratchCartesian.x = subTrianglePositions[i1 * 3]; 
                scratchCartesian.y = subTrianglePositions[i1 * 3 + 1];
                scratchCartesian.z = subTrianglePositions[i1 * 3 + 2];
            }

            cartographic = Cesium.Cartographic.fromCartesian(scratchCartesian);
            height2 = this.viewer.scene.globe.getHeight(cartographic);

            bottomP2 = Cesium.Cartesian3.fromRadians(cartographic.longitude,cartographic.latitude,0);

            maxHeight = Math.max(maxHeight,height2 as number);
            minHeight = Math.min(minHeight,height2 as number);

            if(subTrianglePositions){
                scratchCartesian.x = subTrianglePositions[i2 * 3]; 
                scratchCartesian.y = subTrianglePositions[i2 * 3 + 1];
                scratchCartesian.z = subTrianglePositions[i2 * 3 + 2];
            }

            cartographic = Cesium.Cartographic.fromCartesian(scratchCartesian);
            height3 = this.viewer.scene.globe.getHeight(cartographic);

            bottomP3 = Cesium.Cartesian3.fromRadians(cartographic.longitude,cartographic.latitude,0);

            maxHeight = Math.max(maxHeight,height3 as number);
            minHeight = Math.min(minHeight,height3 as number);


            bottomArea = this.computeAreaOfTriangle(bottomP1,bottomP2,bottomP3);
            totalBottomArea += bottomArea;

            //计算三角体的平均高度
            const avgCubeHeight = ((height1 as number) + (height2 as number) + (height3 as number))/3;

            //判断是 填方还是挖方

            //如果三角体低于基准面，则需要填方
            if(avgCubeHeight <= this.baseHeight){
                totalFillVolume += bottomArea * (this.baseHeight - avgCubeHeight);
            }else { //否则需要挖方
                totalCutVolume += bottomArea * (avgCubeHeight - this.baseHeight);
            }
            //totalCutVolume += bottomArea * ((height1 as number) - minHeight + (height2 as number) - minHeight + (height3 as number) - minHeight)/3;
        }
        
        result.minHeight = minHeight;
        result.maxHeight = maxHeight;
        result.cutVolume = totalCutVolume;
        result.fillVolume = totalFillVolume;
        result.baseArea = totalBottomArea;

        return result;
    }

export class CutAndFillResult {
    minHeight:number = 0.0;
    maxHeight:number = 0.0;
    cutVolume:number = 0.0;
    fillVolume:number = 0.0;
    baseArea:number = 0.0;

}