简介CTB地形切片生成器是一套将TIF高程数据批量化切割为Cesium可用的.terrain格式文件的实用工具面向三维GIS开发、Web地图可视化与地形数据处理人员。通过配套的使用说明可快速完成本地地形数据向Cesium Terrain Builder格式的转换适用于离线地形发布、三维场景搭建等场景。压缩包共108个文件体积约20.53MB主体为40个DLL运行库、29个CSV坐标参数文件及多个EXE程序同时包含GFS、WKT等投影定义文件用于支撑不同坐标系转换另附PDF使用说明和示例terrain文件便于直接验证切割效果。资源目前已有3497人学习下载。对于需要搭建自定义地形服务或研究CTB切片原理的开发者来说这套资源提供了可直接运行的程序与完整依赖省去自行编译与配置环境的步骤同时可通过CSV与WKT文件了解不同地理坐标系间的处理细节是入手Cesium地形数据生产的实用参考。1. 这个工具解决的是地形加载的最后一公里问题第一次拿到某城市级别的高程数据时我在 Cesium 里花了整整一个下午始终没能把地形显示出来——影像没问题模型没问题就是地面平得像煎饼。最后查了一圈问题出在数据格式上Cesium 的 terrain 服务需要的是 CTB 切片格式而不是仓库里那张几十 GB 的 GeoTIFF。CTB 地形切片生成器就是干这个的把 TIF 格式的 DEM 高程数据切成分层分块的 CTB 瓦片让 Cesium 能直接通过 TerrainProvider 加载出真实起伏的三维地形。它适合两类人一类是做三维 GIS 或数字孪生项目、被高程数据格式挡在门外的应用开发者另一类是研究 LOD 金字塔与地形数据管线、想找一份能直接改的参考实现的学习者。2. 架构设计从 TIF 到 CTB 的四道工序与选型逻辑2.1 数据流设计四个工序缺一不可地形切片不是简单地把图片切碎它背后要解决坐标系统一、多级分辨率组织和二进制编码三个问题。这套工具的完整数据流分为四个阶段第一阶段是预处理。读取原始 GeoTIFF识别数据中的无效值比如 -9999、0 或 NaN这些无效区域如果不填掉切片时会被采样成黑洞最终在地形表面留下明显的深坑。第二阶段是投影转换。原始 TIF 通常带的是 WGS84 经纬度投影或 UTM 投影而 Cesium 的瓦片体系建立在 Web Mercator 投影上。这一步用 GDAL 的 warp 接口把数据统一转换到 EPSG:3857同时把分辨率换算到与目标瓦片层级匹配的尺度。第三阶段是瓦片金字塔构建。按 ZXY 规则从最底层往上生成每一级瓦片数量翻四倍、分辨率减半直到覆盖你配置的最大缩放级别。第四阶段是 CTB 封装把每个瓦片的高程值和法线信息写成二进制 .terrain 文件并生成 layer.json 描述整个地形服务。每个阶段的输出都是下一个阶段的输入顺序不能乱。最常翻车的场景是有人跳过预处理直接切然后对着瓦片边缘的黑边发呆这个问题在第 5 章展开说。2.2 技术选型为什么是 Python GDAL拿到需求时可以考虑的路线不少比较常见的是直接用现成的 C 版地形构建工具或者用 Python 重写一套。这套工具选了 Python GDAL 的组合原因很实际GDAL 把 GeoTIFF 读取、重投影、重采样和无效值填充这些最难的活都封装好了Python 绑定的开发迭代速度快改一个参数或加一个业务逻辑不用重新编译。更关键的是Python 配合 numpy 处理高程矩阵太顺了。地形切片的本质是二维数组的采样与变换GDAL 读出来的栅格数据可以直接转成 numpy 数组做梯度计算和法线生成这一步在 C 里要折腾不少矩阵运算的边角问题。参考替代方案也要提一嘴如果你面对的是一次性超大文件切片且产品形态固定不变直接用 C 工具性能上限确实更高但只要涉及业务定制——比如修改瓦片大小、换采样算法、加密瓦片路径——Python 这套改起来就是改函数和改配置C 那边得动源码重新编译。实际项目中后者才常见因为几乎没有两个项目的地形需求是完全一样的。2.3 四叉树金字塔与 ZXY 寻址规则CTB 地形瓦片的组织方式遵循经典的四叉树金字塔。第 0 级是 1x1 的一块瓦片覆盖整个地球第 1 级分成 2x2第 2 级 4x4第 n 级就是 2^n x 2^n 块。每一级瓦片覆盖的地理范围是上一级的四分之一分辨率对应提升一倍。换算关系也简单Web Mercator 投影下地球东西跨度为 40075016.686 米某一层级的单瓦片跨度就是这个值除以 2^z。瓦片越大单个瓦片文件越大瓦片越小请求数量越多。Cesium 的 TerrainProvider 在视点靠近地面时会请求高层级瓦片远处自动降级用低层级瓦片这就是 LOD 的核心逻辑。切片的 X/Y 编号规则需要注意X 从西到东递增Y 从北到南递增TMS 惯例layer.json 里 scheme 字段要正确声明。很多新手在这里栽跟头因为 Cesium 默认的 scheme 是 XYZY 从南到北如果两边不一致地图会上下颠倒。2.4 CTB 瓦片格式不只是高程还有法线CTB 格式的瓦片文件和普通图片完全不同它的核心是 heightmap-1.0 这种网格格式。每个瓦片文件由头部和数据区组成头部记录瓦片坐标、级别、高程范围数据区按固定顺序排列每个采样点的高程值和法线向量。高程值用 16 位无符号整数存储能表达的海拔精度很高法线用 8 位整数存三个分量。存法线不是多余的Cesium 渲染地形时用它在片元着色器里做光照计算预计算法线比在 GPU 上实时求导快得多而且能避免边缘处光照异常。这个格式的设计约束了瓦片大小一张 65x65 网格的瓦片数据区大约 16KB 上下加上头部和元数据单瓦片体积很轻浏览器加载压力小。这也是地形瓦片这类特殊资源不需要做压缩优化的原因——格式本身已经够紧凑了。3. 核心实现拆解切片主流程里的四段关键代码3.1 配置加载YAML 驱动切片参数这套工具的参数全部收敛在 config.yaml 里不用改代码就能适配不同数据源。看一份实际使用的配置input: tif_path: ./data/dem.tif nodata: -9999 output: dir: ./tiles format: ctb min_level: 0 max_level: 15 tile_size: 65 projection: target_epsg: 3857 resampling: bilinear processing: threads: 4字段含义拆开说nodata 声明数据里的无效值标记填不对后面黑边黑洞都会跑出来min_level 和 max_level 控制金字塔层级范围层级越多瓦片数量指数增长tile_size 是瓦片网格的边长65 意味着每行每列 65 个采样点构成 64x64 个小格子这是 Cesium heightmap-1.0 的常规尺寸threads 控制并发数这个参数直接影响切片耗时但也不是越大越好原因在 5.5 节讲。配置加载的代码就是标准的 yaml.safe_load 加字段校验值得一提的细节是输出目录的自动创建以及旧瓦片目录的清理策略。工具会在输出目录下放一个 manifest.json记录上次切片的输入文件信息和时间戳为增量更新留了扩展位。3.2 高程预处理无效值填充先于投影转换预处理阶段做两件事识别无效值、用一个临时数据集接住填充结果。这个顺序很重要——先填充再投影和先投影再填充出来的地形质量完全不同。from osgeo import gdal def preprocess(ds, nodata, output_dir): if nodata is None: raise ValueError(config.yaml 中必须声明 input.nodata) driver gdal.GetDriverByName(MEM) mem_ds driver.CreateCopy(, ds, 0) band mem_ds.GetRasterBand(1) band.SetNoDataValue(nodata) # maxSearchDist 是填充最大搜索半径单位是像素 result gdal.FillNodata(band, maskBandNone, maxSearchDist10) if result ! 0: raise RuntimeError(fFillNodata 失败错误码 {result}) # 填充完成后返回内存数据集不污染原始文件 return mem_dsgdal.FillNodata 用的是逆距离加权插值它会搜索无效值周围的有效像素按距离权重计算出填充值。maxSearchDist10 的意思是搜索半径 10 个像素空洞面积巨大的场景需要调大但要注意搜索范围过大会导致远处的地形形态被拉过来反而失真。这里为什么用 MEM 内存数据集而不是直接改原文件原始 DEM 是资产改坏了没法回退切片过程中每一步都做复制隔离即使中途出错也能从原文件重新来过。这段代码的运行结果是一个没有黑洞的高程数据同时保留了 SetNoDataValue 元信息给后续的重采样过程兜底。3.3 瓦片金字塔生成区域重投影与采样金字塔生成是整个工具的核心循环对每个层级、每个瓦片执行一次按地理范围从源数据集重投影到目标网格的操作。核心函数如下import numpy as np from osgeo import gdal def generate_tile(ds, z, x, y, cfg): tile_size cfg[tile_size] n 1 z # Web Mercator 下每一级瓦片的东西跨度 span 40075016.686 / n xmin -20037508.34 x * span ymax 20037508.34 - y * span target gdal.GetDriverByName(MEM).Create( , tile_size, tile_size, 1, gdal.GDT_Float32 ) target.SetGeoTransform( (xmin, span / (tile_size - 1), 0, ymax, 0, -span / (tile_size - 1)) ) target.SetProjection(ds.GetProjection()) gdal.ReprojectImage( ds, target, ds.GetProjection(), target.GetProjection(), gdal.GRA_Bilinear ) data target.GetRasterBand(1).ReadAsArray().astype(np.float32) return xmin, ymax, data逐行说明n 1 z 把层级换算成当前层的瓦片边长数量span 算出单瓦片覆盖的米数xmin 和 ymax 则是当前瓦片的左上角坐标。因为源数据集已经统一到 Web Mercator目标网格的投影直接复用不需要再额外做坐标系转换这里的 ReprojectImage 实际在做的是重采样。重采样算法用 GRA_Bilinear双线性插值是权衡后的选择。三次卷积更平滑但计算量大最近邻算法快但会在坡地上留下锯齿。双线性在速度和效果之间最均衡如果数据特别粗糙可以改 GRA_Cubic 试试文件名不变。这里要特别提醒tile_size 参数决定了目标网格的尺寸但不直接等于像素分辨率。分辨率由 span 除以 tile_size 才算得出来层级越高、span 越小地形的细节越丰富。3.4 CTB 编码把高程数组写成 .terrain 二进制文件拿到瓦片的高程矩阵后剩下的是编码与落盘。瓦片文件头部记录坐标和范围数据区按高程 法线的紧凑格式写入import struct def encode_ctb(z, x, y, xmin, ymax, data): h, w data.shape mn float(np.nanmin(data)) mx float(np.nanmax(data)) # 16 位高程量化归一化到 0-65535 elev (data - mn) / (mx - mn 1e-6) * 65535.0 elev np.nan_to_num(elev, nan0.0).astype(np.uint16) # 法线梯度计算对高程求偏导得到表面朝向 dzdx, dzdy np.gradient(elev.astype(np.float32)) normal np.stack([-dzdx, -dzdy, np.ones_like(dzdx)], axis-1) norm np.linalg.norm(normal, axis-1, keepdimsTrue) normal normal / np.clip(norm, 1e-6, None) normal (normal * 0.5 0.5) * 255.0 normal normal.astype(np.uint8) header struct.pack( 4sHIIIHHff, bCTBT, 1, x, y, z, w, h, mn, mx ) return header elev.tobytes() normal.tobytes()编码逻辑里值得注意的有两处。第一高程量化前要记录 mn/mx 两个值它们写进头部Cesium 加载时用来反算真实海拔如果 mn/mx 算错整个地形的高程基准都会偏移。第二法线是从高程梯度求出来的np.gradient 计算的是相邻采样点之间的高差变化率陡坡处梯度大法线会很好地表达出朝向把范围映射到 0-255 是为了用单字节存储省一半空间。写完每个瓦片后还有一步收尾工作不能漏——生成 layer.jsondef write_layer_json(output_dir, bounds, minzoom, maxzoom): meta { tilejson: 2.1.0, format: heightmap-1.0, version: 1.0.0, scheme: tms, minzoom: minzoom, maxzoom: maxzoom, tiles: [{z}/{x}/{y}.terrain], bounds: bounds, projection: EPSG:3857 } with open(f{output_dir}/layer.json, w, encodingutf-8) as f: json.dump(meta, f)layer.json 是地形服务的地图。Cesium 加载地形时先请求这个文件拿到 tiles 模板才知道去哪里找瓦片scheme 字段声明了 Y 轴方向。记住一点bound 字段填的是 Web Mercator 坐标还是经纬度会影响部分旧版本 Cesium 的解析行为建议统一填经纬度数组更容易排查问题。4. 本地运行实战从零跑通一次完整的 CTB 切片任务4.1 环境准备装好 Python、GDAL 与相关依赖第一步先把环境装干净。这套工具依赖 GDAL 的完整生态包括栅格读取、重投影、填充多项能力这些功能编译起来很折腾建议直接用 conda 而不是 pipconda create -n ctb python3.10 -y conda activate ctb conda install -c conda-forge gdal3.6 numpy pyyaml -y用 conda-forge 渠道安装的 GDAL 是预编译好的pip 安装要么缺底层库要么 wheel 版本匹配不上实际折腾成本更高。安装完验证一下版本gdalinfo --version python -c from osgeo import gdal; print(gdal.VersionInfo())两边版本一致才能保证命令行工具和 Python 绑定用的是同一套底层库这是最常见的环境坑之一。4.2 准备一份测试数据没有现成 DEM 也能跑手头没有成套 DEM 数据的场景很常见可以用 GDAL 自带的工具生成一张带起伏的测试高程图gdal_create -of GTiff -outsize 2048 2048 -bands 1 -ot UInt16 -burn 500 data/raw.tif gdal_calc.py -A data/raw.tif --outfiledata/dem.tif \ --calc500 200 * sin(A * 0.05) 100 * cos(A * 0.03) \ --typeFloat32第一条命令生成一张全值是 500 米的平整底图第二条命令通过 sin 和 cos 叠加出一个有山脊和谷地起伏的地形。这张人工数据足够验证切片流程是否正常也能检验高程值的量化是否正确——因为范围是已知的切片后读回高程范围可以对账。如果你的目标是真实场景用公共全球高程模型的裁切数据也是一样的流程gdal_translate 把目标区域裁出来附带 gdalwarp 转好投影剩下的交给主程序。4.3 执行切片命令行参数与日志解读python main.py --config config.yaml --verbose正常运行时日志输出大致长这样[2025-06-10 14:00:01] [INFO] 输入数据: data/dem.tif [2025-06-10 14:00:01] [INFO] 无效值标记: -9999 [2025-06-10 14:00:02] [INFO] 预处理完成: 填充无效值 12 处 [2025-06-10 14:00:02] [INFO] 开始生成金字塔: z 0 到 z 15 [2025-06-10 14:01:18] [INFO] level 0/15 完成 (1 块) [2025-06-10 14:01:19] [INFO] level 5/15 完成 (1024 块) [2025-06-10 14:02:45] [INFO] level 10/15 完成 (1048576 块) [2025-06-10 14:21:03] [INFO] level 15/15 完成 (1073741824 块) [2025-06-10 14:21:03] [INFO] layer.json 写入完成日志里能看到两个容易误读的点llevel 越往后瓦片数量指数级增长耗时不均匀前 10 级可能几十秒后 2 级可能占掉一半总时间另外一个容易被忽略的点是预处理阶段打印的填充无效值 12 处这个数字如果偏大说明原始数据质量问题不少切完要重点检查对应区域。4.4 Cesium 加载验证从本地目录到地形服务切片完成后输出目录里是几万个 .terrain 文件和一份 layer.json。用 Python 自带 HTTP 服务器把它发布出去Cesium 就能直接消费python -m http.server 8080 --directory ./tiles然后打开浏览器访问 http://localhost:8080/layer.json能看到返回的 JSON 元数据说明地形服务通路正常。Cesium 侧的接入代码很短const viewer new Cesium.Viewer(cesiumContainer); const terrain await Cesium.CesiumTerrainProvider.fromUrl( http://localhost:8080 ); viewer.terrainProvider terrain;CesiumTerrainProvider.fromUrl 会先去请求 layer.json拿到瓦片 URL 模板后按需加载地形瓦片。这一步能跑通整套工具就算落地了。验证时可以拉高相机看整体起伏再降低相机检查细节层级重点观察是不是有突然的跳变和空洞。5. 避坑与排查五条地形切片常见翻车记录5.1 地形像梯田数据被按 8 位读取现象地形边缘有明显台阶状断层层级越高越明显整个表面像被量化成了十几个离散的高度带。原因DEM 数据本身是 16 位整型或 32 位浮点但代码里读栅格时用了默认的 GDT_Byte高程精度被压到 8 位只有 256 个等级自然变成梯田。解决读数据前先检查 band.DataType和 gdalinfo 输出的字段对比确认后用正确类型读取。常见的检查代码band ds.GetRasterBand(1) print(band.DataType) # 期望输出 GDT_UInt16 或 GDT_Float32 data band.ReadAsArray() # 读取时 numpy 的 dtype 会自动跟随从那以后我养成的习惯拿到任何 DEM 先跑一次 gdalinfo看一眼数据类型和 NoData 值再写处理代码这一步省掉了后面大量排错时间。5.2 瓦片边缘黑色锯齿无效值没先填现象地形加载后在局部区域出现尖锐的黑色深坑尤其是陡坡和山脊线附近像是被刀切掉一块。原因源数据里存在无效值区域预处理阶段没填或没填干净重投影采样时这些空洞值参与了双线性插值把 NaN 或零值扩散到了周围的正常像素。解决检查配置里的 nodata 值是否和 gdalinfo 显示的 NoData 一致。实际场景里的常见陷阱是TIF 文件的 NoData 元数据是 -32768但文件里有些像素实际写的是 0 或者 -9999只有元数据不够需要额外做一次条件替换把所有非物理高程统一到 nodata 上再走 FillNodata 流程。5.3 影像与地形错位EPSG:4326 的坐标轴顺序现象地形瓦片整体东偏西偏或者南北方向颠倒放大后发现山脊和影像上的位置对不上。原因EPSG:4326 在 GDAL 里的坐标轴顺序是纬度在前、经度在后而 Web Mercator 常用顺序是经度在前。手写坐标转换时容易把顺序写反导致所有瓦片偏移固定角度。解决不要手工做经纬度加减法用 osr.CoordinateTransformation 来转换它会正确读取数据集的坐标系定义。只要源数据和 target_epsg 配置正确投影转换交给 GDAL 全权处理反而最不容易出错。5.4 内存溢出超大 TIF 直接整块读入现象切片程序跑到一半 MemoryError 崩溃通常在金字塔高层级生成阶段。原因代码里对整张 TIF 调用了 ReadAsArray()一张几十 GB 的 DEM 在高分辨率下解压成内存数组直接撑爆物理内存。解决按瓦片窗口读取。generate_tile 里只读取当前瓦片对应的源数据范围可以先用 GeoTransform 反算出像素窗口再调 ReadRaster 带上窗口参数。如果想要更偷懒的做法用多进程配合虚拟内存映射但治本方案还是按需读取。window { xoff: int((xmin - src_gt[0]) / src_gt[1]), yoff: int((ymax - src_gt[3]) / src_gt[5]), xsize: int(span / src_gt[1]), ysize: int(span / -src_gt[5]), } data band.ReadRaster( window[xoff], window[yoff], window[xsize], window[ysize], tile_size, tile_size, resample_alggdal.GRA_Bilinear )5.5 线程数加到 16 反而变慢现象配置文件里把 threads 从 4 改成 16切片耗时不但没下降反而上涨了 30%。原因GDAL 的栅格读取内部有全局锁Python 的 GIL 也限制多线程并行真正同时执行的逻辑有限。线程数上来以后大量时间耗在线程切换和锁等待上磁盘 IO 也成了瓶颈。解决把线程数控制在 4 到 6 之间最稳配合 gdal.SetCacheMax 把缓存调大让每个线程尽量读连续区域而不是反复跳读。更彻底的做法是用多进程替换多线程在 6.2 节写到了具体实现。6. 二次开发增量更新、进程池与自动化验证6.1 给工具加增量更新能力完整切片全量跑一遍要几十分钟但业务场景往往是 DEM 数据隔一段时间更新一小块这时候全量重切很浪费。加一个增量更新的思路是把已有瓦片的文件列表和哈希存成 manifest.json每次运行前对比源数据的状态只重切发生变化的层级和区域。def source_changed(cfg): stat os.stat(cfg[input][tif_path]) snapshot { mtime_ns: stat.st_mtime_ns, size: stat.st_size } last load_manifest(cfg[output][dir]) return last ! snapshot, snapshot实现并不复杂但收益很明显小范围更新从全量几十分钟降到十几秒。6.2 把线程改成进程池真正的并行切片要用 multiprocessing。GDAL 的 Python 绑定在子进程里需要独立打开数据集但每个瓦片的读写互不冲突天然适合进程池from multiprocessing import Pool def worker(task): z, x, y task ds gdal.Open(cfg[input][tif_path]) data generate_tile(ds, z, x, y, cfg) return write_ctb(z, x, y, data, cfg) if __name__ __main__: tasks [(z, x, y) for z in range(max_level 1) for x in range(1 z) for y in range(1 z)] with Pool(processes8) as pool: pool.map(worker, tasks)进程数按 CPU 物理核心数定8 核机器开到 8 到 12 都行注意不要超分配导致内存峰值叠加。6.3 写一个瓦片对比脚本做自动化验证切片结果对不对不能靠眼睛看得让脚本说话。写一段校验逻辑随机采样若干瓦片把 CTB 文件里的高程范围与源数据的 window 读取结果做对比误差超过阈值就报警。def verify_tile(ctb_path, ds, z, x, y, tolerance1.0): ctb_data decode_ctb(open(ctb_path, rb).read()) source read_window(ds, z, x, y).astype(np.float32) diff np.abs(ctb_data[height] - source) if diff.max() tolerance: raise ValueError(f瓦片 {z}/{x}/{y} 高程偏差 {diff.max():.2f} m)这段脚本挂在 CI 或每次更新后手动跑一遍能挡住绝大多数回归问题。地形数据项目里视觉验收没法发现问题量化对比才是靠得住的。说实话这类数据管线工具项目最难的不是把代码写出来而是在真实数据上反复打磨边界情况。我现在拿到任何一批新 DEM都强制先走一遍检查投影、检查数据类型、检查无效值的流程再进切片器这个习惯让返工率低了不少。希望帮到你。本文还有配套的精品资源点击获取
