已知影像四角坐标，如何用Python获取像素对应的经纬度？

问题概述 #

拥有Iceye GRD卫星SAR影像，已知四个图像角点的1-based像素坐标与WGS84地理坐标（纬度在前，经度在后）的对应关系：

image pixel    -> georef. location
[1, 1]         -> [1.6033675891108876, 104.05738201292307]
[16892, 1]     -> [1.5291001811499518, 104.42831446776263]
[1, 26855]     -> [1.0079111030191064, 103.93758132038617]
[16892, 26855] -> [0.9334645684110325, 104.30991666976476]

影像采用WGS84地理参考系统，但无有效内置地理元数据，使用rasterio读取时返回全1变换矩阵，无法完成坐标转换，需手动实现像素坐标到地理坐标的转换函数。

最优实现方式 #

由于SAR影像存在透视畸变，**单应性变换（Homography Transformation）**是最优选择——通过4个已知对应点计算单应矩阵，可准确描述像素坐标到地理坐标的透视映射关系；若影像近似线性变形，也可使用最小二乘法拟合仿射变换矩阵，计算更快但精度略低。

轻量Python模块推荐 #

仅需numpy即可实现核心计算，无需重型GIS库；若追求更简洁的API，可搭配scipy的transform模块。

代码实现 #

方法1：单应性变换（推荐，适配SAR透视畸变） #

import numpy as np

def pixel_to_geo(pixel_coord, pixel_points, geo_points):
    # 转换为numpy数组，保持1-based像素坐标规则
    pixel_points = np.array(pixel_points, dtype=np.float32)
    geo_points = np.array(geo_points, dtype=np.float32)
    
    # 构建单应性矩阵求解方程组
    A = []
    for (x, y), (lat, lon) in zip(pixel_points, geo_points):
        A.append([x, y, 1, 0, 0, 0, -x*lat, -y*lat])
        A.append([0, 0, 0, x, y, 1, -x*lon, -y*lon])
    A = np.array(A)
    B = geo_points.flatten()
    
    # 求解单应矩阵H
    H, _, _, _ = np.linalg.lstsq(A, B, rcond=None)
    H = np.append(H, 1).reshape(3, 3)
    
    # 转换输入像素坐标为齐次坐标并计算地理坐标
    x, y = pixel_coord
    pixel_hom = np.array([x, y, 1], dtype=np.float32).reshape(3, 1)
    geo_hom = H @ pixel_hom
    lat = geo_hom[0][0] / geo_hom[2][0]
    lon = geo_hom[1][0] / geo_hom[2][0]
    
    return (lat, lon)

# 初始化已知角点数据
pixel_corners = [
    [1, 1],
    [16892, 1],
    [1, 26855],
    [16892, 26855]
]
geo_corners = [
    [1.6033675891108876, 104.05738201292307],
    [1.5291001811499518, 104.42831446776263],
    [1.0079111030191064, 103.93758132038617],
    [0.9334645684110325, 104.30991666976476]
]

# 测试输入：(10000, 20000)
test_pixel = (10000, 20000)
result = pixel_to_geo(test_pixel, pixel_corners, geo_corners)
print(f"像素坐标{test_pixel}对应的地理坐标：{result}")

方法2：仿射变换（适合线性变形场景） #

import numpy as np

def pixel_to_geo_affine(pixel_coord, pixel_points, geo_points):
    # 转换为numpy数组，统一为0-based像素坐标
    pixel_points = np.array(pixel_points, dtype=np.float32) - 1
    geo_points = np.array(geo_points, dtype=np.float32)
    
    # 构建矩阵并拟合仿射变换参数
    A = np.concatenate([pixel_points, np.ones((4, 1))], axis=1)
    lat_params, _, _, _ = np.linalg.lstsq(A, geo_points[:, 0], rcond=None)
    lon_params, _, _, _ = np.linalg.lstsq(A, geo_points[:, 1], rcond=None)
    
    # 计算目标地理坐标
    x, y = pixel_coord
    x -= 1
    y -= 1
    lat = lat_params[0] * x + lat_params[1] * y + lat_params[2]
    lon = lon_params[0] * x + lon_params[1] * y + lon_params[2]
    
    return (lat, lon)

# 测试调用
test_pixel = (10000, 20000)
result_affine = pixel_to_geo_affine(test_pixel, pixel_corners, geo_corners)
print(f"仿射变换结果：{result_affine}")

注意事项 #

• 像素坐标索引：需保持与已知角点的索引规则一致（1-based或0-based），代码中已做对应适配。
• 地理坐标顺序：已知点为纬度在前，经度在后，若需输出(lon, lat)格式可自行调整返回值。
• 精度验证：可将已知角点代入函数验证结果，若误差过大需检查坐标索引或变换类型是否匹配。

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