Matplotlib-Cartopy全局Orthographic投影Streamplot触发QhullError

我完全懂你遇到的麻烦——在Orthographic这类半球投影上用streamplot处理全球数据时，确实很容易因为投影转换和插值的问题触发Qhull错误。咱们一步步拆解原因和解决办法：

问题根源 #

Orthographic是正射投影，它只能展示半个地球，当你传入全球范围的经纬度数据时：

• 投影转换会把背面半球的点映射到无效/重叠的坐标，导致transform_vectors抛出警告
• scipy的griddata在插值这些包含奇异点（比如极点附近）或共面的转换后坐标时，Qhull无法处理，直接触发错误

可行的解决方案 #

方案1：裁剪数据到投影可见的半球范围 #

既然Orthographic只能显示一个半球，咱们先把数据裁剪到当前投影可见的区域，避免处理无效点：

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

# 生成测试数据
fakelon = np.linspace(-180, 180, 288)
fakelat = np.linspace(-90, 90, 192)
u = xr.DataArray(np.random.rand(len(fakelat), len(fakelon)), coords=[fakelat, fakelon], dims=['lat', 'lon'])
v = xr.DataArray(np.random.rand(len(fakelat), len(fakelon)), coords=[fakelat, fakelon], dims=['lat', 'lon'])

# 定义正射投影（这里以北极中心为例，可替换为其他中心）
proj = ccrs.Orthographic(central_longitude=0, central_latitude=90)
x, y = np.meshgrid(u['lon'], u['lat'])

# 转换坐标到投影空间，筛选可见点（z坐标>0表示在可见半球内）
proj_coords = proj.transform_points(ccrs.PlateCarree(), x, y)
mask = proj_coords[..., 2] > 0

# 裁剪数据
x_vis = x[mask]
y_vis = y[mask]
u_vis = u.values[mask]
v_vis = v.values[mask]

# 绘图
fig, ax = plt.subplots(subplot_kw={'projection': proj})
ax.set_global()
ax.coastlines()

# 使用可见数据绘制流线，regrid_shape控制插值分辨率
ax.streamplot(x_vis, y_vis, u_vis, v_vis, 
              transform=ccrs.PlateCarree(),
              regrid_shape=200)

plt.show()

方案2：手动转换矢量到投影坐标后绘图 #

先把u、v转换到Orthographic投影的坐标系统下，直接在投影坐标上做streamplot，绕开cartopy自动插值的问题：

# 接上方数据生成部分
proj = ccrs.Orthographic(central_longitude=0, central_latitude=0)
# 转换矢量到投影坐标
x_proj, y_proj, u_proj, v_proj = proj.transform_vectors(ccrs.PlateCarree(), x, y, u.values, v.values)

# 筛选有效数据（排除NaN和无效点）
valid_mask = ~np.isnan(u_proj) & ~np.isnan(v_proj)
x_proj_valid = x_proj[valid_mask]
y_proj_valid = y_proj[valid_mask]
u_proj_valid = u_proj[valid_mask]
v_proj_valid = v_proj[valid_mask]

# 在投影坐标下绘制流线，无需指定transform参数
fig, ax = plt.subplots(subplot_kw={'projection': proj})
ax.set_global()
ax.coastlines()
ax.streamplot(x_proj_valid, y_proj_valid, u_proj_valid, v_proj_valid)

plt.show()

方案3：调整插值方法 #

如果一定要用全球数据，可以修改streamplot的插值方法，把默认的linear换成nearest，避免Qhull的共面问题：

# 接原始代码部分
fig, ax = plt.subplots(subplot_kw={'projection':ccrs.Orthographic()})
ax.set_global()
ax.coastlines()
# 修改插值方法为nearest
ax.streamplot(x, y, u.values, v.values, transform=ccrs.PlateCarree(),
              interpolation_method='nearest')
plt.show()

额外提示 #

• 如果你需要展示全球，Orthographic本身不适合，建议考虑Robinson、PlateCarree这类全球投影
• 不同的Orthographic中心（北极/南极/赤道）可能需要调整数据裁剪的逻辑
• 你的cartopy版本（0.16.0）比较旧，升级到新版本可能会修复一些矢量转换的bug，但上面的方案在旧版本也能正常运行

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