简介亚洲DEM文件东亚DEM是一份面向GIS初学者与地理信息从业人员的数字高程模型数据包覆盖中国、日本、韩国、蒙古等东亚区域可用于ArcGIS、MapGIS环境下的地形分析、坡度坡向计算、洪水模拟与地表建模等场景帮助读者在缺少实测地形数据时快速搭建分析底图。压缩包共8个文件约49.45MB其中7个为adf格式的ArcGIS栅格交换文件承担高程格网、投影信息、边界与属性表等核心数据存储另有1个xml元数据文件记录数据来源与坐标系等说明便于软件识别与后续管理。资料已有916人学习下载适合作为DEM加载、坐标系转换、重采样与剖面绘制等操作的练习素材。数据按区域分块存放单块体量适中便于分批加载与局部区域研究也可与土地利用、气候等图层叠加做综合分析。1. 亚洲DEM文件到手之后为什么“能打开”和“能用”是两回事把一份亚洲DEM文件拖进桌面 GIS颜色映射一开青藏高原的轮廓、中亚的盆地、沿海的平原都出来了看着挺像那么回事。可真拿它算坡度、做晕渲、切地形瓦片问题就来了山脊上一圈圈锯齿两块分幅的交界处横着一条台阶海面渲染出来是一片灰白噪点剖面线一拉像楼梯。这不是软件的问题是“能打开”和“能用”之间的差距——分幅没拼、NoData 没声明、投影没换、重采样算法没选对任何一个环节漏掉后面的成果都不可信。这篇把自己处理亚洲DEM文件的整条链路摊开从分辨率和格式怎么挑到上千个瓦片怎么拼成一张图再到跨带投影怎么换、坑在哪最后落到一条能交付的地形底图。2. 分辨率与格式的取舍亚洲DEM文件该选90米、30米还是12米选数据这件事很多人是反过来做的——先下载到哪个用哪个等做到一半发现分辨率不够、空洞太多再回头换源前面的拼接和切片全白干。亚洲这块地方的地形跨度太大既有平均海拔四千米以上的高原也有大面积低缓的冲积平原和群岛同一档分辨率在不同子区域的可用性差别很明显。先把分辨率账算清楚再谈格式。2.1 三档分辨率在亚洲的实际可用性对比公开可分发的全球/区域DEM主流就是90米级、30米级、12米级三档。它们不是简单的“越精细越好”而是各有各的失效区。下表是我自己在选源时习惯先填一遍的对照表数字取自各产品公开的标称指标实际使用要以你手里那份的具体元数据为准档位典型覆盖标称垂直精度单幅体量1°×1°适合做的事主要失效场景90米级合成近全球含高纬约 ±16 m20~30 MB大流域分析、初版底图、宏观选址山脊被抹平窄谷消失30米级雷达测绘南北纬60°以内约 ±10~16 m20~25 MB区域地形分析、坡度、汇流高山阴影区成片空洞30米级光学立体陆地大部分约 ±5~10 m40~60 MB精细地形、工程选线云覆盖、水体匹配差12米级光学立体陆地约 ±3~5 m100 MB 以上小流域、城郊尺度体量大、高频噪声多这张表最该看的是最后一列。90米级在山地区的“抹平”效果非常明显一条几百米宽的峡谷在90米网格里可能只有两三个像元汇流算法跑出来的河网会直接断掉。而30米级雷达测绘产品在高纬度和陡峭阴坡上的空洞如果不做填补坡度计算会在空洞边缘产生极大的伪梯度。我一般会先用90米级跑一遍全局确认研究范围里哪些区块需要更高精度再针对性补30米或12米而不是一上来就全亚洲铺12米——那个体量对绝大多数团队都是灾难。还有一个容易被忽略的点不同源的数据不能直接混拼。90米和30米拼在一起边缘会因为重采样产生系统性抬升或下沉接缝处形成一条肉眼可见的“堤”。如果确实要混用正确做法是先把低精度源重采样到高精度源的网格上标记为“低置信度图层”在后续分析里按掩膜降权而不是让它和高精度数据在同一张栅格上平起平坐。2.2 拿到文件先做元数据体检五个必看字段亚洲DEM文件最容易翻车的地方是元数据缺失或者不可信。很多分幅文件从压缩包里解出来坐标系是空的、NoData 是 None、像元大小是 0.000277777 这种奇怪的值。直接往下做错误会被一路放大。所以我习惯在拼接之前先跑一遍体检脚本把每个文件的五个字段打出来对一遍。import glob import rasterio import numpy as np # 体检输出CRS / 像元大小 / NoData / 数据类型 / 有效值范围 for path in sorted(glob.glob(raw_asia/*.tif))[:20]: with rasterio.open(path) as ds: # maskedTrue 让 NoData 参与掩膜避免被当成 0 米海拔 arr ds.read(1, maskedTrue) print( f{path}\n f CRS : {ds.crs}\n f RES : {ds.res}\n f NODATA : {ds.nodata}\n f DTYPE : {ds.dtypes[0]}\n f RANGE : {arr.min():.1f} ~ {arr.max():.1f} f有效像元占比 {arr.count() / arr.size:.2%} )这段代码的重点不在打印而在三个判断。第一ds.crs如果是None说明这个文件没有嵌入投影信息必须靠文件名或者配套的说明文件补上否则后面gdalbuildvrt会直接报错或者按错误坐标摆放。第二ds.nodata如果是None而数据里确实存在空洞maskedTrue也救不了你因为它不知道该把哪个值当无效值——这时得靠RANGE里的最小值去反推常见的填充值是 -32768、-9999、-99999看到这些就要警觉。第三有效像元占比低于 95% 的分幅说明空洞面积不小要不要用它得先看空洞分布别一股脑全拼进去。体检脚本建议全量跑不要只跑前 20 个。亚洲这么大范围分幅动辄上千混进一两个坐标系标错的文件拼完之后定位非常痛苦。2.3 从 HGT、BIL 到 Cloud-Optimized GeoTIFF格式统一与压缩原始分发格式五花八门HGT 是纯裸数据、无文件头坐标只能从文件名反推BIL/BSQ 是带一个 .hdr 说明的二进制IMG 是老式 Erdas 格式。这些格式在拼接阶段都能被 GDAL 读但读得慢、不支持分块读取、不支持压缩一旦数据量上去就是灾难。统一转成带内部瓦片和压缩的 GeoTIFF是后面所有流程的前提。# 单个分幅转换声明 NoData、加内部瓦片、DEFLATE 压缩、允许 BigTIFF gdal_translate raw_asia/N30E100.hgt clean/N30E100.tif \ -a_srs EPSG:4326 \ -a_nodata -32768 \ -co TILEDYES \ -co BLOCKXSIZE512 -co BLOCKYSIZE512 \ -co COMPRESSDEFLATE -co ZLEVEL6 \ -co BIGTIFFIF_SAFER # 批量用 shell 循环失败的文件名单独记下来 mkdir -p clean failed for f in raw_asia/*.hgt; do base$(basename $f .hgt) gdal_translate $f clean/${base}.tif \ -a_srs EPSG:4326 -a_nodata -32768 \ -co TILEDYES -co COMPRESSDEFLATE -co BIGTIFFIF_SAFER \ || echo $f failed/list.txt done参数逐个说清楚。-a_srs EPSG:4326是给无投影的裸数据“贴”上经纬度坐标HGT 必须加否则 GDAL 只会把它当无地理参考的图片。-a_nodata -32768把无效值显式写进文件头这一步是整个链路里最省后悔药的操作后面gdalbuildvrt、gdalwarp、gdaldem全都依赖它。TILEDYES配 512×512 的块大小是为了后续按窗口读取时只解压需要的那一块不设的话默认按行存储读一个 512×512 的小窗口也要扫过整行。COMPRESSDEFLATE在 DEM 这种连续表面上通常能压到原体积的 40%~60%ZLEVEL 取 6 是速度和压缩率的平衡点取 9 收益很小但慢很多。BIGTIFFIF_SAFER是保险单文件超过 4GB 时普通 TIFF 会直接损坏加了这个参数 GDAL 会自动切换到 BigTIFF代价只是极少数老软件读不了。转换完不要急着删原始文件。至少保留到拼接结果验收通过。我吃过一次亏批量转换时有一个分幅因为磁盘写满被截断gdal_translate返回了非零但循环里没中断结果拼出来的图上有一块 512×512 的矩形空缺找了半天才定位到是转换阶段的问题。所以上面那段脚本最后加了|| echo ... failed/list.txt失败清单必须留。3. 分幅拼接把上千个瓦片拼成一张亚洲图拼接是整个流程里最“体力活”的一步也是最容易埋雷的一步。表面上看把一堆 GeoTIFF 摞在一起就行实际上要处理三件事重叠区的取值优先级、不同分辨率混拼时的重采样、以及拼接过程的峰值内存。亚洲范围按 1°×1° 分幅经度方向一百多幅、纬度方向近百幅总计上万幅不夸张硬拼必然爆内存。3.1 瓦片重叠与接边高差优先级规则怎么定同一档分辨率的相邻分幅理论上边界应该严丝合缝实际上往往有一到两个像元的重叠而且重叠区的高程值不完全相等——不同航次、不同处理版本、不同滤波参数都会造成几十厘米到几米的差异。直接在重叠区做平均会在接缝两侧形成一条平滑但位置固定的“带”直接后来者覆盖前者接缝就变成了硬边坡度一算就是一条直线伪影。我的做法是分层优先级高精度源优先于低精度源同精度源里新版本优先于旧版本同版本里按分幅质量标记有效像元占比、空洞面积排序占比高的优先。实现上不靠手工而是在拼接脚本里维护一个按优先级排序的文件列表先铺低优先级的再用高优先级的覆盖——也就是“先到先得、后来覆盖”但排序完全由规则决定。import glob import os def priority_key(path): 优先级数值越小越先铺后被高优先级覆盖。 规则分辨率越细越优先(0)同分辨率下有效像元占比越高越优先(1)。 with rasterio.open(path) as ds: res abs(ds.res[0]) arr ds.read(1, maskedTrue) valid_ratio arr.count() / arr.size # 分辨率分档0.0004度(约12m) - 00.001度(约30m) - 1其余 - 2 tier 0 if res 0.0004 else (1 if res 0.001 else 2) return (tier, -valid_ratio) tiles sorted(glob.glob(clean/*.tif), keypriority_key) # 低优先级在前高优先级在后后续按顺序覆盖 tiles.reverse()这段的关键是reverse()。排序结果里 tier 小的是高精度reverse()之后高精度排到列表末尾在“后来覆盖”的逻辑里就拥有最高的覆盖权。有效像元占比取负值再排序是为了让占比高的排在前面——同样 reverse 之后它反而靠后同样拥有更高覆盖权。这个顺序一定要在拼接前打印出来人工扫一眼尤其是混了多个数据源的场景顺序错了整张图的接缝都会偏移。3.2 gdalbuildvrt 加 gdal_translate命令行最快的一条路如果全部分幅分辨率一致、坐标系一致最快的方式是先用gdalbuildvrt建一个虚拟镶嵌VRT再用gdal_translate落成实体文件。VRT 本身不复制数据只是一个 XML 索引几秒钟就能建完方便你先看看整体范围对不对再决定要不要落盘。# 1) 生成文件列表顺序就是覆盖顺序低优先级在前 ls clean/*.tif | sort tiles.txt # 2) 建虚拟镶嵌以最高分辨率对齐显式声明源 NoData 和输出 NoData gdalbuildvrt -input_file_list tiles.txt \ -resolution highest \ -srcnodata -32768 \ -vrtnodata -32768 \ -overwrite \ asia_dem.vrt # 3) 落盘分块、压缩、BigTIFF线程拉满 gdal_translate asia_dem.vrt asia_dem.tif \ -co TILEDYES -co BLOCKXSIZE512 -co BLOCKYSIZE512 \ -co COMPRESSDEFLATE -co ZLEVEL6 \ -co BIGTIFFYES \ --config GDAL_NUM_THREADS ALL_CPUS-resolution highest表示按输入里最细的像元大小对齐混拼时低精度源会被最近邻放大——这里要注意gdalbuildvrt不做重采样算法选择它默认最近邻对 DEM 来说最近邻会产生块状台阶。所以这个命令只适合“全部同分辨率”的场景一旦混拼必须换成下一节的 Python 分块方案或者先用gdalwarp把每个源统一到目标网格再拼。-srcnodata和-vrtnodata同时给是为了让 VRT 在读取时把源里的 -32768 当作无效值并输出为 -32768避免它在重叠区参与“有效值”判断。--config GDAL_NUM_THREADS ALL_CPUS在 DEFLATE 压缩阶段能明显提速但要注意它吃内存机器内存小于 16GB 时建议改成 4。还有一个现实问题亚洲全图按 30 米落盘实体文件轻松到几十 GBgdal_translate单进程跑几个小时很正常中途失败就得重来。所以实际项目里我很少直接落一张全图而是按纬度带切成几块每块单独落盘后面用的时候再按需读。这样即使某一块失败也只重跑那一块。3.3 rasterio 分块拼接把峰值内存压在 8GB 以内先说清楚为什么要分块。亚洲范围经度大约 25°E 到 180°E、纬度大约 -11°S 到 82°N跨度约 155°×93°。按 0.001°约 100 米网格算输出尺寸是 155000×93000约 144 亿个像元int16 存储就是 28.8 GB 的裸数组——这还没算中间计算产生的临时副本。任何“先把全图读进内存再处理”的写法都会直接 OOM。正确写法是预先分配输出网格然后逐块读源、按窗口写入目标。关键在于每一步都只碰当前块。import glob import numpy as np import rasterio from rasterio.transform import from_origin from rasterio.windows import from_bounds from rasterio.enums import Resampling # 目标网格覆盖亚洲范围0.001 度分辨率 RES 0.001 LEFT, BOTTOM, RIGHT, TOP 25.0, -11.0, 180.0, 82.0 WIDTH int(round((RIGHT - LEFT) / RES)) HEIGHT int(round((TOP - BOTTOM) / RES)) DST_TRANSFORM from_origin(LEFT, TOP, RES, RES) FILL np.int16(-32768) # 按全局网格分块写每块 4096 x 4096裸数组约 32MB BLK 4096 profile { driver: GTiff, dtype: int16, nodata: int(FILL), width: WIDTH, height: HEIGHT, crs: EPSG:4326, transform: DST_TRANSFORM, tiled: True, blockxsize: 512, blockysize: 512, compress: deflate, BIGTIFF: YES, } tiles sorted(glob.glob(clean/*.tif), keypriority_key) # 低优先级在前 with rasterio.open(asia_dem.tif, w, **profile) as dst: for row0 in range(0, HEIGHT, BLK): for col0 in range(0, WIDTH, BLK): bh min(BLK, HEIGHT - row0) bw min(BLK, WIDTH - col0) block np.full((bh, bw), FILL, dtypenp.int16) # 当前块对应的地理范围 bl LEFT col0 * RES bt TOP - row0 * RES br bl bw * RES bb bt - bh * RES for path in tiles: # 低优先级先铺 with rasterio.open(path) as src: l max(bl, src.bounds.left) r min(br, src.bounds.right) b max(bb, src.bounds.bottom) t min(bt, src.bounds.top) if l r or b t: continue # 无重叠跳过省掉一次解压 win from_bounds(l, b, r, t, src.transform) # 读进块内子范围out_shape 控制目标像元数双线性重采样 data src.read( 1, windowwin, boundlessTrue, fill_valueFILL, out_shape(int(round((t - b) / RES)), int(round((r - l) / RES))), resamplingResampling.bilinear, ) r0 int(round((bt - t) / RES)) c0 int(round((l - bl) / RES)) sub block[r0:r0 data.shape[0], c0:c0 data.shape[1]] # 先到先得只填还没被高优先级覆盖的位置 hole sub FILL sub[hole] data[hole] dst.write(block, 1, windowrasterio.windows.Window(col0, row0, bw, bh))逻辑说明分三点。第一外层双循环遍历的是输出网格的块不是输入文件这样峰值内存只跟BLK有关与总数据量无关把BLK设成 4096 时单块裸数组约 32 MB加上 rasterio 内部的解压缓冲全程峰值通常压在 1~2 GB。第二内层遍历输入分幅先用src.bounds做范围相交判断不相交直接continue——这一步能省掉大量无谓的文件打开和解压实测在千级分幅上能把总耗时砍掉一半以上。第三out_shape让 rasterio 在读取阶段就完成重采样配合Resampling.bilinear保证低精度源放大时不出现块状台阶boundlessTrue允许窗口越界越界部分自动填FILL省掉了大量边界特判。参数上最需要调的是BLK和RES。BLK太小比如 1024会导致频繁打开文件、系统调用开销盖过计算太大比如 16384单块 512 MB多进程时容易撞内存上限。RES要和你后续分析的尺度匹配做流域分析用 0.001° 够用做城市尺度地形必须换到 0.0003° 以下但那样上面那套数字全部要重算输出尺寸会变成天文数字必须改成按子区域分块落盘、最后再建 VRT 索引。3.4 拼完之后必须做的三项校验拼完就交付是新手最容易犯的错。我一般固定跑三项校验任何一项不过就不往下走。第一项是统计量对比。随机抽十个输入分幅计算各自的 min/max/mean再看拼接结果在对应范围内的统计值。差异超过 5% 就要查通常是 NoData 处理不一致或者重采样算法把极值削平了。第二项是接缝抽样。沿几条已知的分幅边界比如整度经纬线取一条剖面打印高程序列看有没有单像元的跳变。正常情况应该是连续过渡如果每隔固定间距出现 ±30 米以上的锯齿说明重叠区覆盖顺序有问题。第三项是无效值分布。统计FILL值在整张图里的占比和空间分布把它叠到范围边界上看。如果空洞大面积出现在研究区内部而不是海洋和境外说明某些分幅没拼上或者被错误地判成无效。# 快速统计最小/最大/均值/标准差以及 NoData 像元数 gdalinfo -stats asia_dem.tif | grep -E Minimum|Maximum|Mean|NoData-stats会触发一次全图扫描并写入.aux.xml几十 GB 的文件跑十几分钟到半小时不等建议放在后台。跑完把统计值和抽样结果记在一个文本文件里和处理参数放在一起——后面出问题回溯时这份记录比任何日志都有用。4. 投影与重采样亚洲DEM为什么不能一路用经纬度算坡度拼接完成得到的是一张经纬度网格的 DEM很多教程到这一步就直接上gdaldem slope了。在赤道附近的小范围还能看铺到亚洲尺度上必然翻车。原因不复杂经纬度网格的像元在东西方向对应的实际距离随纬度剧烈变化而gdaldem默认按像元尺寸当等距处理算出来的坡度在南北方向上系统偏大、东西方向上系统偏小越往高纬偏得越离谱。4.1 经纬度网格在亚洲的三处形变第一处是面积。同样 0.001°×0.001° 的像元在赤道附近约 111 m × 111 m面积约 1.23 万平方米到了北纬 60° 附近东西方向的实际边长缩到约 55 m面积只剩一半多点。用像元个数去估算流域面积、植被覆盖面积高纬地区会被系统性高估。第二处是坡度。坡度本质上是高程对水平距离的偏导。经纬度网格里行方向的像元间距南北向基本恒定列方向东西向随纬度收缩。gdaldem如果不做纬度校正等价于把所有像元当成等距正方形结果是东西向坡度被低估、南北向坡度相对被高估合成的坡向也随之偏转。第三处是距离。任何基于像元邻域的距离量算——比如汇流路径长度、山脊线长度、缓冲半径——在经纬度网格上都不是真实地面距离。范围越大误差累积越明显。解决方式只有一条在做坡度、汇流、距离量算之前先把 DEM 投影到等距或等积的米制坐标系。亚洲跨多个投影带逐带投影再拼接会引入新的接缝所以我一般选一个覆盖整块区域的全局投影牺牲一点局部形变换取一致性。4.2 gdalwarp 的四个关键参数与重采样算法选择全局投影的选择上我用得最多的是 EPSG:6933等积圆柱投影它的优点是面积保真适合做面积统计和汇流分析缺点是角度和距离在高纬仍有形变如果项目重点是坡度坡向的绝对精度就得换到局部横轴墨卡托分带做做完再转回来。这两种方案的取舍没有标准答案取决于你的分析目标是“面积准”还是“局部形状准”。gdalwarp -overwrite \ -t_srs EPSG:6933 \ -tr 250 250 \ -r bilinear \ -srcnodata -32768 -dstnodata -32768 \ -multi -wo NUM_THREADSALL_CPUS \ -co TILEDYES -co BLOCKXSIZE512 -co BLOCKYSIZE512 \ -co COMPRESSDEFLATE -co BIGTIFFYES \ asia_dem.vrt asia_dem_6933.tif-t_srs EPSG:6933是目标坐标系单位是米。-tr 250 250显式指定输出像元大小为 250 米——这里千万不能省省了 GDAL 会按源分辨率和投影变换自动推算一个奇怪的值而且不同版本推算结果不一致导致同一批数据两次处理结果对不上。250 米是个折中30 米源降到 250 米会损失细节但如果目标是全国或大区尺度的地形底图250 米在视觉和计算上都是舒服的档次要做小流域分析就用 30 米甚至更细代价是文件体量翻几十倍。-r bilinear是重采样算法。DEM 是连续表面双线性在大多数场景下是安全选择。但有两个例外一是数据里存在明显的台阶或人工地物比如梯田、矿坑边帮双线性会把硬边糊成斜坡这时候要用near二是做汇流分析前cubic或cubicspline能更好地保持地形曲率减少伪洼地。我的习惯是先用bilinear出一版看整体如果发现山脊变圆、谷底变浅再换cubic重跑对比。-srcnodata和-dstnodata必须成对给。只给前者输出文件里无效值可能变成 0只给后者源里的无效值会被当成有效高程参与插值输出边界出现一圈异常抬升。-multi配NUM_THREADS是并行化对重投影这种逐块独立的操作提速明显但内存占用会翻倍机器内存不足时改成-multi 2。4.3 WarpedVRT 分块重投影不爆内存的写法gdalwarp处理亚洲全图时同样是单进程流式写中途失败无法续跑的问题依旧存在。如果需要在 Python 流程里做重投影并且要和后续的分块处理串起来用WarpedVRT更灵活——它把重投影包装成一个“虚拟数据集”你按窗口读它它在后台自动完成坐标变换和重采样读多少算多少。import rasterio from rasterio.vrt import WarpedVRT from rasterio.enums import Resampling with rasterio.open(asia_dem.vrt) as src: # WarpedVRT 不落盘只描述投影后应该长什么样 with WarpedVRT( src, crsEPSG:6933, resolution(250, 250), resamplingResampling.bilinear, src_nodata-32768, nodata-32768, warp_mem_limit512, # 单次 warp 操作的内存上限单位 MB ) as vrt: profile vrt.profile.copy() profile.update( tiledTrue, blockxsize512, blockysize512, compressdeflate, BIGTIFFIF_SAFER, ) with rasterio.open(asia_dem_6933.tif, w, **profile) as dst: # 按输出块遍历vrt.read 内部只 warp 当前窗口 for _, win in dst.block_windows(1): dst.write(vrt.read(1, windowwin), 1, windowwin)这段的分块逻辑比上一节更省事dst.block_windows(1)直接按输出文件的内部瓦片遍历vrt.read(windowwin)只对当前窗口做坐标变换和重采样内存峰值约等于一个 512×512 块加上warp_mem_limit指定的 512 MB 缓冲。src_nodata和nodata同样必须成对给否则边界处会出现条带状伪影。一个实测经验warp_mem_limit设得太小比如 64会让 GDAL 频繁把中间结果刷到临时文件速度掉一半以上设得太大比如 4096在低内存机器上直接 OOM。512 到 1024 之间是大多数 16GB 机器的舒适区。另外resolution参数只在目标 CRS 是米制时才有直观意义如果换成别的 CRS 要留意单位。5. 避坑亚洲DEM文件处理里最容易翻车的五处下面五条都是我自己或者身边同行真实踩过的按“现象 → 原因 → 解决”记下来。每一条单独看都不复杂但它们有个共同特点错误不会当场报而是等到出图或者出分析结果时才暴露回溯成本极高。5.1 海面变成海拔 0 米NoData 没设引发的连锁反应现象晕渲图上海洋区域是一片均匀的灰白坡度图上沿海出现一圈异常高值带汇流分析里海水被算成了巨大的平地汇水区。原因DEM 文件里的海洋和空洞用某个特定值填充常见 -32768但文件头里的 NoData 字段是空的。读取时如果不加maskedTrue这个值会被当成真实高程如果填的是 0 而没声明问题更隐蔽——海平面正好是 0 米看起来完全合理直到你发现整个海域参与坡度计算时梯度全是 0海岸线附近出现巨大的伪梯度。解决转换阶段就用gdal_translate -a_nodata把无效值写进文件头读取阶段一律maskedTrue或者显式比较arr nodatagdaldem之前用gdalwarp把无效值统一。三道关卡任缺一道这条坑就会以不同形式复现。5.2 接缝处的“堤坝”混拼时没统一网格现象拼好的图上沿着整度经纬线有一条笔直的亮线或暗线宽度一到两个像元坡度图上这条线的值稳定在某个固定角度。原因直接对 90 米和 30 米的分幅做gdalbuildvrt -resolution highest低精度源被最近邻放大到 30 米网格边缘像元对齐不上重叠区取值时一边是原始值、一边是被放大的阶梯值接缝就成了一堵薄墙。解决混拼前必须先把低精度源用gdalwarp -tr显式重采样到目标网格并且用-r bilinear而不是默认的最近邻。更彻底的办法是分层处理——低精度源单独拼一张高精度源单独拼一张分析时按置信度掩膜叠加永远不让它们在同一个数组里混着算。5.3 坡度单位错得离谱经纬度直接上 gdaldem现象坡度结果的最大值动辄七八十度平均坡度比实际大好几倍坡向分布明显偏向南北。原因DEM 是 EPSG:4326像元单位是度。gdaldem slope默认按“像元尺寸即水平距离”计算把 0.001 度当成 0.001 米梯度被放大了十万倍量级加上东西向像元实际距离随纬度收缩南北向和东西向的比例失衡。解决gdaldem之前先用gdalwarp转到米制投影并设定明确的-tr。如果因为某些原因必须在经纬度网格上算至少加-s 111120把度换算成米这个近似仅适用于小范围低纬并且清楚这只是权宜之计高纬地区误差依然很大。5.4 OOM把全图读进内存现象脚本跑到某个分幅时进程被系统杀掉或者报MemoryError日志里没有任何异常堆栈。原因用src.read()不带 window 读整幅或者用numpy做全图运算。亚洲全图 int16 就是 28.8 GB 起加上 float32 转换、中间掩膜、重采样缓冲轻松突破 60 GB。解决任何涉及全图的操作都必须按窗口或按块遍历。用dst.block_windows()或者手写行块循环rasterio.merge只适合小范围numpy运算尽量用out参数写入预分配数组避免产生临时副本。判断标准很简单代码里出现ds.read()不带 window就要警觉。5.5 HGT 文件的坐标反了南北半球搞混现象拼出来的图上本该在北半球的区域跑到了南边或者整张图上下颠倒但单个分幅打开又是正常的。原因HGT 的文件名编码了分幅左下角的经纬度格式是N30E100这种。解析时如果按“先纬度后经度”的固定顺序切字符串遇到S开头或者经度三位数E100就会错位更常见的是把文件名里的纬度当成了分幅的中心而不是左下角导致整体偏移半度。解决解析用正则而不是字符串切片^([NS])(\d{2})([EW])(\d{3})$匹配不上直接报错跳过而不是猜构造 transform 时明确左下角原点from_origin(left, top, res, res)里的 top 是lat 11°×1° 分幅。拼完之后第一件事就是叠一张海岸线或者已知地标点做目视校验这一步花五分钟能省掉后面几小时的排查。6. 从DEM到地形底图晕渲、坡度与一条可交付的最小链路数据和中间成果都齐了最后落一个能交付的东西。地形底图最常见的两种产物是晕渲图和坡度分级图两者都能用gdaldem一条命令出来但参数细节决定成品是“能看”还是“专业”。6.1 晕渲的两个必调参数# 山体阴影光源方位角 315°、高度角 45°垂直夸张 1.5 倍 gdaldem hillshade asia_dem_6933.tif asia_hillshade.tif \ -az 315 -alt 45 -z 1.5 \ -compute_edges \ -co TILEDYES -co COMPRESSDEFLATE -co BIGTIFFYES-az 315是西北方向打光这是制图惯例因为人眼习惯把光源假设在左上反过来会看成凹陷。-alt 45是光源高度角低角度能强化微地形但会把缓坡压成全黑45° 是通用折中。-z 1.5是垂直夸张系数亚洲这种大尺度范围如果按 1:1 渲染山地几乎看不出起伏1.5 到 2.0 之间比较合适再高就会出现明显的“塑料感”。-compute_edges必须加否则瓦片边缘会出现一圈黑边。6.2 坡度分级与配色# 坡度输出单位度-alg Horn 对复杂地形更稳 gdaldem slope asia_dem_6933.tif asia_slope.tif \ -alg Horn -compute_edges \ -co TILEDYES -co COMPRESSDEFLATE -co BIGTIFFYES # 分级配色0-2 平地、2-5 缓坡、5-15 中等、15-30 陡坡、30 极陡 gdaldem color-relief asia_slope.tif slope_colors.txt asia_slope_rgb.tif \ -co TILEDYES -co COMPRESSDEFLATE -co BIGTIFFYESslope_colors.txt每行三到四个值坡度值、R、G、B最后加一行nv指定 NoData 颜色。-alg Horn比默认的ZevenbergenThorne在起伏剧烈的山区更稳定后者在平缓区表现更好选哪个看你研究区的地形类型——我一般两块都出叠着看差异。6.3 我自己的习惯一条完整的链路跑通之后我会把参数固化成脚本但每次都保留一份“参数快照”分辨率、重采样算法、投影、NoData、压缩方式、GDAL 版本号全部写进一个纯文本文件和数据放一起。血泪经验是同一份原始数据在不同参数组合下出来的坡度图差异能到 20% 以上半年后回头查一张图的来源没有这份快照就只能重跑而重跑又未必能复现。另外一个习惯是所有中间文件都不删只做冷归档——磁盘便宜重新下载和处理亚洲全图的时间成本才是真正贵的。处理亚洲DEM文件这件事难点从来不在某一个命令而在于链路上每一步都为下一步留了正确的接口。分辨率选对了、NoData 声明了、网格统一了、投影换过了、内存控住了剩下的都是常规操作。希望帮到你。本文还有配套的精品资源点击获取