简介CTB地形切片生成器面向GIS开发与三维可视化场景可将TIF格式的高程或影像数据批量切割为Cesium可直接加载的terrain地形数据免去手工处理大文件切片与格式转换的麻烦。压缩包共包含108个文件大小约20.53MB主体为40个DLL动态库和4个EXE可执行程序并配有29个CSV坐标投影参数表、8个GFS与6个WKT地理信息定义文件以及XML、Json等配置能支撑从读取原始影像到输出地形切片的全流程。目前已有3497人学习下载。借助随包附带的投影参数表、示例地形文件与说明文档可以快速校验输出结果并接入Cesium场景适合需要自建三维地球或地形底图的开发人员参考使用。1. CTB地形切片生成器把单个GeoTIFF变成Cesium秒开的三维地形做过WebGIS的人都懂手里拿着一块十几GB的GeoTIFF高程数据想用Cesium直接加载预览时卡到帧率个位数浏览器直接内存溢出。传统方案是先切DOM再用CesiumTerrainProvider加载但市面上的切片工具要么收费要么配置复杂到需要专门写文档。这个CTB地形切片生成器就是来解决这个痛点的输入一张带地理坐标的TIF高程栅格输出一套按四叉树层级组织的terrain瓦片和一份layer.jsonCesium加载后地形LOD自动切换加载速度从几分钟降到几百毫秒。适合正在做三维地球、智慧城市、水利淹没分析这类项目的开发者也适合刚接触Cesium地形、想把DEM数据规范化处理的新手。整个转换链路不依赖在线服务离线跑完直接放入静态目录即可。2. 整体架构从TIF到terrain tiles的转换链路CTB地形切片生成器不是一个复杂的服务更像一个离线批处理管道。它的输入是GeoTIFF输出是Cesium地形服务要求的目录结构。理解整条链路之前先要把quantized-mesh和LOD金字塔这两个概念弄清楚否则后面看代码会一头雾水。2.1 关键概念quantized-mesh与LOD金字塔Cesium的地形瓦片格式叫quantized-mesh这是一种紧凑的二进制网格格式每个瓦片只存储顶点坐标、索引、法线、水位信息等必要数据。坐标被量化为16位整数所以体积小、解析快浏览器端直接用Uint16Array读取。它和普通图片瓦片最大的不同在于图片瓦片是栅格图地形瓦片是带拓扑结构的网格Cesium拿到网格顶点后可以做法线贴图、动态LOD、挖方填方这些计算。LOD金字塔是另一块基石。地形切片把全球或者某个局部区域按四叉树划分第0层是一张瓦片覆盖整个范围第1层分成4张第2层分成16张以此类推。每张瓦片代表一个矩形范围当相机靠近某个区域时Cesium请求更细层级的瓦片用更密的网格替代粗网格。生成器要做的事情就是在原始TIF里对每个瓦片范围做采样抽出一张合适分辨率的小网格。2.2 目录结构与模块划分这个生成器的输出目录标准到可以直接扔给Nginx当静态资源服务。典型结构如下terrain_output/ ├── layer.json └── tiles/ ├── 0/ │ └── 0/ │ └── 0.terrain ├── 1/ │ ├── 0/ │ │ └── 0.terrain │ │ └── 1.terrain │ └── 1/ │ └── 0.terrain │ └── 1.terrain └── 2/ └── ...层级目录名就是LOD的level下一层是行号再下一层是列号文件后缀统一为.terrain。layer.json是地形服务的入口描述文件Cesium发起地形请求时先读它再根据相机视锥动态计算需要哪一层哪一行哪一列。代码层面生成器主要分成五个模块TIF读取器、坐标系归一化器、瓦片范围计算器、高程采样器、quantized-mesh编码器。读取器依赖GDAL解析GeoTIFF的GeoTransform和栅格波段范围计算器负责把TIF的地理范围映射到瓦片行列号采样器处理边缘像素越界和插值编码器把二维高程数组变成二进制网格。模块之间通过一个中间高程缓存衔接避免重复读盘。2.3 工作流水线GDAL读取、金字塔构建、瓦片编码整条流水线我按三个阶段来理解。第一阶段是预处理用GDAL打开TIF读取仿射变换参数、投影信息、栅格尺寸然后按目标投影做重投影。大部分TIF是WGS84经纬度坐标Cesium地形服务也默认用WGS84所以这一步大多数时候是直接透传。如果输入是投影坐标比如UTM就需要用gdalwarp或GDAL的Python绑定转成EPSG:4326。第二阶段是金字塔构建。拿到TIF的覆盖范围和像素分辨率后先算出根节点应该在哪一层然后逐层把每个瓦片的地理范围算出来。这一步不直接读像素只做纯数学计算生成一个任务队列。第三阶段是编码输出从任务队列取一个瓦片到原始TIF对应区域采样高程值组成小网格量化为quantized-mesh格式写入目标文件。生成完所有瓦片后最后写layer.json。这三个阶段之所以拆开是因为第二阶段计算量小但决定了整个金字塔形态第三阶段计算量大但逻辑单一适合并发。实际执行时第一、第二阶段通常秒级完成第三阶段占90%以上的耗时。3. 核心模块拆解一个地形瓦片从采样到落盘的完整过程这一章是全文最值得抄的部分。我按一条真实的数据流讲拿到一个TIF先解析出它的范围和分辨率再逐层切瓦片每个瓦片内部做高程采样和二进制编码。每一步我都给出可以直接运行的代码骨架并解释参数含义。3.1 TIF解析与坐标归一化第一步拿GDAL读TIF的头信息。GDAL是地理数据的事实标准处理GeoTIFF、IMG、GBDX这些格式都靠它。核心代码如下from osgeo import gdal import math ds gdal.Open(src_tif) geo ds.GetGeoTransform() width ds.RasterXSize height ds.RasterYSize band ds.GetRasterBand(1) nodata band.GetNoDataValue() # geo六个参数含义: 左上角x, x方向分辨率, x旋转, 左上角y, y方向分辨率, y旋转 x_min geo[0] cell_w geo[1] y_max geo[3] cell_h abs(geo[5]) # 高程栅格的cell_h通常是负数取绝对值 # 转成经纬度范围 x_max x_min width * cell_w y_min y_max - height * cell_h这段代码的关键在geo数组。geo[0]和geo[3]是左上角坐标不是中心点geo[1]是每个像素的宽度单位跟随投影geo[5]是y方向的分辨率正常TIF里是负值表示从上往下递增所以取绝对值。很多人在这一步翻车直接把负分辨率当正数算结果整个地形上下颠倒。拿到范围后我还习惯做一次单位校验。如果TIF的投影是投影坐标x、y单位是米需要先用gdal.ReprojectImage把它转成EPSG:4326再做后续计算。不转的话瓦片行列号会错乱生成出来的地形加载到Cesium会平移几十公里。常见做法是先用gdalinfo看一下数据集投影再用gdalwarp转换gdalwarp -t_srs EPSG:4326 -r bilinear input_utm.tif output_wgs84.tif重采样算法选bilinear而不是near因为高程是连续场邻近采样会在地形边缘出现明显锯齿。3.2 LOD层级计算与瓦片范围划分接下来要决定金字塔从哪层开始、到哪层结束。CTB的做法是根据TIF的原始分辨率反推根节点层级根节点瓦片应该对比原始TIF“足够模糊”但又不能模糊到丢失地形特征。# 全球0层瓦片大约对应156543.03米/像素 root_res 156543.033 # 计算TIF分辨率对应的层级 lod math.ceil(math.log2(root_res / cell_w)) # 如果TIF分辨率是10米lod约为14 max_lod lod 8 # 默认再往下切8级这里有个容易困惑的点lod数字越大瓦片越精细。0层是整个世界一张瓦片分辨率极低14层已经能看清一个城市的山体轮廓。cell_w的单位是度时不能直接用这个公式需要先转成米或者用Cesium源码里那套按纬度修正的算法。我一般偷懒用度转米的近似值1度约111320米在低纬度误差可接受。瓦片范围划分用经典的墨卡托或经纬度四叉树。Cesium地形默认用经纬度做瓦片坐标轴所以逻辑比较简单根节点覆盖TIF的地理范围每往下切一层把上一层的矩形均匀分成4块。def tile_bounds(level, row, col, min_x, min_y, max_x, max_y): n 2 ** level w (max_x - min_x) / n h (max_y - min_y) / n return ( min_x col * w, min_y (n - 1 - row) * h, min_x (col 1) * w, min_y (n - row) * h )注意行号的方向Cesium瓦片行号从南到北递增而TIF的y是从北到南所以计算时要取n - 1 - row。这个符号方向是经典踩坑点生成出来的瓦片行列号反过来会导致地形在南北方向镜像。3.3 高程采样与Quantized-Mesh编码拿到瓦片范围后要从原始TIF里抠出对应像素块。为了性能我通常先算好这个瓦片在该层级应该有的网格大小再通过缩放采样到固定尺寸。Cesium官方推荐每个瓦片网格为65x65或129x129个顶点太小网格粗糙太大文件体积翻倍。CTB默认按65x65采样。import numpy as np # 从原始TIF读取指定范围的栅格 def read_elevation(ds, geo, x_min, y_min, x_max, y_max, target_size65): # 把地理坐标转为像素坐标 px0 int((x_min - geo[0]) / geo[1]) py0 int((y_max - geo[3]) / geo[5]) px1 int((x_max - geo[0]) / geo[1]) py1 int((y_max - geo[3]) / geo[5]) # 防止越界 px0, py0 max(0, px0), max(0, py0) px1 min(ds.RasterXSize, px1) py1 min(ds.RasterYSize, py1) # 波段读入为numpy数组 arr ds.GetRasterBand(1).ReadAsArray(px0, py0, px1 - px0, py1 - py0) # 双线性缩放到target_size x target_size from scipy.ndimage import zoom scale_y target_size / arr.shape[0] scale_x target_size / arr.shape[1] resampled zoom(arr, (scale_y, scale_x), order1) return resampled这里ReadAsArray的参数是像素行列号并不是经纬度。把地理坐标转像素坐标时要注意y方向geo[3]是北边界py0应该用y_max减而不是y_min加。很多人在这一步把上下颠倒生成出来的地形是反的。Quantized-Mesh编码是另一块核心。格式不复杂由头信息、顶点数据、索引数据、法线数据四部分组成。头信息固定74字节包含中心点坐标、最小高度、最大高度、边界范围、误差值其中误差值用来告诉Cesium这个瓦片在什么距离下可以被下一层级替换填0会导致无限加载下一层。import struct def encode_quantized_mesh(vertices, indices, center, min_h, max_h): # vertices是Nx3的numpy数组, 坐标已归一化 header struct.pack(3d3d3d2d3d, *center, # 中心点经纬度 *center, # 边界最小角 *center, # 边界最大角 min_h, max_h, # 最小/最大高度 *[0.0]*3) # 误差 body vertices.astype(u2).tobytes() idx indices.astype(H).tobytes() return header body idx这段代码省略了法线和水位部分但主体结构是对的。实际项目里法线需要由顶点坐标差分计算可以用numpy向量化实现不必逐三角形循环。indices要按行列顺序生成规则网格索引相邻两行三角形交替绕序否则地形会出现周期性的裂缝和穿模。3.4 layer.json索引生成的细节所有瓦片落盘后最后写一份layer.json。这个文件是Cesium的入口内容非常固定但有几个字段必须算准。{ tilejson: 2.2.0, format: quantized-mesh-1.0, version: 1.0.0, bounds: [-180, -90, 180, 90], available: [[0, 0, 0, 0], [1, 0, 0, 1, 1, 1]], projection: EPSG:4326 }bounds是地形覆盖范围必须跟实际生成的瓦片范围一致。如果TIF只覆盖中国某省bounds填全球Cesium在请求全球其他区域的瓦片时会得到404地形直接缺一块。available数组描述每层可用的行号和列号范围格式是[层级, 最小行, 最小列, 最大行, 最大列]。生成时遍历一遍实际落盘的瓦片文件统计出来不能手写手写必错。layer.json还有一个容易忽略的字段attribution如果TIF数据来自某测绘机构这里要写上来源。这既是合规要求也是习惯问题我在所有生成器里都默认带上。4. 运行与调参让工程化的地形生成跑起来代码讲清楚之后得让它真的跑起来。这一章覆盖从零开始的环境准备、命令行参数含义和一次完整的批量处理。我用的是Linux环境加Python 3.9Windows下逻辑相同但GDAL库的安装方式略有区别。4.1 环境准备与依赖安装生成器核心依赖只有三个GDAL、numpy、scipy。GDAL负责TIF读写numpy处理栅格数组scipy的双线性插值用于网格缩放。安装顺序有讲究先装libgdal系统库再装Python绑定否则容易碰到版本不匹配。sudo apt install libgdal-dev gdal-bin pip install GDAL$(gdal-config --version) numpy scipy注意pip装GDAL时版本号必须和系统gdal-config输出的版本一致不一致时编译报错。Windows用户建议直接装预编译的whl包或者用conda环境conda install -c conda-forge gdal numpy scipy装完验证一下python -c from osgeo import gdal; print(gdal.VersionInfo())能打印版本号就说明绑定成功。我见过不少人在这一步卡住原因大多是系统里残留了多个GDAL版本Python导入时把动态库解析到了旧版本上。4.2 命令行参数详解生成器的主体是一个CLI程序参数不算多但每个都直接影响产物质量。我把常用参数整理成一张表实际使用按需组合。参数类型默认值说明-i / --input路径无输入的GeoTIFF文件必填-o / --output路径./terrain输出目录自动创建-l / --max-level整数12最大LOD层级不宜超过18-s / --size整数65每张瓦片的网格顶点数-e / --error浮点0.1瓦片误差影响LOD分级-f / --format字符串quantized-mesh输出格式暂支持qm-p / --projection字符串EPSG:4326输出投影通常不用改-c / --compression开关关闭是否对瓦片做zlib压缩max-level参数最影响生成时间和占用空间。每降1级瓦片数量约为上一级的4倍所以从12级提到14级生成时间直接差一个数量级。我建议先用小范围TIF试跑确认效果后再调高层级。size参数也别乱改65x65是Cesium加载效率和渲染质量的平衡点改到129会让单瓦片重4倍但视觉效果提升有限。4.3 一个完整的批量处理场景假设手上有某个区域的DEM数据是UTM投影、5米分辨率的GeoTIFF。第一步先重投影再切片gdalwarp -t_srs EPSG:4326 -r bilinear -co TILEDYES dem_utm.tif dem_wgs84.tif python ctb_slicer.py -i dem_wgs84.tif -o ./terrain_dem -l 15 -s 65跑完后检查输出目录的瓦片数量是否对称。一个正方形的TIF第15层应该有2^15乘2^15规模的行列结构但由于TIF不是全球范围实际生成的瓦片只占一部分不是全覆盖。这时用官方layer.json里available字段能明确记录实际覆盖。批量处理多块TIF时我习惯先写一个拼接脚本用gdal_merge.py把所有TIF合成一张大的再做一次切片。这样能减少接缝处可能出现的重复采样。如果TIF之间位置有重叠合并时要指定优先级否则重叠区域的高程会取到平均值地形细节被抹平。gdal_merge.py -o merged.tif tile_a.tif tile_b.tif gdalwarp -t_srs EPSG:4326 -r bilinear merged.tif merged_wgs84.tif python ctb_slicer.py -i merged_wgs84.tif -o ./terrain_all -l 15合并时注意nodata值必须统一否则边界区域会出现黑洞。通常在地形数据里无值区域用-32768或0我建议在处理前统一赋一个小于所有有效海拔的值例如-10000并在layer.json对应的采样时把它剔除后续加载才不会有异常尖刺。5. CTB避坑指南坐标系、空洞与层级跳变的排查切片工具本身逻辑不复杂但地形数据太容易出各种隐藏问题。这一章我把实际踩过的高频坑整理出来每条都按现象、原因、解决的顺序写遇到类似问题可以照着排查。5.1 地形被整体拉伸或压扁椭球体基准不一致现象Cesium里地形看起来“变形”山体被拉高或者整个区域下沉几百米和真实地形对不上。原因原始TIF的高程基准和Cesium默认的WGS84椭球不一致。国内不少DEM数据用的是1985国家高程基准和WGS84椭球面之间有几十米的系统性偏差直接切片后Cesium把它当作椭球高渲染自然整体偏移。解决切片前先确认高程基准。用gdalinfo查看TIF的垂直坐标系如果发现高程基准不是WGS84就需要做垂直基准转换。GDAL本身不带完整的垂直变换参数我一般用两种方式兜底一是查到该区域的固定偏移值在采样时减去二是用EGM96网格做修正。生产环境里最省事的方案是获取该区域的已知控制点高程差推算出常量偏移量在生成器中加一个global_offset参数。5.2 瓦片边缘出现裂缝采样边界没对齐现象Cesium加载后相邻瓦片之间有一条条黑色缝隙放大后能看到三角形网格断裂。原因每个瓦片独立采样时边缘像素取了四舍五入后不同的结果。瓦片A右边界取到了第n列瓦片B左边界取到了第n1列两者共享的顶点值不一致网格就在接缝处开裂。解决采样时强制让相邻瓦片共用同一行像素。我通常的做法是给每个瓦片向外多读1个像素overlap在编码时只写入内部网格顶点交界处取左、上瓦片写入的值。更工程化的做法是采样时按“包含关系”计算像素范围即右边界使用左闭右开区间保证相邻瓦片不会重复取到同一个像素。调整完后再跑一遍对比裂缝会明显减少。5.3 内存暴涨或OOM大TIF直接全量读入现象处理十几GB的TIF时跑到一半进程被系统杀掉或者内存占用飙到几十GB。原因代码里用了ReadAsArray(0, 0, width, height)把整个栅格一次性读入内存再开始切瓦片。TIF像素数达到数亿级别时float32数组就是几个GB起步再加上金字塔所有的采样临时数组内存直接打爆。解决改成按瓦片范围局部读取即上一章代码里的方式每个瓦片只读它覆盖的那一小块。这样内存占用稳定在几百MB以内。如果TIF本身太大还可以用GDAL的VRT虚拟栅格先建金字塔读取时GDAL会自动按需调度。另外要注意numpy数组的类型读出来的int16转float32再计算不要全程用float64能省一半内存。5.4 加载后地形空洞layer.json的bounds写错现象Cesium地形加载后某些区域始终是平的甚至缺一块刷新也出不来。原因layer.json里的bounds范围写错了或者available数组把实际没有生成的瓦片也算进去了。Cesium发起请求后拿到404会一直重试或者直接放弃该区域的地形渲染。解决从输出目录里重新统计瓦片层级和行列号反向生成available。我在这块吃过亏后来干脆写了一个小函数遍历tiles目录把存在文件的层、行、列取最大最小值自动生成数组。这样就算某次切片中途失败导致部分瓦片缺失layer.json也能准确反映实际产物Cesium不会请求不存在的瓦片。5.5 LOD切换卡顿或闪烁最大层级设置过高现象相机拉近时地形瓦片加载特别费劲拉远又闪动像在“呼吸”。原因max-level设得太高瓦片总数成指数级膨胀Cesium并发请求几十个瓦片本地服务器响应不过来。同时误差值填0会让Cesium认为下一层永远需要加载导致无限细分。解决max-level根据原始TIF分辨率来控制。5米分辨率的TIF切到15级左右已经接近原始精度再往下没有新增细节还白白增加瓦片数。误差值不要填0按TIF实际分辨率给一个合理值比如0.5米。这样Cesium会在地面分辨率低于该值时停止请求下一层LOD切换会顺滑很多。6. 二次开发从“能跑”到“能改”的落地技巧工具能用是一回事要真正融入自己的项目得学会改。我在实际项目中做了三处改造每一处都对应一个常见的业务需求生成速度、瓦片存储、服务对接。6.1 改造一多线程切瓦片单线程切瓦片在大范围TIF上耗时很长一个10GB的TIF跑十几分钟很正常。因为瓦片之间完全独立适合用多线程并发。我用Python的concurrent.futures实现了一个简单的线程池把瓦片任务按层级分组提交每层内部用多个worker并行处理。from concurrent.futures import ThreadPoolExecutor def process_level(level, rows, cols, ds, output_dir): with ThreadPoolExecutor(max_workers8) as pool: for row in range(rows): for col in range(cols): pool.submit(generate_tile, level, row, col, ds, output_dir) def generate_tile(level, row, col, ds, output_dir): bounds tile_bounds(level, row, col, ...) elev read_elevation(ds, bounds, target_size65) mesh encode_quantized_mesh(elev, ...) write_tile(output_dir, level, row, col, mesh)注意线程数不是越多越好GDAL底层读文件时会有锁竞争8个线程是性价比比较高的点另外生成瓦片写文件时线程安全没问题但如果用了shared的numpy数组做缓存就需要加锁。改造后同样的TIF10分钟能压到3分钟左右瓶颈反而变成了IO。6.2 改造二瓦片缓存与增量更新真实项目中DEM数据会定期更新每次都全量重新切片太浪费。我加了一个增量模式先对比TIF的修改时间和瓦片文件的生成时间跳过没有变化的区域。这个逻辑看起来简单但实现时要能定位出“某块区域对应的瓦片层级”和“该瓦片对应的源数据像素范围”两个条件才能判断是否需要重新生成。增量更新后我还加了一层瓦片缓存策略输出目录里保留一份version.json记录每个瓦片的生成用到的TIF的md5值。下次跑任务时先比对md5一致就直接跳过不一致再重新生成。这个习惯帮我省了大量重复计算时间尤其是源数据几乎没变但程序参数微调的场景。6.3 改造三对接自有地形服务项目里没有用Cesium自带的在线地形服务而是把生成的瓦片部署到内部静态文件服务器上。这时要修改的是Cesium端的加载方式与生成器本身无关但我在生成器里加了一个输出选项把layer.json的路径写成相对路径并生成一份s3m格式兼容的目录说明方便运维配置CDN。我一般会在生成完瓦片后用CesiumTerrainProvider先本地验证一次const viewer new Cesium.Viewer(cesiumContainer); const terrainProvider await Cesium.CesiumTerrainProvider.fromUrl(/terrain/layer.json); viewer.terrainProvider terrainProvider;注意Cesium新版本里terrainProvider是个属性老版本用setTerrain API。验证时我会故意飞到某个已知坐标的高山上对比实际海拔确认输出没有垂直偏移。从那以后我每次跑完切片任务都强制自己先做一遍这个“飞高点”验证再交付避免把坐标系问题带到下游项目里。如果你也被地形数据折磨过希望这套思路能帮到你少走几步弯路。本文还有配套的精品资源点击获取
