简介：面向GIS从业者和初学者的全球DEM地形数据获取指引文档，针对全球高程数据下载渠道分散、速度不稳等问题，整理了可自行选择的百度网盘与天翼云盘入口，覆盖全球陆地区域，适用于地形分析、坡度计算、水文与城市规划等场景。资源包共1个docx文件，约12KB，篇幅精简但信息集中，除DEM数据链接外，还汇总了免费GIS软件与GIS学习教程入口，便于按需查阅。已有6002人学习或下载，说明其在相关人群中具备一定参考价值。读者可据此获取全球地形数据，并衔接QGIS、GRASS GIS、MapInfo等免费工具进行读取、处理与分析，也可借助教程了解数据加载、空间处理和制图流程，减少四处搜集资料的时间。

1. 全球 DEM 数据包：一份能直接开工的高程底图

做水文模拟或者选址分析的同行，十有八九在同一个地方卡过壳：全球 DEM 高程数据去哪儿弄。真正麻烦的是找一圈之后发现，要么分辨率对不上，要么下回来是一堆需要自己拼的零散瓦片，要么干脆是坐标系不匹配的半成品。手上这份资源包按三块分好了：全球地形数据、免费 GIS 软件、一套 GIS 教程，云盘直取，不用注册也不用做任务。它的价值不在“有数据”，而在于替你省掉了在几个数据站点之间反复试错的那两天。适合三类人：刚接手地形分析任务、手上还没有底图的 GIS 从业者；预算有限只能用开源工具的独立开发者；以及想在 ArcGIS 里快速验证一个坡度或汇流思路、不想先花半天搞定数据源的人。

2. 全球 DEM 的分辨率与基准：选错一档，后面的地形数据全得重来

拿到一个几十 GB 的地形数据包，最容易犯的错是“先解压再说”。我见过太多人把 90 米数据当 30 米用，做完坡度分析发现沟道全是锯齿；也见过拼完一整块区域，跟实测水准点一比差了四十多米。这些问题都不是软件操作层面的，全出在选型和基准上。这一章先把数据本身的属性捋清楚，再谈怎么把文件安稳落到本地。

2.1 三种主流分辨率，各自的适用边界

全球 DEM 产品按地面采样间隔大致分三档，每档的能力上限差别很大。下面这张表是我自己选数据时的对照依据，不是产品说明书，是踩过坑之后总结的取舍线。

选型的关键不是“越精细越好”，而是你的分析尺度必须比像元大一个量级以上。做一条长度 5 公里的沟道汇流，30 米数据够用；做某个县全域的坡度分级统计，90 米完全能交差，还能把内存占用降一个数量级。反过来，如果你要算的是几十公顷范围的土方量，30 米数据的每个像元代表 900 平方米，误差会大到没有参考意义，这时候该换的是实测点云或者无人机成果，再高的全球产品也救不回来。

还有一层容易被忽略：同样是 30 米，不同产品的高程来源不一样。雷达干涉测量出来的产品在植被茂密区会把树冠高度算进去，光学立体像对在云覆盖区会有空洞。资源包里如果同时给了几种产品，别急着只留一个，先用同一小块区域叠一下，看差异分布，这个动作花不了二十分钟，但能避免后面返工。

2.2 下载前的三个确认项

动手取数据之前，有三件事必须确认，顺序不能颠倒。

第一是范围。全球数据是分幅存放的，命名一般带经纬度，比如 N30E110 表示北纬 30 度、东经 110 度这一幅。你要做的是先在地图上框出研究区的外接矩形，算出覆盖的纬度带和经度带，再去核对文件清单。直接整包下载再筛，下 200 GB 用 20 GB，时间和磁盘都是浪费。

第二是垂直基准。这是最隐蔽的一个坑。部分产品的高程是相对于 EGM96 大地水准面的正高，另一部分是相对于 WGS84 椭球的大地高，两者在中国境内能差二三十米，而且这个差值随位置变化。如果你的分析要跟实测水准点、水文站记录对比，不统一基准，结果就是系统性偏移，而且偏移量还不是常数，没法简单加个固定值修掉。资源包里如果只有一种产品，问题不大；混用两种以上，必须先把基准统一到同一个参考面。

