rasterio仿射变换深度解析：像素坐标与地理坐标互转及窗口裁剪避坑

先前我用仿射变换处理一批无人机影像的时候，犯过一个特别蠢的错误：把行列号当成经纬度的输入顺序，程序跑得顺顺利利，结果一核对地面实测点，整体偏出去好几百米，而且不报任何错。后来我才意识到，rasterio 的 Transform 根本不是"六个数字"那么简单，它背后藏着一整套像素坐标与地理坐标的约定。

这篇文章想把这些约定一次说透：什么是仿射变换、六个参数各自代表什么、行列号与坐标互转时到底该谁在前，以及窗口裁剪后为什么必须换变换矩阵。适合刚开始用 rasterio 处理遥感影像、或者被坐标偏移问题折磨过的开发者参考。看懂之后，你在栅格数据定位这件事上基本不会再生坑。

1. 仿射变换的 6 个数字到底在描述什么

1.1 一张 GeoTIFF 里的 6 个浮点数

栅格影像本身只是一个二维数组，数组里的每个元素只有"第几行第几列"这种像素坐标，并没有真实的地理坐标。GeoTIFF 之所以能放进 GIS 软件里跟矢量数据完美叠合，靠的就是文件头里记录的一组仿射变换参数。

rasterio 中这组参数被封装成affine.Affine对象，也就是src.transform。它本质上对应这样一个线性映射：

x = a * col + b * row + c y = d * col + e * row + f

这里col是像素列号，row是像素行号，(x, y)是地理坐标。六个参数分别是：

• a：列方向的分辨率，也就是一个像素在东西方向上占多少距离，通常为正。
• e：行方向的分辨率，一个像素在南北方向上占多少距离，标准北朝上影像中通常为负。
• b、d：旋转和倾斜项。数据未校正、或经过多项式纠正时会出现非零值。
• c、f：影像左上角的地理坐标，对应左上角像素的左上角。

很多人第一次接触时会想：这跟平时用的经纬度有什么关系？其实不管数据是经纬度还是投影坐标（比如 UTM），都只是把上面的(x, y)换一种单位而已。矩阵的映射逻辑完全一样。

1.2 为什么 y 方向分辨率一定是负的

这是新手最容易困惑的点：为什么e参数通常是负数？

原因在数组坐标和地理坐标的轴向定义上。影像数组的第 0 行在顶部，行号向下增长；而地理坐标的 Y 轴（北方向）是向上增长的。为了让"行号 +1 时，地理位置向下移动"，而向下对应地理坐标 Y 值减小，所以行方向的分辨率必须带负号。

举个例子，10 米分辨率的正射影像，a通常是10.0，e则是-10.0。如果你手工构造变换矩阵时把e写成了正数，结果就是影像在南北方向发生镜像翻转，但程序依然不会报错。这种问题非常隐蔽，后面我会专门讲怎么自检。

1.3 Affine 对象和 GDAL GeoTransform 的对应关系

做过 GDAL 开发的人应该熟悉GetGeoTransform()返回的六元组。rasterio 也提供了src.geotransform属性，返回同样的六元组。两者顺序不一样，对照关系如下：

所以把 GDAL 六元组转成 Affine 时，不是简单平移下标，而是交叉对应的。rasterio 内部已经做好了封装，日常开发我建议直接用src.transform这个 Affine 对象，别自己去拆六元组，减少一层出错概率。

2. Transform 从哪来：读取、构造函数与 GCP 数据

2.1 直接读取：src.transform与src.geotransform

处理标准 GeoTIFF 时，最省心的办法就是打开文件后直接读取。

import rasterio with rasterio.open("示例影像.tif") as src: print(src.crs) # 坐标系统 print(src.transform) # Affine 对象 print(src.geotransform) # GDAL 风格六元组

正常生产的数据都会写入这些信息。拿到src.transform之后，可以配合src.width、src.height和src.bounds快速了解影像覆盖范围。这里我养成了一个习惯：任何一批新数据到手，第一件事就是打印crs、transform、bounds三个值。坐标系是否是我预期的、像元大小是否合理、范围是否落在目标区域，一眼就能看出问题。

2.2 从角点和边界手工构造：from_origin与from_bounds

如果你的数据不是从文件读出来的，而是内存里一个 numpy 数组需要写成带坐标的 GeoTIFF，那就要手工构造仿射变换。

已知左上角坐标和分辨率时，用from_origin：

