第09章 - 地形数据处理

第09章：地形数据处理

9.1 地形概述

9.1.1 地形系统架构

CesiumJS 的地形系统提供了高精度的全球高程数据支持：

┌─────────────────────────────────────────────────────────────────┐
│                      地形系统架构                                 │
├─────────────────────────────────────────────────────────────────┤
│                                                                  │
│  TerrainProvider (地形提供者)                                    │
│  ├── 请求地形瓦片                                               │
│  ├── 解码高程数据                                               │
│  └── 提供元数据                                                 │
│                                                                  │
│  支持的格式：                                                    │
│  ├── Quantized Mesh   - Cesium 原生格式，高效压缩               │
│  ├── Heightmap        - 高程图格式                              │
│  └── Google 3D Tiles  - Google 真实感地形                       │
│                                                                  │
│  特性：                                                          │
│  ├── LOD 层次细节                                               │
│  ├── 水体遮罩 (Water Mask)                                      │
│  ├── 顶点法线 (Vertex Normals)                                  │
│  └── 地形遮挡 (Terrain Occlusion)                               │
│                                                                  │
└─────────────────────────────────────────────────────────────────┘

9.1.2 地形提供者类型

类型	Provider 类	描述
Cesium World Terrain	CesiumTerrainProvider	Cesium 官方全球地形
ArcGIS Terrain	ArcGISTiledElevationTerrainProvider	Esri 高程服务
Google Terrain	GoogleEarthEnterpriseTerrainProvider	Google 企业版地形
VR World Terrain	VRTheWorldTerrainProvider	VR 世界地形
Ellipsoid	EllipsoidTerrainProvider	无地形（椭球面）
Custom	CustomHeightmapTerrainProvider	自定义高程图

9.2 内置地形服务

9.2.1 Cesium World Terrain

// 基本使用
viewer.terrainProvider = Cesium.createWorldTerrain();

// 完整配置
viewer.terrainProvider = Cesium.createWorldTerrain({
    requestWaterMask: true,      // 水体遮罩（海洋、湖泊）
    requestVertexNormals: true   // 顶点法线（光照效果）
});

// 使用 async/await
const worldTerrain = await Cesium.createWorldTerrainAsync({
    requestWaterMask: true,
    requestVertexNormals: true
});
viewer.terrainProvider = worldTerrain;

// 从 Ion 资源创建
const ionTerrain = await Cesium.CesiumTerrainProvider.fromIonAssetId(1, {
    requestWaterMask: true,
    requestVertexNormals: true
});
viewer.terrainProvider = ionTerrain;

9.2.2 ArcGIS 高程服务

// ArcGIS World Elevation
const arcgisTerrain = await Cesium.ArcGISTiledElevationTerrainProvider.fromUrl(
    'https://elevation3d.arcgis.com/arcgis/rest/services/WorldElevation3D/Terrain3D/ImageServer'
);
viewer.terrainProvider = arcgisTerrain;

// 自定义 ArcGIS 高程服务
const customArcgisTerrain = await Cesium.ArcGISTiledElevationTerrainProvider.fromUrl(
    'https://your-arcgis-server/arcgis/rest/services/Elevation/ImageServer',
    {
        token: 'YOUR_ARCGIS_TOKEN'  // 如果需要认证
    }
);

9.2.3 椭球体地形（无高程）

// 无地形（平面地球）
viewer.terrainProvider = new Cesium.EllipsoidTerrainProvider();

// 切换地形
function toggleTerrain(enabled) {
    if (enabled) {
        viewer.terrainProvider = Cesium.createWorldTerrain({
            requestWaterMask: true,
            requestVertexNormals: true
        });
    } else {
        viewer.terrainProvider = new Cesium.EllipsoidTerrainProvider();
    }
}

9.3 自定义地形服务

9.3.1 CesiumTerrainProvider 配置

// 自定义 Quantized Mesh 地形服务
const customTerrain = await Cesium.CesiumTerrainProvider.fromUrl(
    'https://your-terrain-server/terrain/',
    {
        requestWaterMask: true,
        requestVertexNormals: true,
        requestMetadata: true,
        credit: 'Custom Terrain Data'
    }
);
viewer.terrainProvider = customTerrain;

9.3.2 自定义高程图地形

