如何为NWC-SAF卫星影像精准叠加Basemap/Cartopy海陆边界？

我有NWC-SAF提供的某区域卫星影像，以及存储每个像素经纬度值的两个ndarray。尝试用mpl_toolkits.basemap为影像添加精准海陆边界，脚本如下：

def add_basemap(foreground_dimensions: tuple[int, int], low_left_corner_lat, low_left_corner_lon, upper_right_corner_lat, upper_right_corner_lon, dpi):
    """
    Given a satellite image and its coordinates (corners' latitudes & longitudes), Add borders and coastlines
    """
    plt.clf()
    width, height = foreground_dimensions[0], foreground_dimensions[1]
    fig = plt.figure(figsize=(width / (dpi), height / (dpi)), dpi=dpi)
    map_figure = Basemap(projection='cyl',
                         llcrnrlat=low_left_corner_lat,
                         llcrnrlon=low_left_corner_lon,
                         urcrnrlat=upper_right_corner_lat,
                         urcrnrlon=upper_right_corner_lon,
                         resolution='i',
                         )

    fig.subplots_adjust(left=0, right=1, bottom=0, top=1)
    plt.axis('tight')
    plt.axis('off')
    plt.margins(0)
    map_figure.drawcoastlines(linewidth=dpi / 8)

其中foreground_dimensions是卫星影像尺寸，但结果地图与影像匹配不准确。推测是投影差异导致——NWC-SAF用极立体投影（polar stereographic），尝试改用Basemap的npstere投影：

map_figure = Basemap(projection='npstere', boundinglat=max_lat, lon_0=middle_lon, resolution='l')

max_lat是影像最大纬度，middle_lon是影像中心像素经度，但生成的地图完全不符合预期。请问如何用Basemap或Cartopy生成精准匹配的地图？

一、使用Basemap实现精准匹配 #

NWC-SAF的极立体投影有特定参数，不能仅用boundinglat和lon_0，需要严格对齐投影参数：

• 确定基准参数：NWC-SAF北半球极立体投影通常以70°N为标准纬线，中央经线设为影像中心经度，同时指定WGS84地球椭球参数。
• 转换坐标：利用已有的像素经纬度ndarray，通过Basemap方法将经纬度转换为投影坐标，以此定位影像和海岸线。
• 完整示例代码：

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np

def add_basemap_with_polar_stereographic(img_data, lats, lons, dpi=100):
    # 获取影像尺寸
    height, width = img_data.shape[:2]
    fig = plt.figure(figsize=(width/dpi, height/dpi), dpi=dpi)
    
    # 初始化极立体投影（北半球），匹配NWC-SAF参数
    map_figure = Basemap(
        projection='npstere',
        boundinglat=np.min(lats),  # 用影像最小纬度作为边界纬度，避免范围偏移
        lon_0=np.mean(lons),       # 影像中心经度作为中央经线
        lat_ts=70,                 # NWC-SAF标准纬线70°N
        ellps='WGS84',             # 地球椭球参数
        resolution='i'
    )
    
    # 将像素经纬度转换为投影坐标
    x, y = map_figure(lons, lats)
    
    # 绘制卫星影像
    map_figure.imshow(img_data, extent=(np.min(x), np.max(x), np.min(y), np.max(y)), origin='upper')
    
    # 绘制海岸线
    map_figure.drawcoastlines(linewidth=dpi/8, color='white')
    
    # 调整布局，去除边距
    fig.subplots_adjust(left=0, right=1, bottom=0, top=1)
    plt.axis('off')
    plt.margins(0)
    return fig

注：若为南半球，将投影改为spstere，lat_ts设为-70°。

二、使用Cartopy实现精准匹配（推荐，Basemap已停止维护） #

Cartopy对投影的支持更灵活，能直接对齐NWC-SAF的极立体投影：

• 定义极立体投影：使用cartopy.crs.NorthPolarStereo（北半球），指定标准纬线、中央经线和椭球参数。
• 定位影像：将影像的经纬度坐标转换为Cartopy的投影坐标系，确保影像和海岸线使用同一投影。
• 完整示例代码：

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

def add_basemap_cartopy(img_data, lats, lons, dpi=100):
    height, width = img_data.shape[:2]
    fig = plt.figure(figsize=(width/dpi, height/dpi), dpi=dpi)
    
    # 定义NWC-SAF北半球极立体投影
    proj = ccrs.NorthPolarStereo(
        central_longitude=np.mean(lons),
        true_scale_latitude=70,
        globe=ccrs.Globe(ellipse='WGS84')
    )
    
    ax = fig.add_subplot(1, 1, 1, projection=proj)
    
    # 转换经纬度为投影坐标，生成网格
    x, y = proj.transform_points(ccrs.PlateCarree(), lons, lats)[..., :2].T
    
    # 绘制影像
    ax.imshow(img_data, extent=(np.min(x), np.max(x), np.min(y), np.max(y)), origin='upper', transform=proj)
    
    # 添加海岸线
    ax.coastlines(resolution='50m', color='white', linewidth=dpi/8)
    
    # 调整布局
    ax.set_extent([np.min(x), np.max(x), np.min(y), np.max(y)], crs=proj)
    fig.subplots_adjust(left=0, right=1, bottom=0, top=1)
    plt.axis('off')
    return fig

注：南半球使用ccrs.SouthPolarStereo，true_scale_latitude设为-70°。

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