第三是 NoData 值。不同产品的空洞标记五花八门，有 -32768、-9999、0，也有直接留 0 但没标 NoData 的。这个值要是不提前查清楚，后面做镶嵌、算坡度的时候，这些格子会被当成真实高程参与运算，出来一堆诡异的悬崖。

2.3 从取件到本地校验：一个能省下三小时的流程

资源给了两个取件入口。一个是覆盖面最广的公共网盘，好处是什么设备都能开，坏处是限速明显，几十 GB 的东西挂一晚上可能还没过半；另一个是企业级网盘，实测能拉满带宽，同样的包一两个小时就下来了。我的建议很直接：只要链接能打开，优先走快的那个，别为了省事在慢的那个上耗一整天。大文件下载中断的概率跟时长正相关，越慢越容易断，断了还得从头来。

下完之后别急着解压，先做完整性校验。这一步能挡掉大部分“文件损坏”类的玄学问题。

# 1. 先看清单和体积，确认包的数量对得上 ls -lh ./dem_pack/ du -sh ./dem_pack/ # 2. 压缩包先做完整性测试，不产生解压文件 for f in ./dem_pack/*.zip; do unzip -t "$f" >/dev/null 2>&1 && echo "OK $f" || echo "FAIL $f" done # 3. 解压到独立目录，别和已有数据混在一起 mkdir -p ./dem_raw unzip -q ./dem_pack/global_dem_30m.zip -d ./dem_raw # 4. 统计瓦片数量，核对命名规律 find ./dem_raw -name "*.tif" | wc -l find ./dem_raw -name "*.tif" | sort | head -5 find ./dem_raw -name "*.tif" | sort | tail -5

unzip -t会逐条读取压缩包的校验和，不做实际写盘，几十 GB 的包跑几分钟就有结论，比解压到一半报错强得多。第四步的排序输出很关键：全球分幅数据的文件名是按经纬度顺序排的，看一眼头尾就知道纬度带是不是从南到北连续、经度带有没有缺号，缺的那几幅就是下载时漏掉的。

目录结构上，我习惯分成dem_raw（原始瓦片，只读不改）、dem_work（镶嵌、裁剪的中间成果）、dem_out（最终交付）三层。原始数据这一层永远不动，后面任何处理出问题，删掉工作目录重来就行，不用重新下载。这个习惯值不了多少钱，但能救你几次心态。

3. QGIS 跑通 GIS 链路：拼接、重投影、坡度出图

数据落到本地之后，第一件事是把它变成“一块完整的地形数据”，而不是几百个碎瓦片。开源工具链里，QGIS 负责看和交互，真正干重活的是底层的 GDAL。很多人习惯在 QGIS 菜单里点“栅格 - 杂项 - 合并”，几十幅还行，上百幅就容易卡死。命令行这条路更快也更可控。

3.1 QGIS 装好之后先改三个默认设置

资源包里的免费 GIS 软件装完之后，别急着加载数据，先花五分钟改三个设置，后面能少很多麻烦。

第一个是坐标系显示。默认情况下 QGIS 会做“即时重投影”，也就是把不同坐标系的数据动态投到同一个显示坐标系上。看图没问题，但如果你直接拿显示结果去量算距离或面积，得到的是投影后的值，跟数据本身的坐标单位可能对不上。我的做法是把工程坐标系固定成研究区所在的投影带，让所有数据都显式转换，不依赖动态投影。

第二个是临时文件目录。GDAL 处理大栅格时会在临时目录写缓存，默认位置经常在系统盘，处理全球数据时几十 GB 的缓存能把 C 盘写满。在设置里把临时目录指到空间充足的机械盘或者数据盘，这个问题就消失了。

第三个是栅格金字塔的构建策略。加载大影像时 QGIS 会问要不要建金字塔，选“是”并勾选“不再询问”。金字塔是数据的一个缩略副本，能让缩放和平移流畅很多；不建的话，每次拖动视图都在重采样全分辨率数据，几 GB 的影像能把界面卡到无响应。

