When attempting to plot wind barbs using matplotlib on a cartopy map, I get a QhullError. I've never seen this error before and my code hasn't changed since I last used it. I've also made sure the packages are up to date and the grib2 file being used is valid by printing the xarray variables. Below is the code:
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')
When running the script through terminal, the below errors show:
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
This error is most likely caused by your input latitude coordinate containing the pole. The following code will reproduce the error you are seeing:
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()
But when the pole is removed there is no error:
lats = np.array([30, 40, 50])
Since your example is not runnable I cannot suggest the exact fix for your case. However, you can likely exclude such points from your input data prior to plotting to avoid this problem and still use your desired projection.