// 使用自定义高程图
class CustomHeightmapProvider {
    constructor(options) {
        this._tilingScheme = new Cesium.GeographicTilingScheme({
            numberOfLevelZeroTilesX: 2,
            numberOfLevelZeroTilesY: 1
        });
        this._levelZeroMaximumGeometricError = 
            Cesium.TerrainProvider.getEstimatedLevelZeroGeometricErrorForAHeightmap(
                this._tilingScheme.ellipsoid,
                256,
                this._tilingScheme.getNumberOfXTilesAtLevel(0)
            );
        this._ready = true;
        this._errorEvent = new Cesium.Event();
        this._heightmapWidth = 256;
        this._heightmapHeight = 256;
        this._baseUrl = options.url;
    }
    
    get tilingScheme() { return this._tilingScheme; }
    get ready() { return this._ready; }
    get errorEvent() { return this._errorEvent; }
    get hasWaterMask() { return false; }
    get hasVertexNormals() { return false; }
    get availability() { return undefined; }
    
    getLevelMaximumGeometricError(level) {
        return this._levelZeroMaximumGeometricError / (1 << level);
    }
    
    requestTileGeometry(x, y, level, request) {
        const url = `${this._baseUrl}/${level}/${x}/${y}.terrain`;
        
        return Cesium.Resource.fetchArrayBuffer({ url: url })
            .then(buffer => {
                return new Cesium.HeightmapTerrainData({
                    buffer: new Uint16Array(buffer),
                    width: this._heightmapWidth,
                    height: this._heightmapHeight,
                    childTileMask: 15,  // 所有子瓦片都存在
                    structure: {
                        heightScale: 1.0,
                        heightOffset: 0.0,
                        elementsPerHeight: 1,
                        stride: 1,
                        elementMultiplier: 256.0,
                        isBigEndian: false
                    }
                });
            });
    }
    
    getTileDataAvailable(x, y, level) {
        return level <= 15;  // 最大级别
    }
}

// 使用自定义高程图提供者
const customHeightmap = new CustomHeightmapProvider({
    url: 'https://your-server/heightmaps'
});
viewer.terrainProvider = customHeightmap;

9.4 地形配置与控制

9.4.1 Globe 地形设置

const globe = viewer.scene.globe;

// ===== 地形深度测试 =====
globe.depthTestAgainstTerrain = true;   // 启用地形遮挡

// ===== 地形夸张 =====
globe.terrainExaggeration = 2.0;         // 地形高度夸张倍数
globe.terrainExaggerationRelativeHeight = 0;  // 相对高度

// ===== 地形瓦片设置 =====
globe.maximumScreenSpaceError = 2;       // 最大屏幕空间误差（越小越精细）
globe.tileCacheSize = 100;               // 瓦片缓存大小
globe.preloadSiblings = false;           // 预加载相邻瓦片
globe.preloadAncestors = true;           // 预加载祖先瓦片

// ===== 地形加载设置 =====
globe.loadingDescendantLimit = 20;       // 加载后代限制
globe.fillHighlightColor = undefined;    // 填充高亮颜色

// ===== 光照设置 =====
globe.enableLighting = true;             // 启用光照
globe.dynamicAtmosphereLighting = true;  // 动态大气光照
globe.dynamicAtmosphereLightingFromSun = false;

// ===== 水体效果 =====
globe.showWaterEffect = true;            // 显示水体效果（需要水体遮罩）

// ===== 地球颜色 =====
globe.baseColor = Cesium.Color.BLUE;     // 基础颜色（无数据区域）
globe.showGroundAtmosphere = true;       // 显示地面大气

// ===== 地下可见性 =====
globe.undergroundColor = Cesium.Color.BLACK;
globe.undergroundColorAlphaByDistance = new Cesium.NearFarScalar(1000, 0.0, 10000, 1.0);

9.4.2 地形阴影

// 启用阴影
viewer.shadows = true;
viewer.terrainShadows = Cesium.ShadowMode.ENABLED;

// 阴影模式选项
// Cesium.ShadowMode.DISABLED      - 禁用
// Cesium.ShadowMode.ENABLED       - 启用（投射和接收）
// Cesium.ShadowMode.CAST_ONLY     - 仅投射阴影
// Cesium.ShadowMode.RECEIVE_ONLY  - 仅接收阴影

// 配置阴影贴图
const shadowMap = viewer.scene.shadowMap;
shadowMap.enabled = true;
shadowMap.size = 2048;           // 阴影贴图大小
shadowMap.softShadows = true;    // 软阴影
shadowMap.darkness = 0.3;        // 阴影暗度
shadowMap.maximumDistance = 10000;  // 最大阴影距离

9.5 地形采样与分析

9.5.1 高程采样

