1. 坐标查高程这件事，卡住新手的从来不是 GDAL 的 API

用 GDAL 拿地理坐标去 DEM 里查高程值，这个需求听起来简单得像调用一个函数：给个经纬度，返回一个海拔数字。我最早接触这个需求的时候也是这么想的，结果第一版代码跑出来，同一座山头在相邻两个点位上差了将近 40 米，还有一批点直接返回 0。后来花了整整两个下午才把整条链路捋清楚——问题根本不在 GDAL 会不会用，而在于对"坐标"这两个字的理解有偏差。

这篇文章讲的就是这条链路：DEM 数据源怎么勘察、地理坐标怎么换算成栅格的行列号、单点和批量查询分别怎么写、插值该选哪一种、以及那些查询结果看着正常其实是错的情况。适合已经能跑通 GDAL 打开文件、但对"为什么查出来的高程不准"心里没底的人；也适合完全没接触过栅格数据、只想照着抄一份可用代码的人。代码基于 Python 的 GDAL 绑定，换成 C++ 或 Java 绑定，逻辑完全一致。

先把最核心的一句话放在这里：用地理坐标查 DEM 高程，本质是一次带投影换算的像素寻址，再加一次数值采样。前半句决定你找没找对位置，后半句决定你取得准不准。绝大多数翻车都发生在前半句，而前半句恰恰是 GDAL 的 API 文档里最容易被跳过的部分。

1.1 一次查询实际发生了什么：从经纬度到像素下标的完整链路

很多人以为band.ReadAsArray(x, y, 1, 1)里的 x、y 就是经纬度，这是个非常典型的误解。栅格数据的存储结构是二维数组，索引是行列号，原点在左上角，行号向下增大。你手里的经纬度是一套球面坐标，两者之间隔着至少两层转换：第一层是坐标参考系转换，把经纬度投影到 DEM 自己使用的平面坐标系；第二层是仿射变换，把平面坐标映射成浮点形式的行列号。

完整的链路是这样的：输入是 WGS84 经纬度 (lon, lat)；先做坐标转换，得到 DEM 坐标系下的 (X, Y)，单位可能是米也可能是度，取决于 DEM 的投影；再用 GeoTransform 反算出浮点行列号 (col_f, row_f)；然后根据采样策略决定是四舍五入取整（最近邻），还是取周围四个像素做加权（双线性）；最后取到的值还要判断是不是 NoData，是的话就得换策略。

这条链路里每一步都可能出问题。第一层如果投影搞反了，查询位置会直接跳出图幅，返回空值；第二层如果 GeoTransform 的符号搞错了，行号方向会整个反过来，你查山顶拿到的是山脚；第三层插值选错，在平坦区域看不出差别，一到陡坡就露馅；第四层 NoData 没判，你会得到一堆看似合理的 0 米或 -32767 米。我在实际项目中见过最离谱的一次，是有人把行号当成了 x、列号当成了 y，代码跑得通、返回值也不报错，但整条路线的海拔剖面看上去完全不像地形——因为他是沿着图幅的对角线在采样。

所以别急着写查询函数。先花十分钟把下面两件事确认清楚，后面能省掉几个小时的调试。

1.2 三种坐标身份不清，后面全是无用功

在 DEM 这个场景里，"坐标"至少有三种身份，混用是新手最常见的错误来源。

第一种是地理坐标，也就是经纬度，单位是度，EPSG:4326 是它的常见代号。手机定位、地图上点选、GNSS 设备输出的都是这一种。它的特点是数值范围小（经度 -180 到 180，纬度 -90 到 90），但你绝对不能拿它直接去和以米为单位的 GeoTransform 做除法。

第二种是投影坐标，单位通常是米，比如各种 UTM 分带或者高斯投影。DEM 产品为了控制面积变形，绝大多数都用投影坐标存储。这一种才是和 GeoTransform 同一个量纲的，反算行列号必须在这一层做。

第三种是像素坐标，也就是浮点行列号，或者是"像素中心／像素角点"这种细分概念。GDAL 的 GeoTransform 描述的是像素角点（左上角）的位置，而实际采样一般按像素中心来算，这两者之间差了半个像素。这个半像素在 30 米分辨率的 DEM 上就是 15 米的位置偏差，在陡坡上足以造成几米的高程误差。

三种坐标的对应关系我整理成一张表，调试的时候对着看会快很多：

