3

我无法从 GFS 数据中绘制来自 metpy.calc 的 Q 向量；ax.set_extent应用and时，我无法正确绘制向量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)

好像要排队...

但是，当我查看 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）我在绘图中缺少什么？

4

1 回答 1

3

所以我认为你正确地计算了 Q 向量，但是有一个更简单的解决方案。发生错误是因为您传递了 和 的 2D 数组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())

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

于 2019-04-08T20:12:56.270 回答