3.2 拼接：gdalbuildvrt 比 Merge 更稳的原因

拼接这件事，直接上手gdal_merge.py是常见的误用。它会把所有输入读进内存再写出，几幅没问题，上百幅 30 米瓦片就是内存溢出的标准剧本，而且跑到一半被系统杀掉，前面的时间全白费。

正确做法是用gdalbuildvrt先建一个虚拟镶嵌。VRT 文件本身只有几十 KB，里面记的是每幅瓦片的位置索引，不复制任何像素。后续的所有操作都读这个索引，GDAL 会按需分块读取，内存占用跟数据总量基本脱钩。

# 先生成文件清单，用 sort 保证顺序稳定 find ./dem_raw -name "*.tif" | sort > tile_list.txt wc -l tile_list.txt # 建虚拟镶嵌：分辨率取最高，统一 NoData gdalbuildvrt \ -input_file_list tile_list.txt \ -resolution highest \ -srcnodata -32768 \ -vrtnodata -32768 \ mosaic.vrt # 检查镶嵌结果的范围、波段数、统计值 gdalinfo -stats mosaic.vrt | head -40

-input_file_list是清单模式，比在命令行里列上百个路径干净得多，也避开了 shell 的参数长度限制。-resolution highest的意思是当输入瓦片分辨率不一致时取最精细的那一档，如果你的数据本来就是同一产品，这个参数不影响结果，但保留它是个保险。-srcnodata和-vrtnodata要设成同一个值，前者告诉 GDAL 输入里的空洞标记是什么，后者决定输出里继续用哪个值标记空洞。两处不一致，接缝处就会冒出黑边或白边。

gdalinfo -stats会扫描一遍数据算最小最大值和均值，这一步在大数据上要花点时间，但值得。看最小值是不是那个负的 NoData 值、最大值是不是在合理的高程区间内（地球陆地最高约 8848 米），基本就能判断镶嵌有没有出岔子。

虚拟镶嵌只适合中间处理，不适合交付。要落到实体文件，用gdal_translate转一次，同时加上分块和压缩：

gdal_translate \ -co COMPRESS=DEFLATE \ -co PREDICTOR=2 \ -co TILED=YES \ -co BIGTIFF=IF_SAFER \ mosaic.vrt mosaic.tif

PREDICTOR=2对高程这种连续型数据很有效，通常能再压下去两到三成。TILED=YES让文件内部按 256×256 分块存储，后面做局部读取时性能提升明显。BIGTIFF=IF_SAFER是防止文件超过 4 GB 时格式报错，全球数据必加。

3.3 重投影与裁剪：gdalwarp 的四个关键参数

拿到镶嵌好的全球或大区数据，下一步是裁出研究区并转到投影坐标系。地理坐标系（经纬度）下的像元大小单位是度，直接拿去算坡度、算面积都会出问题，必须转成米制投影。

gdalwarp -overwrite \ -t_srs EPSG:32650 \ -tr 30 30 \ -r bilinear \ -te 500000 3000000 800000 3300000 \ -te_srs EPSG:32650 \ -multi -wo NUM_THREADS=ALL_CPUS \ -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=IF_SAFER \ mosaic.vrt study_area.tif

-t_srs是目标坐标系。EPSG 编码按研究区所在的 UTM 带选，北半球 326 开头，南半球 327 开头，后面的数字是带号。选错带的直接后果是形状畸变，面积量算偏差能到百分之几。

-tr 30 30指定输出像元大小，单位跟随目标坐标系，这里是米。从地理坐标转过来时，如果不指定，GDAL 会按输入推算一个近似的米制分辨率，结果往往不是整数，后面做多源数据对齐很麻烦。显式给 30 米，输出就是规整的网格。

-r bilinear是重采样算法。高程数据我从不用最邻近，因为它会把地形弄成台阶状；双线性在平滑和保真之间平衡得比较好。如果是分类栅格（比如土地利用），那必须用最邻近，不能插值出新类别。这个区别值得记一下。

-te加-te_srs是裁剪范围。四个数字的顺序是xmin ymin xmax ymax，配合-te_srs明确这个范围的坐标系，避免 GDAL 猜错。范围可以直接从研究区矢量数据的边界框读出来，或者用ogrinfo查。

