ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

DEM高程数据碎图镶嵌合并全攻略:ArcGIS与QGIS操作详解

DEM高程数据碎图镶嵌合并全攻略:ArcGIS与QGIS操作详解 1. 这些高程碎图是怎么来的不合并会怎样1.1 分幅数据的常见来源与分块逻辑干GIS这行几乎没有人能绕开高程数据。做淹没分析、汇水区划分、坡度坡向计算、天际线模拟、土方量估算第一步永远是拿DEM或DSM。问题是你很难一次性下载到覆盖整个项目区的完整高程影像。国内很多公开的数据平台、国外的USGS、OpenTopography都倾向于把数据按标准分幅或者按经纬度网格切成小块提供下载。这种切法有它的道理文件体积可控单幅数据大小适合网络传输服务器端也方便做瓦片化缓存。但到了用户手里就变成了一堆让人头疼的碎图。举个例子。你在平台上框选了一个覆盖三个县的项目区域系统按1:1万标准图幅切分一查下载列表可能跳出来六七十个TIF文件。每个文件单独看看范围是一条窄窄的条带或者一块小方格单拿出来根本没法用。你得把它们先拼成一张完整的大范围DEM后面所有分析才好统一处理。这个拼接动作在专业术语里叫镶嵌Mosaic或者叫镶嵌合并。简单说就是把多张有相邻或重叠关系的高程栅格按一定规则合并成一张连续无缝、坐标和像元大小统一的新栅格。1.2 碎图的三个日常麻烦加载、配色与接缝有段时间我在做一个片区的地形分析项目范围不大但正好骑在四条分幅带上。我没合并就直接把四张图拖进ArcMap里用结果给自己挖了一堆坑。第一个麻烦是加载和显示效率低。每次打开工程四张图要同时参与渲染缩放大范围时要同时处理四个文件卡顿不说配色还各管各的。同一片山体落在不同分幅里的颜色深浅不一样接缝处一眼就能看出来放到汇报截图里非常难看。第二个麻烦是后续提取数据的步骤要重复四遍。我需要在特定范围内提取高程断面线如果数据不合并就得先裁剪四张图再分别提取再把结果手动拼接中间坐标系、范围、分辨率稍微不一致结果就对不上。第三个麻烦最隐蔽——接缝处的高程值可能存在微小不连续。不做镶嵌直接分析时坡度图在山谷接缝位置容易出现一条虚假的“刀痕”明明地面很平缓算出来却像有一道台阶原因就是相邻分幅边缘的像元值存在细微差异。这三个麻烦叠在一起基本上等于逼着你学会镶嵌合并。早一点把数据拼好后面省出来的时间足够你在项目交付前多跑两遍质检。1.3 镶嵌合并的价值与适用边界镶嵌合并的价值不只是“把图拼在一起”这么简单。合并后的单一大栅格在空间分析中的意义主要体现在三个方面。第一分析范围的完整性。很多工具比如水文分析里的填洼、流向计算输入数据如果是一堆分离的碎图边缘会出现大量无效结果因为这些工具假设数据是连续覆盖的。只有合并成完整DEM分析结果才可信。第二处理效率提升。一张大TIF在后续的裁剪、重投影、统计、可视化中流程上更顺也更容易写脚本批量处理。第三数据管理的规范性。项目交付时给甲方提供一份合并后的标准DEM比给十几份零散分幅文件专业得多也方便归档。当然镶嵌也不是万能的它解决的是“空间范围连续”的问题无法改变原始数据的质量和精度。如果源数据本身存在明显的错误高程值、投影混乱、范围不重叠等问题镶嵌以后这些问题只会被放大。所以动手镶嵌之前有几个准备工作比打开工具更重要这部分内容放到下一节细说。2. 动手前的三道检查比开工具更重要很多新手拿到碎图就急着打开镶嵌工具一顿操作之后发现结果是歪的或者接缝处全是黑边白边然后开始怀疑软件坏了。实际上大部分镶嵌失败案例根源都在准备工作没做好。我自己的习惯是“先查三件事再打开工具”。2.1 坐标与投影一致性检查第一件事确认所有碎图的坐标系是否一致。你可以把所有分幅文件加载到ArcMap或者QGIS里打开图层属性看坐标系也可以用工具批量读取。最稳妥的办法是在文件管理层面检查元数据或者利用GIS软件自带的属性查看功能。对于从同一个数据平台、同一次下载任务里拿到的分幅数据坐标系通常是一致的比如都是CGCS2000 3 Degree GK CM 111E或者都是WGS 1984 UTM Zone 50N。这种情况下直接镶嵌没有问题。但如果你从多个渠道收集数据比如一个项目跨越了不同省份的公开数据接口那么不同来源的DEM很可能混用了不同的投影参数甚至有的是地理坐标系有的是投影坐标系。这时候要进行镶嵌必须先统一投影否则拼接出来的结果会出现明显的位移和错位表面上看是“接不上”实际上是坐标基准就不对。统一投影的操作本身不难ArcGIS里用Project Raster工具QGIS里用Raster Projections选好目标坐标系后可以把全部数据批量重投影。需要注意的一点是重投影会改变像元值和像元大小最好同时设置好重采样方式对高程数据来说最常用的重采样方式是双线性插值Bilinear因为它比最近邻法Nearest Neighbor更平滑又不像三次卷积Cubic那样计算量大。这一步完成后把所有碎图加载到同一个地图文档里肉眼检查一下边界是否严丝合缝再进入下一步。2.2 NoData值确认黑边白边的最常见源头第二件事也是最容易出问题的一环——NoData值的确认。高程栅格和普通影像不一样它里面的有效数据是地形表面但边缘往往有大量的无效区域。这些无效区域在文件里通常用一个大负数或者特殊值填充比如-9999、-3.402823e038也有用0的但用0的情况不多因为0表示海平面容易被误判为真实高程。镶嵌的时候如果工具没有正确识别NoData值这些无效区域就会被当作真实的高程值参与计算。结果就是合并产物的边缘出现黑边或白边有的甚至在数据内部出现一条条“灰色撕裂带”非常显眼。更麻烦的是这类错误会直接影响高程统计比如你用镶嵌结果计算最大高程不小心把NoData的-9999也算进去了最大值直接就废了。所以动手之前一定要确认每个碎图的NoData值到底是什么。ArcGIS的Raster Properties里可以查看QGIS的图层属性里也能看到。如果不同分幅的NoData值不一致最好先统一设置。ArcGIS里可以用Copy Raster工具在NoData Value参数里设定统一值QGIS里可以用GDAL的Gdal_translate命令行或者在建虚拟栅格时统一下设置。这一步做扎实了后面能有八成把握不出现黑边白边问题。2.3 重叠区域情况摸底数据之间是严丝合缝还是彼此交叠第三件事摸清碎图之间的空间关系。相邻分幅之间可能边界刚好接上也可能存在一小条重叠带甚至可能出现缝隙。如果你的分幅数据是按标准图幅切的通常边缘是刚好相接的或有一两行像元的压边重叠。这种重叠在镶嵌时影响不大选好融合规则即可。如果重叠带很宽比如无人机数据处理得到的航测成果或某些采集方式产生的大面积重叠那就需要格外注意融合规则的选择否则重叠区域会出现重影和模糊。如果分幅之间存在缝隙也就是两幅图之间有几行像元谁都没覆盖到镶嵌后的结果里会留下一道空白带需要额外处理比如用插值工具补洞或者重新找数据源补覆盖但这属于比较极端的情况。怎么摸底最简单的方法是把所有碎图叠加到地图视图里把图层符号设置为空心、加粗边框然后缩放到全图范围眼睛扫一遍空间关系。碎图数量多的时候也可以借助ArcGIS的Scanned Map或QGIS的Overlay分析来快速统计重叠范围。这个步骤看起来不起眼但能让你在后面选择融合规则和镶嵌范围定义时心里有数。3. 两条最快路径ArcGIS镶嵌与QGIS虚拟拼接准备工作做完就可以正式动手了。我平时最常用的工具是ArcGIS和QGIS两条路径各有特点。这里把两条路都讲清楚你按自己的软件环境选一条走通就行。3.1 ArcGIS路径Mosaic To New Raster参数逐项说明ArcGIS桌面环境里做一次性镶嵌最顺手的工具是Mosaic To New Raster。它在ArcToolbox里的位置是Data Management Tools → Raster → Raster Dataset → Mosaic To New Raster。这个工具的功能就是把多个输入栅格合并生成一个新的栅格数据集而不是修改原始数据适合我们这种“拼一张图出来”的需求。工具对话框打开以后有几个参数值得逐项仔细过一遍。Input Rasters把你要合并的碎图文件一次性添加进去。可以多选也可以直接拖入整个文件夹里的多个文件。这里我建议不要用通配符全选文件夹里的所有文件除非你确认里面没有混入其他东西否则很容易把一些辅助文件或者文档一起选中导致报错。Output Location输出结果要存放的位置建议指向一个文件地理数据库File Geodatabase或者某文件夹最好不要直接放在根目录下方便管理。Raster Dataset Name with Extension给输出栅格起名字。注意如果你选择存到File Geodatabase里后缀名写不写都行如果存到文件夹里需要带上.tif扩展名。Coordinate System for Raster这里如果留空工具会默认用第一个输入文件的坐标系。如果你在前面检查投影时发现各碎图坐标系一致这里不用特别设置如果不一致我建议先重投影再镶嵌而不是在这里选坐标系因为在这里选坐标系并不能自动帮你把数据重投影它只是给输出结果打上一个坐标标记数据本身的像元值并不会重新采样到新坐标系下。这一点非常容易踩坑切记。Pixel Type像元类型也就是输出数据的数据深度。常见的高程DEM是16位有符号整型16 bit signed integer如果源数据是浮点型DEM那就选32 bit floating point。不知道怎么选的时候看一下某一张源图的属性参考它的Bit Depth和类型来选否则可能导致高程值精度损失或溢出。Cell Size输出像元大小默认取第一个输入文件一般情况下保持默认即可。如果碎图分辨率不一致这里就需要手动指定一个目标分辨率建议使用所有碎图中最小最精细的那个分辨率以免高分辨率数据被强行降采样导致细节丢失。Number of Bands波段数高程数据通常为单波段选1。Mosaic Operator和Mosaic Colormap Mode是两个容易让人犯迷糊的参数它们的本质是决定重叠区域“听谁的”和“颜色怎么过渡”放到下一节专门展开讲。所有参数设置完点OK工具会开始跑。输出结果像元大小、范围和坐标系都统一后会在指定位置生成一个新栅格。处理碎图数量大的时候这一步可能需要几分钟到十几分钟耐心等待即可。3.2 QGIS路径先用VRT预览再导出成正式TIF如果你是QGIS用户或者工作中偏爱开源工具那么GDAL的镶嵌功能是最稳的选择。QGIS里最常规的玩法是先用Build Virtual Raster构建虚拟栅格生成一个VRT文件这个文件本质上是个XML目录它并不拷贝像元数据只是“引用”了各碎图极其范围和坐标信息。VRT的生成几乎瞬间完成适合快速预览镶嵌效果。具体操作路径菜单Raster → Miscellaneous → Build Virtual Raster。在弹出的对话框里把碎图文件添加进Input layersResolution参数可以选First layer的默认值也可以手动指定Separate band设为不勾选最重要的一项——Place each input file into a separate band这个千万别勾勾了以后生成的VRT会把每个文件放在单独波段里看起来就像一张多波段影像而不是我们想要的单波段镶嵌结果。设置好Output VRT路径点Run几秒钟后就能看到一张虚拟的“镶嵌图”了。VRT的好处是快作用相当于预览。但虚拟栅格有一个明显的局限它还依赖原始碎图文件存在且路径不变如果你把碎图移动了VRT就失效。所以正式交付或作为后续分析数据时还是要把VRT“实体化”转成真正的TIF文件。方法也很简单菜单Raster → Conversion → Translate把VRT作为输入选择输出格式为GTiff目标路径填一个.tif文件名点Run。或者在命令行里用GDAL自带命令gdal_translate -of GTiff -co COMPRESSLZW -co TILEDYES input.vrt output.tif这里的COMPRESSLZW表示无损压缩TILEDYES表示输出为金字塔分块存储的TIF后续加载显示会快很多。如果你希望同时生成金字塔统计信息还可以追加一段命令gdaladdo -r average output.tif 2 4 8 16 32这段命令会生成不同层级的概览金字塔在大型TIF缩放显示时明显提升速度。QGIS路径的优势在于透明、灵活、可脚本化尤其是碎图数量较多时用命令行批量处理比鼠标点选高效得多。3.3 工具选型建议什么时候该用Mosaic Dataset前面两条路径都是“一次性镶嵌成一张新图”。但还有一种使用场景比如你管理的是全省或整个流域的大范围高程库数据量可能有几千上万幅碎图每次全量镶嵌一次不仅耗时还会生成一个几十GB的单一大文件操作起来反而不方便。这时候更合适的是ArcGIS里的Mosaic Dataset镶嵌数据集而不是Mosaic To New Raster。Mosaic Dataset本身不复制栅格数据它只是创建了一个“目录式”的管理结构把碎图的路径、坐标、轮廓、属性统一登记起来并在需要显示或分析时动态调用对应分幅数据。它和VRT有点像但功能更强支持按需裁剪、实时投影变换、自定义服务发布等功能而且能更好地处理海量数据。Mosaic Dataset的缺点是概念和操作相对复杂对于只有几十幅碎图的临时项目“杀鸡用牛刀”反而增加工作量。我的建议是碎图数量在几十幅以内、项目是一次性任务用Mosaic To New Raster或QGIS的VRT导出就够如果是长期维护的数据基础设施比如要做地形数据服务对外发布那就老老实实建Mosaic Dataset。4. 重叠区融合规则的本质选对Operator等于选对接缝策略第一次做镶嵌的人看到Mosaic Operator这个下拉列表里面写着First、Last、Blend、Maximum、Minimum这些选项大概率会随便选一个完事。实际上这个参数直接决定了重叠区域像元值的最终走向。选对了接缝平滑无痕选错了数据里会留下明显的条带状伪影。4.1 三种常用融合规则的原理与利弊我用一个生活化的方式来解释这些规则。想象你面前有两块窗帘边缘重叠了十厘米现在要把它们缝成一块更大的窗帘。First直接用第一块窗帘的布料覆盖重叠区第二块在重叠区的部分被完全舍弃。好处是处理速度快逻辑简单坏处是在接缝位置如果两块布料的纹理高程值不一致会留下一条明显的“台阶”。Last相反重叠区直接用最后一块窗帘。逻辑同样简单但同样会有接缝问题。Blend这是我最常用的规则。两块窗帘在重叠区按权重渐变融合——靠近第一块的地方主要显示第一块的纹理靠近第二块的地方主要显示第二块的纹理中间部分两边各占一半。这样接缝被平滑地摊开了视觉上几乎看不出交界线。Maximum/Minimum重叠区取两侧像元的较大值或较小值。这种规则的逻辑是明确的比如合并不同测区的高程数据时如果只关心最高点才会用Maximum但它很少用于常规的DEM镶嵌因为地形表面不是“取最大”就合理的。Mean把重叠区两侧像元求平均。结果比First平滑一些但如果两侧数据存在系统性的高程偏移比如两个数据源的椭球面基准不一致求平均后重叠区会呈现出中间渐变的“过渡带”看起来有点像两道坡之间的缓坡。因此对于绝大多数高程碎图镶嵌任务我的首选是Blend。它的计算量稍大但换来了重叠区最自然的过渡接缝痕迹最不明显。如果数据重叠带很窄甚至刚好无缝拼接那么First也能用基本看不出区别。4.2 从一次失败案例理解BLEND权重渐变在很多年前的一次项目里我拿到两块相邻的山地DEM拼接片重叠区大约有5个像元的宽度。第一次操作时我图省事Mosaic Operator选了First。结果拼接完成后我放大到接缝附近发现坡面线在这里出现了一条细小的“折痕”原本平滑的山坡在重叠区突然像是被削了一刀。用ArcScene做了三维显示更是触目惊心整片山体在拼接处有一道微弱的亮线仿佛地形裂开了。后来我改用Blend重新镶嵌只做了这一个参数的调整再检查剖面线折痕消失了坡面恢复连续。这件事给我的教训很深刻对于地形数据来说接缝处的微小高程突变在平面图上可能肉眼看不太出来但一旦参与坡度、曲率等二阶分析就会被成倍放大。不要小看一个参数的选择它对结果质量的影响经常是决定性的。Blend权重的计算逻辑本身并不复杂。假设重叠区从碎图A到碎图B跨越N个像元那么靠近A方向的像元值主要由A决定权重线性过渡到0另一侧则主要由B决定。GDAL和ArcGIS在实现时都采用的是类似的线性权重或者三角权重策略。正是这种渐变过渡让接缝实现了“软着陆”。4.3 采样方式与像元类型容易被忽视的隐藏参数除了Operator还有几个容易被忽视的参数会影响最终的镶嵌结果。一是重采样方式Resampling。如果你的碎图分辨率完全一致像元网格能严格对齐那这个参数无关紧要但如果分幅之间分辨率有轻微差异比如有的地方是10米有的是12.5米那么镶嵌时必然牵扯到重采样。对高程数据来说最近邻法Nearest虽然快但会把原来的值生硬地搬到新网格上地形细节有损失。双线性插值Bilinear兼顾了平滑和处理速度是DEM镶嵌最常用的选项。三次卷积Cubic更平滑但计算量大且容易在地形陡峭处产生过冲反而出现不合理的负值一般不建议使用。二是像元类型。如果源数据是16位整型输出设置了8位整型高程值会被截断山的细节直接丢失如果源数据是浮点型但输出设成了整型高程小数部分全部抹掉平地和高山的分辨能力严重下降。我见过有人把整型DEM和浮点型DEM放在一起镶嵌输出选整型结果后续做填洼分析时微小洼地全部识别不出来。所以在设置输出像元类型时务必先了解所有源数据的像元类型再选择一个能容纳全部数据精度的类型。5. 镶嵌完成不等于结束接缝检查与异常修复镶嵌工具跑完不代表能直接拿去用。凡是做过几次镶嵌的人都会同意真正的功夫在于镶嵌后的质量检查与异常修复。这一节分享我自己的质检流程和踩过的坑。5.1 我自己常用的质检三连拉伸显示、直方图、地形剖面我给自己定了一套“质检三连”每次镶嵌完都强制走一遍基本能抓住九成以上的问题。第一拉伸显示检查。把镶嵌结果加载到ArcMap或QGIS里符号化方式设为拉伸Stretch色带选一个地形色带然后缩放到全图范围。先整体扫一遍重点看接缝区域有没有明显的颜色突变。接着放大到每一个原始图幅的边界附近沿着边界走一圈。正常的镶嵌结果从大局到细节都应该是连续过渡的。如果某处颜色明显变化成一个台阶带十有八九是接缝问题。第二直方图检查。打开栅格图层的属性查看直方图分布。正常DEM的直方图应该像一个单峰或双峰的平滑曲线集中在某个高程区间。如果直方图两端出现极值毛刺比如-9999或者一个巨大的正值说明NoData没有处理好或者源数据存在异常值。这时可以进一步用栅格计算器筛选出这些异常像元的空间分布定位具体位置。直方图检查还能发现一个常见问题——镶嵌后某些碎图高程整体比其他碎图高一截直方图会出现“双峰”一个峰是这部分数据另一个峰是另一部分数据结合空间分布能判断是否存在不同测区数据混用。第三地形剖面检查。在ArcScene或ArcGIS Pro的3D环境下用剖面线工具在接缝附近画一条横向穿越剖面线观察高程折线是否光滑。如果剖面线在接缝位置出现台阶或锯齿说明重叠区域融合没做好。QGIS里也有类似的地形剖面工具插件可以使用。剖面检查是发现“隐蔽折痕”的利器比肉眼盯色带可靠得多。5.2 黑边和折痕的常见成因与处理办法说到修问题最常见的两类问题是黑边和折痕。黑边白边同理的直接成因就是NoData值没设置好。解决思路是老生常谈的“先设置再镶嵌”但如果已经生成了结果怎么办有两个补救方法。 一个是在ArcGIS里用Copy Raster工具把结果重新拷贝一份拷贝时指定NoData Value另一个是用栅格计算器把异常值所在的区域重新赋值为NoData。例如Con(rast -9999, -9999, rast)这句会保留有效的高程值只把等于-9999的像元继续标记为NoData。需要注意后续处理中要把NoData值设置为正确识别否则问题依旧。折痕问题稍微复杂。它可能来自重叠区融合规则选得不对也可能来自源数据接缝处本来就存在高程不连续。如果砖已经烧好了非要在结果上修折痕最现实的办法是利用局部插值方法来处理接缝带比如取折痕两侧一定缓冲区内的高程值构建一个局部地表面把折痕区域重新插值平滑掉。但这个办法有风险处理不当会引入新的地形假象所以我在实际工作中更倾向于回到原始碎图检查接缝两侧是否需要先做水平校或垂直校再重新镶嵌。很多时候接缝处的微小偏移其实是坐标系转换引起的把源数据放到同一基准下重新镶嵌问题自然就消失了。5.3 大范围数据的内存问题和分段拼法最后聊聊大范围镶嵌时的硬件瓶颈。几十幅碎图一次性镶嵌像元数量动辄上亿内存和磁盘IO压力都很大。我曾有一次拼接一个地级市范围的高程数据二十多幅1米分辨率DEM直接在Mosaic To New Raster里全选运行跑了一个小时输出到一半报“999999Error executing function”一看是临时磁盘空间不足气到无语只能清磁盘重来。如果你也遇到类似问题我的建议是分段拼接。比如把碎图按行或按列分成若干组每组先镶嵌成一个中间成果最后再把中间成果镶嵌成最终文件。这样每一步处理的像元数量可控内存占用低得多出错后也方便定位是哪个环节的问题。同时为了保证每个分段的接缝规则一致所有中间成果的Mosaic Operator都要设置成统一的规则不要一半用Blend一半用First。另一个实用小技巧是先输出一个低分辨率版本的预拼接结果比如设置输出像元大小是最终目标的两倍用来快速检查空间关系和颜色过渡。等确认无误后再用全分辨率正式跑一遍。这个做法能帮你省下大量因为返工而浪费的时间。最后再分享一个我个人的小习惯处理完高程碎图镶嵌后把源文件的投影信息、NoData值、分辨率、输出参数记在一个TXT文档里和数据放在同一个文件夹。这个项目做完过几个月甲方要追加范围重新下载一批相邻碎图这时候你只需要翻出当时的参数记录照着重跑一遍新老数据就能严丝合缝地接上。否则单靠记忆去回忆当初设的像元类型和融合规则大概率会踩回之前踩过的坑。
返回列表