Geostrophic wind appears off for NARR

[yeechianlow]

Hi,

The geostrophic wind function does not seem to produce accurate results on the NARR; the winds are not parallel to the height contours like they should be (top image). It seems to do better on the CFSR (bottom image). I am running Python 2.7 on a Linux server. Any help or insight on this issue would be appreciated!

Code (for NARR):

# Import the necessary libraries
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
from netCDF4 import num2date
import numpy as np
import xarray as xr
import scipy.ndimage as ndimage

# Load in one month of NARR height data to a dataset
hgt_ds=xr.open_dataset('/nas/repository/hgt.199601.nc')

# Get lat/lon data from file
lat = hgt_ds.lat.data
lon = hgt_ds.lon.data
time = hgt_ds.time

# Make all longitudes positive for plotting purposes
lon_360 = lon
lon_360[lon<0]=lon_360[lon<0]+360

# Use helper function defined above to calculate distance between lat/lon grid points
dx_lat_lon, dy_lat_lon = mpcalc.lat_lon_grid_deltas(lon_360, lat)

#Calculate Coriolis parameter
f = mpcalc.coriolis_parameter(np.deg2rad(lat)).to(units('1/sec'))

# Extract 500-hPa heights and use them to calculate 500-hPa geostrophic winds
hgt_500 = np.asarray(hgt_ds['hgt'].metpy.sel(level=500, time='1996-01-25 06:00'))
(uwnd,vwnd) = mpcalc.geostrophic_wind(units('m')*hgt_500,f,dx_lat_lon, dy_lat_lon, dim_order='yx')
uwnd = np.asarray(uwnd)
vwnd = np.asarray(vwnd)

# Set projection of data and grab data for plotting state boundaries
dataproj = ccrs.PlateCarree()
states_provinces = cfeature.NaturalEarthFeature(category='cultural',name='admin_1_states_provinces_lakes',scale='50m',facecolor='none')

# Make figure and set title
fig = plt.figure(figsize=(17., 11.))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
plt.title('500-hPa heights (m) and geostrophic wind vectors, Jan 25, 1996 06z', {'fontsize':18})

# Set extent and plot map lines
ax.set_extent([-90., -60, 35., 55.], ccrs.PlateCarree())
ax.coastlines('50m', edgecolor='gray', linewidth=0.75)
ax.add_feature(states_provinces, edgecolor='gray', linewidth=0.5)

# Plot 500-hPa heights
clevs_500_hght = np.arange(0, 8000, 60)
cs = ax.contour(lon_360, lat, hgt_500, clevs_500_hght, colors='black', linewidths=1.5, linestyles='solid', transform=dataproj)
plt.clabel(cs, fmt='%d')

# Plot Wind Barbs
ax.quiver(lon_360[0::4,0::4], lat[0::4,0::4], uwnd[0::4,0::4], vwnd[0::4,0::4], scale=2e3, width=0.0015, pivot='middle', transform=ccrs.PlateCarree())
plt.savefig('../figures/NARR_geos_winds2.png')
plt.close()

The only difference in the code for CFSR is to replace

hgt_ds=xr.open_dataset('/nas/repository/hgt.199601.nc')
lat = hgt_ds.lat.data
lon = hgt_ds.lon.data
time = hgt_ds.time
hgt_500 = np.asarray(hgt_ds['hgt'].metpy.sel(level=500, time='1996-01-25 06:00'))

with

hgt_ds=xr.open_dataset('../data/z500.gdas.199601.nc')
lat = np.repeat(np.array([hgt_ds.lat_0.data]),720,axis=0).T
lon = np.repeat(np.array([hgt_ds.lon_0.data]),361,axis=0)
time = hgt_ds.initial_time0_hours
hgt_500 = np.asarray(hgt_ds['HGT_P0_L100_GLL0'].metpy.sel(initial_time0_hours='1996-01-25 06:00'))

Yeechian

[jthielen]

Would you be able to post links to the dataset files used in your two examples? I cannot find anything immediately wrong with your geostrophic wind calculation, so it may be something incorrect in your coordinate/projection metadata in the files, or in the plotting after the calculation.

In either case, I think it would be very good to try figuring out the issue soon, since there will be work on the kinematics calculations (like geostrophic wind here) in the next few days in preparation for the v1.0 release.

[jrleeman]

I agree with @jthielen that it feels like a projection issue, but would you also try with barbs instead of quiver just to make sure it's not something specific to quiver?

[yeechianlow]

@jthielen The NARR file is ftp://ftp.cdc.noaa.gov/Datasets/NARR/pressure/hgt.199601.nc and the CFSR file is https://rda.ucar.edu/data/ds093.1/1996/z500.gdas.199601.grb2 (login required).

@jrleeman I tried using barbs and I get the same problem. I had also tried printing some values from the uwnd and vwnd arrays and confirmed that they were somewhat off, so it's not the plotting that's the problem. A couple of other students in my group had the same problem.

[jthielen]

@yeechianlow Thank you for passing along links to the data files!