3.4 坡度、坡向、山影：gdaldem 的算法选择

地形分析的三个基础衍生量都在gdaldem里，关键是算法和单位这两个参数别选错。

# 投影坐标系下算坡度，输出百分比 gdaldem slope study_area.tif slope_pct.tif \ -alg Horn -p -compute_edges # 经纬度原始数据必须加 -s，否则坡度全是 89.9 gdaldem slope mosaic.vrt slope_deg.tif \ -alg Horn -s 111120 # 坡向：-zero_for_flat 让平坦区输出 0 而不是 -9999 gdaldem aspect study_area.tif aspect.tif \ -zero_for_flat -compute_edges # 山影：方位角 315、高度角 45 是制图惯例 gdaldem hillshade study_area.tif hillshade.tif \ -az 315 -alt 45 -z 1.5 -compute_edges

-alg有两个选项，Horn 和 ZevenbergenThorne。Horn 是 3×3 窗口的加权差分，对噪声更宽容，是多数场景的默认选择；ZevenbergenThorne 只用四邻域，对陡崖和线状地物的响应更锐利，适合地形起伏剧烈的区域。选哪个没有绝对对错，但对同一批数据要保持一致，不然两块区域的坡度没法放在一起统计。

-p表示输出百分比坡度而不是度数。度数在 0 到 90 之间线性分布，百分比在陡坡段会迅速放大，视觉上更容易看出危险区。做工程边坡分级时我用百分比，做日照和辐射分析时用度数，因为后续公式里的三角函数吃的是角度。

-s这个参数是经纬度数据的专属参数。地理坐标系下，一个像元的横向距离随纬度变化，GDAL 需要一个“垂直单位与水平单位的比值”才能算出真实坡度。全球平均取 111120（一度纬度约 111.12 公里）。不加这个参数，算出来的坡度会趋近 90 度，整个图面一片白，这属于最经典的坑之一。

-compute_edges会让边缘像元也参与计算，不加的话数据外框会有一圈 NoData。做镶嵌拼接时这圈空洞很烦人，所以只要后续还要拼接，就加上。

4. ArcGIS 里的地形分析：填洼、流向、坡度的正确顺序

QGIS 那条链路适合批量处理和数据准备，ArcGIS 这边的优势在于水文分析工具的完整性和可视化表达。但工具全不代表顺序可以乱，同一个数据，四个工具换个先后，结果能差出十倍。这一段把顺序和参数讲透。

4.1 填洼为什么必须排在流向之前

DEM 是栅格化的地表采样，每个像元取一个代表值。这种离散化会在局部造出一批“假洼地”：四周都比它高，但它并不是真实的湖。这类像元通常来自雷达数据的噪声、或者两条等高线之间的平滑插值误差。

如果不填洼直接算流向，水流到这些假洼地就断了，算出来的汇流累积量在洼地出口处会突然归零，下游整条河的网络全断。更麻烦的是这种断裂在上万平方公里的区域里非常零散，人工看不出规律，只会觉得“结果不对但不知道哪不对”。

填洼的代价是改变地形。所以填的时候要留个心眼：填洼工具一般有个“允许的最大填充深度”参数，设得太小填不干净，设得太大可能把真实的小型洼地也抹掉。我的习惯是先不设限跑一遍，看填了多少像元，如果填掉的面积占比不到千分之一，说明数据质量不错；超过百分之一，就该回头查数据源是不是有系统性问题。

流向算法有两个主流选择：D8 和 D-infinity。D8 把每个像元的水流指派给八个邻域中坡降最大的那一个，简单快速，产生的河网清晰锐利，是提取河道的主流做法。D-infinity 允许水流按角度分配给两个相邻像元，更适合计算坡面漫流和土壤侵蚀，但河网会显得模糊。做河网提取我用 D8，做泥沙输移用 D-infinity。

4.2 坡度工具的 Z 因子与单位

ArcGIS 的坡度工具参数不多，但每一个都能让结果翻车。

