企业官网建设流程全解析

简介：亚洲DEM高程数据集面向从事地形分析与空间建模的GIS从业者、测绘与地理专业师生，提供可直接加载的东亚地区数字高程模型，解决区域地形数据获取与底图准备的起步问题。压缩包内共8个文件，以7个adf与1个xml为主：adf为ESRI Grid二进制栅格的分块与索引文件，承载高程栅格本体及投影、统计信息，xml则记录元数据说明；整体约49.45MB，可在ArcGIS与MapGIS中直接打开。已有916人学习下载。借助这套数据，读者可完成坡度坡向计算、地形剖面绘制、洪水淹没模拟、日照与地质灾害评估等练习，也可将其与土地利用、气候数据叠加，观察东亚地形起伏与环境要素的关系；分块存放的子集结构便于按区域单独取用，适合作为教学实验与项目预研的基础底图。

1. 亚洲DEM文件：同一座山两次取高程差 30 米，问题出在哪

亚洲DEM文件这四个字，落到硬盘上是一堆N27E086.hgt这样的分幅包，落到屏幕上却常常变成麻烦：同一座山，两个工具读出来的高程差了三十米。我见过最典型的一次出问题，是有人把 SRTM 分幅直接拼起来算坡度，出来一片整齐的条纹，还以为是数据源精度不行——实际上那批瓦片在垂直基准、NoData 约定、像素对齐三项约定里，有两项没处理。

亚洲DEM文件真正的难点在于：来自不同源、不同基准、不同分幅规则的栅格，要统一成一套能算坡度、算汇流、算淹没的数值场；下载本身只是最简单的环节。它适合三类人：做水文与洪涝模拟的、做路桥选线与场地平整的、做大范围地形可视化与切片服务的。下面按“选源—拼接—验收—排错”的顺序推进，参数和判据都写清楚，能照着跑，也能看懂为什么这么设。

2. 亚洲DEM文件从哪来：SRTM、ALOS 与 Copernicus 的精度对照与选型

选源这一步决定后面 80% 的麻烦。亚洲这块地方的特殊性太强：青藏高原与喜马拉雅南坡是全球地形起伏最剧烈的地带，雷达侧视成像在这里叠掩和阴影成片出现，空洞天然集中；东南亚热带雨林的地表被十几米到二十几米的树冠盖着，雷达和光学拿到的都是冠层顶，跟真实地面差一大截；中亚荒漠裸露地表反而是精度最好的区域。所以“哪个数据源更准”这个问题，在亚洲没有统一答案，得按你研究区的地表类型分。

2.1 分辨率相同不等于能互换：四种常见源的关键参数

网格间距都写着 30 米的产品，互相之间可能差出十几米。下面这张表是我在选型时会先摆出来的对照，标称值来自各产品公开的技术说明，实际精度还得靠自己抽检查点验证，别拿标称值当验收结论。

选型的判断顺序我一般这么走：先看研究区是否落在 SRTM 覆盖范围外（亚洲北部高纬会出界），出界就只能选 ALOS 或 Copernicus；再看是否要做水文，要做水文就必须接受“所有这些都是 DSM”这个事实，植被区要么换专门的 DTM，要么先做降采样把冠层噪声平均掉；最后看垂直基准，如果下游要和 GNSS 椭球高、或者和已有 EGM2008 的工程数据混算，优先选 Copernicus，省一次基准转换。

提示：同一个项目里混用两种源是最容易出问题的地方。两块拼接区域的高程差会直接变成一条直线状的地形陡坎，在晕渲图上非常显眼，但在数值统计里往往被平均掉，看不出来。

2.2 分幅命名与范围推导：从 N27E086 到一份不重不漏的瓦片清单

1°×1° 分幅的命名规则是N27E086：N 后面两位是分幅西南角的纬度整数，E 后面三位是西南角的经度整数。它覆盖的是北纬 27° 到 28°、东经 86° 到 87° 这一格，正好落在喜马拉雅中段。HGT 格式本身还有两条约定必须记住：数值是 16 位有符号大端整数，所以读之前要转字节序；行列顺序是从北往南，第一行对应分幅的北边界，这一点跟绝大多数影像库“第一行在南”的直觉相反，直接imshow出来是上下颠倒的。

