将EUMETSAT卫星GEOS投影数据重采样为墨卡托投影的问题

问题背景 #

我有一个EUMETSAT卫星的scene对象，感兴趣区域经纬度范围为**[40, -10, 105, 40]**（西经40°、南纬10°、东经105°、北纬40°），原始数据采用GEOS投影，区域定义如下：

Area ID: msg_seviri_unknown_3km
Description: MSG SEVIRI unknown area definition with 3 km resolution
Projection: {'a': '6378169', 'h': '35785831', 'lon_0': '45.5', 'no_defs': 'None', 'proj': 'geos', 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 1848
Number of rows: 1665
Area extent: (-604581.2379, -1096647.3571, 4940163.8125, 3899023.914)

尝试转墨卡托投影时遇到两个问题：

1. 使用原始场景的区域参数创建墨卡托投影定义，重采样后投影正常但区域被裁剪：

my_area = create_area_def('my_area', {'proj': 'merc', 'lon_0': 45.5},
                        width=1848, height=1665,
                        area_extent= (-604581.2379, -1096647.3571, 4940163.8125, 3899023.914))

new_scn = scn.resample(my_area)
new_scn.show('natural_color')

2. 显式指定经纬度范围作为area_extent，图像显示完整但水平方向压缩变形（推测是手动指定宽高导致），且原始场景右上角坐标转经纬度得到(inf, inf)（应为卫星视野边缘的太空区域NaN值），新墨卡托区域右上角经纬度与预期不符。

需求：实现GEOS转墨卡托投影时，无需手动指定宽高，生成保持真实比例（类似谷歌地图的最佳实践）的图像。

解决方案：基于目标经纬度范围自动计算匹配比例的墨卡托区域 #

核心思路是：先将感兴趣的经纬度范围转换为墨卡托投影坐标，再结合原始数据分辨率自动推导宽高，避免手动指定导致的比例失调。

代码实现 #

import pyproj
from pyresample import create_area_def

# 假设已加载原始scene对象：scn = Scene(...)

# 1. 定义感兴趣的经纬度范围 [西, 南, 东, 北]
lon_min, lat_min, lon_max, lat_max = 40, -10, 105, 40

# 2. 初始化投影转换器：经纬度 -> 目标墨卡托投影
merc_proj = pyproj.CRS({'proj': 'merc', 'lon_0': 45.5, 'ellps': 'WGS84'})
lonlat_to_merc = pyproj.Transformer.from_crs("EPSG:4326", merc_proj, always_xy=True)

# 转换区域四个角点的坐标
x1, y1 = lonlat_to_merc.transform(lon_min, lat_min)
x2, y2 = lonlat_to_merc.transform(lon_max, lat_max)

# 3. 自动计算宽高：沿用原始3km分辨率，保证比例与清晰度
resolution = 3000  # 单位：米
area_extent_merc = (x1, y1, x2, y2)

width = int((x2 - x1) / resolution)
height = int((y2 - y1) / resolution)

# 4. 创建墨卡托区域定义并重采样
my_area = create_area_def(
    'my_merc_area',
    merc_proj.to_dict(),
    width=width,
    height=height,
    area_extent=area_extent_merc
)

new_scn = scn.resample(my_area)
new_scn.show('natural_color')

关键说明 #

• 比例自动匹配：通过分辨率推导宽高，确保墨卡托投影区域的像素比例与真实地理范围一致，杜绝拉伸/压缩变形
• 边缘NaN处理：仅针对感兴趣经纬度范围创建区域，原始GEOS投影边缘的太空NaN区域会被自动排除，无需额外处理
• 通用基准对齐：采用WGS84椭球体的墨卡托投影，与谷歌地图等通用平面地图标准一致，符合最佳实践

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