// 单点高程采样
async function sampleHeight(longitude, latitude) {
    const positions = [Cesium.Cartographic.fromDegrees(longitude, latitude)];
    
    const updatedPositions = await Cesium.sampleTerrainMostDetailed(
        viewer.terrainProvider,
        positions
    );
    
    return updatedPositions[0].height;
}

// 使用
const height = await sampleHeight(116.4, 39.9);
console.log('高程:', height, '米');

// 多点高程采样
async function sampleHeights(coordinates) {
    const positions = coordinates.map(coord => 
        Cesium.Cartographic.fromDegrees(coord[0], coord[1])
    );
    
    const updatedPositions = await Cesium.sampleTerrainMostDetailed(
        viewer.terrainProvider,
        positions
    );
    
    return updatedPositions.map(pos => pos.height);
}

// 使用
const heights = await sampleHeights([
    [116.0, 39.0],
    [117.0, 39.5],
    [118.0, 40.0]
]);
console.log('高程数组:', heights);

// 指定级别采样（更快但精度较低）
async function sampleHeightAtLevel(longitude, latitude, level = 11) {
    const positions = [Cesium.Cartographic.fromDegrees(longitude, latitude)];
    
    const updatedPositions = await Cesium.sampleTerrain(
        viewer.terrainProvider,
        level,
        positions
    );
    
    return updatedPositions[0].height;
}

9.5.2 地形剖面分析

// 地形剖面线采样
async function getTerrainProfile(startLon, startLat, endLon, endLat, samples = 100) {
    // 生成采样点
    const positions = [];
    for (let i = 0; i <= samples; i++) {
        const fraction = i / samples;
        const lon = startLon + (endLon - startLon) * fraction;
        const lat = startLat + (endLat - startLat) * fraction;
        positions.push(Cesium.Cartographic.fromDegrees(lon, lat));
    }
    
    // 采样高程
    const sampledPositions = await Cesium.sampleTerrainMostDetailed(
        viewer.terrainProvider,
        positions
    );
    
    // 计算距离
    const startCartesian = Cesium.Cartesian3.fromDegrees(startLon, startLat);
    const endCartesian = Cesium.Cartesian3.fromDegrees(endLon, endLat);
    const totalDistance = Cesium.Cartesian3.distance(startCartesian, endCartesian);
    
    // 返回剖面数据
    return sampledPositions.map((pos, index) => ({
        distance: (index / samples) * totalDistance,
        height: pos.height,
        longitude: Cesium.Math.toDegrees(pos.longitude),
        latitude: Cesium.Math.toDegrees(pos.latitude)
    }));
}

// 使用
const profile = await getTerrainProfile(116.0, 39.0, 117.0, 40.0, 50);
console.log('剖面数据:', profile);

// 可视化剖面线
function visualizeProfile(viewer, profileData) {
    const positions = profileData.map(p => 
        Cesium.Cartesian3.fromDegrees(p.longitude, p.latitude, p.height)
    );
    
    viewer.entities.add({
        name: '地形剖面线',
        polyline: {
            positions: positions,
            width: 3,
            material: Cesium.Color.RED,
            clampToGround: false
        }
    });
    
    // 标注高低点
    const minHeight = Math.min(...profileData.map(p => p.height));
    const maxHeight = Math.max(...profileData.map(p => p.height));
    const minPoint = profileData.find(p => p.height === minHeight);
    const maxPoint = profileData.find(p => p.height === maxHeight);
    
    viewer.entities.add({
        name: '最低点',
        position: Cesium.Cartesian3.fromDegrees(minPoint.longitude, minPoint.latitude, minPoint.height),
        point: { pixelSize: 10, color: Cesium.Color.BLUE },
        label: { text: `最低: ${minHeight.toFixed(0)}m`, font: '12px sans-serif' }
    });
    
    viewer.entities.add({
        name: '最高点',
        position: Cesium.Cartesian3.fromDegrees(maxPoint.longitude, maxPoint.latitude, maxPoint.height),
        point: { pixelSize: 10, color: Cesium.Color.RED },
        label: { text: `最高: ${maxHeight.toFixed(0)}m`, font: '12px sans-serif' }
    });
}

9.5.3 地形可视域分析

