绘制风刺时出现 QhullError

Question

当尝试在 cartopy 地图上使用 matplotlib 绘制风刺时，出现 QhullError。我以前从未见过这个错误，自从我上次使用它以来我的代码没有改变。我还通过打印 xarray 变量确保包是最新的并且正在使用的 grib2 文件是有效的。下面是代码：

file = xr.open_dataset('/Users/victoralvarez/prog2/grib/&var_UGRD=on&var_VGRD=on.grb',
                       engine='cfgrib')

# Mask the barbs where winds < 50.
masknum = int(input('BARB THRESHOLD: '))

# Extract the lon/lat.
x = file.variables['longitude'].values
y = file.variables['latitude'].values

# Extract the desired data.
u_wind = file.variables['u'].values * units('m/s')
v_wind = file.variables['v'].values * units('m/s')

# Calculate the wind speed.
wndspd = mpcalc.wind_speed(u_wind, v_wind).to('kt')
wnds_f = wndspd.astype(float)

mask = np.ma.masked_less_equal(wnds_f, masknum).mask
u_wind[mask] = np.nan
v_wind[mask] = np.nan

fig = plt.figure(1, figsize=(15,15))
ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-100,
                                               central_latitude=35,
                                               standard_parallels=(30, 60)))

ax.set_extent([-121, -75, 25, 50], ccrs.PlateCarree())

ax.add_feature(cfeature.OCEAN.with_scale('50m'),  facecolor='#626262',
                                                  edgecolor='black',
                                                  zorder=0,
                                                  linewidth=.5)
ax.add_feature(cfeature.LAND.with_scale('50m'),   edgecolor='black',
                                                  facecolor='#626262',
                                                  zorder=1)
ax.add_feature(cfeature.STATES.with_scale('50m'), linewidth=.5,
                                                  edgecolor='black',
                                                  zorder=5)

b1 = ax.barbs(x, y, u_wind.to('kt').m, v_wind.to('kt').m,
            color='black', length=4.5, regrid_shape=20, pivot='middle',
            linewidth=1.5, zorder=103, transform=ccrs.PlateCarree())
b2 = ax.barbs(x, y, u_wind.to('kt').m, v_wind.to('kt').m,
            color='white', length=4.5, regrid_shape=10, pivot='middle',
            linewidth=0.5, zorder=104, transform=ccrs.PlateCarree())

plt.savefig('img.png', dpi=300, bbox_inches='tight')

当运行脚本通过终端时，出现以下错误：

Traceback (most recent call last):
  File "winds_barb.py", line 63, in <module>
    linewidth=0.5, zorder=104, transform=ccrs.PlateCarree())
  File "/anaconda3/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py", line 1826, in barbs
    target_extent=target_extent)
  File "/anaconda3/lib/python3.7/site-packages/cartopy/vector_transform.py", line 146, in vector_scalar_to_grid
    return _interpolate_to_grid(nx, ny, x, y, u, v, *scalars, **kwargs)
  File "/anaconda3/lib/python3.7/site-packages/cartopy/vector_transform.py", line 68, in _interpolate_to_grid
    method='linear'),)
  File "/anaconda3/lib/python3.7/site-packages/scipy/interpolate/ndgriddata.py", line 222, in griddata
    rescale=rescale)
  File "interpnd.pyx", line 248, in scipy.interpolate.interpnd.LinearNDInterpolator.__init__
  File "qhull.pyx", line 1828, in scipy.spatial.qhull.Delaunay.__init__
  File "qhull.pyx", line 354, in scipy.spatial.qhull._Qhull.__init__
scipy.spatial.qhull.QhullError: QH6019 qhull input error: can not scale last coordinate.  Input is cocircular
   or cospherical.   Use option 'Qz' to add a point at infinity.

While executing:  | qhull d Qz Qt Qbb Q12 Qc
Options selected for Qhull 2015.2.r 2016/01/18:
  run-id 1133843321  delaunay  Qz-infinity-point  Qtriangulate  Qbbound-last
  Q12-no-wide-dup  Qcoplanar-keep  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood

Answer 1

似乎对此的解决方案（大概）是使用与我最初使用的 LambertConformal 投影不同的投影。不确定哪里出了问题，所以这只是对原始问题的规避。

Answer 2

这个错误很可能是由于您输入的包含极点的纬度坐标造成的。以下代码将重现您看到的错误：

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

lons = np.array([-110, -100, -90])
lats = np.array([-90, 30, 40, 50])
u = np.ones([len(lats), len(lons)]) * 10
v = np.ones_like(u) * 10

p = ccrs.LambertConformal(
    central_longitude=-100,
    central_latitude=35,
    standard_parallels=(30, 60),
)

ax = plt.axes(projection=p)
ax.coastlines()
ax.set_extent([-121, -75, 25, 50], crs=ccrs.PlateCarree())
ax.barbs(lons, lats, u, v, regrid_shape=3, transform=ccrs.PlateCarree())

plt.show()

但是拆杆的时候没有报错:

lats = np.array([30, 40, 50])

由于您的示例不可运行，因此我无法针对您的情况建议确切的修复方法。但是，您可以在绘图之前从输入数据中排除这些点以避免此问题并仍然使用您想要的投影。

绘制风刺时出现 QhullError

QhullError When Plotting Wind Barbs

python

matplotlib

cartopy