After digging into this a bit more with your data files, it luckily is just an issue with specifying the projections, and there doesn't seem to be anything wrong with the geostrophic wind calculation itself. When plotting a vector quantity with Cartopy, you need to have your coordinates be the data's native dimension coordinates. So, for CFSR, which has lat/lon dimensions, your code works perfectly, however, it fails for the y/x dimensions of the NARR data.

See this notebook for an updated implementation of your code (compatible with MetPy 0.11 0.12...sorry had wrong version number), which should solve the projection issue.

[yeechianlow]

@jthielen Interesting, thanks for figuring it out! What does parse_cf do? I get AttributeError: 'MetPyDataArrayAccessor' object has no attribute 'longitude' at the line dx, dy = mpcalc.lat_lon_grid_deltas(hgt.metpy.longitude, hgt.metpy.latitude), though I got around that by directly accessing x and y from the dataset without parsing.

If it's the plotting projection that's the problem, why does accessing the wind arrays with the same indices as the height dataset also give the erroneous values and why doesn't the error also occur on the height field? To use the geostrophic wind for further calculations (such as vorticity or temperature advection), is there any special procedure I need to do to make sure the calculations are currently done on the same grid or projection?

[jthielen]

What does parse_cf do?

This provides the crs coordinate from the CF metadata in the dataset, which is what allows the data_proj Cartopy CRS object to be defined.

I get AttributeError: 'MetPyDataArrayAccessor' object has no attribute 'longitude' at the line dx, dy = mpcalc.lat_lon_grid_deltas(hgt.metpy.longitude, hgt.metpy.latitude), though I got around that by directly accessing x and y from the dataset without parsing.

It sounds like you're on version 0.11 of MetPy. The .longitude and .latitude was just added in 0.12 (I'm sorry, I had the wrong version number in my previous comment), so you could try updating to that.

If it's the plotting projection that's the problem, why does accessing the wind arrays with the same indices as the height dataset also give the erroneous values and why doesn't the error also occur on the height field?

This is the difference between a vector and a scalar field when using Cartopy. While I'm not sure if this is fully technically correct with the implementation, here's how I think of it: a vector field has two components tied to the grid of the data, so the coordinates have to be consistently defined with the components themselves (see https://scitools.org.uk/cartopy/docs/v0.15/matplotlib/geoaxes.html#cartopy.mpl.geoaxes.GeoAxes.quiver). However, with a scalar field, it's essentially just a collection of point values, so you can specify the coordinates for the points however you wish as long as you have the right transform specified.

To use the geostrophic wind for further calculations (such as vorticity or temperature advection), is there any special procedure I need to do to make sure the calculations are currently done on the same grid or projection?

With the calculations themselves, it is just based on having consistent grid spacings and dimension ordering. So, having dx and dy like you've done (and not transposing or reordering anything along the way) should get you the correct results. The errors you encountered appear to be only in the plotting of the data, rather than their calculation.

Also, in not too long, MetPy v1.0 will be released, which is planned to have physically consistent kinematics calculations that automatically "do the right thing" with the coordinate data available from xarray (see #893). So, your calculation will be able to be simplified to just

uwnd, vwnd = mpcalc.geostrophic_wind(hgt_500)

(since f, dx, and dy can all be automatically calculated from the coordinates and projection information).

[dopplershift]

To follow onto @jthielen 's excellent explanation, with the wind vectors (u, v), you can have them defined in one of two ways:

1. Movement along the x/y direction as defined by the grid ('grid-relative wind')
2. Movement north/south ('earth-relative wind')

It's not always readily clear from the data what you're given. Regardless, when you call barbs and quiver in CartoPy, it assumes that the u/v components are in the same coordinate system as the x/y values you give it. So if you're giving it lat/lon coordinates, and passing a transform of something like 'PlateCarree()`, you need to be passing earth-relative wind components. If you are passing x/y values and a native transform (e.g. Lambert Conformal), your wind components need to be grid-relative.

Based on what worked in @jthielen's notebook, I think the NARR is giving you grid-relative winds...which I think is what we've noticed in the past. @kgoebber can you confirm or correct?

[kgoebber]

@dopplershift is correct that the NARR wind components are grid relative, not earth-relative, despite what any sparse documentation exists for the NARR grib files. This same issue arrises with NAM output as well, and they happen to be the same exact grid. There is an example in the Python Gallery that addresses this by plotting the wind barbs in their native x, y coordinate and using the information from the file about its projection to set up the proper transformation.

500-hPa Example: https://unidata.github.io/python-training/gallery/xarray_500hpa_map/

[jthielen]

@dopplershift and @kgoebber Note as well that this example doesn't even use the NARR wind components, just calculated geostrophic wind, so it's automatically grid-relative.

Also, having this discussion while working on #893 brought up some concerns that may be worth discussing...see my comment there.

[dopplershift]

I think @kgoebber meant: https://unidata.github.io/python-training/gallery/500hpa_absolute_vorticity_winds/

[yeechianlow]

@jthielen @dopplershift @kgoebber Thanks for your great explanations! That resolved some headaches for me.

[dopplershift]

No problem. It’s really a sticky problem, and unfortunately the metadata we have available isn’t sufficient to do the right thing automatically. If you happen to think up anything that would have made your life easier, please let us know.