第一个是输出单位。DEGREE输出 0 到 90 的度数，PERCENT_RISE输出百分比。这两个不是简单换算关系，百分比在 45 度时正好是 100%，超过 45 度后继续放大，90 度时会趋向无穷。选哪个取决于下游怎么用，别混着来。

第二个是 Z 因子。这个参数的意思是“垂直方向的一个单位等于水平方向的多少个单位”。如果数据是投影坐标系、单位是米，高程单位也是米，Z 因子填 1。如果数据是地理坐标系、单位是度，那垂直的米和水平的度根本不是一个量纲，必须换算。纬度为 φ 的地方，一个度大约对应 111320×cos(φ) 米，所以 Z 因子约等于 1/(111320×cos(φ))。

这个换算因子随纬度变化，在北纬 30 度附近约 0.0000104，在北纬 60 度附近要翻倍。有人在赤道附近调好的参数，拿到高纬度直接用，坡度就系统性偏大。这也是我不推荐在经纬度数据上直接做地形分析的原因：与其记一堆随位置变化的系数，不如先投影一次，把问题一次性消掉。

4.3 山影的方位角、高度角与垂直夸张

山影不是分析量，是表达量，但它的参数设置直接影响读图的人对地形的判断。

方位角默认 315 度，也就是光从西北方向照过来。这个惯例来自人类对自然光照的直觉：北半球的人习惯光从左上角打过来，倒过来会让人觉得地形凹凸反转，山脊看成沟谷。除非有特殊表达需求，别改这个值。

高度角默认 45 度，模拟太阳在天空中的高度。角度太低，阴影拉得很长，山脊线会很夸张但谷底一片黑；角度太高，阴影对比弱，地形起伏看不出来。做区域概览图我用 45 度，做峡谷、陡崖这类微地形我用 30 度左右。

垂直夸张系数 Z 默认是 1。在地形起伏平缓的平原区，把 Z 调到 2 到 3 能让微地形显出来，但那已经不叫真实比例了，图例上必须注明。我见过一份成果图把 Z 设成 5，读者以为那里是山区，实际高差不到 50 米。

4.4 用 arcpy 把四步串成批处理

研究区大的时候，一个个点工具不现实。arcpy 的写法比模型构建器更可控，出错也能定位到具体哪一步。

import arcpy from arcpy.sa import Fill, FlowDirection, FlowAccumulation, Slope, Hillshade arcpy.env.workspace = r"D:\gis\dem\study" arcpy.env.overwriteOutput = True arcpy.env.parallelProcessingFactor = "100%" arcpy.CheckOutExtension("Spatial") dem = "study_area.tif" # 1. 填洼，先把假洼地处理掉 filled = Fill(dem) filled.save("dem_fill.tif") # 2. 流向，D8 算法，输出编码 1-128 fdir = FlowDirection("dem_fill.tif", "NORMAL") fdir.save("flow_dir.tif") # 3. 汇流累积，用浮点避免大流域溢出 facc = FlowAccumulation("flow_dir.tif", None, "FLOAT") facc.save("flow_acc.tif") # 4. 坡度，投影数据 Z 因子为 1 slope = Slope("dem_fill.tif", "PERCENT_RISE", 1) slope.save("slope_pct.tif") # 5. 山影，方位角 315、高度角 45、垂直夸张 1.5 hs = Hillshade("dem_fill.tif", 315, 45, 1.5) hs.save("hillshade.tif") arcpy.CheckInExtension("Spatial") print("done")

几个参数值得单独说。FlowDirection的第二个参数是强制流向选项，NORMAL表示按最陡下降计算，FORCE会强行把所有边缘像元的水流导向外部，只在研究区边界确实有出水口时用。默认的NORMAL更安全。

FlowAccumulation的第三个参数是数据类型，FLOAT支持大数值，INTEGER更省空间但流域面积一大就溢出成负数，那个负值会一直传下去，最后统计出来的河道长度全是错的。宁可多占点磁盘。

坡度用的是dem_fill.tif而不是原始 DEM，这是有意的。填洼只改了极少数像元，但对坡度的连续性有好处，尤其是洼地边缘那一圈。如果研究重点是微地形，也可以直接用原始数据对比一次，看差异是否可接受。