from rasterio.transform import from_origin transform = from_origin( west=500000, # 左上角 X north=3000000, # 左上角 Y xsize=10, # 列方向分辨率 ysize=10 # 行方向分辨率 )

注意ysize传正值即可，函数内部会把它转成负的e参数，你不需要自己加负号。

已知整个范围框和目标网格大小时，用from_bounds：

from rasterio.transform import from_bounds transform = from_bounds( west=499900, south=2999900, east=501000, north=3001000, width=110, height=110 )

from_bounds会按照范围框的宽度和高度，自动算出列方向和行方向的分辨率。它不保证像元一定是正方形，如果你的范围框宽高比和网格行列数不匹配，得到的像元就是长方形。这种情况只要在合理范围内，一般不影响使用。

2.3 没有仿射变换的数据：GCP 定位的识别

还有一种情况挺常见：文件本身没有可靠的 GeoTransform，只带了地面控制点。比如一些老的扫描地形图、未经严格几何校正的航片，会以 GCP 的形式保存少量控制信息。

rasterio 中可以通过src.gcps获取：

gcps, gcp_crs = src.gcps if len(gcps) > 0: print("这份数据主要靠地面控制点定位")

遇到这种情况，别急着用src.transform去转坐标，因为此时它很可能已经退化成一个没有实际意义的默认值。你要先检查 GCP 的分布、数量、残差，判断这份数据的定位精度到底够不够用。如果必须把 GCP 拟合成一个仿射变换，通常的做法是用 GDAL 的 GCP 多项式拟合能力，但这属于近似计算，拟合残差一定要单独评估，不能当准确值用。

3. 行列号与地理坐标互转：官方约定和容易搞混的参数顺序

3.1 正确的乘法姿势：列号在前、行号在后

rasterio 的transform是仿射矩阵，直接用乘法就可以把像素坐标变成地理坐标：

col = 50 row = 30 x, y = transform * (col, row)

这里第一个位置放的是列号col，第二个才是行号row。原因在于仿射矩阵的数学约定：第一输入轴对应东西方向（列方向），第二输入轴对应南北方向（行方向）。

很多人看到经纬度习惯先写纬度、再写经度，于是顺手写成transform * (row, col)，这在标准北朝上影像里就会造成坐标分量对调。假设影像左上角是(1000, 2000)，分辨率 10 米，正确应该输入(col=50, row=30)，得到：

x = 10 * 50 + 0 * 30 + 1000 = 1500 y = 0 * 50 + (-10) * 30 + 2000 = 1700

如果误写成(30, 50)，得到的是(1300, 1500)，数值差得不多，但对于 10 米分辨率来说，等于错了 28 个像元。旋转参数非零的数据，这种情况会更离谱。

3.2 像素中心 vs 像素角点：0.5 个像素的偏差

如果你直接用transform * (col, row)计算，得到的是像素左上角的地理坐标，而不是像素中心。因为 rasterio / GDAL 的像素坐标约定里，像元覆盖的是[col, col+1] × [row, row+1]这个范围，整数(col, row)对应的是格网交点。

想取像素中心，需要自己加 0.5：

x_center = transform * (col + 0.5, row + 0.5)

当然，日常开发我更推荐用src.xy这个封装方法，它默认返回的就是像元中心坐标：

row = 30 col = 50 x_center, y_center = src.xy(row, col)

注意src.xy的参数顺序又是行号在前、列号在后，跟transform * (col, row)正好相反。这个别扭很容易把人绕晕，我的建议是：代码里注释写清楚每一步参数是什么，别省这个注释。

低分辨率数据上，中心点和角点之间的偏差会变得非常大。比如 1 公里分辨率的全球产品，0.5 个像素就是 500 米。这种数据你要是用角点坐标去跟地面站点观测值做匹配，误差直接就是致命的。

3.3 逆变换与取整：用src.index还是~transform

反过来，已知地理坐标想求它落在哪个像元，可以用逆仿射变换。Affine 对象支持用~运算符取逆矩阵：

col_float, row_float = ~src.transform * (x, y)

这里返回的是浮点行列号，需要自己取整得到像元索引。但注意两点：

• ~transform * (x, y)返回的顺序是(col, row)，和正变换保持一致。
• 取整方式会影响最终结果，通常用四舍五入round()。

实际上 rasterio 已经封装了一个更符合直觉的 API，直接传地理坐标，返回整数行列号：

row_idx, col_idx = src.index(x, y)

src.index返回的顺序是(row, col)，内部完成逆变换和四舍五入。而且它会对超出影像范围的点返回一个相对位置，你在取像元值之前最好检查行列号是否落在[0, height-1]和[0, width-1]范围内。

这里有个很容易混乱的地方需要反复强调：src.index(x, y)返回(row, col)，而transform * (col, row)的输入是(col, row)。两者方向不同但都遵循"第一参数是 X 方向、第二参数是 Y 方向"的底层逻辑。我见过不少人把src.index的解包写成col_idx, row_idx = src.index(x, y)，然后数组索引就全乱了。

4. 窗口裁剪之后的 Transform：一个非常隐蔽的偏移坑

4.1 为什么窗口裁剪后原来的 transform 不能继续用

处理大影像时，经常只读取一个窗口区域。假设影像起点地理坐标是(500000, 3000000)，分辨率 10 米，我用Window(col_off=100, row_off=200, width=256, height=256)切了一块。

窗口左上角对应的地理坐标是：

x = 500000 + 100 * 10 = 501000 y = 3000000 - 200 * 10 = 2998000

如果裁剪之后，子数组的第一个元素仍然沿用整幅影像的src.transform，那子图的(0, 0)会被错误地当成(500000, 3000000)。东西偏 1000 米，南北偏 2000 米，这个偏差对于后续叠加矢量、提取像元值来说，完全不可接受。

4.2window_transform的应用示例

正确的做法是使用src.window_transform(window)，它会基于窗口偏移量计算出一个只适用于该子图的局部变换矩阵：

from rasterio.windows import Window window = Window(col_off=100, row_off=200, width=256, height=256) with rasterio.open("示例影像.tif") as src: local_transform = src.window_transform(window) sub = src.read(1, window=window) print(local_transform)

local_transform的a、e分辨率参数和全局保持一致，但c、f已经置换成了窗口左上角的地理坐标。这样sub[0, 0]就能正确对应到(501000, 2998000)这个位置。

4.3 裁剪成 GeoTIFF 时如何保护正确的定位信息

把裁剪结果保存成新文件时，同样要用local_transform覆盖 profile，而不是用原来的src.transform：

profile = src.profile.copy() profile.update(width=256, height=256, transform=local_transform) with rasterio.open("裁剪子图.tif", "w", **profile) as dst: dst.write(sub, 1)

很多新手都会在这里踩坑：裁剪出来的影像单独打开看不出任何问题，但放进 GIS 软件里跟矢量叠加，整幅图都会偏移。而且因为裁剪区域往往是整幅影像的一小部分，这种偏移不会让影像"消失"，只会让地物对不上，排查起来特别费劲。

如果你需要裁剪之后的影像范围和原来的某个矢量完全对齐，我建议保存前再验证一次：用local_transform去算子图四个角的坐标，跟矢量边界对比一下，没问题再落盘。

5. 常见翻车现场与 5 分钟自检方法

5.1 参数顺序颠倒但"看起来正常"

这是最典型的错误：程序不报错，坐标输出范围也大体合理，但跟地图一比就错位。原因就是这章前面反复提到的参数顺序问题。

判断方法：取影像最中心的一个像元，用src.xy得到中心坐标，再用你手写的transform * (col, row)去算，如果发现两者结果相差恰好约 0.5 个像元，说明你的参数顺序没错，只是把角点当成了中心；如果相差好几个像元且方向不稳定，多半是row和col对调了。

5.2 负分辨率丢失导致的镜像问题

手工构造 Affine 时，e参数写成正数是最容易犯的低级错误。一旦e为正，整个影像的行方向坐标会反向增长，结果是南北方向发生镜像。在 GIS 里打开时，影像会和周围地物颠倒，像照镜子一样。

排查方式很简单：打印transform，检查e是不是负数。另外可以对比src.bounds的top和bottom，正常北朝上影像top > bottom，如果反了，基本就是e出了问题。

5.3 跨 CRS 直接套用矩阵的危险

仿射变换的六个参数是相对于某个特定坐标系统的。如果影像的 CRS 是 UTM，而你传入的x、y是经纬度，运算结果会离谱到天际线。API 不报错，是因为它根本不知道你给的数值代表什么含义。

所以任何坐标转换前，先确认影像的src.crs和你手上的坐标是否一致。如果不一致，用rasterio.warp.transform_bounds先把范围框转换过去，再配合src.index使用：

from rasterio.warp import transform_bounds # 把影像范围转换到 WGS84 经纬度 wgs84_bounds = transform_bounds(src.crs, "EPSG:4326", *src.bounds)

严格来说，这里是把整个边界框做了转换，得到的范围和逐点投影会有细微差异，但对大多数判断场景足够用。

5.4 快速自检矩阵是否正确的两个办法

第一个办法是中心校验法。打开影像后取正中间的像元，打印它的中心坐标，再和src.bounds对比：

with rasterio.open("示例影像.tif") as src: row, col = src.height // 2, src.width // 2 x_center, y_center = src.xy(row, col) print(x_center, y_center) print(src.bounds)

中心点的x应该落在bounds的左右区间之内，y也应该在上下区间之内，并且大致在中间。如果中心点跑到了范围外，矩阵八成有问题。

第二个办法是已知控制点法。如果你有一个独立来源的地面控制点，比如实测的农田角点、房角坐标，直接用src.index(x, y)转成行列号，再读取该位置的像元值，看是否符合预期地物。这个方法最直接，能同时验证 transform、CRS、以及影像本身的质量。

我在实际处理流程里定了一条铁律：处理流水线一旦涉及坐标换算，必须保留一份日志，记录每批数据的crs、transform、bounds和抽样验证结果。这样即使后面发现问题，也能快速定位是哪一步、哪一批数据出了偏差，而不是从头重新翻找原始文件。

最后再分享一条个人经验。rasterio Transform 本身不复杂，真正的复杂度来自历史包袱和隐含约定。我被坑得最惨的一次，就是高分辨率影像上混淆了角点和中心坐标，半像素偏差在那种数据上直接影响点位精度。从那以后，我的原则特别简单：能用src.xy和src.index解决的需求，绝不手写矩阵；必须手写时，先把六个参数对照 GDAL 六元组和 Affine 字段核对一遍。做到这三点，基本不会再在栅格定位这种事上浪费一整天。