还有第四个隐形的坑：GDAL 在 EPSG:4326 下的轴顺序。很多传统 GIS 库按"经度在前、纬度在后"处理，但 OGC 规范里 EPSG:4326 的定义是"纬度在前、经度在后"。GDAL 3.x 之后默认遵循规范，于是TransformPoint(lon, lat)会得到一个经纬度对调的结果。解决办法是在构造 SpatialReference 之后调用SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)，强制回到"经度在前"的习惯。这一行代码看起来无关紧要，但少了它，你查到的可能是地球另一侧某个点的海拔。

2. DEM 数据源的预处理勘察：先看清再动手

拿到一景 DEM，我的习惯是先跑一遍勘察，把关键元数据抄在纸上或者写进注释，再动手写查询逻辑。这一步花不了几分钟，但能挡掉后面八成的诡异问题。

2.1 gdalinfo 看不到但必须知道的字段

命令行敲一句gdalinfo your_dem.tif，输出信息很长，但真正要盯住的就那么几行。第一是Size is X, Y，也就是列数和行数，这决定了你反算出来的行列号有没有越界。第二是Coordinate System is后面那段 WKT，重点看末尾的 EPSG 代号或者投影名称，判断是地理坐标系还是投影坐标系。第三是Origin，这就是 GeoTransform 的前两个数。第四是Pixel Size，对应 GeoTransform 的第 2 和第 6 个数，注意第二个数通常是正的、第六个数通常是负的，这个符号关系就是"行号向下增大"的数学表达。第五是NoData Value，有些 DEM 产品这一项是空的，那就得靠值域统计去推断。

除了这些，还有两个信息 gdalinfo 不一定直观展示，但非常关键。一个是数据类型，Int16 和 Float32 在处理上是不同的，前者遇到 0.5 的插值结果会被截断，如果你做双线性插值却忘了先转成浮点，精度会白白损失。另一个是波段数，绝大多数 DEM 只有一个波段，但偶尔会遇到带掩膜波段或者多波段的产品，GetRasterBand(1)拿的不一定是高程。

用 Python 直接把这几项打出来，比看命令行输出更省事：

from osgeo import gdal, osr ds = gdal.Open("dem.tif", gdal.GA_ReadOnly) band = ds.GetRasterBand(1) print("尺寸:", ds.RasterXSize, ds.RasterYSize) print("波段数:", ds.RasterCount) print("数据类型:", gdal.GetDataTypeName(band.DataType)) print("GeoTransform:", ds.GetGeoTransform()) print("NoData:", band.GetNoDataValue()) print("最小值/最大值:", band.ComputeRasterMinMax()) print("投影:", ds.GetProjection()[:120]) srs = osr.SpatialReference() srs.ImportFromWkt(ds.GetProjection()) print("是否地理坐标系:", srs.IsGeographic()) print("EPSG:", srs.GetAuthorityCode(None))

这段代码最重要的一行其实是srs.IsGeographic()。它返回 True 意味着 GeoTransform 里的 Pixel Size 单位是度，返回 False 则单位大概率是米。这个判断直接决定了你后面坐标转换的目标坐标系怎么设。我见过有人在投影坐标的 DEM 上做经纬度转换，转出来的数值是几千几万，拿去反算行列号直接越界；也见过反过来的，把米当成度去算，行列号永远是 0。

2.2 GeoTransform 六个数的物理含义与行列号反算公式

GeoTransform 是一个长度为 6 的元组，GDAL 里的定义是：

gt = (X0, dX, dY_rot, Y0, dX_rot, dY)

• gt[0]：栅格左上角像素左上角点的 X 坐标
• gt[1]：单个像素在 X 方向的宽度（列步长），一般是正数
• gt[2]：X 方向的旋转项，正北向的 DEM 这一项是 0
• gt[3]：栅格左上角像素左上角点的 Y 坐标
• gt[4]：Y 方向的旋转项，正北向的 DEM 这一项是 0
• gt[5]：单个像素在 Y 方向的高度（行步长），一般是负数

最容易被忽略的是gt[5]是负数这件事。因为栅格原点在左上角，行号越大位置越靠南（Y 值越小），所以行步长必须为负。如果你在公式里忘了这个符号，直接用(Y - gt[3]) / abs(gt[5])，得到的行号是镜像的——查北边的点，拿到的是南边的值。

正算公式，已知行列号求坐标（取像素中心）：

X = gt[0] + (col + 0.5) * gt[1] + (row + 0.5) * gt[2] Y = gt[3] + (col + 0.5) * gt[4] + (row + 0.5) * gt[5]

