气象卫星局部PNG图像重投影与配准问题求助

我的问题和气象卫星图像重采样相关，但手里只有局部气象卫星PNG图像，而非完整图像，目前投影校正遇到了障碍。

已尝试的代码 #

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

def transform_extent_pts(extent_pts, map_proj, pt_crs):
    
    xul, yul = map_proj.transform_point(
        x = extent_pts[0],
        y = extent_pts[3],
        src_crs = pt_crs)
    
    xlr, ylr = map_proj.transform_point(
        x = extent_pts[1],
        y = extent_pts[2],
        src_crs = pt_crs)

    return [xul, xlr, ylr, yul]


sat_image1 = ROOT
df = plt.imread(sat_image1)

# 重投影到墨卡托投影
map_proj = ccrs.Mercator()

# 图像范围（静止卫星坐标系）:
data_crs = ccrs.Geostationary(central_longitude=0.0)
ax2 = plt.axes(projection=data_crs)
img_extent_sat = ax2.get_extent(crs=data_crs)    
img_extent_sat = [1.03*x for x in img_extent_sat]
# img_extent_sat=[-32.150361957, 30.150361957, 7.150361956, 42.150361956]
# 转换到墨卡托投影
img_extent_merc = transform_extent_pts(img_extent_sat, map_proj, ccrs.Geodetic())
plt.close() 

fig = plt.figure(figsize=(10,10))
ax = plt.axes(projection=map_proj)
ax.coastlines(color='blue')
ax.gridlines(color='black', alpha=0.5, linestyle='--', linewidth=0.75, draw_labels=True)

# 地图范围（经纬度坐标系）:
map_extent_deg = (50., -20., -40., 40.) # 非洲大陆
map_extent_deg = (-31, 38.009232, 2.880476, 42)
# 转换到墨卡托投影
map_extent_merc = transform_extent_pts(map_extent_deg, map_proj, ccrs.Geodetic())
ax.set_extent(map_extent_merc, map_proj)
    
plt.imshow(df, origin='upper', transform=data_crs, extent=img_extent_sat)

核心问题 #

目前数据投影不正确，图像存在轻微偏移。因为是PNG格式的局部图像，用静止卫星投影的transform函数似乎不起作用，想确认：

• 是否有其他修改投影的可行方案？
• 能否在图像加载后事后修改投影？

已知参考点信息 #

我已掌握两个地理坐标点（对应代码中的sta_lon、sta_lat），它们应当和图像中的两个白块精准对齐，相关代码如下：

plt.plot(sta_lon, sta_lat, marker='o', color='red', markersize=8,
        alpha=0.7, transform=ccrs.Geodetic())
plt.text(sta_lon, sta_lat+0.25, sta_name, ha='center', fontsize=18,
         color='red', transform=ccrs.Geodetic()) 

解决方案建议 #

1. 利用已知参考点做地理配准
   既然有对应地理坐标的参考点，这是最可靠的校正方式，步骤如下：

   • 记录图像中两个白块的像素坐标（如(px1, py1)、(px2, py2)）
   • 使用rasterio或affine库计算像素坐标到地理坐标的转换矩阵
   • 给图像赋予地理参考后再进行重投影
     示例代码（基于rasterio）：

   import rasterio
   from rasterio.transform import from_gcps
   from rasterio.warp import reproject, Resampling
   
   # 假设已知像素坐标与对应地理坐标
   gcps = [
       rasterio.control.GroundControlPoint(px1, py1, sta_lon1, sta_lat1),
       rasterio.control.GroundControlPoint(px2, py2, sta_lon2, sta_lat2)
   ]
   transform = from_gcps(gcps)
   
   # 定义输入输出投影
   src_crs = ccrs.Geodetic().to_wkt()
   dst_crs = ccrs.Mercator().to_wkt()
   
   # 临时保存带地理参考的图像
   with rasterio.open('temp_sat.tif', 'w', 
                      driver='GTiff',
                      height=df.shape[0],
                      width=df.shape[1],
                      count=3,  # 根据PNG实际通道数调整
                      dtype=df.dtype,
                      crs=src_crs,
                      transform=transform) as dst:
       dst.write(df.transpose(2, 0, 1))  # 调整通道顺序为(C, H, W)
   
   # 重投影到目标投影
   with rasterio.open('temp_sat.tif') as src:
       dst_transform, dst_width, dst_height = rasterio.warp.calculate_default_transform(
           src.crs, dst_crs, src.width, src.height, *src.bounds)
       kwargs = src.meta.copy()
       kwargs.update({
           'crs': dst_crs,
           'transform': dst_transform,
           'width': dst_width,
           'height': dst_height
       })
   
       with rasterio.open('sat_mercator.tif', 'w', **kwargs) as dst:
           for i in range(1, src.count + 1):
               reproject(
                   source=rasterio.band(src, i),
                   destination=rasterio.band(dst, i),
                   src_transform=src.transform,
                   src_crs=src.crs,
                   dst_transform=dst_transform,
                   dst_crs=dst_crs,
                   resampling=Resampling.bilinear)
   

2. 修正Cartopy的范围与投影参数
   你的代码中使用默认静止卫星全局范围，和局部图像不匹配是偏移的核心原因：

   • 不要用ax2.get_extent()，直接根据图像对应的实际地理范围设置img_extent_sat
   • 确保data_crs的参数和卫星实际参数一致，比如satellite_height、sweep_axis等，这些参数会直接影响投影精度

3. 事后修改投影的可行性
   完全可以事后修改投影，但前提是图像已经具备正确的地理参考信息（坐标系统、范围或控制点）。如果没有，需要先通过参考点或已知范围给图像赋予地理参考，再执行重投影操作。

内容的提问来源于stack exchange，提问作者S.Kociok