简介面向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 产品按地面采样间隔大致分三档每档的能力上限差别很大。下面这张表是我自己选数据时的对照依据不是产品说明书是踩过坑之后总结的取舍线。分辨率典型产品单幅覆盖全球体量能干什么干不了什么1 弧秒约 30 米SRTM 1 arc-second、ASTER GDEM、ALOS AW3D30、Copernicus DEM 30m1°×1°百 GB 以上小流域水文、坡度坡向、工程选址、淹没模拟城市级精细建模、单体建筑周边分析3 弧秒约 90 米SRTM 3 arc-second、GMTED2010 7.5 弧秒1°×1° 或 5°×5°十几 GB省域尺度地形概览、全国性制图小于 1 平方公里的汇流路径30 弧秒约 1 公里GMTED2010 30 弧秒、GTOPO30全球整幅几百 MB全球气候模式、宏观水系提取任何需要局部精度的分析选型的关键不是“越精细越好”而是你的分析尺度必须比像元大一个量级以上。做一条长度 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 21 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 -5unzip -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 COMPRESSDEFLATE \ -co PREDICTOR2 \ -co TILEDYES \ -co BIGTIFFIF_SAFER \ mosaic.vrt mosaic.tifPREDICTOR2对高程这种连续型数据很有效通常能再压下去两到三成。TILEDYES让文件内部按 256×256 分块存储后面做局部读取时性能提升明显。BIGTIFFIF_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_THREADSALL_CPUS \ -co TILEDYES -co COMPRESSDEFLATE -co BIGTIFFIF_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 rD:\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 TILEDYES并且分块处理——先按经纬度切成几块每块单独处理完再拼把内存峰值控在单块的量级。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, encodingutf-8) as f: writer csv.DictWriter(f, fieldnamesrows[0].keys()) writer.writeheader() writer.writerows(rows) print(findexed {len(rows)} tiles)这段代码的关键点在ds.bounds它直接给出每幅瓦片的外接矩形比解析文件名可靠得多——文件名里的经纬度有时是瓦片左上角有时是中心点不同产品约定不一样靠字符串截取早晚要出错。ds.res给出像元大小能顺手筛出分辨率混进来的异常文件。跑完之后同一目录下再出现重复下载或者不同分辨率混放的情况打开 CSV 按res_x排一下序就能发现。索引建好之后按范围取数的逻辑就三行import csv def pick_tiles(lon_min, lon_max, lat_min, lat_max, indextile_index.csv): with open(index, encodingutf-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 没设。从那以后我每次拿到新的地形数据包第一件事就是建索引、跑巡检、抽三个点核高程三步走完再动分析工具。希望帮到你。本文还有配套的精品资源点击获取
简介面向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 产品按地面采样间隔大致分三档每档的能力上限差别很大。下面这张表是我自己选数据时的对照依据不是产品说明书是踩过坑之后总结的取舍线。分辨率典型产品单幅覆盖全球体量能干什么干不了什么1 弧秒约 30 米SRTM 1 arc-second、ASTER GDEM、ALOS AW3D30、Copernicus DEM 30m1°×1°百 GB 以上小流域水文、坡度坡向、工程选址、淹没模拟城市级精细建模、单体建筑周边分析3 弧秒约 90 米SRTM 3 arc-second、GMTED2010 7.5 弧秒1°×1° 或 5°×5°十几 GB省域尺度地形概览、全国性制图小于 1 平方公里的汇流路径30 弧秒约 1 公里GMTED2010 30 弧秒、GTOPO30全球整幅几百 MB全球气候模式、宏观水系提取任何需要局部精度的分析选型的关键不是“越精细越好”而是你的分析尺度必须比像元大一个量级以上。做一条长度 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 21 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 -5unzip -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 COMPRESSDEFLATE \ -co PREDICTOR2 \ -co TILEDYES \ -co BIGTIFFIF_SAFER \ mosaic.vrt mosaic.tifPREDICTOR2对高程这种连续型数据很有效通常能再压下去两到三成。TILEDYES让文件内部按 256×256 分块存储后面做局部读取时性能提升明显。BIGTIFFIF_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_THREADSALL_CPUS \ -co TILEDYES -co COMPRESSDEFLATE -co BIGTIFFIF_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 rD:\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 TILEDYES并且分块处理——先按经纬度切成几块每块单独处理完再拼把内存峰值控在单块的量级。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, encodingutf-8) as f: writer csv.DictWriter(f, fieldnamesrows[0].keys()) writer.writeheader() writer.writerows(rows) print(findexed {len(rows)} tiles)这段代码的关键点在ds.bounds它直接给出每幅瓦片的外接矩形比解析文件名可靠得多——文件名里的经纬度有时是瓦片左上角有时是中心点不同产品约定不一样靠字符串截取早晚要出错。ds.res给出像元大小能顺手筛出分辨率混进来的异常文件。跑完之后同一目录下再出现重复下载或者不同分辨率混放的情况打开 CSV 按res_x排一下序就能发现。索引建好之后按范围取数的逻辑就三行import csv def pick_tiles(lon_min, lon_max, lat_min, lat_max, indextile_index.csv): with open(index, encodingutf-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 没设。从那以后我每次拿到新的地形数据包第一件事就是建索引、跑巡检、抽三个点核高程三步走完再动分析工具。希望帮到你。本文还有配套的精品资源点击获取