parallelProcessingFactor设成 100% 会吃满所有核，跑得快但机器会卡，本地作业时按需调整。批处理脚本里我通常设成 50%，留点余量给其他程序。

5. 排查清单：DEM 处理里最容易翻车的五个地方

下面这五条基本都是我在实际项目里撞过的，每条按“看到什么现象、根因在哪、怎么修”的顺序展开。遇到问题时按顺序排查，能省掉大量瞎试的时间。

5.1 拼完跟实测点差了二三十米

现象：整个研究区的高程值系统性地偏高或偏低，而且这个差值在不同位置还不一样，没法用一个固定偏移量修掉。

原因：不同来源的 DEM 用了不同的垂直基准。一部分产品输出的是正高（相对大地水准面），另一部分输出的是大地高（相对参考椭球），两者之间的差距就是高程异常，在一个国家范围内可以从负几十米变到正几十米，随位置连续变化。

解决：先在元数据或gdalinfo的输出里确认每个产品的垂直基准标识。如果是混用，用gdalwarp配合高程基准转换网格做一次统一，或者在 GIS 软件里用基准转换工具处理。实在查不到元数据，就找几个已知高程的控制点做比对，拟合出一个趋势面再从数据里减掉，虽然不够严谨，但比带着系统性偏差往下做要好。

5.2 拼接或算坡度时进程被杀

现象：命令行跑到一半直接退出，没有报错信息；或者系统提示内存不足，整个桌面卡住。

原因：用了会全量加载的工具，比如把上百幅瓦片直接喂给合并工具。30 米分辨率的全球数据，铺满地球陆地是百 GB 量级，解压成内存里的浮点数组还要再翻几倍，任何机器都扛不住。

解决：换用gdalbuildvrt建虚拟镶嵌，后续操作全部基于 VRT 做分块读写。如果必须生成实体文件，用gdal_translate加-co TILED=YES，并且分块处理——先按经纬度切成几块，每块单独处理完再拼，把内存峰值控在单块的量级。QGIS 里的操作同理，能用命令行就别用图形界面处理大数据。

5.3 坡度图一片白或者全是 0

现象：算出来的坡度值集中在 89 度以上，或者整幅图接近 0，但地形明明有起伏。

原因：单位没对齐。地理坐标系下，水平单位是度、垂直单位是米，两者量级差了五个数量级。不加比例因子，GDAL 会认为水平方向几乎没有距离变化，垂直方向却有几百米，算出来的坡降趋于无穷，坡度自然逼近 90 度。反过来，如果把投影数据的像元大小误填成 0.0008 而不是 30，坡度就会趋近 0。

解决：先看gdalinfo输出里的坐标系类型。地理坐标系必须加-s 111120（或者按研究区纬度精确算），投影坐标系不加。ArcGIS 里对应的是 Z 因子参数，逻辑一样。改完参数后拿一小块已知坡度的区域验算一下，比如一段 45 度的边坡，算出来应该在 100% 左右。

5.4 数据边缘出现黑边、白边或者一圈零值

现象：镶嵌结果里，瓦片接缝处有明显的一条暗线或亮线；裁剪后的成果图外框有一圈突兀的 0 值。

原因：NoData 没统一。输入瓦片用 -32768 标空洞，输出没指定，GDAL 就按默认的 0 处理，0 米在沿海地区是合法高程，在高原上就是明显的异常值。另一种情况是计算坡度时没加-compute_edges，边缘像元因为邻域不全被判为无效。

解决：gdalbuildvrt和gdalwarp都显式指定-srcnodata和-dstnodata，两处设成同一个值。gdaldem系列加上-compute_edges。已经生成的数据可以用gdal_translate -a_nodata -32768补一次声明，不用重算。显示层面的黑边有时只是渲染问题，在符号化里把 NoData 设为透明即可，别急着改数据。

5.5 下载到九成断掉，解压报校验错误

现象：大文件下载进度卡住不动，最后提示失败；或者解压到某个文件时报 CRC 错误，剩下的都解不出来。

