Plotting MetPy Q-vectors

Question

I'm having trouble plotting Q-vectors from metpy.calc from GFS data; I can't get the vectors to plot correctly when applying the ax.set_extent and ax.quiver

Calculation code:

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)

Now I was trying to run the dx and dy through the q-vector function:

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

But I was getting an error I believe was related to the dx and dy dimensions:

IndexError: too many indices for array

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

Ok, so that was clearly a problem; I need a 1d array for each, so I probed the shape of the temperature array:

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

So I tried taking a subset of the dx and dy:

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

Then I figured this should line up with the temp and press arrays and I tried the calculation again based off these subsets:

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

No errors, feeling good now. Checking the dimensions of Q which I assumed to be the x and y-components:

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

Seems to line up...

However, when I look at the dimensions of lats and lons:

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

It seems to be backwards from the shapes of dx and dy?

Am I missing something or did I just go about calculating the Q-vectors completely wrong? This was my first attempt and I'm not sure if what I was doing was correct in the first place.

The real problem comes when I try and plot them with any projection with ax.quiver

Plotting code:

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

Resulting map:

When I leave out the ax.set_extent it seems to plot the Q-vectors, just with no map background now...

So I guess my two questions are:

1) Did I calculate the Q-vectors appropriately from the GFS data?

2) What am I missing for the plotting?

Answer 1

So I think you calculate the Q-vectors properly, but there was a simpler solution. The error occurred because you were passing 2D arrays of dx and dy, but your fields temperature and pressure_for_calc are 3D. NumPy doesn't know that it should repeat the dx and dy for each height level. You can accomplish this without subsetting by:

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

What this does is insert a size 1 dimension as the first dimension for dx and dy, leaving the remaining dimensions unaffected. This allows everything to line up properly with the other arrays.

As far as plotting is concerned, this is a classic CartoPy gotcha. Your call to quiver should look like:

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

Note the change to pass transform=ccrs.PlateCarree(). This is how to tell CartoPy that the data you are passing to quiver are in a longitude/latitude coordinate system. This also assumes that the vectors you are plotting are properly referenced in this coordinate system--they should be since you passed dx, dy calculated from mpcalc.lat_lon_grid_deltas(). Note in this case, since CartoPy is going to be doing some reprojection of the vectors, you need to use .m to drop the units.