哥白尼30米DEM全国分省打包处理全流程解析

GIS圈里折腾地形数据的朋友这两年基本绕不开一个名字哥白尼30米DEM。先别急着划走这次分享的不是那种随手甩个原始链接就完事的资源帖而是我自己把数据从原始瓦片一直处理到“全国分省、解压即用”的全过程。这个30米分辨率的FathomDEM数据属于哥白尼计划的再优化版本在原始GLO-30基础上把平坦地区的水文走向、陡坡伪影都做了修正对做水文分析、道路选线、山体阴影图的朋友来说比直接下原始HGT瓦片要省心得多。这次我按省级边界重新裁剪打包顺便把处理过程中的坑都记录下来适合刚接触DEM、需要全国尺度数据但不想自己折腾拼接投影的新手也适合那些已经下过原始数据、被各种黑边缝隙气到头疼的老手。1. 哥白尼30米DEM和FathomDEM到底强在哪1.1 先从数据源说起哥白尼DEM是目前全球公开可下载的高程数据里性价比最高的那一档。它由两颗雷达卫星不同视角对地表进行干涉测量处理后生成了一套全球数字高程模型。官方公开版按分辨率分为GLO-10和GLO-30两档大家口中常说的“哥白尼30米”对应的是GLO-30版本。与老牌SRTM 30米相比它的优势非常明显纬度覆盖到北纬60度与南纬56度之间范围更完整山区和林区的雷达盲区大幅减少空洞率低在平坦区域的竖向精度实测能到1到2米级别拿来当项目前期地形分析的底图完全撑得住。但这个数据有个让所有人头疼的特点官方发布的原始文件是按1度乘1度的经纬度网格切好的全国范围横跨好几十个网格下载下来是几百个小碎片。更麻烦的是这些WGS84地理坐标系下的分片直接扔进ArcGIS或者QGIS里做坡度、做填洼在高纬度或大跨度区域计算出来的结果会变形失真。所以“拿到手先拼接、再投影、后裁剪”是整个流程里绕不开的三个步骤也是我这次分省打包的核心工作。1.2 FathomDEM在哪些地方被优化了FathomDEM是根据哥白尼GLO-30数据进一步优化出来的水文修正版本。之前有朋友问我说直接下原始GLO-30不就行了多这一道优化有什么意义这里需要解释清楚官方GLO-30在做全球拼图时个别区域存在雷达信号失相干导致的尖峰、凹陷以及水面上的起伏噪声。这些噪声渲染山体阴影图时会表现为一条条放射状的虚假纹理做水文分析时又会让河网提取结果出现断头河和伪沟谷。FathomDEM的优化思路说白了就是“专门为水流走向服务”在保持真实地形骨架不变的前提下把河道位置压低、把平原区的微小起伏做合理平滑同时清理掉卫星数据里明显的飞点噪声。我实测下来的感受是如果只做宏观地形展示、坡度分级底图原始GLO-30完全够用但换成和流水方向强相关的场景比如洪水淹没模拟、排水路径分析、流域划分优化后的版本能少处理很多莫名其妙的异常小沟。当然它也不是万能药在西南高山峡谷地区部分陡崖的细节仍被平滑掉了一些这个后面使用建议里我会专门提醒。1.3 为什么30米分辨率是当前的最优平衡点可能有人会问既然要分享为什么不直接上分辨率更高的数据道理很简单高分辨率数据要么收费昂贵要么公开覆盖有限。5米、2米级的商业DEM在县城尺度好使放到全国范围数据量直接爆炸后期存储和处理成本一般人吃不消。而10米公开数据目前在全球多数区域覆盖不全真正能稳定获取的全球趋近一致产品就是SRTM 30米和哥白尼30米这两个档位。SRTM虽然同为30米时间却老了一代空洞和噪声偏多哥白尼30米无论从时相、覆盖度还是细节上看都已经远远甩开了前辈。对全国尺度的分析来说30米格网在计算量和地形细节表现之间踩在了一个很舒服的位置流域面积计算、坡度带划分、通视分析都不至于因为格网太粗而失真又不会因为数据量太大让普通台式机跑不起来。这一版全国分省打包正是为了让更多人能直接跨过“数据获取和处理”的坎专心做自己的业务分析。2. 全国分省打包的整体思路和实操流程2.1 拿到数据后先做整体规划动手之前先想清楚最终要交付什么。我的目标是全国按省级边界分好的GeoTIFF文件每个省一个文件坐标系采用CGCS2000地理坐标系方便国内项目直接叠加各种矢量边界使用。整个流程分成六步第一步下载原始瓦片按区域筛选第二步用虚拟栅格VRT做逻辑拼接不急着物化大文件第三步把拼接结果重投影到目标坐标系并重采样到30米分辨率第四步用省级边界矢量做裁剪第五步对裁剪结果统一处理NoData和异常值顺手构建金字塔第六步分省压缩打包附上说明文档。这里有一个关键思路我要强调不要先拼全国再裁剪除非你的电脑是服务器级别。全国范围的原始瓦片加起来几十GB桌面GIS软件直接merge卡死是常态。按省处理的话每个省最多涉及几十个小瓦片内存占用完全可控。如果某个大省跨了好几个经度带后面单独补投影细节即可难度低得多。2.2 拼接、重投影和裁剪的具体操作我这次用QGIS配合GDAL命令行完成QGIS选用长期支持版GDAL版本在3.6以上这样对COG格式的支持比较完善后续分省文件还能保持云优化结构。第一步拼接。官方原始瓦片格式一般是TIFF或HGT推荐先建一个虚拟拼接索引gdalbuildvrt -overwrite -srcnodata -32767 -vrtnodata -32767 E:/raw_tiles/*.tif E:/working/province_A.vrt注意命令里的NoData值一定要统一。哥白尼DEM的无效值在多数瓦片里是-32767但也有部分版本用0表示所以先执行gdalinfo逐一检查把每个瓦片的NoData确认好再建VRT。gdalbuildvrt只是建立逻辑索引几秒钟就完事不会生成占用几十GB的大文件速度极快。第二步重投影。因为后续裁剪用的省级边界是CGCS2000地理坐标系而栅格原始是WGS84虽然两者在平面表现上差异很小但严格项目不能混用。这里更要紧的是重采样到“真30米”在地理坐标系里没法直接说“30米”必须切到以米为单位的投影坐标系。我选用了中国范围常用的Albers等积投影参数中央经线105度两条标准纬线25度和47度用PROJ字符串直接传给GDALgdalwarp -overwrite \ -t_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumCGCS2000 unitsm no_defs \ -tr 30 30 -r bilinear \ -srcnodata -32767 -dstnodata -32767 \ E:/working/province_A.vrt E:/working/province_A_aea.tif-tr 30 30的含义是重采样后像元为30米乘30米重采样算法选双线性。为什么不选最邻近最邻近会把地形切出锯齿状台阶三次卷积精度好但容易把NoData边缘扩散一圈出去反而增加后处理负担。双线性是DEM重采样最稳妥的默认选择。第三步裁剪直接用带-cutline参数的gdalwarp一步完成省得先生成裁剪再用工具提取gdalwarp -overwrite \ -cutline E:/boundary/province_A.shp -crop_to_cutline \ -dstnodata -32767 \ E:/working/province_A_aea.tif E:/output/province_A_cgcs2000.tif裁剪完记得检查输出栅格的范围防止出现空白输出或者只有一小条的情况。常见原因是省界矢量坐标系和栅格不一致叠加前先统一。2.3 关键参数选择背后的用意很多人第一次处理DEM最容易纠结的就是“我要不要填洼”。我在这套分省数据里默认不做填洼处理原因很简单填洼是破坏原始信息的一步操作是否执行完全取决于下游应用。做流域提取填洼能避免平地上的水卡在原地做凹凸坡分析、做地形粗糙度填洼反而会把真实地貌抹掉。所以正确姿势是“分省原样数据用户按需填洼”而不是我在打包阶段替用户做决定。另一个参数是数据格式。我最终选用COGCloud Optimized GeoTIFF输出配合DEFLATE压缩。COG格式有一个很实用的特性数据内部带金字塔和块状结构软件按需读取不用把整个省的文件全部加载进内存。压缩率方面对于地形数据这种邻居值高度相关的栅格DEFLATE通常能把文件压到原始大小的一半以下而且绝大多数新版GIS软件都能直接读取兼容性也好。内部块大小我选了256乘256既不浪费IO又能让局部读取很快。这些细节不影响高程数值但直接影响你加载省数据时的等待时长。3. 分省打包的规范和发布说明3.1 分省文件的命名与目录结构整理好的分省数据目录设计要让人一眼就懂我采用如下的结构FathomDEM_China_30m/ ├── 00_README.txt ├── 01_数据说明与坐标系.txt ├── data/ │ ├── 华东某省.tif │ ├── 东北某省.tif │ ├── 西南某省.tif │ └── ...各省份文件 └── preview/ └── 全国山体阴影示意图.jpgdata目录下按行政区划习惯排序文件名直接用省级名称方便用户在业务系统里批量挂载。preview目录里放一张全国范围的缩略山体阴影图一张图就能直观看出数据有没有黑斑、有没有接边断裂。很多分享资源根本不提供预览图用户下载半天打开发现数据有问题非常浪费时间所以我强烈建议打包者一定顺手出一张预览图。3.2 打包压缩方式与校验信息分省文件单独压缩我还是用了ZIP格式压缩级别选“最快”或“标准”即可。为什么不选7z或RAR一是ZIP在Windows、Linux、macOS上都是原生支持任何人拿到都能直接解压零门槛二是DEM文件的压缩率受限于内部已经做过的DEFLATE外层再压7z提升有限反而增加解压成本。每个省压缩包内部除了主TIFF文件建议同时放一份与该省同名的QMD或AUX辅助文件吗不推荐辅助文件太多容易丢失GeoTIFF自带的坐标系信息已经足够。压缩完成后我习惯在包根目录放一个SHA256校验文件列出每个压缩包的校验值。下载数据这个场景文件中途损坏非常常见。有了校验文件用户解压前可以先跑一遍sha256sum避免用到坏文件后排查半天结果发现是源数据坏了。3.3 README里必须写清楚的内容打包数据最忌讳的是“只给数据、不给说明”。我写了这几项全部踩过坑之后才改出来的经验首先是坐标系和投影参数明确写上“CGCS2000 / Albers等积投影”并附上中央经线、标准纬线的具体数值。其次是NoData值统一为-32767单位是米。如果有人加载后看到大片黑色往往就是没看这一行把NoData当0米处理了。第三是分辨率、像素深度和压缩方式方便高级用户判断能不能直接做定量分析。第四是数据来源和处理日期注明原始数据文件版本、处理软件和GDAL版本这既是对自己工作的记录也能在数据异常时快速定位问题。最后一定要写清楚“建议使用场景和禁用场景”明确说明这套数据用于流域分析、区域规划、专题制图没问题但不适合作为厘米级工程设计的基础地形。4. 处理过程中的实战踩坑记录4.1 下载阶段怎么选瓦片才能不漏不重原始数据按经纬度网格分片一个省往往横跨好几个网格凭肉眼数会数漏。我的办法是先用QGIS创建经纬度网格矢量再把省界叠加上去用“按位置选择”选中所有相交网格然后批量生成下载清单。下载器方面Windows下我推荐命令行工具配合断点续传参数后台挂着慢慢下比如用带-c参数的wget或curl。如果网络不稳定中途断掉的线程不会从头再来实测能节省大量重试时间。全部下完以后必须比对每个文件的大小是否和源站一致不一致的直接删掉重新下千万不要带着残缺瓦片进入拼接环节否则会出现明明范围重叠却偏偏缺了一块的情况。4.2 拼接阶段内存溢出和VRT解法第一次按省拼接的时候我直接用QGIS的“合并”工具结果一个平原大省就把16GB内存吃满软件直接崩溃。后来改用gdalbuildvrt问题瞬间解决。VRT本质上是一个XML索引只记录每个瓦片的空间位置和路径不复制任何像素。等你重投影或裁剪时GDAL再按需读取瓦片并实时计算。所以整个处理路线的核心原则是能懒则懒能索引则不物化把真正的像素读取推迟到最后一步。如果实在因为业务需要必须物化一张临时大图建议用带-wm和-co TILEDYES参数的gdal_merge.py限制单线程内存池大小即便几十GB的全国数据也能慢慢跑完只是时间会拉长不影响结果正确性。4.3 裁剪后的黑边、缝隙和横条纹裁剪之后最容易出现的三类视觉问题我一个个说。黑边几乎都是NoData值不统一造成的瓦片A的无效值被记成-32767瓦片B里却是0拼接时0成了一片似真似假的高程裁剪后边缘就成了黑洞。解决方法就是我在2.2节里说的建VRT前先统一NoData。缝隙是瓦片之间在某一行列重叠不足多数出现在数据源快照不一致的拼接处。对这种细小缝隙可以用GDAL的gdal_fillnodata插值填补它会把周围有效高程的加权平均值填进空洞。注意这个方法只适合填补狭长的小缝隙如果是整片大面积空洞填出来的地形完全是“想象值”宁可保留NoData交给用户自行决定也不要给出错误的地形走向。横条纹在部分山区比较明显这是雷达侧视成像的固有现象表现为条带状的细微高程起伏。处理思路是用低通滤波减弱但滤波器半径不能太大否则把山谷和山脊的真实落差也给糊掉了。我实际用3乘3的均值窗口能把横条纹抑制到肉眼几乎不可见又不会伤到地形主结构。4.4 坐标偏移和对不齐有一次裁剪出来的结果和道路矢量叠加后明显错位我查了半天发现省界矢量是CGCS2000栅格还是WGS84。这两个坐标系在国内大陆区域差值大约在几十厘米到一两米之间放大到制图比例时视觉上完全能看出重影。解决办法是第一时间用gdalinfo确认每个文件的CRS描述最后强制输出到目标坐标系。图层一旦多起来务必用QGIS“项目属性-坐标参考系统”把整个工程设成同一个坐标系让软件自动重投影。另外还有一个容易被忽视的点省界边界数据本身的精度。如果手里的边界是多年前的勘界版本和最新行政区划有出入裁剪出来的面积和范围就会和官方统计对不上。做省级数据时尽量找最新的标准边界这一步偷懒后续所有分析结果都会跟着偏。4.5 高程异常值排查处理完的省文件用gdalinfo -stats看一遍最大最小值诸如“最低负一万多”这种明显不正常的结果基本就是NoData数值被当成有效高程统计进去了。ArcGIS里有时候拉伸显示后地面出现许多极亮或极暗的点则是飞点噪声。针对离散飞点我会再用条件判断找出超出合理范围的像元把其替换为邻域中值。注意这个操作只针对少数异常像元绝不能全图滤波否则高山区真实裸岩也会被误伤。5. 拿到数据之后怎么用才不算暴殄天物5.1 加载与可视化技巧分省数据下载后第一件事先看预览图确认没有大面积黑边然后在QGIS里加载tif文件。如果显示灰蒙蒙一片说明符号化还是默认的线性拉伸。先右键图层属性把渲染类型改成“山体阴影”或者用“单波段灰度色带”的组合底色用从绿到棕到灰白的渐变表达低到高再叠加20%透明度的山体阴影立体感马上出来。ArcGIS里则可以在符号系统里选择“拉伸”拉伸类型选“直方图均衡化”再叠加“山体阴影”作底图视觉冲击力很强特别适合做汇报前的区域地形展示。5.2 水文分析和坡度计算的正确姿势这套数据最擅长的就是水文分析。按照“填洼→流向→汇流累积→阈值提取”的标准流程走下来效果比旧版SRTM好很多但有几个细节不能说没踩过坑。填洼前先想清楚分析尺度在全国范围提取大河流填洼后提取结果基本合理到了小流域尺度30米格网上一条两三米宽的天然沟渠可能就三四个像元提取结果会断断续续。如果发现河网在平原区疯长到处都是伪河道先别急着怀疑数据大概率是汇流累积阈值设得太低。阈值推荐从1000或2000开始试在结果和视觉噪声之间找平衡。坡度分析方面如果按度为单位输出平缓区域会出现明显的阶梯状色块这是学位换算造成的量化效应。把输出单位改成百分比或者使用方向统计量视觉上的阶梯感会明显减弱。另外在处理分省数据时边界像元由于缺少相邻省份的信息坡度值会系统性偏低跨省项目建议先拼接再分析分析完再裁剪。5.3 与其他DEM产品的对比取舍这套FathomDEM数据发布于近年时相较新比ASTER GDEM和SRTM在细节上优势明显。但如果手头要做的事对某个山谷的局部形态非常敏感比如微型水电站选址我还是建议用更高分辨率的无人机航测或激光雷达点云去复核。任何30米全球数据都只能提供一个宏观骨架不能承载精细工程的全部信息。另外提供一点实战心得把FathomDEM和旧的SRTM叠加对比时两者在高山区绝对高程差可达十几米这不是谁的错而是不同传感器、不同成像年代带来的系统性差异。跨年度的地形变化分析最忌讳混用两套基准不同的DEM要么都用同一套数据做相减要么先做校正对齐否则得出的“地表变化量”很可能是数据源差异造成的假象。这版全国分省FathomDEM数据我前前后后处理了两周最深的体会是数据处理的功夫七分花在数据源整理上三分才是参数调优。把瓦片下载、坐标系统一、NoData清洗做到位后面所有分析都能安心跑。最后再分享一个小技巧如果你要拿分省数据做全国拼接相邻省份之间在边界上会有很小的高程跳变用QGIS里“虚栅格目录”加载所有省份文件再把渲染模式设为“使用栅格目录内嵌金字塔”就可以无缝浏览全国数据不必真的合出一张几十GB的大图来。希望大家拿到这版数据后都能少熬夜、多出活。