绘制风刺时QhullError

问题描述 投票:0回答:2

[尝试在mattoplotlib上在地形图上绘制风钩时,出现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
python matplotlib cartopy
2个回答
0
投票

看来,解决方案(排序)是简单地使用与我最初使用的LambertConformal投影不同的投影。不确定原因到底是什么,因此这只是对原始问题的一种规避。


0
投票

此错误很可能是由您输入的包含极点的纬度坐标引起的。以下代码将重现您看到的错误:

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])

由于您的示例无法运行,因此我无法为您的案例提供确切的解决方案。但是,您可以在绘制之前从输入数据中排除这些点,以避免出现此问题并仍然使用所需的投影。

© www.soinside.com 2019 - 2024. All rights reserved.