从研究区外接矩形推瓦片清单，看着简单，但负经纬度会出问题。亚洲范围跨到南半球（印尼南部到 -11°），也跨到西经（亚洲最西端在高加索一带其实还是东经，但如果研究区含北非边缘就会碰到西经）。用int()截断对负数是向零取整，int(-11.2)得 -11，可分幅应该落到 -12 那一格，结果就漏了一整行。必须用math.floor。

import math def tile_names(west, south, east, north): """按 1°x1° 分幅规则生成覆盖研究区的瓦片名，不重不漏。 西/南边界向下取整，东/北边界向上取整，range 右开区间自动处理闭合。""" lat0, lat1 = math.floor(south), math.ceil(north) lon0, lon1 = math.floor(west), math.ceil(east) names = [] for lat in range(lat0, lat1): for lon in range(lon0, lon1): ns = "N" if lat >= 0 else "S" ew = "E" if lon >= 0 else "W" names.append(f"{ns}{abs(lat):02d}{ew}{abs(lon):03d}") return names # 中亚到东亚的一块：东经 73-135，北纬 18-54 tiles = tile_names(73, 18, 135, 54) print(len(tiles), tiles[0], tiles[-1]) # 2232 N18E073 N53E134

这段的关键在三处。math.floor和math.ceil保证边界格一定被包含，宁可多取一格也不能漏；abs()加格式化宽度02d／03d保证S、W前缀下的位数和北半球一致，很多数据源的文件名就是这么排的；range(lat0, lat1)右开，当north恰好是整数 54 时ceil得 54，最后一格是 53，正好覆盖 53°–54°，不会多取一格空瓦片。

清单出来以后，先跟手上的文件列表做一次差集，缺哪几格一目了然。亚洲陆域整块按 1° 分幅是四位数到五位数级别的瓦片量，靠人工核对必然出错，这一步一定要脚本化。

2.3 用 GDAL 建虚拟栅格：先看覆盖率，再谈处理

拿到一堆分幅后不要急着gdal_merge.py。那个脚本会把所有数据读进内存再写，亚洲级别的数据量会直接把机器拖垮。正确做法是先建 VRT——它本质上只是一份 XML 索引，记录每个瓦片的位置和读取偏移，不复制任何像素。

# 1) 先确认单个瓦片的真实情况：尺寸、类型、NoData、范围 gdalinfo N27E086.hgt | head -30 # 2) 把所有瓦片路径写进清单，一行一个 ls -1 *.hgt > tiles.txt # 3) 建虚拟栅格。-resolution highest 表示按最精细的那层走 gdalbuildvrt -resolution highest \ -srcnodata -32768 -vrtnodata -32768 \ -input_file_list tiles.txt \ asia_raw.vrt # 4) 看整体覆盖：范围对不对、有没有缺块、统计量是否合理 gdalinfo -stats asia_raw.vrt | tail -20

gdalinfo那一步重点看四行：Size is 3601, 3601（1 弧秒分幅的标准尺寸，3 弧秒是 1201）、Type=Int16、NoData Value=-32768、Origin的经纬度。如果某个瓦片的尺寸不是 3601，说明它被重采样过或者根本不是 1 弧秒产品，混进去会让gdalbuildvrt的接边错半个像素。

-srcnodata和-vrtnodata同时给，是为了在 VRT 层就把空洞标记出来。只给-srcnodata的话，VRT 输出端的 NoData 是未定义的，后面gdalwarp做重采样时会把 -32768 当成真实高程参与加权平均，产生一堆负一万多米的“深坑”，这就是第 1 章里那个 30 米误差的另一个变种，只不过量级更夸张。

-stats会让 GDAL 扫一遍全区算最小最大值和均值。亚洲整块扫下来要几分钟到几十分钟，取决于磁盘。看到STATISTICS_MINIMUM是 -32768 就说明 NoData 没被正确排除；看到最小值是 -500 左右，说明有部分空洞被当成了数据；看到最小值在 -430 附近（死海一带），那才是正常的。