// 简单的可视域分析
async function viewshedAnalysis(observerLon, observerLat, observerHeight, radius, resolution = 36) {
    // 观察点位置
    const observerCartographic = Cesium.Cartographic.fromDegrees(
        observerLon, observerLat
    );
    
    // 采样观察点高程
    await Cesium.sampleTerrainMostDetailed(viewer.terrainProvider, [observerCartographic]);
    const observerElevation = observerCartographic.height + observerHeight;
    
    const visiblePoints = [];
    const blockedPoints = [];
    
    // 按方位角射线检测
    for (let angle = 0; angle < 360; angle += (360 / resolution)) {
        const radAngle = Cesium.Math.toRadians(angle);
        
        // 沿射线采样
        const rayPoints = [];
        for (let dist = 100; dist <= radius; dist += 100) {
            const offsetLon = observerLon + (dist / 111000) * Math.sin(radAngle);
            const offsetLat = observerLat + (dist / 111000) * Math.cos(radAngle);
            rayPoints.push(Cesium.Cartographic.fromDegrees(offsetLon, offsetLat));
        }
        
        // 采样射线点高程
        await Cesium.sampleTerrainMostDetailed(viewer.terrainProvider, rayPoints);
        
        // 检查可视性
        let blocked = false;
        let maxAngle = -Infinity;
        
        for (const point of rayPoints) {
            const distance = Cesium.Cartesian3.distance(
                Cesium.Cartesian3.fromRadians(observerCartographic.longitude, observerCartographic.latitude, observerElevation),
                Cesium.Cartesian3.fromRadians(point.longitude, point.latitude, point.height)
            );
            
            const elevationAngle = Math.atan2(
                point.height - observerElevation,
                distance
            );
            
            if (elevationAngle > maxAngle) {
                maxAngle = elevationAngle;
                if (!blocked) {
                    visiblePoints.push({
                        longitude: Cesium.Math.toDegrees(point.longitude),
                        latitude: Cesium.Math.toDegrees(point.latitude),
                        height: point.height
                    });
                }
            } else {
                blocked = true;
                blockedPoints.push({
                    longitude: Cesium.Math.toDegrees(point.longitude),
                    latitude: Cesium.Math.toDegrees(point.latitude),
                    height: point.height
                });
            }
        }
    }
    
    return { visiblePoints, blockedPoints };
}

9.6 地形裁剪

9.6.1 裁剪平面

// 地形裁剪（挖掘效果）
const clippingPlanes = new Cesium.ClippingPlaneCollection({
    planes: [
        new Cesium.ClippingPlane(new Cesium.Cartesian3(1, 0, 0), 0),
        new Cesium.ClippingPlane(new Cesium.Cartesian3(-1, 0, 0), 0),
        new Cesium.ClippingPlane(new Cesium.Cartesian3(0, 1, 0), 0),
        new Cesium.ClippingPlane(new Cesium.Cartesian3(0, -1, 0), 0)
    ],
    edgeWidth: 1.0,
    edgeColor: Cesium.Color.WHITE,
    enabled: true
});

viewer.scene.globe.clippingPlanes = clippingPlanes;

// 设置裁剪区域
function setClippingRegion(centerLon, centerLat, width, height) {
    const center = Cesium.Cartesian3.fromDegrees(centerLon, centerLat);
    const transform = Cesium.Transforms.eastNorthUpToFixedFrame(center);
    
    clippingPlanes.modelMatrix = transform;
    
    const halfWidth = width / 2;
    const halfHeight = height / 2;
    
    clippingPlanes.removeAll();
    clippingPlanes.add(new Cesium.ClippingPlane(new Cesium.Cartesian3(1, 0, 0), halfWidth));
    clippingPlanes.add(new Cesium.ClippingPlane(new Cesium.Cartesian3(-1, 0, 0), halfWidth));
    clippingPlanes.add(new Cesium.ClippingPlane(new Cesium.Cartesian3(0, 1, 0), halfHeight));
    clippingPlanes.add(new Cesium.ClippingPlane(new Cesium.Cartesian3(0, -1, 0), halfHeight));
}

// 使用
setClippingRegion(116.4, 39.9, 5000, 5000);  // 5km x 5km 的挖掘区域

9.7 本章小结

本章详细介绍了地形数据处理：

1. 地形架构：TerrainProvider 系统结构
2. 内置服务：Cesium World Terrain、ArcGIS 高程
3. 自定义地形：自定义服务配置、高程图
4. 地形配置：深度测试、夸张、阴影、光照
5. 地形采样：单点/多点高程采样
6. 地形分析：剖面分析、可视域分析
7. 地形裁剪：裁剪平面应用

在下一章中，我们将详细介绍 3D Tiles 大规模数据处理。

9.8 思考与练习

1. 实现地形高度夸张的动态调整功能。
2. 开发地形剖面分析工具，支持绘制剖面图表。
3. 实现基于地形的量测功能（坡度、坡向）。
4. 创建地形裁剪效果，模拟地下管线展示。
5. 对比不同地形源的精度和性能差异。