反算公式，已知坐标求浮点行列号：

col_f = (X - gt[0]) / gt[1] row_f = (Y - gt[3]) / gt[5]

注意反算出来的是浮点值，且这个浮点值对应的是像素的左上角点，不是中心。要得到以像素中心为基准的浮点行列号，两个方向各减 0.5：

col_c = col_f - 0.5 row_c = row_f - 0.5

这个 0.5 的偏移是双线性插值正确性的前提。如果你的采样策略是最近邻，那取整的时候用round还是int差别不大（在像素边界附近会有分歧）；但只要你打算做插值，就必须用中心化之后的浮点行列号，否则所有结果都会系统性地偏向西北方向半个像素。

对于gt[2]和gt[4]不为 0 的斜仿射数据（比如某些航摄影像衍生的 DEM），上面的简化反算公式是不成立的，需要用矩阵求逆。工程上遇到这种数据，直接调gdal.InvGeoTransform()更稳妥：

inv_gt, ok = gdal.InvGeoTransform(ds.GetGeoTransform()) # 返回的是 (col_f, row_f)，仍是左上角基准 col_f, row_f = gdal.ApplyGeoTransform(inv_gt, X, Y)

我建议不管数据有没有旋转，都统一用InvGeoTransform+ApplyGeoTransform这套组合。代码多两行，但省去了自己判断旋转项的心智负担，也不会在换了数据源之后突然失效。

2.3 投影不一致时的处理顺序与陷阱

真实的业务输入几乎不会是 DEM 自己的投影坐标。用户在手机上点一下给的是 WGS84 经纬度，GNSS 轨迹是 WGS84，某些行业系统内部用的是 CGCS2000 或者地方坐标系。所以坐标转换是绕不过去的一步，而这一步的顺序有讲究。

正确的做法是：先把输入坐标转到 DEM 的坐标系，再做行列号反算。顺序反过来是不行的，因为行列号只在特定投影下才有意义。转换用osr.CoordinateTransformation，源坐标系是输入坐标的坐标系，目标坐标系从 DEM 的 WKT 里导入。

def build_transformer(ds, src_epsg=4326): src = osr.SpatialReference() src.ImportFromEPSG(src_epsg) src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) dst = osr.SpatialReference() dst.ImportFromWkt(ds.GetProjection()) dst.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) return osr.CoordinateTransformation(src, dst)

这里有两个细节值得展开。第一是SetAxisMappingStrategy必须对源和目标都设置，只设一边在某些 GDAL 版本上依然会出现轴序混乱。第二是TransformPoint返回的是三元组(X, Y, Z)，第三个值是高程，但它是坐标转换过程中带出来的椭球高或者直接是 0，绝对不是你要的 DEM 高程。这一点坑过不少人，因为返回 0 看起来毫无破绽。

还有一个更隐蔽的坑：如果你的 DEM 本身是地理坐标系（Pixel Size 单位是度），那么转换后的 X、Y 也是度。这个时候反算出来的行列号在数学上是对的，但你在做双线性插值时，权重计算所依据的"距离"其实是角度而非米。在低纬度地区，经度方向 1 度约 111 公里，纬度方向也差不多；但到了高纬度，经度方向的 1 度会缩短到几十公里甚至更少。如果 DEM 覆盖范围不大、纬度跨度小，这个差异可以忽略；但如果你的数据横跨十几个纬度，插值权重的各向异性就会带来额外误差。稳妥的做法是：插值在像素坐标系里做，不要在投影坐标系里做，这样就不受单位影响。

3. Python 实现：单点查询的三层写法

代码部分我分三层递进：最短能跑的版本、工程上更推荐的版本、以及一个用来交叉验证的第三方写法。三层各有适用场景，不是谁替代谁的关系。

3.1 最近邻版本：最短代码与它的适用边界

最近邻就是四舍五入取最近的那个像素，不做加权。代码只有十几行，跑得也快，很多人第一版就是这么写的。

