metpy - 绘制 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-vector 函数运行 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,)

然后我认为这应该与 temp 和 press 数组对齐，我根据这些子集再次尝试计算:

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)

好像在排队...

但是，当我查看经纬度的尺寸时:

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

生成的 map :

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

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

1) 我是否根据 GFS 数据适本地计算了 Q 向量？

2) 我在绘图中缺少什么？

Answer 1

所以我认为您正确计算了 Q 向量，但有一个更简单的解决方案。发生错误是因为您正在传递 dx 和 dy 的二维数组，但您的字段 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 来删除单位。