3. 拼成一块亚洲 DEM：NoData、重采样与垂直基准的参数怎么定

拼接这一步的每一个参数都会在成果里留下痕迹，而且大多是不可逆的。我一般把它拆成四件事按顺序做：统一 NoData、定重采样算法、定目标投影、统一垂直基准。顺序不能乱，尤其是垂直基准要放在最后，因为它只改数值不改几何，放前面会让后面所有重采样的中间结果都带着一个偏移量，排查起来很麻烦。

3.1 -32768 不是高程：NoData 的三层约定与统计自检

NoData 在 DEM 处理里有三层含义，很多人只处理了第一层。第一层是文件级的：HGT 用 -32768，GeoTIFF 常见 -9999 或 -32767，也有产品用 NaN。第二层是语义级的：0 到底代表海平面还是代表“无数据”？Copernicus DEM 的海面就是 0，而且这是有效值；但有些产品把水体也压成 0 当空洞处理。第三层是运算级的：任何一次重采样、任何一次滤波，都必须让 NoData 参与掩膜而不是参与计算。

第二层最容易出错。我见过有人写了一句arr[arr == 0] = np.nan，本意是清理空洞，结果把里海、咸海、贝加尔湖连同所有沿海低地全挖成了坑，后面算汇流时整个流域的水系断成几十段。判断办法是先看产品的技术说明里 0 的定义，再用海岸线矢量做一次叠加抽样：抽一百个海岸线上的点，看高程是不是稳定在 0 附近。

import numpy as np from osgeo import gdal ds = gdal.Open("asia_raw.vrt") band = ds.GetRasterBand(1) print("NoData =", band.GetNoDataValue()) # None 就说明没设，必须补 # 只读一小块做统计，别整幅 ReadAsArray arr = band.ReadAsArray(0, 0, 4000, 4000).astype("float32") print("分位数", np.percentile(arr, [0, 0.1, 50, 99.9, 100])) nd = band.GetNoDataValue() if nd is not None: valid = arr[arr != nd] print("有效像元占比 %.2f%%" % (100 * valid.size / arr.size)) print("低于 -500 m 的有效像元数", int((valid < -500).sum()))

分位数那一行是核心判据。正常的亚洲 DEM，0.1% 分位大概在 -400 上下（死海、吐鲁番盆地一带），99.9% 分位在 6000 出头（喜马拉雅、喀喇昆仑）。如果 0.1% 分位跑到 -30000，就是 NoData 没排除；如果 99.9% 分位超过 8000，多半是重采样时cubic在珠峰附近产生了过冲。

ReadAsArray的四个参数是xoff, yoff, xsize, ysize，注意它按栅格坐标而不是经纬度，所以想读某个经纬度范围得先换算。整幅读入是内存爆掉的头号原因，亚洲 30 米整块浮点数组大概在百 GB 量级，别试。

3.2 重采样算法与目标投影：坡度和水系最敏感的两个参数

重采样算法的选择有一条硬规则：连续高程场绝不能用最近邻。near会在等值面上产生台阶，算出来的坡度是块状噪声，而且这种错误在数值统计里几乎看不出来，只有画晕渲图才会暴露。剩下几个算法按用途分：

• average：降采样首选。30 米降到 250 米，取 8×8 窗口均值，能顺带把树冠噪声和雷达噪声压下去，代价是山脊被削平一点。
• bilinear：升降采样都稳，不会过冲。做坡度、做汇流、做淹没分析，我都用它。
• cubic/cubicspline：视觉上更锐，但会在陡崖和雪线附近产生过冲，出现负高程或超出物理上限的峰值。只用于出图，不用于分析。
• mode/max/min：特定用途，比如提取建筑高度包络。别用在通用流程里。

目标投影的选择在亚洲比在欧洲重要得多。直接在 EPSG:4326 上算坡度是错的，因为经纬度网格的东西向像素宽度随纬度按 cos(lat) 收缩，北纬 45° 处一个 1 弧秒的像素在东西向只有南北向的 70% 宽，gdaldem却按正方形像素算，坡度会整体偏陡。跨越 70°E 到 135°E、纬度从 -11° 到 55° 的范围，这种偏差在不同纬度还不一样，没法用一个系数修正。