原因：公共网盘对非会员限速厉害，几十 GB 的包要挂很久，链路一抖动就断。断了之后很多客户端不支持续传，只能从头来。压缩包不完整，解压时自然在断点处报错。

解决：优先走限速宽松的那个入口，同样的包能快一个数量级。下载前先确认客户端的断点续传是开着的。拿到包之后先跑unzip -t做完整性测试，别等解压到一半才发现。如果只是缺了几幅瓦片，根据命名规律单独补下那几幅就行，不用整包重来——这也是为什么第 2 章要你先看文件清单的头尾。

6. 建一张瓦片索引表：把全球高程数据变成随取随用的底图库

数据攒到一定量之后，真正的效率瓶颈就从“怎么下”变成了“怎么找”。几百上千幅瓦片躺在同一个目录里，每次新项目都要重新翻一遍，靠文件名猜覆盖范围，这个动作重复几十次就该自动化了。我的做法是给整个数据目录建一张索引表，把所有瓦片的范围、分辨率、NoData 值抽出来存成 CSV，以后按研究区范围查一下就能定位到具体文件。

import os, csv import rasterio rows = [] src_dir = "dem_raw" for name in sorted(os.listdir(src_dir)): if not name.lower().endswith((".tif", ".hgt")): continue path = os.path.join(src_dir, name) with rasterio.open(path) as ds: b = ds.bounds rows.append({ "file": name, "lon_min": round(b.left, 4), "lon_max": round(b.right, 4), "lat_min": round(b.bottom, 4), "lat_max": round(b.top, 4), "res_x": round(ds.res[0], 8), "res_y": round(ds.res[1], 8), "nodata": ds.nodata, "crs": str(ds.crs), }) with open("tile_index.csv", "w", newline="", encoding="utf-8") as f: writer = csv.DictWriter(f, fieldnames=rows[0].keys()) writer.writeheader() writer.writerows(rows) print(f"indexed {len(rows)} tiles")

这段代码的关键点在ds.bounds，它直接给出每幅瓦片的外接矩形，比解析文件名可靠得多——文件名里的经纬度有时是瓦片左上角，有时是中心点，不同产品约定不一样，靠字符串截取早晚要出错。ds.res给出像元大小，能顺手筛出分辨率混进来的异常文件。跑完之后，同一目录下再出现重复下载或者不同分辨率混放的情况，打开 CSV 按res_x排一下序就能发现。

索引建好之后，按范围取数的逻辑就三行：

import csv def pick_tiles(lon_min, lon_max, lat_min, lat_max, index="tile_index.csv"): with open(index, encoding="utf-8") as f: rows = list(csv.DictReader(f)) return [ r["file"] for r in rows if float(r["lon_max"]) > lon_min and float(r["lon_min"]) < lon_max and float(r["lat_max"]) > lat_min and float(r["lat_min"]) < lat_max ] files = pick_tiles(108.5, 112.0, 28.0, 31.5) print(len(files), "tiles selected")

选出来的文件名列表可以直接喂给gdalbuildvrt -input_file_list，整条链路就接上了：查索引、出清单、建虚拟镶嵌、裁剪、出坡度，全程不需要打开图形界面。

索引表还有个附带好处是能做数据质量巡检。按nodata字段分组，看看有没有瓦片用了跟其他不一样的空洞标记；按crs分组，检查坐标系是否统一。这两类问题在单幅数据上看不出来，一旦拼起来就是接缝处莫名其妙的暗线。

验证环节我习惯抽三个点做交叉核对：从成果栅格上取一个像元值，跟原始瓦片同位置的像元值比，确认镶嵌没有引入偏移；再找一个公开的已知高程点比一次，确认基准没问题；最后看一眼坡度直方图，正常地形应该是单峰偏态分布，如果出现两个峰，多半是两块区域的数据源不一样，混进同一张图里了。

血泪经验是这套流程一定要在项目开始时就跑一遍，别等做到一半才发现基准不统一、NoData 没设。从那以后我每次拿到新的地形数据包，第一件事就是建索引、跑巡检、抽三个点核高程，三步走完再动分析工具。希望帮到你。

本文还有配套的精品资源，点击获取