from osgeo import gdal, osr def get_elevation_nn(dem_path, lon, lat, src_epsg=4326): ds = gdal.Open(dem_path, gdal.GA_ReadOnly) if ds is None: raise IOError("DEM 打开失败: %s" % dem_path) band = ds.GetRasterBand(1) gt = ds.GetGeoTransform() nodata = band.GetNoDataValue() src = osr.SpatialReference() src.ImportFromEPSG(src_epsg) src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) dst = osr.SpatialReference() dst.ImportFromWkt(ds.GetProjection()) dst.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) ct = osr.CoordinateTransformation(src, dst) X, Y, _ = ct.TransformPoint(lon, lat) col = int((X - gt[0]) / gt[1]) row = int((Y - gt[3]) / gt[5]) if col < 0 or row < 0 or col >= ds.RasterXSize or row >= ds.RasterYSize: return None value = band.ReadAsArray(col, row, 1, 1)[0, 0] if nodata is not None and value == nodata: return None return float(value)

这段代码能跑，但它有三个明确的局限。

第一，int()截断不是四舍五入。对于负数行列号（虽然已经越界判断了），以及 0.9 这种本该进位的情况，int的行为和round不一样。虽然越界检查兜住了大部分场景，但在图幅边缘一两个像素内会有系统性偏移。

第二，每次都gdal.Open一遍。单点调用无所谓，但如果你要查一万个点，这个开销会非常可观——打开文件、解析头信息、建立内部缓存，这些操作重复一万次的代价远超你的想象。

第三，最近邻的精度在陡坡区域很差。30 米分辨率的 DEM，一个像素内的高差在山区可以达到二三十米。你查的那个点恰好在像素边缘还是中心，结果可能差出十几米。对于只做粗略展示的场景够用，但凡涉及坡度计算、剖面分析、淹没模拟，最近邻都不能用。

3.2 双线性插值版本：DEM 高程取值的常态选择

双线性插值取周围四个像素做加权平均，权重由目标点在像素内的相对位置决定。它假设高程在像素之间是线性变化的，这个假设在平坦和缓坡区域相当合理，在陡崖、山脊线上会低估峰值、高估谷值，但整体误差远小于最近邻。

实现的关键是先算出以像素中心为基准的浮点行列号，然后分离整数部分和小数部分：

def get_elevation_bilinear(ds, band, lon, lat, ct, gt, nodata=None): X, Y, _ = ct.TransformPoint(lon, lat) # 以像素左上角为基准的浮点行列号 col_f = (X - gt[0]) / gt[1] row_f = (Y - gt[3]) / gt[5] # 转成以像素中心为基准 col_c = col_f - 0.5 row_c = row_f - 0.5 c0 = int(col_c // 1) r0 = int(row_c // 1) dc = col_c - c0 dr = row_c - r0 if c0 < 0 or r0 < 0 or c0 + 1 >= ds.RasterXSize or r0 + 1 >= ds.RasterYSize: return None win = band.ReadAsArray(c0, r0, 2, 2).astype("float64") if nodata is not None and (win == nodata).any(): return None v00 = win[0, 0] v01 = win[0, 1] v10 = win[1, 0] v11 = win[1, 1] top = v00 * (1 - dc) + v01 * dc bottom = v10 * (1 - dc) + v11 * dc return float(top * (1 - dr) + bottom * dr)

这里有几点必须说清楚。ReadAsArray(c0, r0, 2, 2)一次读四个像素，比读四次单像素要快得多，因为 GDAL 的块读取有内部缓存。astype("float64")是必须的，如果 DEM 是 Int16，不转浮点的话加权计算会被截断成整数，精度损失直接体现为阶梯状的高程值。

//运算符在负数上的行为是向负无穷取整，这一点比int()更符合我们的需要——当目标点落在第一个像素的左半边时，col_c会是 -0.3 这种值，//得到 -1，然后越界检查会正确拦下；而int(-0.3)得到 0，会被误判为合法。

NoData 的判断我用了"四邻域中任意一个是 NoData 就整体返回 None"的保守策略。也有做法是排除 NoData 像素、只用有效像素重新归一化权重，这在图幅边缘能多救回一些点，但代价是插值结果的可靠性下降。我的建议是：除非业务明确要求边缘尽量出值，否则直接返回 None，让上层决定怎么补，比在下层悄悄造一个可能失真的值要好。

还有一个容易忽略的问题：gt[1]和gt[5]如果为 0（损坏的 GeoTransform），上面所有除法都会抛异常或者返回 inf。生产代码里加一句if gt[1] == 0 or gt[5] == 0: raise ValueError(...)，能让你在出问题的时候少走弯路。

3.3 用 rasterio 对照验证，快速排错

GDAL 的 Python 绑定写起来偏底层，一旦结果不对，很难判断是逻辑错了还是 API 用错了。这时候我会用 rasterio 写一个五行的对照版本，两者结果对得上，说明逻辑没问题，问题在输入数据；对不上，说明我的实现有 bug。

