如何将Polar Pcolormesh图投影叠加到Cartopy地图指定经纬度？

问题核心 #

需要把生成的极坐标Pcolormesh图叠加到Cartopy地图的任意经纬度位置，直接将theta/phi替换为经纬度不可行，需完成局部极坐标到全局地理坐标的球面转换。

现有基础代码 #

极坐标图生成代码 #

import numpy as np
import matplotlib.pyplot as plt

phis = np.linspace(1e-5,10,10) # 从nadir向上测量的半锥角
thetas = np.linspace(0,2*np.pi,361) # 方位角，0与速度矢量重合

X,Y = np.meshgrid(thetas,phis)
Z   = np.sin(X)**10 + np.cos(10 + Y*X) * np.cos(X)

fig, ax = plt.subplots(figsize=(4,4),subplot_kw=dict(projection='polar'))
im = ax.pcolormesh(X,Y,Z, cmap=plt.cm.jet_r,shading='auto') 

ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2.0)
ax.grid(True)
plt.show()

Cartopy地图代码 #

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import gc

flatMap = ccrs.PlateCarree()
resolution = '110m'
fig = plt.figure(figsize=(12,6), dpi=96)
ax = fig.add_subplot(111, projection=flatMap)

ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -90, 90])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, edgecolor='black', facecolor=cfeature.COLORS['land']))
gc.collect()
plt.show()

解决思路与实现步骤 #

核心逻辑 #

把以目标经纬度为中心的局部极坐标（theta=方位角，phi=从nadir向上的半锥角），转换为全局的经纬度网格，再用Cartopy绘制pcolormesh。

具体步骤 #

1. 定义投影中心：指定要放置极坐标图的经纬度点（比如北京116.4°E, 39.9°N）。
2. 极坐标参数转换：将phi（从nadir向上的角度）转换为与当地天顶的夹角zenith_angle = np.pi - phi，用于计算球面距离。
3. 球面坐标转换：
   • 小范围场景：用平面近似计算经纬度增量，得到全局经纬度网格。
   • 大范围场景：用Cartopy的Geodesic类做精确球面计算，避免平面近似误差。
4. 叠加到地图：将转换后的经纬度网格传入ax.pcolormesh，并指定正确的坐标转换参数。

完整实现代码（精确球面转换版） #

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.geodesic import Geodesic
import gc

# 1. 生成极坐标数据
phis = np.linspace(1e-5, 10, 10)  # 从nadir向上的半锥角（弧度）
thetas = np.linspace(0, 2*np.pi, 361)  # 方位角，0与速度矢量重合
X, Y = np.meshgrid(thetas, phis)
Z = np.sin(X)**10 + np.cos(10 + Y*X) * np.cos(X)

# 2. 指定投影中心经纬度
target_lon, target_lat = 116.4, 39.9
earth_radius = 6371000  # 地球平均半径（米）

# 3. 转换极坐标参数：phi转天顶角，theta转地理方位角（正北为0，顺时针）
zenith_angle = np.pi - Y  # 从当地天顶向下的角度（弧度）
distance = earth_radius * zenith_angle  # 球面距离（米）
azimuth = np.rad2deg(X)  # 转换为度数，Cartopy方位角为正北顺时针

# 4. 精确计算全局经纬度网格
geod = Geodesic()
lon_grid, lat_grid, _ = geod.direct(
    np.full_like(X, target_lon),
    np.full_like(Y, target_lat),
    azimuth,
    distance
)

# 5. 绘制Cartopy地图并叠加投影图
flatMap = ccrs.PlateCarree()
resolution = '110m'
fig = plt.figure(figsize=(12,6), dpi=96)
ax = fig.add_subplot(111, projection=flatMap)

# 绘制背景
ax.imshow(np.tile(np.array([[cfeature.COLORS['water'] * 255]], dtype=np.uint8), [2, 2, 1]), 
          origin='upper', transform=ccrs.PlateCarree(), 
          extent=[-180, 180, -90, 90])
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', resolution, 
                                            edgecolor='black', facecolor=cfeature.COLORS['land']))

# 叠加投影后的pcolormesh
im = ax.pcolormesh(lon_grid, lat_grid, Z, cmap=plt.cm.jet_r, shading='auto', transform=flatMap)
fig.colorbar(im, ax=ax, shrink=0.8)

# 放大到目标区域，方便查看
ax.set_extent([target_lon-20, target_lon+20, target_lat-20, target_lat+20])
ax.gridlines(draw_labels=True, linestyle='--')

gc.collect()
plt.show()

关键说明 #

• 方位角一致性：用户代码中theta的0方向与速度矢量重合，需确保与Cartopy的方位角定义（正北为0，顺时针）匹配，若不匹配需调整theta的偏移量。
• phi的定义转换：因phi是从nadir（天底）向上测量，需转换为天顶角才能正确计算球面距离，若phi的定义变更，此转换逻辑需同步调整。
• 投影适配：若使用非PlateCarree的投影（比如Mercator），只需修改projection参数，pcolormesh的transform仍用ccrs.PlateCarree()即可，Cartopy会自动完成投影转换。

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