常见做法是投到等积或等距的米制投影。全球范围我一般用 EPSG:6933（等积圆柱），亚洲区域用 Albers 等积圆锥更贴合，参数大致是中央经线 105°E、标准纬线 25°N 和 47°N，这个组合在中亚到东亚的形变比较均匀。

gdalwarp -overwrite \ -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +datum=WGS84 +units=m +no_defs" \ -tr 250 250 -tap \ -r bilinear \ -srcnodata -32768 -dstnodata -9999 \ -multi -wo NUM_THREADS=8 \ -co TILED=YES -co COMPRESS=DEFLATE -co PREDICTOR=2 -co BIGTIFF=IF_SAFER \ asia_raw.vrt asia_250m.tif

逐个说参数。-tr 250 250是目标分辨率，单位是目标投影的米，所以这里是 250 米格网；如果想保持像元数不变，用-ts width height指定行列数。-tap强制输出格网对齐到分辨率的整数倍，这是分块处理的前提——不加上它，两次gdalwarp的成果在拼接时会错半个像素，接缝处出现一条细缝。

-multi加-wo NUM_THREADS=8是并行读写，对磁盘 IO 密集的大范围重投影提速明显，线程数别超过物理核数。-co PREDICTOR=2只对整型和浮点有效，配合DEFLATE能把高程数据的压缩率提高三到五成。BIGTIFF=IF_SAFER是防止输出超过 4 GB 时静默失败，超过 4 GB 的 GeoTIFF 在老版本 GDAL 下会直接写坏，这个后悔药必须提前吃。

-dstnodata -9999是刻意换了个值。用 -32768 的好处是和源一致，坏处是它在 Int16 里是边界值，某些软件做符号扩展时会出问题；用 -9999 更安全，代价是必须用 Int16 或 Float32 存，不能用无符号类型。

3.3 垂直基准统一：EGM96 与 EGM2008 差的那一米怎么补

高程数值必须说清楚“相对谁”。SRTM 和 ALOS World 3D 给的是相对 EGM96 大地水准面的正高，Copernicus DEM 给的是相对 EGM2008 的正高，GNSS 测出来的是相对 WGS84 椭球的大地高。三者之间的差值不是常数：椭球高减正高等于该点的大地水准面差距，这个差距在亚洲从 -30 米到 +50 米都有；而 EGM96 与 EGM2008 两个模型之间的差，在亚洲大部分地区是 0.5 到 2 米。

一两米对水文模拟确实无关紧要，但对土方计算、对淹没水深只有 0.5 米的湿地生态分析、对和实测水准点做比对，就是结论级的误差。统一基准在 GDAL 里可以靠复合坐标系一次做完：

# EGM96 正高 -> EGM2008 正高，中间自动经过 WGS84 椭球高中转 gdalwarp -overwrite \ -s_srs "EPSG:4326+5773" \ -t_srs "EPSG:4326+3855" \ -tr 0.0002777778 0.0002777778 -tap \ -r bilinear \ -srcnodata -32768 -dstnodata -32768 \ srtm_egm96.tif srtm_egm08.tif

EPSG:4326+5773是“WGS84 水平 + EGM96 高程”的复合坐标系，EPSG:4326+3855是“WGS84 水平 + EGM2008 高程”。-tr这里写的是度，0.0002777778 就是 1 弧秒，因为输出还是地理坐标。整个变换分两步：先把 EGM96 正高加上该点的 EGM96 差距变成椭球高，再减去 EGM2008 差距变成新的正高。

这里有个问题值得单独说：这一步依赖 PROJ 的大地水准面格网文件。本地没装格网时，gdalwarp要么直接报错说找不到 grid，要么在某些版本下静默退化成“只做水平变换、不做垂直改正”——你会得到一份看起来完全正常、实际上偏移了一两米的成果。判断办法很简单：找一个已知两种基准高程差的点，转换前后各读一次，看差值是不是落在预期区间。别信“命令没报错就是对的”。

4. 亚洲DEM文件能不能用：检查点 RMSE、坡度接缝与水系自洽

