关于使用Rasterio的index()方法获取栅格像素坐标及仿射变换的技术问询

1. index()方法是否为最优方案？ #

是的，直接使用data.index(longitude, latitude)是推荐的最优方案之一。rasterio的index()方法已经封装了地理坐标到像素行列的转换逻辑，内部会自动调用栅格的仿射变换矩阵完成计算，无需手动编写转换代码，既简洁又能避免手动计算可能出现的错误。

但要注意一个关键前提：输入的地理坐标必须和栅格的CRS（坐标参考系）一致。比如你的代码里打印了data.crs，如果输出是UTM类投影坐标，而你输入的是WGS84经纬度（13.3886, 52.5174），直接调用index()会得到错误结果。这种情况下需要先把经纬度转换为栅格对应的CRS，示例代码如下：

import rasterio
from rasterio.warp import transform

with rasterio.open("LC08_L2SP_190037_20190619_20200827_02_T1_ST_B10.TIF") as data:
    print(data.crs)
    # 输入坐标的CRS是WGS84 (EPSG:4326)
    src_crs = "EPSG:4326"
    dst_crs = data.crs
    # 转换坐标到栅格的CRS
    x, y = transform(src_crs, dst_crs, [13.3886], [52.5174])
    row, col = data.index(x[0], y[0])
    print("Pixel coordinates of point ({}, {}): ({}, {})".format(13.3886, 52.5174, col, row))

2. 仿射变换的使用场景与方法 #

什么是仿射变换？ #

栅格的仿射变换矩阵（data.transform）是一组6个参数的数组，定义了地理坐标与像素行列坐标之间的线性转换关系，核心公式为：

x = transform[0] + col * transform[1] + row * transform[2]
y = transform[3] + col * transform[4] + row * transform[5]

其中(x,y)是地理坐标，(row,col)是像素行列索引。

使用场景 #

• 批量转换多个点：如果需要转换大量地理坐标到像素坐标，循环调用index()效率偏低，可手动用rasterio.transform.rowcol()方法，传入变换矩阵和批量坐标一次性完成转换：

  from rasterio.transform import rowcol
  
  coords = [(13.3886, 52.5174), (13.4, 52.52)]
  # 先确保坐标和栅格CRS一致（此处假设已完成转换）
  xs = [x for x, y in coords]
  ys = [y for x, y in coords]
  rows, cols = rowcol(data.transform, xs, ys)
  

• 反向转换（像素转地理坐标）：当需要从像素行列索引获取对应的地理坐标时，直接调用data.transform * (col, row)即可得到结果：

  col, row = 100, 200
  x, y = data.transform * (col, row)
  

• 自定义转换逻辑：如果需要调整转换规则（比如修改四舍五入方式、增加坐标偏移等），可以基于变换矩阵手动编写转换代码。

3. 获取像素坐标前是否需要执行仿射变换？ #

不需要手动执行。data.index()方法内部已经自动调用栅格的仿射变换矩阵，完成了地理坐标到像素行列的转换计算。你只需要确保输入坐标的CRS与栅格CRS一致即可；如果不一致，需要先做坐标转换，而非手动执行仿射变换。

内容的提问来源于stack exchange，提问作者Salah Amani