import rasterio from rasterio.warp import transform as rio_transform def get_elevation_rasterio(dem_path, lon, lat): with rasterio.open(dem_path) as src: xs, ys = rio_transform({"init": "epsg:4326"}, src.crs, [lon], [lat]) vals = list(src.sample([(xs[0], ys[0])], indexes=1)) return vals[0][0]

src.sample()内部用的是最近邻采样，所以它只适合和我的最近邻版本对照。如果要比双线性，就得手动算行列号再用src.read(1, window=...)。

这个对照方法在我排查"高程差 40 米"那次问题的时候立了功：rasterio 和我的 GDAL 版本结果一致，说明代码没问题，最后定位到是输入的经纬度顺序颠倒了。如果只用一个库，我大概率会在代码里反复折腾半天。

4. 批量查询的性能工程：从逐点循环到窗口分块

单点查询怎么写都行，真正拉开差距的是批量场景。一万个点的轨迹要提取高程剖面，或者一整个区域的网格点要做地形分析，逐点循环加逐次打开文件，跑起来可能就是几十分钟和几秒钟的区别。

4.1 数据集句柄复用与内存窗口读取

最直接的优化是把gdal.Open提到循环外面，只打开一次，所有查询共用同一个 dataset 和 band 对象。这一条能砍掉大部分时间。

def batch_query(dem_path, points, src_epsg=4326): ds = gdal.Open(dem_path, gdal.GA_ReadOnly) band = ds.GetRasterBand(1) gt = ds.GetGeoTransform() nodata = band.GetNoDataValue() ct = build_transformer(ds, src_epsg) results = [] for lon, lat in points: results.append(get_elevation_bilinear(ds, band, lon, lat, ct, gt, nodata)) return results

第二层优化是减少ReadAsArray的调用次数。每次调用都有固定开销，哪怕只读一个像素。如果查询点在图幅上的分布比较集中，可以算出这批点的行列号包围盒，一次性把整个窗口读进内存，之后所有查询都在内存数组里做索引。这个思路对于"沿一条路线提取高程剖面"特别有效，因为路线的包围盒通常是个细长条，数据量不大。

def query_in_window(ds, band, gt, ct, points, pad=2): cols, rows = [], [] for lon, lat in points: X, Y, _ = ct.TransformPoint(lon, lat) cols.append((X - gt[0]) / gt[1]) rows.append((Y - gt[3]) / gt[5]) c_min = max(int(min(cols)) - pad, 0) r_min = max(int(min(rows)) - pad, 0) c_max = min(int(max(cols)) + pad, ds.RasterXSize) r_max = min(int(max(rows)) + pad, ds.RasterYSize) win = band.ReadAsArray(c_min, r_min, c_max - c_min, r_max - r_min) return win, c_min, r_min

注意pad参数。双线性插值需要当前像素和右下邻域，如果窗口刚好卡在边界上，插值就会缺角。留 2 个像素的余量是最低要求，留多了浪费内存，留少了边界点取不到值。窗口读取之后，坐标反算要减去c_min和r_min的偏移量，这一步经常被漏掉，表现就是所有结果都整体平移了。

内存占用也要心里有数。一景 10000×10000 的 Float32 DEM 大概占 400MB，如果窗口覆盖了整幅图，那就是直接把全图读进来了。所以在做窗口读取之前，先算一下包围盒面积，超过阈值就分块处理。

4.2 多幅 DEM 拼接成 VRT 后的统一查询

单幅 DEM 覆盖不了业务范围是常态，尤其是跨省、跨流域的项目。手写多幅拼接逻辑很痛苦：判断点落在哪一幅、坐标系是否一致、接边处是否有重叠。GDAL 的 VRT 就是为这个场景准备的。

