绘制MetPy Q向量

Question

我很难从GFS数据中从metpy.calc绘制q向量；当应用ax.set_extent和ax.quiver时，无法正确地绘制这些向量。

计算代码：

import metpy.calc as mpcalc

query.variables('Temperature_isobaric', 'Geopotential_height_isobaric',
                'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')

data = subset_access.get_data(query)

lat = data.variables['lat'][:]
lon = data.variables['lon'][:]

press = data.variables['isobaric'][:] * units.Pa

# Make the pressure same dimensions as the temperature and winds
pressure_for_calc = press[:, None, None]

temperature = data.variables['Temperature_isobaric'][0] * units.kelvin

u = data.variables['u-component_of_wind_isobaric'][0] * units.meter / 
    units.second

v = data.variables['v-component_of_wind_isobaric'][0] * units.meter / 
    units.second

dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)

现在，我试图通过q向量函数运行dx和dy：

Q = mpcalc.q_vector(u,v,temperature,pressure_for_calc,dx,dy)

但我发现了一个错误，我认为这与dx和dy维数有关：

IndexError: too many indices for array

dx.shape, dy.shape
>>> ((101, 160), (100, 161))

好的，这显然是个问题，我需要为每个一维数组，所以我探索了温度数组的形状：

print(temperature.shape)
>>> (31, 101, 161)

所以我尝试取dx和dy的一个子集：

print(dx[:,0].shape, dy[0,:].shape)
>>> (101,) (161,)

然后，我认为这应该与临时数组和按下数组一致，我再次尝试了基于这些子集的计算：

Q = mpcalc.q_vector(u,v,temperature,pressure_for_calc,dx[0,:],dy[:,0])

没有错误，现在感觉很好。检查我假设为x和y分量的q的维数：

 print(Q[0].shape, Q[1].shape)
 >>> (31, 101, 161)
 >>> (31, 101, 161)

好像排好队..。

然而，当我看lats和lons的尺寸时：

lat.shape, lon.shape
>>> ((101,), (161,))

好像是从dx和dy的形状倒过来的？

是我遗漏了什么，还是我只是完全错误地计算了q-向量？这是我的第一次尝试，我不确定我所做的一开始是否正确。

真正的问题是当我尝试用ax.quiver的任何投影来绘制它们时

绘图代码：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# Set Projection of Data
datacrs = ccrs.PlateCarree()

# Set Projection of Plot
plotcrs = ccrs.LambertConformal()

# Add Map Features
states_provinces = cfeature.NaturalEarthFeature(category='cultural',
    name='admin_1_states_provinces_lakes',scale='50m', facecolor='none')

country_borders = cfeature.NaturalEarthFeature(category='cultural',
    name='admin_0_countries',scale='50m', facecolor='none')

# Lat/Lon Extents [lon0,lon1,lat0,lat1]
extent = [-130., -70, 20., 60.]

# Create a plot
fig = plt.figure(figsize=(17., 11.))

# Add the map
ax = plt.subplot(111, projection=plotcrs)

# Add state boundaries to plot
ax.add_feature(states_provinces, edgecolor='k', linewidth=1)

# Add country borders to plot
ax.add_feature(country_borders, edgecolor='black', linewidth=1)

lon_slice = slice(None, None, 8)
lat_slice = slice(None, None, 8)

ax.quiver(lon[lon_slice],lat[lat_slice],Q[0][0,lon_slice,lat_slice], Q[1][0,lon_slice,lat_slice],
color='k',transform=plotcrs)

ax.set_extent(extent, datacrs)

plt.show()

产生的地图：

​

​

当我省略ax.set_extent时，它似乎在绘制Q向量，只是现在没有地图背景.

​

​

所以我想我的两个问题是：

1)我是否从GFS数据中适当地计算了q向量？

( 2)我在密谋中遗漏了什么？

Answer 1

所以我认为你正确地计算了q向量，但是有一个更简单的解。发生此错误是因为您传递的是dx和dy的2D数组，但是您的字段temperature和pressure_for_calc是3D的。NumPy不知道它应该对每个高度重复dx和dy。您可以通过以下方法来完成这一任务，而不需要添加以下内容：

Q = mpcalc.q_vector(u, v, temperature,pressure_for_calc, dx[None, :], dy[None, :])

这样做是插入尺寸1维作为dx和dy的第一个维度，使其余维度不受影响。这使得一切都可以与其他数组正确地排列在一起。

就密谋而言，这是一个经典的CartoPy抓取。您对quiver的调用应该如下所示：

ax.quiver(lon[lon_slice], lat[lat_slice],
          Q[0][0,lon_slice,lat_slice].m, Q[1][0,lon_slice,lat_slice].m,
          color='k', transform=ccrs.PlateCarree())

注意要传递transform=ccrs.PlateCarree()的更改。这是如何告诉CartoPy，您要传递给quiver的数据位于经纬度坐标系中。这还假设您正在绘制的向量在这个坐标系中被正确引用--应该是因为您通过了dx，dy是从mpcalc.lat_lon_grid_deltas()计算出来的。注意，在这种情况下，由于CartoPy将对向量进行重新投影，您需要使用.m来删除这些单元。