成品出来以后，最危险的心态是“图看着挺像那么回事”。晕渲图对高程误差极不敏感，偏移 20 米在视觉上几乎看不出来，但足以让一条分水岭挪位置。验收必须用数值判据，我一般走三条线：检查点误差、坡度接缝、水系自洽。

4.1 检查点抽多少、怎么分层、RMSE 和偏差要分开看

检查点的来源有三类：实测水准点或 GNSS 静态观测点，精度最高但数量少；更高精度的 DEM 产品做交叉验证，数量管够但存在相关性，只能做相对判据；大比例尺地形图上的等高线交点，亚洲大部分地区都能找到，精度中等但空间分布均匀。

抽样策略比样本量更重要。整块亚洲只报一个 RMSE 是没有意义的，因为误差跟地表类型强相关。我一般按地表分层：平原农业区、丘陵、高山峡谷、热带雨林、荒漠戈壁，每层至少 30 个点，层内按经纬度均匀撒。总量控制在几百个就够，多了边际收益很低。

import numpy as np def rmse_bias(pred, truth): """pred 是待验收 DEM 的高程，truth 是检查点真值，单位一致""" pred = np.asarray(pred, dtype="float64") truth = np.asarray(truth, dtype="float64") m = np.isfinite(pred) & np.isfinite(truth) d = pred[m] - truth[m] # 正 = DEM 偏高 return { "n": int(m.sum()), "bias": float(d.mean()), # 系统性偏移，先看这个 "rmse": float(np.sqrt((d ** 2).mean())), "le90": float(np.percentile(np.abs(d), 90)), "p95": float(np.percentile(np.abs(d), 95)), } print("高山峡谷层", rmse_bias(pred_hi, truth_hi)) print("平原层", rmse_bias(pred_pl, truth_pl))

判读顺序是先bias后rmse。bias明显不为零（比如整体偏高 15 米）说明有系统性因素，通常是植被冠层、垂直基准没统一、或者重采样算法带偏，这种误差可以整体减掉，但要先搞清楚来源再减。bias接近零而rmse很大，说明是随机误差或者局部错误，得回去查空洞填补和接边。

le90是 90% 分位的绝对误差，很多产品规范用这个口径，方便和标称值对比。p95是我自己加的，因为亚洲的高山峡谷层误差分布尾部很重，le90看着漂亮，p95可能翻一倍。

4.2 坡度接缝与水系自洽：两个能一眼看穿拼接错误的检查

接边错位在原始高程上几乎看不出来，但在坡度上会被放大成一条亮线。原理很简单：坡度是高程的一阶差分，半个像素的错位造成的高程跳变除以像素间距，就是一条人为的陡坎。

# 已经是米制投影，-s 用默认 1.0；Horn 算法对噪声比 ZevenbergenThorne 稳健 gdaldem slope -alg Horn -compute_edges -s 1.0 asia_250m.tif asia_slope.tif # 沿一条经线抽剖面，找孤立的高梯度位置 python3 - <<'PY' import numpy as np from osgeo import gdal ds = gdal.Open("asia_250m.tif") b = ds.GetRasterBand(1) w = ds.RasterXSize col = w // 2 prof = b.ReadAsArray(col, 0, 1, ds.RasterYSize).astype("float64").ravel() d = np.abs(np.diff(prof)) idx = np.argsort(d)[-10:][::-1] print("最大跳变(米):", d[idx].round(1)) print("位置(行):", idx) PY

-alg Horn是 3×3 窗口的经典算法，比ZevenbergenThorne更抗噪，代价是山脊被略微平滑。-compute_edges让边缘一圈也能算出坡度，不加的话输出会比输入小一圈，后续做掩膜时对不上。-s 1.0是水平单位与垂直单位的比值，在米制投影下就是 1；如果误在经纬度数据上用，-s得填 111120 之类的近似值，但那种做法在不同纬度都是错的，别用。

剖面检查的判读标准：自然地形的相邻像元高程差在 250 米格网上很少超过 100 米，除了极少数垂直崖壁。如果最大跳变出现在某个整度数的行号上（比如正好是 1° 对应的行数整数倍），那基本可以确定是接边错位。

