如何使用Rasterio对卫星图像子集的经纬度坐标进行插值？

对于你提到的80-100公里范围的卫星图像，直接用线性插值计算子集的经纬度坐标完全可行，不需要额外考虑地球曲率的影响，原因如下：

• 这个尺度下，地球曲率带来的经纬度偏差极小，远小于常规卫星图像的地理定位误差，对结果的影响可以忽略。
• Rasterio本身的地理坐标映射逻辑就是基于线性插值实现的，完全适配这类小尺度场景。

具体实现方法 #

方法1：利用Rasterio自带的地理变换（推荐）

Rasterio加载图像时会自动读取其地理变换矩阵，直接用这个矩阵就能快速转换像素坐标与经纬度：

import rasterio

# 加载图像
with rasterio.open("your_sat_image.tif") as src:
    # 获取地理变换矩阵，结构为 (左上角经度, 经度方向分辨率, 0, 左上角纬度, 0, -纬度方向分辨率)
    transform = src.transform
    # 假设你截取的子集像素范围是：起始行row_start、结束行row_end，起始列col_start、结束列col_end
    subset_bbox = (
        transform * (col_start, row_start),    # 子集左上角经纬度
        transform * (col_end, row_end)        # 子集右下角经纬度
    )
    print(f"子集BoundingBox：{subset_bbox}")

如果需要批量计算子集中所有像素的坐标，可以用rasterio.transform.xy函数：

from rasterio.transform import xy

# 获取子集内所有像素中心的经纬度（offset可选'ul'表示左上角）
lons, lats = xy(
    transform,
    rows=range(row_start, row_end),
    cols=range(col_start, col_end),
    offset='center'
)

方法2：手动构建线性映射（如果你只有四个角的坐标）

如果没有图像的元数据，只有四个角的经纬度(X1,Y1)（左上）、(X2,Y2)（右上）、(X3,Y3)（左下）、(X4,Y4)（右下），可以手动计算分辨率后做线性插值：

# 假设图像宽W像素，高H像素
W = src.width
H = src.height

# 计算经度、纬度方向的分辨率
lon_res = (X2 - X1) / W
lat_res = (Y1 - Y3) / H

# 计算子集左上角像素(col_start, row_start)的经纬度
subset_lon1 = X1 + col_start * lon_res
subset_lat1 = Y1 - row_start * lat_res

# 计算子集右下角像素(col_end, row_end)的经纬度
subset_lon2 = X1 + col_end * lon_res
subset_lat2 = Y1 - row_end * lat_res

关键总结 #

80-100公里属于典型的小尺度地理范围，线性插值的精度完全满足需求，无需引入复杂的曲率校正模型。直接用Rasterio的内置工具或手动线性映射即可高效完成坐标转换。

内容的提问来源于stack exchange，提问作者Becky W