gdalbuildvrt mosaic.vrt tiles/*.tif

一行命令，把目录下所有分幅 DEM 拼成一个虚拟数据集。VRT 本身不复制数据，只是一份索引，查询开销几乎可以忽略。打开 VRT 的方式和普通 TIFF 完全一样，ReadAsArray会自动代理到对应的底层文件。

这里有两个必须提前确认的前提。第一，所有分幅的坐标系必须一致，如果不一致，gdalbuildvrt默认会报错或者给出警告，你需要先用gdalwarp统一到同一个坐标系再拼。第二，NoData 值要显式指定，否则接边处的空缺会被当成 0 参与计算：

gdalbuildvrt -srcnodata -32767 -vrtnodata -32767 mosaic.vrt tiles/*.tif

我踩过一个坑：不同批次的 DEM 产品 NoData 值不一样，有的是 -9999，有的是 -32767，拼接的时候没统一，结果在接缝附近出现了一条明显的高程异常带。排查了很久才想到去查每幅图的 NoData 字段。现在的习惯是拼接前先跑一遍脚本，把所有分幅的 NoData 和数据类型打印出来核对。

还有一点，VRT 的GetGeoTransform返回的是整个虚拟图幅的变换参数，行列号反算按它来就行，不需要关心具体落在哪一幅。但如果你的查询点非常集中在某一幅的范围内，直接打开那一幅会比走 VRT 略快，因为少了一层代理。这个差异在少量查询时可以忽略，千万级查询时值得考虑做一次预分区。

4.3 向量化查询的落地方案

再往上走一步，就是把整个批量查询彻底向量化，不写 Python 层的循环。GDAL 本身不提供现成的向量化接口，但可以借用 NumPy 把计算密集的部分搬到数组层面。

思路是这样：先用ct.TransformPoints()（注意是复数形式）一次性把所有点转换到 DEM 坐标系，得到一个 N×3 的数组；然后用 NumPy 的广播机制一次性算出所有浮点行列号；接着筛掉越界点；最后对合法点做批量取值。

import numpy as np def batch_vectorized(ds, band, ct, gt, lon_arr, lat_arr, nodata=None): pts = ct.TransformPoints(list(zip(lon_arr, lat_arr))) arr = np.asarray(pts, dtype="float64") X, Y = arr[:, 0], arr[:, 1] col_f = (X - gt[0]) / gt[1] - 0.5 row_f = (Y - gt[3]) / gt[5] - 0.5 c0 = np.floor(col_f).astype(int) r0 = np.floor(row_f).astype(int) valid = (c0 >= 0) & (r0 >= 0) & \ (c0 + 1 < ds.RasterXSize) & (r0 + 1 < ds.RasterYSize) result = np.full(len(lon_arr), np.nan) if not valid.any(): return result # 一次性读全图（小图适用），或按包围盒分块读 data = band.ReadAsArray().astype("float64") dc = col_f[valid] - c0[valid] dr = row_f[valid] - r0[valid] v00 = data[r0[valid], c0[valid]] v01 = data[r0[valid], c0[valid] + 1] v10 = data[r0[valid] + 1, c0[valid]] v11 = data[r0[valid] + 1, c0[valid] + 1] vals = (v00 * (1 - dc) + v01 * dc) * (1 - dr) + \ (v10 * (1 - dc) + v11 * dc) * dr if nodata is not None: vals = np.where(vals == nodata, np.nan, vals) result[valid] = vals return result

TransformPoints比循环调TransformPoint快一个数量级，因为它把投影转换的初始化开销摊薄到了整批点。ReadAsArray()全图读取在小图上完全没问题，大图就得改成按包围盒分块，把上面的逻辑包在一个循环里，逐块处理。

这里有个隐藏的精度问题：np.floor对于恰好落在整数上的浮点值，行为和//一致；但由于浮点误差，一个本该是 3.0 的值可能算出来是 2.9999999，floor之后变成 2。这种误差在像素边界上会造成一个像素的跳变。缓解办法是加一个极小的容差，比如np.floor(col_f + 1e-9)，代价几乎为零。我在做长路线剖面时遇到过这个问题，表现是剖面曲线偶尔出现单个点的尖刺，一开始以为是 DEM 本身的噪声，查到最后就是浮点边界问题。

5. 踩坑记录：那些查出来的高程是错的

前面讲的是"怎么写对"，这一章讲"怎么发现写错了"。因为用地理坐标查 DEM 高程这件事，最麻烦的地方在于错误的结果往往看起来是对的——没有异常抛出，返回值在合理范围内，曲线也画得出来，只是数值不对。

5.1 NoData 伪装成 0 米的值域陷阱

这是杀伤力最大的一类问题。有些 DEM 的 NoData 值是 0，有些是 -9999，有些是 -32767，还有一些根本没有声明 NoData，边界用 0 填充。如果你不判 NoData，这些占位值会和真实高程混在一起。

-32767 这种极端值比较容易发现，画个直方图一眼就看到了。麻烦的是 0。沿海地区的真实高程可能是 0 到 5 米，内陆盆地的湖面高程也可能是几百米，你很难从数值上判断一个 0 到底是"海平面附近的真实值"还是"无数据占位"。

我的排查顺序是这样的：先ComputeRasterMinMax()看全局极值，如果最小值是一个非常整齐的负数（-9999、-32767、-3.4e38），基本可以确认是 NoData；再看GetNoDataValue()有没有返回值；两者都对不上，就把图幅四周一圈的像素值打出来看看，因为 NoData 通常集中在边界。

还有一种情况是 NoData 通过掩膜波段表达，GetNoDataValue()返回 None，但实际数据里有一块是无效的。这时候要用GetMaskBand()：

mask = band.GetMaskBand() mask_arr = mask.ReadAsArray(c0, r0, 2, 2) if (mask_arr == 0).any(): return None

掩膜为 0 表示该像素无效。这个判断在数据源不统一的项目里非常必要，因为不同厂商的产品对无效区域的处理方式完全不同。

5.2 经纬度顺序颠倒与"整条路线漂移"

这个坑我在不同项目里见过至少三次，每次的表现都很像：查询结果不是全空，也不是明显异常，而是所有点的高程都偏低或者偏高，剖面曲线的形状大致还在，但整体平移了。

原因就是TransformPoint(lon, lat)的参数顺序。GDAL 3.x 默认按 EPSG:4326 的官方轴序（纬度在前）解释，如果你传的是(lat, lon)，转换结果的 X 和 Y 就互换了。在中国境内，经度大约 73 到 135，纬度大约 3 到 54，互换之后 X 变成了 3 到 54，Y 变成了 73 到 135，这个坐标在地球上对应的是完全不同的位置——可能落在非洲或者北冰洋。但由于后续的行列号反算依然能算出一个数（只是越界或者落在图幅的某个角落），如果你恰好没做严格的越界检查，就会查到一堆"看起来合理"的高程值。

防线有三道。第一道，SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)必须加，而且要加在源和目标两个 SpatialReference 上。第二道，越界检查一定要严格，col和row都必须落在[0, RasterXSize)和[0, RasterYSize)内。第三道，转换后的坐标做一次量纲校验：如果 DEM 是投影坐标系，转换结果的绝对值应该在几万到几百万之间（米），如果你的输入是中国境内的经纬度，转出来的 X 却只有几十，那一定是轴序错了。

我还建议在开发阶段加一句断言，检查转换前后的坐标是否落在同一个大致范围内：

X, Y, _ = ct.TransformPoint(lon, lat) assert -180 <= X <= 180 or abs(X) > 1000, "坐标量纲异常，检查轴序"

这种土办法比事后调试有效得多。

5.3 越界坐标的静默返回

ReadAsArray有一个非常坑的行为：当请求的窗口超出栅格范围时，它不会抛异常，而是返回一个只包含有效部分的较小数组，或者在某些参数组合下直接返回 None。如果你拿这个返回值直接取[0, 0]，得到的是窗口左上角那个像素的值，而不是你请求的那个位置的值。

value = band.ReadAsArray(col, row, 1, 1)[0, 0]

这一行在col或row越界时的行为是不确定的，可能报 IndexError，可能返回别的值，取决于 GDAL 版本和数据源类型。稳妥的做法是在读之前就做越界判断，而不是依赖读操作的报错。我的习惯是封装一个统一的取值函数，越界一律返回 None，让上层的业务逻辑去决定怎么处理。

对于批量查询，越界的比例也是个体检指标。正常情况下，如果你的点位都在 DEM 覆盖范围内，越界率应该是 0；如果超过了 5%，要么是坐标系统一没做对，要么是数据覆盖范围有大片空缺，这时候应该先停下来查数据，而不是继续跑。

5.4 接边处的半像素错位

多幅 DEM 拼接之后，接边处的高程经常出现台阶状跳变。很多人第一反应是数据质量问题，实际上更常见的原因是相邻两幅的 GeoTransform 没有对齐到同一个像素网格。

假设第一幅的原点 X 是 500000，像素宽 30 米；第二幅的原点 X 是 500015，像素宽也是 30 米。两幅拼在一起，栅格采样点在物理空间上是错开半个像素的。VRT 会按第一幅的网格重采样第二幅，这个重采样过程会引入误差，在陡坡处表现为接缝两侧高程不一致。

排查方法是把接缝附近的坐标打出来，分别用两幅图各自查询，对比结果。如果差异在合理插值误差范围内（比如 1 到 2 米），那是正常的；如果差出十几米，基本可以确认是网格错位。解决办法是在拼接前用gdalwarp把所有分幅重采样到统一的网格（-te指定范围，-tr指定分辨率），代价是一次性的预处理时间，但换来的是接边处的干净。

6. 精度取舍与工程封装建议

写到最后，我想聊两个更偏工程的问题：插值方式到底带来多大误差，以及这段逻辑该怎么封装才能在项目里长期用下去。

6.1 分辨率和插值方式带来的误差量级

很多人关心"用 GDAL 查出来的高程误差有多大"，这个问题没法脱离数据本身回答。误差来源主要是三个：DEM 本身的垂直精度、水平分辨率导致的代表性误差、以及插值方法引入的误差。

以常见的 30 米分辨率 DEM 为例，标称垂直精度通常在 7 到 16 米之间（具体看产品和地形）。这个量级是底噪，任何查询方法都突破不了。在此之上，水平分辨率带来的代表性误差更值得关注：一个 30 米像素覆盖的实际地表是一个斜坡，像素存储的值是该范围的一个统计代表值。你查询的那个精确位置的真实高程，和像素值之间本来就有一个差值。在坡度 20 度的山区，30 米像素内的自然高差可以到 10 米以上。

插值方法的差异体现在这个基础上。最近邻直接取像素值，误差就是上面说的那个代表性误差。双线性插值用四个像素加权，相当于对地表做了一个局部线性拟合，在平滑地形上能把误差压低一些，但在山脊、陡崖这类地形突变处会低估峰值。我做过一次对比：在同一景 30 米 DEM 上取一百个采样点，双线性相对最近邻的差值，平坦区域基本在 1 米以内，缓坡区在 1 到 3 米，陡坡区能到 5 米以上。

如果你的业务对高程精度有硬性要求（比如淹没分析、坡度分级），那不能只靠插值优化，得从数据源下手：换更高分辨率的 DEM，或者用实测点做局部校正。GDAL 只是取值的工具，它不会凭空提升数据的精度。

6.2 封装成工具函数时该定哪些接口

这段逻辑在一个项目里往往会被调用几十次，封装质量直接决定了后期维护成本。我总结了几个自己踩过坑之后定下来的约定。

第一，数据源只开一次。用一个类或者一个带缓存的工厂函数管理数据集句柄，不要在每次查询时gdal.Open。缓存 key 用文件路径加修改时间，避免文件被替换后还在用旧句柄。

第二，越界和 NoData 明确返回 None，不返回 0 或 -9999。上层业务用None判断"这个点没查到"是最清晰的，返回哨兵值容易在后续计算里被当成真实数据参与运算。

第三，输入坐标系作为显式参数。不要默认所有输入都是 WGS84，接口上写清楚src_epsg=4326，调用方需要传别的就直接改。遇到过有项目内部用的是本地平面坐标，因为接口写死了 4326，导致整批数据全错。

第四，插值方式做成可选参数。默认双线性，需要极致速度时切最近邻。两种实现共用坐标转换和越界检查的逻辑，只在最后取值那一步分叉。

第五，加一个自检方法。传入一组已知点位的坐标和预期高程，跑一遍验证。这个自检在换 DEM 数据源、升级 GDAL 版本之后特别有用，能在第一时间发现回归问题。

一个简单的接口签名大概是这样：

class DemElevationReader: def __init__(self, dem_path, src_epsg=4326, method="bilinear"): ... def query(self, lon, lat): ... def query_batch(self, points): ... def self_check(self, samples): ...

写到这里，我想起一件事。最早做这类查询的时候，我总觉得 GDAL 的 API 又长又啰嗦，总想找个更简单的库替代。后来发现，那些"更简单"的封装在遇到投影不一致、NoData、接边错位这些问题时，反而因为隐藏了细节而更难排查。GDAL 的接口确实不友好，但它把每一层都摊开给你看，出了问题你知道该往哪儿找。这大概也是它在这个领域活了这么多年的原因。

如果你的 DEM 是分幅的，记得先拼 VRT 再查；如果查询量在十万级以上，记得把坐标转换批量做、把窗口读取用起来；如果结果看着正常但心里没底，找一个已知高程的地标点验证一下，比如机场跑道、水库水面，这些地方的高程是公开且稳定的，拿来当标定基准比自己猜靠谱得多。