水系自洽更直接：填洼、算流向、算汇流累积，然后看河网是不是连续的。接缝错位会在汇流累积图上表现为一条笔直的、横跨整幅图的异常高值带——因为水在接缝处被挡住了，全积在那里。这个检查不需要真值，只要有常识就能判断，是我最喜欢用的一招。

4.3 晕渲图与空洞掩膜：三十秒的目视初筛

数值检查之前先做目视初筛，能省很多时间。晕渲图的参数里，-az方位角建议固定在 315°（西北方向打光），这是 cartographic 惯例，反过来打光会让人产生“山变成谷”的错觉，看久了会误判。-z垂直夸张在亚洲这种起伏剧烈的地方用 1.5 到 2 就够，再大会到处是假陡崖。

gdaldem hillshade -z 2 -az 315 -alt 45 -compute_edges asia_250m.tif asia_hs.tif # 把 NoData 单独导成二值图，空洞分布一目了然 gdal_calc.py -A asia_250m.tif --calc="(A==-9999)*255" \ --outfile=asia_mask.tif --type=Byte --NoDataValue=0

-alt 45是光源高度角，45° 是视觉上最平衡的选择，太高失去立体感，太低阴影拉得过长。空洞掩膜那张图重点看两处：高山区的空洞是不是成条带状（雷达阴影的典型形态），以及水体边缘是不是出现了不该有的空洞（0 值被误判的迹象）。这两类问题的处理方式完全不同，前者要插值填补，后者要改 NoData 判据。

5. 亚洲DEM文件避坑排查：五个反复出现的问题

下面五条是我这些年反复遇到的，按现象、原因、解决三段写。每一条都对应过真实的返工，最短的一次耽误两天，最长的一次整套水文模型重跑。

5.1 数值与 NoData 类

现象一：成果里出现大片 -16000 米左右的“深坑”，集中在山区。

原因：gdalbuildvrt或gdalwarp只设了-srcnodata没设输出端 NoData，重采样时 -32768 被当成真实高程参与加权平均。两个空洞像元和两个有效像元做双线性，结果就是 -16000 这个量级。

解决：输入端和输出端都要显式声明，-srcnodata -32768 -dstnodata -9999。已经错了的成果没法靠后处理救回来，只能重跑——这是典型的没有后悔药的操作，参数要先想清楚再按回车。

现象二：里海、咸海变成了大坑，沿海低地被挖掉一块。

原因：把 0 一刀切当成了 NoData。很多产品里 0 是有效的海平面高程，Copernicus DEM 尤其如此。

解决：先去产品说明里确认 0 的语义，再决定是否掩膜。真要区分水体和空洞，得引入独立的水体掩膜数据做叠加，不能靠高程值猜。已经误掩的情况，用原始数据重新生成一遍比修补快。

5.2 投影与基准类

现象三：坡度整体偏陡，且高纬区域偏得更厉害，能到 30%。

原因：在 EPSG:4326 地理坐标系上直接跑gdaldem slope。经纬度网格的东西向像元宽度是南北向的 cos(lat) 倍，gdaldem按正方形处理，等于把东西向的坡度系统性放大了。

解决：先投到米制等积投影再算坡度，Albers 或 EPSG:6933 都行。注意重投影要在计算坡度之前，不能先算坡度再投影，那样投出来的是错的坡度场。

现象四：和平原区的水准点比对整体差一两米，山区差得更多，但水平位置看着是对的。

原因：垂直基准混用。SRTM 是 EGM96，Copernicus 是 EGM2008，两者在亚洲差 0.5 到 2 米；或者更糟，gdalwarp做基准转换时本地缺大地水准面格网，静默退化成只做水平变换。

解决：整条链路统一到 EGM2008，转换后用已知差异的点验证一次。离线环境要手动把格网文件放到PROJ_DATA指向的目录，别指望它自己下载。这一步我现在的习惯是：转换前后各抽十个点打印高程，差值不为零才继续往下走。

5.3 性能与工程化类

