在Cartopy的正射投影上绘制风向

问题描述

我有一个GRIB文件，其中包含ECMWF的温度和风数据，该数据在三个压力水平上，我想使用以一定经度，纬度为中心的正交投影与cartopy绘制。

附图显示了我当前获得的结果。该脚本效果很好，但是使用quiver()时，两极周围有一些怪异的副作用，其中风矢量以圆形形式出现。

我尝试使用regrid_shape=20关键字选项。然后效果消失了，我得到了我要寻找的情节。但是，这仅在我在lat_0=90中使用ccrs.Orthographic(lon_0,lat_0)时有效。

如果lat_0小于90，则出现以下错误：

Traceback (most recent call last):
  File "xarray_read_grib_pl.py",line 74,in <module>
    regrid_shape=20
  File "/Users/assink/opt/miniconda3/envs/py37/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py",line 310,in wrapper
    return func(self,*args,**kwargs)
  File "/Users/assink/opt/miniconda3/envs/py37/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py",line 1834,in quiver
    target_extent=target_extent)
  File "/Users/assink/opt/miniconda3/envs/py37/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 "/Users/assink/opt/miniconda3/envs/py37/lib/python3.7/site-packages/cartopy/vector_transform.py",line 68,in _interpolate_to_grid
    method='linear'),)
  File "/Users/assink/opt/miniconda3/envs/py37/lib/python3.7/site-packages/scipy/interpolate/ndgriddata.py",line 221,in griddata
    rescale=rescale)
  File "interpnd.pyx",line 248,in scipy.interpolate.interpnd.LinearNDInterpolator.__init__
  File "qhull.pyx",line 1839,in scipy.spatial.qhull.delaunay.__init__
  File "qhull.pyx",line 357,in scipy.spatial.qhull._Qhull.__init__
scipy.spatial.qhull.QhullError: QH6019 qhull input error (qh_scalelast): can not scale last coordinate to [   0,inf].  Input is cocircular or cospherical.   Use option 'Qz' to add a point at infinity.

While executing:  | qhull d Qbb Qz Qt Q12 Qc
Options selected for Qhull 2019.1.r 2019/06/21:
  run-id 1029403189  delaunay  Qbbound-last  Qz-infinity-point  Qtriangulate
  Q12-allow-wide  Qcoplanar-keep  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood  _maxoutside  0

下面是我使用的完整代码。任何帮助将不胜感激！

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

fid = 'haiti_hemi_era5_pl_20100122_00.grib'
ds = []

for shortName in ('z','t','u','v','w','q'):
    ds += [xr.open_dataset(fid,engine='cfgrib',backend_kwargs={'filter_by_keys': {'shortName': shortName}})]
ds = xr.merge(ds,combine_attrs='override')

temp = ds.t - 273.15
temp_lim = {'min': -50.0,'max': 0.0}

windspeed = np.sqrt(ds.u**2+ds.v**2)
wind_lim = {'min': 20.0,'max': 100.0}

lon_0 = -35.0
lat_0 =  70.0
n_levels = len(ds.isobaricInhPa.values)

lon_idx = slice(None,None,10)
lat_idx = slice(30,10)

fig,ax = plt.subplots(2,n_levels,subplot_kw={'projection': ccrs.Orthographic(lon_0,lat_0)},figsize=(12,8)
                       )

for ilev in range(0,n_levels):
    level_hpa = ds.isobaricInhPa.values[ilev]
    print('Plotting level: {} hPa'.format(level_hpa))
    im1 = temp.isel(isobaricInhPa=ilev).plot.imshow(
                                        cmap='Spectral_r',interpolation='bicubic',ax=ax[0,ilev],transform=ccrs.PlateCarree(),vmin=temp_lim['min'],vmax=temp_lim['max'],add_colorbar=False
                                        )
    ax[0,ilev].set_title('{} hPa'.format(level_hpa))
    ax[0,ilev].coastlines()

    v_mag = windspeed.isel(isobaricInhPa=ilev)
    im2 = v_mag.plot.imshow(cmap='inferno_r',ax=ax[1,vmin=wind_lim['min'],vmax=wind_lim['max'],add_colorbar=False
                           )
    
    # superimpose wind direction
    lat_ = ds.u.isel(isobaricInhPa=ilev).latitude.values[lat_idx]
    lon_ = ds.u.isel(isobaricInhPa=ilev).longitude.values[lon_idx]
    u_ = (ds.u.isel(isobaricInhPa=ilev)).where(
        v_mag > wind_lim['min']).values[lat_idx,lon_idx]
    v_ = (ds.v.isel(isobaricInhPa=ilev)).where(
        v_mag > wind_lim['min']).values[lat_idx,lon_idx]
    ax[1,ilev].quiver(lon_,lat_,u_,v_,scale=1000,headlength=5,minlength=2,headwidth=3,width=0.009,edgecolor='w',linewidth=0.75,regrid_shape=20
                )
    ax[1,ilev].set_title('')
    ax[1,ilev].coastlines()

fig.subplots_adjust(bottom=0.2,top=0.8,left=0.1,right=0.9,wspace=0.1,hspace=0.1)

cb_ax = fig.add_axes([0.15,0.925,0.75,0.02])
cbar = fig.colorbar(im1,cax=cb_ax,orientation='horizontal',extend='both')
cbar.set_label(label='Temperature [deg C]',size=12)

cb_ax = fig.add_axes([0.15,0.1,0.02])
cbar = fig.colorbar(im2,extend='both')
cbar.set_label(label='Wind speed [m/s]',size=12)

plt.savefig('temperature-wind-pl.png')

print(ds)

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）