简介面向全球陆地与海底地形联合分析需求这份ETOPO1地形数据以GeoTIFF栅格形式提供约2km分辨率的全球DEM坐标系统为WGS1984可用于GIS制图、地形起伏分析、海陆一体化建模、三维可视化底图及地学教学实验适合具备基础GIS操作能力的研究者与学习者使用。压缩包为rar格式共3个文件、约266.01MB主体为tif栅格数据另附ovr金字塔概览文件与xml辅助元数据便于在常见GIS软件中快速浏览、缩放和读取坐标信息。已有3036人学习或下载说明其在全球尺度地形应用中具有较高参考价值。读者可据此直接开展全球高程与海底地形提取、坡度坡向计算、区域裁剪和专题图制作也可作为气候、生态、水文及海洋相关研究的统一空间底图减少多源地形数据拼接与坐标转换的准备工作。1. 全球地形数据tif格式为什么第一件事不是打开看图把一张 5°×5° 的全球地形数据 tif 直接拖进普通看图软件多半是转圈、全黑或者一片灰白分不出山脉和海沟。不是文件坏了是这类 tif 原本就不是为肉眼浏览准备的——它是带地理坐标的栅格矩阵像素值是高程米数覆盖范围可以很大文件也可能达到几十 GB。全球地形数据 tif 格式要解决的核心问题只有一个把起伏的地表切成规则网格每个格子记一个高程值再用 GeoTIFF 标签把这些格子钉回经纬度坐标上。它支撑地形可视化、水文汇流分析、选址与通视评估这类工作新手容易在“能打开”和“能用”之间卡很久熟手则更在意分幅规则、无效值、垂直基准这些边界参数。本章先讲清这类文件的结构与参数再落到拼接、重采样、派生产品和上线压缩的具体命令。2. 全球地形 tif 里到底装了什么拿到一批瓦片先别急着拼。花十分钟把单幅的元数据读出来比后面折腾几小时返工划算得多。全球 DEM 的分发方式是“一套产品切成上千个矩形块”块与块之间靠文件名的经纬度对齐任何一块的坐标参考、无效值或数据类型和其他块不一致整张拼图就会在接缝处出现台阶、黑边或者数值断层。2.1 GeoTIFF 的六个关键字段普通 tif 只存像素GeoTIFF 在文件头里多塞了一组地理标签。真正影响后续处理的是这六个宽高与波段数、数据类型、坐标参考系CRS、仿射变换六参数、无效值NoData、压缩与分块方式。仿射变换是最容易被忽略的一个它定义了“左上角像素的左上角在地球上的哪个点”以及“每个像素跨多少度”六个参数依次是左上角 x、x 方向像素宽、旋转项、左上角 y、旋转项、y 方向像素高通常是负值因为图像行号向下增长。import rasterio # 读一个瓦片的元数据不做任何计算 with rasterio.open(tile_N00E000.tif) as ds: print(尺寸:, ds.width, x, ds.height, 波段:, ds.count) print(数据类型:, ds.dtypes) # 常见 int16 / float32 print(坐标参考:, ds.crs) # 全球 DEM 多为 EPSG:4326 print(仿射变换:, ds.transform) # 六参数注意 e 为正、f 为负 print(无效值:, ds.nodata) # 没设就是 None后面必然出事 print(分块:, ds.block_shapes[0]) # (256, 256) 之类 print(压缩:, ds.profile.get(compress))这段代码跑完重点看两个地方。第一nodata如果打印出None说明这份数据没有声明无效值海面或空洞区可能是 0、-32768 或者极大负数统计均值、算坡度都会把无效区域算进去。第二transform里的a和e像素宽高如果不是同一个数量级或者a是 0.0083333 这种数字说明它是 30 弧秒级产品不是常见的 1 弧秒单位是“度”不是“米”这一点决定后面所有距离和坡度计算要不要换算。ds.profile这个字典值得留一份写出新文件时直接复制再改几个键比手写一堆参数稳。常见写法是profile src.profile.copy()然后改driver、dtype、compress、tiled、blockxsize、blockysize、nodata。改完再传给rasterio.open(..., w, **profile)元数据基本不会丢。2.2 分辨率不是固定米数1 弧秒、3 弧秒与瓦片命名全球地形产品常说的“30 米”“90 米”是近似说法真实单位是角秒。1 弧秒的纬度跨度在地表约 30.87 米子午线周长约 40007 公里除以 360 度再除以 3600 秒这个值全球基本恒定但 1 弧秒的经度跨度要乘以纬度的余弦赤道处约 30.87 米到北纬 60 度只剩约 15.4 米。所以同一份 1 弧秒数据在高纬地区像素在地面上是“窄长条”不是正方形。这个差异带来的直接后果是任何以“米”为单位的计算都不能直接拿像素个数乘 30。算面积要先投影到等积坐标系算坡度要让 GDAL 知道纬度换算比例算缓冲区要按实际地面距离而不是行列数。很多人第一次做全球尺度分析发现高纬区域的坡度整体偏大、面积整体偏小根子就在这里。产品级别角秒赤道地面分辨率全球瓦片数量级单幅 5°×5° 体积量级粗分辨率全球网格30 弧秒约 900 米数百几 MB中分辨率3 弧秒约 90 米数千数十 MB高分辨率1 弧秒约 30 米上万上百 MB瓦片命名一般是“半球 纬度 半球 经度”的形式比如N00E000、S45W120数字是瓦片左下角或左上角看产品约定的整度坐标。拼接前务必确认三件事命名用的是左下角还是左上角、瓦片边长是 1 度还是 5 度、经度是 0~360 还是 -180~180。这三件事任一搞错拼出来就是错位或者中间一条空白带。2.3 垂直基准与数据类型高程到底相对哪个面高程值本身不带说明它一定相对某个参考面。公开全球 DEM 大多给的是正高相对大地水准面常见 EGM96 模型而 GPS 和很多测绘软件默认输出的是椭球高相对 WGS84 椭球。两者在全球范围内的差值大约在负一百米到正八十多米之间波动具体到某个区域可能是几十米。做洪水淹没、管线纵断面、机场净空这类对绝对高程敏感的分析时混用两套基准就是事故级错误做可视化山体晕渲则影响不大。import numpy as np import rasterio with rasterio.open(tile_N00E000.tif) as ds: arr ds.read(1, maskedTrue) # maskedTrue 自动把 nodata 屏蔽 print(最小值:, float(arr.min()), 最大值:, float(arr.max())) print(无效值占比:, 1 - arr.count() / arr.size) # 快速判断单位全球陆地最高约 8848 米 # 若最大值接近 29000基本可以确定单位是英尺需要乘 0.3048数据类型同样要看清。int16的取值范围是 -32768 到 32767装全球陆地高程绰绰有余海面常用负值表示水深或者干脆用 -32768 当无效值。float32支持小数代价是体积翻倍。如果一份int16数据里出现了 -32768 且nodata已声明它是无效值如果nodata没声明它就会被当成“负三万米的海沟”参与运算晕渲图上会出现一片突兀的黑块。提示拿到任何一份全球地形 tif先跑一遍min / max / count / nodata四个值再决定后面怎么处理。这一步不做后面所有分析结果都不可信。3. 用 GDAL 与 rasterio 把全球瓦片拼成一张可用底图拼接全球地形有两条路物理合并成一个大文件或者用虚拟栅格只建索引。前者适合做最终交付后者适合先看清全局、验证瓦片没错位。我一般的顺序是先用虚拟栅格扫一遍找问题确认无误再决定落盘的范围和分辨率避免上来就 merge 几十 GB 然后中途内存爆掉。3.1 环境与最小读取Python 侧直接用pip install rasterio它自带打包好的 GDAL 二进制省去编译麻烦需要命令行工具就装gdal的 conda 包。两边版本尽量对齐否则rasterio调用的 GDAL 和命令行调用的 GDAL 对同一个 COG 文件的解析行为可能不一致。import rasterio from rasterio.windows import Window # 按窗口读取避免把整幅数据一次性读进内存 with rasterio.open(tile_N00E000.tif) as ds: win Window(col_off0, row_off0, width2048, height2048) block ds.read(1, windowwin) print(block.shape, block.dtype) # 只看某个经纬度点的高程用 index 把坐标转成行列号 with rasterio.open(tile_N00E000.tif) as ds: row, col ds.index(12.5, 45.2) # 注意参数顺序是 (经度, 纬度) value ds.read(1)[row, col] print(该点高程:, value)ds.index()的参数顺序是先经度后纬度和日常口语的“纬度经度”相反这是新手最常弄反的一处。如果拿到的高程明显不合理先检查这里。Window的四个参数是列偏移、行偏移、宽、高全部以像素为单位不是坐标。3.2 gdalbuildvrt先建索引不搬数据虚拟栅格VRT是一个 XML 索引文件只记录“哪块数据在哪个文件的哪个位置”本身几乎不占空间。几十 GB 的全球瓦片用 VRT 几秒钟就能“拼”好之后所有读取、重投影、切片操作都按需从原始瓦片取数。# 先把所有瓦片路径写进一个清单文件一行一个 find ./tiles -name *.tif tiles.txt # 建虚拟栅格统一无效值分辨率取最高的那份 gdalbuildvrt \ -input_file_list tiles.txt \ -resolution highest \ -srcnodata -32768 \ -vrtnodata -32768 \ global.vrt # 查一下拼出来的范围对不对 gdalinfo -stats global.vrt | head -40-resolution highest表示按输入中分辨率最高的那份统一网格average和lowest分别取平均和最粗。如果各瓦片分辨率本来一致这三个参数结果相同混了 1 弧秒和 3 弧秒数据时highest会把粗数据插值放大文件虚拟尺寸变大但信息量没增加这时更该用-tr显式指定目标分辨率。-srcnodata和-vrtnodata成对出现前者告诉工具每份输入里哪个值是无效的后者决定输出里用什么值表示无效两边设成一样最省事。gdalinfo -stats会顺带算一遍全局统计这一步在 VRT 上跑等同于扫一遍全部原始数据几十 GB 的数据可能要几分钟。值得跑因为min / max能直接暴露单位是不是英尺、有没有异常值把范围撑爆。3.3 分块写出与 COG确认 VRT 没问题后落盘成一个自包含的文件。全球范围 1 弧秒全分辨率约合十几万乘几万像素未压缩体积能到 TB 级必须分块、压缩、加金字塔也就是常说的 COG云优化 GeoTIFF。它的特点是内部按 256 或 512 像素分块块内压缩并内置多级缩略图读取端可以只取需要的那一块。# 直接用 COG 驱动输出一步完成分块、压缩、金字塔 gdal_translate \ -of COG \ -co COMPRESSDEFLATE \ -co PREDICTOR2 \ -co BLOCKSIZE512 \ -co OVERVIEWSAUTO \ -co NUM_THREADSALL_CPUS \ global.vrt global_cog.tifPREDICTOR2针对整型数据做行差分预测配合 DEFLATE 通常能再压下去两三成浮点数据要用PREDICTOR3。BLOCKSIZE512是体积和随机读取效率的折中256 更省内存、512 读取次数更少浏览器端渲染一般 256 就够。OVERVIEWSAUTO让工具按 2 的幂次自动生成缩略层级缩略图不生成的话任何缩放浏览都会退化成全分辨率扫描卡到怀疑人生。如果不用 COG 驱动用 rasterio 手动写也可以好处是能顺手做无效值替换和类型转换import rasterio from rasterio.enums import Resampling with rasterio.open(global.vrt) as src: profile src.profile.copy() profile.update( driverGTiff, dtypeint16, nodata-32768, tiledTrue, blockxsize512, blockysize512, compressdeflate, predictor2, bigtiffIF_SAFER, # 超过 4GB 自动切 BigTIFF ) with rasterio.open(global_out.tif, w, **profile) as dst: # 逐块搬运内存占用只跟块大小有关 for _, window in dst.block_windows(1): data src.read(1, windowwindow) dst.write(data, 1, windowwindow) # 写完补金字塔缩放浏览才不会卡 dst.build_overviews([2, 4, 8, 16, 32], Resampling.average) dst.update_tags(nsrio_overview, resamplingaverage)bigtiffIF_SAFER这一行很关键全球范围数据几乎必然超过 4GB普通 TIFF 到 4GB 就写不进去了报错信息还特别含糊。build_overviews的重采样方式选average对高程数据来说比nearest更能保留地形起伏的平均特征如果这份数据要用来做分类或掩膜才改用nearest。3.4 重投影与重采样大部分情况不该转 Web Mercator习惯上把地形数据转成 Web MercatorEPSG:3857方便和在线底图叠加但这个投影在高纬度会把面积和形状放大得离谱格陵兰看起来比非洲还大就是它的“功劳”。做全球尺度的坡度、汇流、面积统计应该保留地理坐标 EPSG:4326或者转等积投影只有在做纯展示、且范围不超过中低纬度时才考虑 Web Mercator。# 保留地理坐标只做重采样把 1 弧秒降到 3 弧秒 gdalwarp \ -tr 0.0008333333 0.0008333333 \ -r average \ -multi -wo NUM_THREADSALL_CPUS \ -co COMPRESSDEFLATE -co TILEDYES -co BIGTIFFYES \ global.vrt global_3arcsec.tif # 需要等积投影做面积统计时换目标坐标系并把重采样改成 bilinear gdalwarp -t_srs projcea lat_ts30 datumWGS84 unitsm \ -tr 1000 1000 -r bilinear -multi \ global.vrt global_cea.tif-tr的单位跟着目标坐标系走地理坐标下是“度”投影坐标下是“米”写错一个数量级就是全图变糊或者文件大到打不开。重采样算法选择上average适合降分辨率避免抽样抽到噪点bilinear和cubic适合投影变换near只用于分类数据和掩膜用在高程上会产生台阶状伪影。-multi开多进程配合NUM_THREADSALL_CPUS能把重投影时间压下来一大截代价是内存占用上升机器内存小于数据量十分之一时别开。4. 从 DEM 派生产品晕渲、坡度、等高线的参数怎么定原始高程栅格能直接看的信息有限实际交付里九成是派生产品山体晕渲图、坡度图、坡向图、等高线。这几个工具的默认参数都不是为全球数据准备的直接跑出来的结果往往“能看但不对”参数差一点视觉效果和数值含义都变。4.1 晕渲图的方位角、高度角与垂直夸张晕渲的本质是给每个像素算一个光照强度三个参数决定效果光照方位角、光照高度角、垂直夸张系数。方位角默认 315 度西北方向打光这符合人眼长期形成的习惯——光从左上角来凸起看起来是凸起凹下去看起来是凹下去改成 135 度东南方向同一片山地会整体反转成“坑”这是视觉上的经典翻车点。# 先削减掉无效值再晕渲避免海域出现黑块 gdal_translate -a_nodata -32768 global.vrt dem_clean.tif gdaldem hillshade \ -az 315 -alt 45 -z 2 \ -compute_edges \ -co COMPRESSDEFLATE -co TILEDYES \ dem_clean.tif hillshade.tif-alt 45是光照高度角低于 30 会让阴影拉得过长、细节糊成一片高于 60 则立体感消失。-z 2是垂直夸张系数全球尺度上地形起伏相对水平距离很小不夸张几乎看不出山局部小区域几十公里用 1 就够夸到 5 以上会显得像针山。-compute_edges让边缘像素也参与计算不加的话瓦片接缝处会有一圈明显的黑线。全球范围直接跑晕渲输出仍是十几万像素宽普通看图软件打不开。常见做法是先降到 3 弧秒或更低再晕渲或者晕渲完再生成金字塔。另外晕渲结果是无符号字节型只有 0~255 的灰度不要再拿它做数值分析。4.2 坡度坡向单位换算与 NoData 传染坡度的单位有两种度0~90和百分比0~无穷。gdaldem slope默认输出度加-p输出百分比。地理坐标系下必须加-s参数告诉工具“一个度对应多少米”否则它会把度当米算算出来的坡度接近 0看上去整张图都是平原——这是使用地理坐标 DEM 时最典型的错误。# 地理坐标EPSG:4326必须给 -s值约为纬度 1 度对应的米数 gdal_dem_slope() { gdaldem slope dem_clean.tif slope.tif \ -s 111120 \ -compute_edges \ -co COMPRESSDEFLATE -co TILEDYES } # 坡向0~360 度正北为 0顺时针增加平地输出 -9999 gdaldem aspect dem_clean.tif aspect.tif -compute_edges-s 111120是赤道处 1 度纬度对应的米数在低中纬度误差可接受做高纬度精细分析时应该先重投影到适合该区域的投影坐标系再算坡度让三个方向的单位统一成米误差最小。另外-s对所有行用同一个值跨度很大的全球数据本身就存在系统性偏差这一点无法靠调参数消除只能靠分带处理。坡向在平地上没有意义工具会输出 -9999。这个值和真正的 NoData 值比如 -32768不是一回事做后续统计时要单独处理否则平均值会被拉偏。还有一个更隐蔽的问题原始 DEM 边缘是 NoData 时gdaldem会在边界内侧算出一圈虚假的陡坡-compute_edges只能缓解不能根除正确做法是裁剪时在 NoData 外围留出一圈缓冲算完再裁回来。4.3 等高线抽稀与平滑gdal_contour从栅格直接生成矢量线参数少但每个都要想清楚。等高距设太小全球数据能生成几千万条线文件打不开设太大地形细节丢失。一般 1 弧秒数据用 20~50 米等高距3 弧秒数据用 100 米做全球示意用 500 米。# 生成 100 米等高距带高程字段 gdal_contour \ -a elev \ -i 100 \ -f GPKG \ dem_clean.tif contours.gpkg-a elev指定把高程写成属性字段不加这个字段生成的线只有几何没有高度值后续没法做标注和分级。生成的线会非常密且锯齿明显用ogr2ogr做一次简化ogr2ogr -simplify 0.0005 contours_smooth.gpkg contours.gpkg-simplify的容差单位跟着坐标系走地理坐标下 0.0005 度约 55 米这个量级对全球图够用局部图要调到 0.00005。简化过度会让等高线自相交生成拓扑错误宁可多留一点锯齿也别过分抽稀。5. 全球地形 tif 处理中的避坑清单下面五条都是实际处理中反复出现的问题每一条都能让结果“看起来正常但数值全错”排查起来比崩溃更费时间。5.1 无效值没声明0 值悄悄混进统计现象统计平均高程得到几百米明显偏高晕渲图上大片海域是纯黑像是被挖空了。原因原始文件的nodata元数据缺失海面用 0 或 -32768 填充。0 被当成“海拔 0 米的陆地”参与均值计算-32768 被当成深不见底的海沟晕渲时算出极端坡度。解决读取时统一用maskedTrue或者在处理前用gdal_translate -a_nodata显式声明再用gdalwarp -dstnodata统一各瓦片的无效值。统计前先确认arr.count()和总像素数的比例是否合理全球陆地占比约 29%偏差太大就说明无效值没处理干净。5.2 PixelIsPoint 与 PixelIsArea 差半个像素现象拼接后接缝处有 1 像素宽的错位或者在和矢量边界叠加时整体偏移一点点放大才看得出来。原因GeoTIFF 里有个GTRasterTypeGeoKey标记PixelIsArea表示坐标对应像素左上角PixelIsPoint表示对应像素中心。两套数据混用时仿射变换会差半个像素全球尺度上就是几百米到几公里的偏移。解决用gdalinfo看输出里的Raster Type一行或者读ds.tags(ns)。混用前统一转成同一种类型gdal_translate -mo GTRasterTypeGeoKey1可以显式改写。这类偏移在地形上几乎看不出来叠加到建筑轮廓或者行政边界上就会暴露。5.3 各瓦片 CRS 不一致拼接后整体错位现象拼图中间出现一条明显的空白带或者数值断层单看每一块都正常。原因部分瓦片来自不同版本的同一系列产品或者经过别人重投影又改回来CRS 名义上都是 EPSG:4326但实际基准或轴顺序有差异。解决拼接前先跑一遍全量检查把每块的 CRS、分辨率、仿射参数、波段数导成一张表逐列比对。用gdalbuildvrt -strict让工具在分辨率不一致时直接报错而不是默默插值。# 批量导出关键元数据快速找出不一致的那几块 for f in ./tiles/*.tif; do gdalinfo $f | grep -E EPSG|Pixel Size|Origin|Size is echo ---- $f done5.4 直接 merge 几十 GB内存当场爆掉现象gdal_merge.py跑到一半被系统杀掉或者机器开始疯狂交换进度条几小时不动。原因合并工具会把输出尺寸的整幅数据在内存里组织一遍全球 1 弧秒全分辨率需要的内存远超普通工作站。解决改用 VRT 做索引再gdal_translate分块写出或者按纬度带分批合并再拼。gdalwarp的-multi和-wo NUM_THREADS要配合内存设置机器内存不足时宁可关掉多线程。写文件时加-co BIGTIFFYES避免 4GB 处突然报错。5.5 没建金字塔渲染端卡到打不开现象数据在桌面 GIS 里能打开但拖动一次要等十几秒Web 端直接超时用代码按窗口读取却很快。原因文件没有内建缩略图overview任何缩放操作都要按全分辨率扫描一遍IO 量随缩放级别暴涨。解决写出时加-co OVERVIEWSAUTO或用gdaladdo -r average补生成压缩方式选 DEFLATE 或 ZSTD别用不压缩的裸 TIFF。已经发布的文件补金字塔要注意先确认下游有没有缓存旧副本否则补了也白补。注意这五条的共同点是“不报错但结果错”。崩溃容易发现数值偏差往往要到出图或者交付评审时才暴露所以每完成一步就存一份中间产物和统计摘要是唯一可靠的后悔药。6. 把全球 DEM 压到能上线的体积块大小、压缩与金字塔的取舍全球范围的数据最终都要落到“别人能打开”这个现实约束上。几个 TB 的原始数据没人愿意用压缩到几十 GB 的 COG 才是可用状态。压缩方式的选择不只看压缩比还要看解码速度和下游工具链的支持面。压缩方式有无损典型压缩比解码速度适用场景DEFLATE无损中等快通用首选兼容性最好LZW无损略低于 DEFLATE快老工具链必须兼容时ZSTD无损高于 DEFLATE很快新版本 GDAL 环境追求吞吐JPEG有损高很快只做展示的底图不做数值分析高程数据绝对不能选有损压缩。JPEG 压缩会把相邻像素做频域变换几百米的高程差可能被抹平坡度算出来全是噪声而且这种损失在文件里看不出来属于典型的“黑匣子”问题。做纯展示的晕渲底图可以用原始高程一定用无损。块大小的取舍在于随机读取效率和文件开销。块越小随机读一个点需要传输的数据越少但块头开销和压缩效率都变差块越大反之。512×512 是通用值Web 端按瓦片请求时 256×256 更贴合。block 尺寸设置完后不要随意改改一次等于整个文件重写。金字塔层级按 2 的幂次生成一直到最长边小于 256 像素为止。层级太少缩小到洲际尺度会退化成读次高层再插值层级太多文件头部的缩略图索引本身也会占空间一般 5~7 级足够覆盖从全球到省级的缩放范围。重采样方式用average保留地形趋势用nearest保留极值和分类边界两者不要混用在同一份文件上否则不同缩放级别的图看起来像两份数据。项目输出分两类一类用于分析保留原始分辨率、无损压缩、带完整元数据另一类用于展示降采样到合适层级、生成晕渲、附带金字塔体积压到前者的一成以内。两类文件都保留下来比试图用一份文件兼顾两边省事得多我踩过的坑就是反复在同一个文件上改来改去最后连哪一版是准的都对不上。现在的习惯是每处理一步就写一份info.txt记录源文件清单、命令、输出尺寸和统计摘要和成品放在一起。过几个月回头查“这个数据为什么和图上不一样”这份摘要就是唯一能救命的线索。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/26 18:27:59 · 阅读 1.2k 全球地形TIF处理指南：GeoTIFF拼接、重采样与COG压缩
简介面向全球陆地与海底地形联合分析需求这份ETOPO1地形数据以GeoTIFF栅格形式提供约2km分辨率的全球DEM坐标系统为WGS1984可用于GIS制图、地形起伏分析、海陆一体化建模、三维可视化底图及地学教学实验适合具备基础GIS操作能力的研究者与学习者使用。压缩包为rar格式共3个文件、约266.01MB主体为tif栅格数据另附ovr金字塔概览文件与xml辅助元数据便于在常见GIS软件中快速浏览、缩放和读取坐标信息。已有3036人学习或下载说明其在全球尺度地形应用中具有较高参考价值。读者可据此直接开展全球高程与海底地形提取、坡度坡向计算、区域裁剪和专题图制作也可作为气候、生态、水文及海洋相关研究的统一空间底图减少多源地形数据拼接与坐标转换的准备工作。1. 全球地形数据tif格式为什么第一件事不是打开看图把一张 5°×5° 的全球地形数据 tif 直接拖进普通看图软件多半是转圈、全黑或者一片灰白分不出山脉和海沟。不是文件坏了是这类 tif 原本就不是为肉眼浏览准备的——它是带地理坐标的栅格矩阵像素值是高程米数覆盖范围可以很大文件也可能达到几十 GB。全球地形数据 tif 格式要解决的核心问题只有一个把起伏的地表切成规则网格每个格子记一个高程值再用 GeoTIFF 标签把这些格子钉回经纬度坐标上。它支撑地形可视化、水文汇流分析、选址与通视评估这类工作新手容易在“能打开”和“能用”之间卡很久熟手则更在意分幅规则、无效值、垂直基准这些边界参数。本章先讲清这类文件的结构与参数再落到拼接、重采样、派生产品和上线压缩的具体命令。2. 全球地形 tif 里到底装了什么拿到一批瓦片先别急着拼。花十分钟把单幅的元数据读出来比后面折腾几小时返工划算得多。全球 DEM 的分发方式是“一套产品切成上千个矩形块”块与块之间靠文件名的经纬度对齐任何一块的坐标参考、无效值或数据类型和其他块不一致整张拼图就会在接缝处出现台阶、黑边或者数值断层。2.1 GeoTIFF 的六个关键字段普通 tif 只存像素GeoTIFF 在文件头里多塞了一组地理标签。真正影响后续处理的是这六个宽高与波段数、数据类型、坐标参考系CRS、仿射变换六参数、无效值NoData、压缩与分块方式。仿射变换是最容易被忽略的一个它定义了“左上角像素的左上角在地球上的哪个点”以及“每个像素跨多少度”六个参数依次是左上角 x、x 方向像素宽、旋转项、左上角 y、旋转项、y 方向像素高通常是负值因为图像行号向下增长。import rasterio # 读一个瓦片的元数据不做任何计算 with rasterio.open(tile_N00E000.tif) as ds: print(尺寸:, ds.width, x, ds.height, 波段:, ds.count) print(数据类型:, ds.dtypes) # 常见 int16 / float32 print(坐标参考:, ds.crs) # 全球 DEM 多为 EPSG:4326 print(仿射变换:, ds.transform) # 六参数注意 e 为正、f 为负 print(无效值:, ds.nodata) # 没设就是 None后面必然出事 print(分块:, ds.block_shapes[0]) # (256, 256) 之类 print(压缩:, ds.profile.get(compress))这段代码跑完重点看两个地方。第一nodata如果打印出None说明这份数据没有声明无效值海面或空洞区可能是 0、-32768 或者极大负数统计均值、算坡度都会把无效区域算进去。第二transform里的a和e像素宽高如果不是同一个数量级或者a是 0.0083333 这种数字说明它是 30 弧秒级产品不是常见的 1 弧秒单位是“度”不是“米”这一点决定后面所有距离和坡度计算要不要换算。ds.profile这个字典值得留一份写出新文件时直接复制再改几个键比手写一堆参数稳。常见写法是profile src.profile.copy()然后改driver、dtype、compress、tiled、blockxsize、blockysize、nodata。改完再传给rasterio.open(..., w, **profile)元数据基本不会丢。2.2 分辨率不是固定米数1 弧秒、3 弧秒与瓦片命名全球地形产品常说的“30 米”“90 米”是近似说法真实单位是角秒。1 弧秒的纬度跨度在地表约 30.87 米子午线周长约 40007 公里除以 360 度再除以 3600 秒这个值全球基本恒定但 1 弧秒的经度跨度要乘以纬度的余弦赤道处约 30.87 米到北纬 60 度只剩约 15.4 米。所以同一份 1 弧秒数据在高纬地区像素在地面上是“窄长条”不是正方形。这个差异带来的直接后果是任何以“米”为单位的计算都不能直接拿像素个数乘 30。算面积要先投影到等积坐标系算坡度要让 GDAL 知道纬度换算比例算缓冲区要按实际地面距离而不是行列数。很多人第一次做全球尺度分析发现高纬区域的坡度整体偏大、面积整体偏小根子就在这里。产品级别角秒赤道地面分辨率全球瓦片数量级单幅 5°×5° 体积量级粗分辨率全球网格30 弧秒约 900 米数百几 MB中分辨率3 弧秒约 90 米数千数十 MB高分辨率1 弧秒约 30 米上万上百 MB瓦片命名一般是“半球 纬度 半球 经度”的形式比如N00E000、S45W120数字是瓦片左下角或左上角看产品约定的整度坐标。拼接前务必确认三件事命名用的是左下角还是左上角、瓦片边长是 1 度还是 5 度、经度是 0~360 还是 -180~180。这三件事任一搞错拼出来就是错位或者中间一条空白带。2.3 垂直基准与数据类型高程到底相对哪个面高程值本身不带说明它一定相对某个参考面。公开全球 DEM 大多给的是正高相对大地水准面常见 EGM96 模型而 GPS 和很多测绘软件默认输出的是椭球高相对 WGS84 椭球。两者在全球范围内的差值大约在负一百米到正八十多米之间波动具体到某个区域可能是几十米。做洪水淹没、管线纵断面、机场净空这类对绝对高程敏感的分析时混用两套基准就是事故级错误做可视化山体晕渲则影响不大。import numpy as np import rasterio with rasterio.open(tile_N00E000.tif) as ds: arr ds.read(1, maskedTrue) # maskedTrue 自动把 nodata 屏蔽 print(最小值:, float(arr.min()), 最大值:, float(arr.max())) print(无效值占比:, 1 - arr.count() / arr.size) # 快速判断单位全球陆地最高约 8848 米 # 若最大值接近 29000基本可以确定单位是英尺需要乘 0.3048数据类型同样要看清。int16的取值范围是 -32768 到 32767装全球陆地高程绰绰有余海面常用负值表示水深或者干脆用 -32768 当无效值。float32支持小数代价是体积翻倍。如果一份int16数据里出现了 -32768 且nodata已声明它是无效值如果nodata没声明它就会被当成“负三万米的海沟”参与运算晕渲图上会出现一片突兀的黑块。提示拿到任何一份全球地形 tif先跑一遍min / max / count / nodata四个值再决定后面怎么处理。这一步不做后面所有分析结果都不可信。3. 用 GDAL 与 rasterio 把全球瓦片拼成一张可用底图拼接全球地形有两条路物理合并成一个大文件或者用虚拟栅格只建索引。前者适合做最终交付后者适合先看清全局、验证瓦片没错位。我一般的顺序是先用虚拟栅格扫一遍找问题确认无误再决定落盘的范围和分辨率避免上来就 merge 几十 GB 然后中途内存爆掉。3.1 环境与最小读取Python 侧直接用pip install rasterio它自带打包好的 GDAL 二进制省去编译麻烦需要命令行工具就装gdal的 conda 包。两边版本尽量对齐否则rasterio调用的 GDAL 和命令行调用的 GDAL 对同一个 COG 文件的解析行为可能不一致。import rasterio from rasterio.windows import Window # 按窗口读取避免把整幅数据一次性读进内存 with rasterio.open(tile_N00E000.tif) as ds: win Window(col_off0, row_off0, width2048, height2048) block ds.read(1, windowwin) print(block.shape, block.dtype) # 只看某个经纬度点的高程用 index 把坐标转成行列号 with rasterio.open(tile_N00E000.tif) as ds: row, col ds.index(12.5, 45.2) # 注意参数顺序是 (经度, 纬度) value ds.read(1)[row, col] print(该点高程:, value)ds.index()的参数顺序是先经度后纬度和日常口语的“纬度经度”相反这是新手最常弄反的一处。如果拿到的高程明显不合理先检查这里。Window的四个参数是列偏移、行偏移、宽、高全部以像素为单位不是坐标。3.2 gdalbuildvrt先建索引不搬数据虚拟栅格VRT是一个 XML 索引文件只记录“哪块数据在哪个文件的哪个位置”本身几乎不占空间。几十 GB 的全球瓦片用 VRT 几秒钟就能“拼”好之后所有读取、重投影、切片操作都按需从原始瓦片取数。# 先把所有瓦片路径写进一个清单文件一行一个 find ./tiles -name *.tif tiles.txt # 建虚拟栅格统一无效值分辨率取最高的那份 gdalbuildvrt \ -input_file_list tiles.txt \ -resolution highest \ -srcnodata -32768 \ -vrtnodata -32768 \ global.vrt # 查一下拼出来的范围对不对 gdalinfo -stats global.vrt | head -40-resolution highest表示按输入中分辨率最高的那份统一网格average和lowest分别取平均和最粗。如果各瓦片分辨率本来一致这三个参数结果相同混了 1 弧秒和 3 弧秒数据时highest会把粗数据插值放大文件虚拟尺寸变大但信息量没增加这时更该用-tr显式指定目标分辨率。-srcnodata和-vrtnodata成对出现前者告诉工具每份输入里哪个值是无效的后者决定输出里用什么值表示无效两边设成一样最省事。gdalinfo -stats会顺带算一遍全局统计这一步在 VRT 上跑等同于扫一遍全部原始数据几十 GB 的数据可能要几分钟。值得跑因为min / max能直接暴露单位是不是英尺、有没有异常值把范围撑爆。3.3 分块写出与 COG确认 VRT 没问题后落盘成一个自包含的文件。全球范围 1 弧秒全分辨率约合十几万乘几万像素未压缩体积能到 TB 级必须分块、压缩、加金字塔也就是常说的 COG云优化 GeoTIFF。它的特点是内部按 256 或 512 像素分块块内压缩并内置多级缩略图读取端可以只取需要的那一块。# 直接用 COG 驱动输出一步完成分块、压缩、金字塔 gdal_translate \ -of COG \ -co COMPRESSDEFLATE \ -co PREDICTOR2 \ -co BLOCKSIZE512 \ -co OVERVIEWSAUTO \ -co NUM_THREADSALL_CPUS \ global.vrt global_cog.tifPREDICTOR2针对整型数据做行差分预测配合 DEFLATE 通常能再压下去两三成浮点数据要用PREDICTOR3。BLOCKSIZE512是体积和随机读取效率的折中256 更省内存、512 读取次数更少浏览器端渲染一般 256 就够。OVERVIEWSAUTO让工具按 2 的幂次自动生成缩略层级缩略图不生成的话任何缩放浏览都会退化成全分辨率扫描卡到怀疑人生。如果不用 COG 驱动用 rasterio 手动写也可以好处是能顺手做无效值替换和类型转换import rasterio from rasterio.enums import Resampling with rasterio.open(global.vrt) as src: profile src.profile.copy() profile.update( driverGTiff, dtypeint16, nodata-32768, tiledTrue, blockxsize512, blockysize512, compressdeflate, predictor2, bigtiffIF_SAFER, # 超过 4GB 自动切 BigTIFF ) with rasterio.open(global_out.tif, w, **profile) as dst: # 逐块搬运内存占用只跟块大小有关 for _, window in dst.block_windows(1): data src.read(1, windowwindow) dst.write(data, 1, windowwindow) # 写完补金字塔缩放浏览才不会卡 dst.build_overviews([2, 4, 8, 16, 32], Resampling.average) dst.update_tags(nsrio_overview, resamplingaverage)bigtiffIF_SAFER这一行很关键全球范围数据几乎必然超过 4GB普通 TIFF 到 4GB 就写不进去了报错信息还特别含糊。build_overviews的重采样方式选average对高程数据来说比nearest更能保留地形起伏的平均特征如果这份数据要用来做分类或掩膜才改用nearest。3.4 重投影与重采样大部分情况不该转 Web Mercator习惯上把地形数据转成 Web MercatorEPSG:3857方便和在线底图叠加但这个投影在高纬度会把面积和形状放大得离谱格陵兰看起来比非洲还大就是它的“功劳”。做全球尺度的坡度、汇流、面积统计应该保留地理坐标 EPSG:4326或者转等积投影只有在做纯展示、且范围不超过中低纬度时才考虑 Web Mercator。# 保留地理坐标只做重采样把 1 弧秒降到 3 弧秒 gdalwarp \ -tr 0.0008333333 0.0008333333 \ -r average \ -multi -wo NUM_THREADSALL_CPUS \ -co COMPRESSDEFLATE -co TILEDYES -co BIGTIFFYES \ global.vrt global_3arcsec.tif # 需要等积投影做面积统计时换目标坐标系并把重采样改成 bilinear gdalwarp -t_srs projcea lat_ts30 datumWGS84 unitsm \ -tr 1000 1000 -r bilinear -multi \ global.vrt global_cea.tif-tr的单位跟着目标坐标系走地理坐标下是“度”投影坐标下是“米”写错一个数量级就是全图变糊或者文件大到打不开。重采样算法选择上average适合降分辨率避免抽样抽到噪点bilinear和cubic适合投影变换near只用于分类数据和掩膜用在高程上会产生台阶状伪影。-multi开多进程配合NUM_THREADSALL_CPUS能把重投影时间压下来一大截代价是内存占用上升机器内存小于数据量十分之一时别开。4. 从 DEM 派生产品晕渲、坡度、等高线的参数怎么定原始高程栅格能直接看的信息有限实际交付里九成是派生产品山体晕渲图、坡度图、坡向图、等高线。这几个工具的默认参数都不是为全球数据准备的直接跑出来的结果往往“能看但不对”参数差一点视觉效果和数值含义都变。4.1 晕渲图的方位角、高度角与垂直夸张晕渲的本质是给每个像素算一个光照强度三个参数决定效果光照方位角、光照高度角、垂直夸张系数。方位角默认 315 度西北方向打光这符合人眼长期形成的习惯——光从左上角来凸起看起来是凸起凹下去看起来是凹下去改成 135 度东南方向同一片山地会整体反转成“坑”这是视觉上的经典翻车点。# 先削减掉无效值再晕渲避免海域出现黑块 gdal_translate -a_nodata -32768 global.vrt dem_clean.tif gdaldem hillshade \ -az 315 -alt 45 -z 2 \ -compute_edges \ -co COMPRESSDEFLATE -co TILEDYES \ dem_clean.tif hillshade.tif-alt 45是光照高度角低于 30 会让阴影拉得过长、细节糊成一片高于 60 则立体感消失。-z 2是垂直夸张系数全球尺度上地形起伏相对水平距离很小不夸张几乎看不出山局部小区域几十公里用 1 就够夸到 5 以上会显得像针山。-compute_edges让边缘像素也参与计算不加的话瓦片接缝处会有一圈明显的黑线。全球范围直接跑晕渲输出仍是十几万像素宽普通看图软件打不开。常见做法是先降到 3 弧秒或更低再晕渲或者晕渲完再生成金字塔。另外晕渲结果是无符号字节型只有 0~255 的灰度不要再拿它做数值分析。4.2 坡度坡向单位换算与 NoData 传染坡度的单位有两种度0~90和百分比0~无穷。gdaldem slope默认输出度加-p输出百分比。地理坐标系下必须加-s参数告诉工具“一个度对应多少米”否则它会把度当米算算出来的坡度接近 0看上去整张图都是平原——这是使用地理坐标 DEM 时最典型的错误。# 地理坐标EPSG:4326必须给 -s值约为纬度 1 度对应的米数 gdal_dem_slope() { gdaldem slope dem_clean.tif slope.tif \ -s 111120 \ -compute_edges \ -co COMPRESSDEFLATE -co TILEDYES } # 坡向0~360 度正北为 0顺时针增加平地输出 -9999 gdaldem aspect dem_clean.tif aspect.tif -compute_edges-s 111120是赤道处 1 度纬度对应的米数在低中纬度误差可接受做高纬度精细分析时应该先重投影到适合该区域的投影坐标系再算坡度让三个方向的单位统一成米误差最小。另外-s对所有行用同一个值跨度很大的全球数据本身就存在系统性偏差这一点无法靠调参数消除只能靠分带处理。坡向在平地上没有意义工具会输出 -9999。这个值和真正的 NoData 值比如 -32768不是一回事做后续统计时要单独处理否则平均值会被拉偏。还有一个更隐蔽的问题原始 DEM 边缘是 NoData 时gdaldem会在边界内侧算出一圈虚假的陡坡-compute_edges只能缓解不能根除正确做法是裁剪时在 NoData 外围留出一圈缓冲算完再裁回来。4.3 等高线抽稀与平滑gdal_contour从栅格直接生成矢量线参数少但每个都要想清楚。等高距设太小全球数据能生成几千万条线文件打不开设太大地形细节丢失。一般 1 弧秒数据用 20~50 米等高距3 弧秒数据用 100 米做全球示意用 500 米。# 生成 100 米等高距带高程字段 gdal_contour \ -a elev \ -i 100 \ -f GPKG \ dem_clean.tif contours.gpkg-a elev指定把高程写成属性字段不加这个字段生成的线只有几何没有高度值后续没法做标注和分级。生成的线会非常密且锯齿明显用ogr2ogr做一次简化ogr2ogr -simplify 0.0005 contours_smooth.gpkg contours.gpkg-simplify的容差单位跟着坐标系走地理坐标下 0.0005 度约 55 米这个量级对全球图够用局部图要调到 0.00005。简化过度会让等高线自相交生成拓扑错误宁可多留一点锯齿也别过分抽稀。5. 全球地形 tif 处理中的避坑清单下面五条都是实际处理中反复出现的问题每一条都能让结果“看起来正常但数值全错”排查起来比崩溃更费时间。5.1 无效值没声明0 值悄悄混进统计现象统计平均高程得到几百米明显偏高晕渲图上大片海域是纯黑像是被挖空了。原因原始文件的nodata元数据缺失海面用 0 或 -32768 填充。0 被当成“海拔 0 米的陆地”参与均值计算-32768 被当成深不见底的海沟晕渲时算出极端坡度。解决读取时统一用maskedTrue或者在处理前用gdal_translate -a_nodata显式声明再用gdalwarp -dstnodata统一各瓦片的无效值。统计前先确认arr.count()和总像素数的比例是否合理全球陆地占比约 29%偏差太大就说明无效值没处理干净。5.2 PixelIsPoint 与 PixelIsArea 差半个像素现象拼接后接缝处有 1 像素宽的错位或者在和矢量边界叠加时整体偏移一点点放大才看得出来。原因GeoTIFF 里有个GTRasterTypeGeoKey标记PixelIsArea表示坐标对应像素左上角PixelIsPoint表示对应像素中心。两套数据混用时仿射变换会差半个像素全球尺度上就是几百米到几公里的偏移。解决用gdalinfo看输出里的Raster Type一行或者读ds.tags(ns)。混用前统一转成同一种类型gdal_translate -mo GTRasterTypeGeoKey1可以显式改写。这类偏移在地形上几乎看不出来叠加到建筑轮廓或者行政边界上就会暴露。5.3 各瓦片 CRS 不一致拼接后整体错位现象拼图中间出现一条明显的空白带或者数值断层单看每一块都正常。原因部分瓦片来自不同版本的同一系列产品或者经过别人重投影又改回来CRS 名义上都是 EPSG:4326但实际基准或轴顺序有差异。解决拼接前先跑一遍全量检查把每块的 CRS、分辨率、仿射参数、波段数导成一张表逐列比对。用gdalbuildvrt -strict让工具在分辨率不一致时直接报错而不是默默插值。# 批量导出关键元数据快速找出不一致的那几块 for f in ./tiles/*.tif; do gdalinfo $f | grep -E EPSG|Pixel Size|Origin|Size is echo ---- $f done5.4 直接 merge 几十 GB内存当场爆掉现象gdal_merge.py跑到一半被系统杀掉或者机器开始疯狂交换进度条几小时不动。原因合并工具会把输出尺寸的整幅数据在内存里组织一遍全球 1 弧秒全分辨率需要的内存远超普通工作站。解决改用 VRT 做索引再gdal_translate分块写出或者按纬度带分批合并再拼。gdalwarp的-multi和-wo NUM_THREADS要配合内存设置机器内存不足时宁可关掉多线程。写文件时加-co BIGTIFFYES避免 4GB 处突然报错。5.5 没建金字塔渲染端卡到打不开现象数据在桌面 GIS 里能打开但拖动一次要等十几秒Web 端直接超时用代码按窗口读取却很快。原因文件没有内建缩略图overview任何缩放操作都要按全分辨率扫描一遍IO 量随缩放级别暴涨。解决写出时加-co OVERVIEWSAUTO或用gdaladdo -r average补生成压缩方式选 DEFLATE 或 ZSTD别用不压缩的裸 TIFF。已经发布的文件补金字塔要注意先确认下游有没有缓存旧副本否则补了也白补。注意这五条的共同点是“不报错但结果错”。崩溃容易发现数值偏差往往要到出图或者交付评审时才暴露所以每完成一步就存一份中间产物和统计摘要是唯一可靠的后悔药。6. 把全球 DEM 压到能上线的体积块大小、压缩与金字塔的取舍全球范围的数据最终都要落到“别人能打开”这个现实约束上。几个 TB 的原始数据没人愿意用压缩到几十 GB 的 COG 才是可用状态。压缩方式的选择不只看压缩比还要看解码速度和下游工具链的支持面。压缩方式有无损典型压缩比解码速度适用场景DEFLATE无损中等快通用首选兼容性最好LZW无损略低于 DEFLATE快老工具链必须兼容时ZSTD无损高于 DEFLATE很快新版本 GDAL 环境追求吞吐JPEG有损高很快只做展示的底图不做数值分析高程数据绝对不能选有损压缩。JPEG 压缩会把相邻像素做频域变换几百米的高程差可能被抹平坡度算出来全是噪声而且这种损失在文件里看不出来属于典型的“黑匣子”问题。做纯展示的晕渲底图可以用原始高程一定用无损。块大小的取舍在于随机读取效率和文件开销。块越小随机读一个点需要传输的数据越少但块头开销和压缩效率都变差块越大反之。512×512 是通用值Web 端按瓦片请求时 256×256 更贴合。block 尺寸设置完后不要随意改改一次等于整个文件重写。金字塔层级按 2 的幂次生成一直到最长边小于 256 像素为止。层级太少缩小到洲际尺度会退化成读次高层再插值层级太多文件头部的缩略图索引本身也会占空间一般 5~7 级足够覆盖从全球到省级的缩放范围。重采样方式用average保留地形趋势用nearest保留极值和分类边界两者不要混用在同一份文件上否则不同缩放级别的图看起来像两份数据。项目输出分两类一类用于分析保留原始分辨率、无损压缩、带完整元数据另一类用于展示降采样到合适层级、生成晕渲、附带金字塔体积压到前者的一成以内。两类文件都保留下来比试图用一份文件兼顾两边省事得多我踩过的坑就是反复在同一个文件上改来改去最后连哪一版是准的都对不上。现在的习惯是每处理一步就写一份info.txt记录源文件清单、命令、输出尺寸和统计摘要和成品放在一起。过几个月回头查“这个数据为什么和图上不一样”这份摘要就是唯一能救命的线索。希望帮到你。本文还有配套的精品资源点击获取