现象五：处理到一半进程被系统杀掉，日志停在一个随机位置，没有任何报错。

原因：整幅ReadAsArray()或者用了gdal_merge.py。亚洲 30 米整块的数据量在这个量级上，任何全量读入都会吃光内存，被 OOM killer 干掉时是没有任何堆栈的。

解决：全部改成窗口读，或者用 VRT +gdalwarp -multi让 GDAL 自己按块调度；输出加-co TILED=YES和BIGTIFF=IF_SAFER。分块大小按 512×512 或 1024×1024 起，太小会让 IO 次数暴涨，太大又失去分块意义。

6. 亚洲DEM文件的进阶用法：COG 加分块，把整块亚洲装进一台笔记本

做完上面几步，你手上大概会有一份几百 GB 的 GeoTIFF 或 VRT。这时候真正的工程问题才浮出来：下游应用怎么读。整块加载不现实，切片服务又需要快速随机访问，解决方案是把成果转成云优化 GeoTIFF（COG），再用窗口读做分块计算。

COG 的核心是把文件按 512×512 或 256×256 的瓦片重排，同时把各层级金字塔内嵌进同一个文件，读取任意一小块时只需要几次寻道，不用扫全文件。

gdal_translate -of COG \ -co COMPRESS=DEFLATE -co PREDICTOR=2 -co LEVEL=9 \ -co BLOCKSIZE=512 -co OVERVIEWS=AUTO -co OVERVIEW_RESAMPLING=AVERAGE \ asia_250m.tif asia_250m_cog.tif

参数里有三个值得说明。OVERVIEWS=AUTO让 GDAL 自己决定建几层金字塔，规则是每层降采样 2 倍直到最小边小于 256，够用；OVERVIEW_RESAMPLING=AVERAGE必须显式指定，默认是最近邻，金字塔层会出现块状噪声，缩略图看着还行，一旦拿来做快速预览的坡度分析就露馅了；PREDICTOR=2配DEFLATE对高程这种空间连续的数据压缩效果最好，LEVEL=9是最高压缩比，代价是转换时间翻倍，但只做一次，值得。

实际用的时候，分块读长这样：

import rasterio from rasterio.windows import Window import numpy as np def block_iter(path, block=2048): """按块遍历，返回 (window, array)，内存占用只跟 block 有关""" with rasterio.open(path) as src: for y in range(0, src.height, block): for x in range(0, src.width, block): w = min(block, src.width - x) h = min(block, src.height - y) win = Window(x, y, w, h) arr = src.read(1, window=win, masked=True) yield win, arr # 逐块算坡度均值，统计里顺手排除 NoData for win, arr in block_iter("asia_250m_cog.tif"): if arr.mask.all(): continue print(win.col_off, win.row_off, float(np.ma.mean(arr)))

masked=True是关键，它让 rasterio 按文件的 NoData 声明自动生成掩膜，后续所有np.ma运算都会跳过无效像元，不用自己写判断。block=2048在 250 米格网上大约是 512 公里见方，单块内存约 16 MB，笔记本随便跑。这个循环最大的好处是内存占用恒定，处理亚洲整块和一景影像的代码完全一样。

顺手提一个分块处理的延伸做法：如果下游是训练模型或者做统计分析，可以在 COG 之上再建一层“瓦片索引”——把每个块的行列号、经纬度范围、有效像元占比、高程均值方差写进一张 Parquet 表。后面要按高程分层抽样、要跳过空洞密集的块、要做并行任务分发，都是查这张表的事，不用再碰原始栅格。这一步做过一次之后，再处理别的大区域数据基本就是复用。

最后说个我自己的习惯。每次交付前，我会把成果丢回最开始那批检查点上跑一遍rmse_bias，并且只在这时候才看bias的符号——如果它和我动手之前实测的符号相反，说明中间某一步引入了新的系统偏移，不管rmse多漂亮都得回头查。这个习惯救过我两次，一次是垂直基准退化，一次是重采样算法被同事改成了cubic（理由是“看起来更清楚”）。放在这条链路上的经验就一句：地形数据的所有参数都在数值里留指纹，图好看不算数，数字对得上才算数。希望帮到你。

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