MetPy geostrophic wind for WRF data

Question

Edit: I'm starting to suspect the problems arising below are due to the metadata, because even after correcting the issues raised regarding units mpcalc.geostrophic_wind(z) still issues warnings about the coordinates and ordering. Maybe the function is unable to identify the coordinates from the file? Perhaps this is because WRF output data is non-CF compliant?

I would like to compute geostrophic and ageostrophic winds from WRF-ARW data using the MetPy function mpcalc.geostrophic_wind.

My attempt results in a bunch of errors and I don't know what I'm doing wrong. Can someone tell me how to modify my code to get rid of these errors?

Here is my attempt so far:

#
import numpy as np
from netCDF4 import Dataset
import metpy.calc as mpcalc

from wrf import getvar

# Open the NetCDF file
filename = "wrfout_d01_2016-10-04_12:00:00"
ncfile = Dataset(filename)

# Extract the geopotential height and wind variables
z = getvar(ncfile, "z", units="m")
ua = getvar(ncfile, "ua", units="m s-1")
va = getvar(ncfile, "va", units="m s-1")

# Smooth height data
z = mpcalc.smooth_gaussian(z, 3)

# Compute the geostrophic wind
geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(z)

# Calculate ageostrophic wind components
ageo_wind_u = ua - geo_wind_u
ageo_wind_v = va - geo_wind_v
#

The computation of the geostrophic wind throws several warnings:

>>> # Compute the geostrophic wind
>>> geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(z)
/mnt/.../.../metpy_en/lib/python3.9/site-packages/metpy/xarray.py:355: UserWarning: More than one time coordinate present for variable.
  warnings.warn('More than one ' + axis + ' coordinate present for variable'
/mnt/.../.../lib/python3.9/site-packages/metpy/xarray.py:1459: UserWarning: Horizontal dimension numbers not found. Defaulting to (..., Y, X) order.
  warnings.warn('Horizontal dimension numbers not found. Defaulting to '
/mnt/.../.../lib/python3.9/site-packages/metpy/xarray.py:355: UserWarning: More than one time coordinate present for variable "XLAT".
  warnings.warn('More than one ' + axis + ' coordinate present for variable'
/mnt/.../.../lib/python3.9/site-packages/metpy/xarray.py:1393: UserWarning: y and x dimensions unable to be identified. Assuming [..., y, x] dimension order.
  warnings.warn('y and x dimensions unable to be identified. Assuming [..., y, x] '
/mnt/.../.../lib/python3.9/site-packages/metpy/calc/basic.py:1274: UserWarning: Input over 1.5707963267948966 radians. Ensure proper units are given.
  warnings.warn('Input over {} radians. '

Can anyone tell me why I'm getting these warnings?

And then trying to compute an ageostrophic wind component results in a bunch of errors:

>>> # Calculate ageostrophic wind components
>>> ageo_wind_u = ua - geo_wind_u
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/mnt/.../lib/python3.9/site-packages/xarray/core/_typed_ops.py", line 209, in __sub__    
return self._binary_op(other, operator.sub)
  File "/mnt/.../lib/python3.9/site-packages/xarray/core/dataarray.py", line 4357, in _binary_op    f(self.variable, other_variable)
  File "/mnt/.../lib/python3.9/site-packages/xarray/core/_typed_ops.py", line 399, in __sub__
    return self._binary_op(other, operator.sub)
  File "/mnt/.../lib/python3.9/site-packages/xarray/core/variable.py", line 2639, in _binary_op
    f(self_data, other_data) if not reflexive else f(other_data, self_data)
  File "/mnt/iusers01/fatpou01/sees01/w34926hb/.conda/envs/metpy_env/lib/python3.9/site-packages/pint/facets/numpy/quantity.py", line 61, in __array_ufunc__
    return numpy_wrap("ufunc", ufunc, inputs, kwargs, types)
  File "/mnt/.../lib/python3.9/site-packages/pint/facets/numpy/numpy_func.py", line 953, in numpy_wrap return handled[name](*args, **kwargs)
  File "/mnt/.../lib/python3.9/site-packages/pint/facets/numpy/numpy_func.py", line 513, in _subtract (x1, x2), output_wrap = unwrap_and_wrap_consistent_units(x1, x2)
  File "/mnt/.../lib/python3.9/site-packages/pint/facets/numpy/numpy_func.py", line 130, in unwrap_and_wrap_consistent_units args, _ = convert_to_consistent_units(*args, pre_calc_units=first_input_units)
  File "/mnt/.../lib/python3.9/site-packages/pint/facets/numpy/numpy_func.py", line 111, in convert_to_consistent_units tuple(convert_arg(arg, pre_calc_units=pre_calc_units) for arg in args),
  File "/mnt/.../lib/python3.9/site-packages/pint/facets/numpy/numpy_func.py", line 111, in <genexpr> tuple(convert_arg(arg, pre_calc_units=pre_calc_units) for arg in args),
  File "/mnt/.../lib/python3.9/site-packages/pint/facets/numpy/numpy_func.py", line 93, in convert_arg raise DimensionalityError("dimensionless", pre_calc_units)
pint.errors.DimensionalityError: Cannot convert from 'dimensionless' to 'meter / second'

Any help would be appreciated.

(By the way, I looked at the script at https://github.com/Unidata/python-training/blob/master/pages/gallery/Ageostrophic_Wind_Example.ipynb and did not find it helpful because I'm not sure which of the data manipulations near the top I need to do for the WRF data.)

Answer 1

The data presented by raw WRF-ARW datasets and by variables extracted via wrf-python do not have metadata that interact well with MetPy's assumptions about unit attributes, coordinate variables, and grid projections (from the CF Conventions). Instead, I would recommend using xwrf, a recently released package for working with WRF data in a more CF-Conventions-friendly way. With xwrf, your example would look like:

import metpy.calc as mpcalc
import xarray as xr
import xwrf

# Open the NetCDF file
filename = "wrfout_d01_2016-10-04_12:00:00"
ds = xr.open_dataset(filename).xwrf.postprocess()

# Extract the geopotential height and wind variables
z = ds['geopotential_height']
ua = ds['wind_east']
va = ds['wind_north']

# Smooth height data
z = mpcalc.smooth_gaussian(z, 3)

# Compute the geostrophic wind
geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(z)

# Calculate ageostrophic wind components
ageo_wind_u = ua - geo_wind_u
ageo_wind_v = va - geo_wind_v

Answer 2

I've made some progress: an updated script and the resulting plot are included below. Part of the problem was that I needed to pass dx, dy, and lat into the function metpy.calc.geostrophic_wind, as they were seemingly not being read automatically from the numpy array.

There are still (at least) two problems:

I've passed x_dim=-2 and y_dim=-1 in an effort to set [X,Y] order. (The documentation here https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.geostrophic_wind.html says the default is x_dim = -1 and y_dim=-2 for [...Y,X] order, but does not say what to set x_dim and y_dim to for [...X,Y] order, so I just guessed.) However, I am still getting ``UserWarning: Horizontal dimension numbers not found. Defaulting to (..., Y, X) order.''

Secondly, as you can see in the plot there is something weird going on with the geostrophic wind component at the coastlines.

u-component of geostrophic wind at 300 mb

Here is my current script:

import numpy as np
from netCDF4 import Dataset
import metpy.calc as mpcalc
from metpy.units import units
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap

from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords

# Open the NetCDF file
filename = "wrfout_d01_2016-10-04_12:00:00"
ncfile = Dataset(filename)

z = getvar(ncfile, "z", units="m") * units.meter

# Smooth height data
z = mpcalc.smooth_gaussian(z, 3)

dx = 4000.0 * units.meter
dy = 4000.0 * units.meter

lat = getvar(ncfile, "lat") * units.degrees

geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(z,dx,dy,lat,x_dim=-2,y_dim=-1)

#####

p = getvar(ncfile, "pressure")
z = getvar(ncfile, "z", units="m")

ht_300 = interplevel(z, p, 300)

#geostrophic wind components on 300 mb level
geo_wind_u_300 = interplevel(geo_wind_u, p, 300)
geo_wind_v_300 = interplevel(geo_wind_v, p, 300)

# Get the lat/lon coordinates
lats, lons = latlon_coords(ht_300)

# Get the basemap object
bm = get_basemap(ht_300)

# Create the figure
fig = plt.figure(figsize=(12,12))
ax = plt.axes()

# Convert the lat/lon coordinates to x/y coordinates in the projection space
x, y = bm(to_np(lons), to_np(lats))

# Add the 300 mb height contours
levels = np.arange(8640., 9690., 40.)
contours = bm.contour(x, y, to_np(ht_300), levels=levels, colors="black")
plt.clabel(contours, inline=1, fontsize=10, fmt="%i")

# Add the wind contours
levels = np.arange(10, 70, 5)
geo_u_contours = bm.contourf(x, y, to_np(geo_wind_u_300), levels=levels, cmap=get_cmap("YlGnBu"))
plt.colorbar(geo_u_contours, ax=ax, orientation="horizontal", pad=.05, shrink=0.75)

# Add the geographic boundaries
bm.drawcoastlines(linewidth=0.25)
bm.drawstates(linewidth=0.25)
bm.drawcountries(linewidth=0.25)

plt.title("300 mb height (m) and u-component of geostrophic wind (m s-1) at 1200 UTC on 04-10-2016", fontsize=12)

plt.savefig('geo_u_300mb_04-10-2016_1200_smoothed.png', bbox_inches='tight')

Answer 3

wrfpython's getvar function, while it takes units as a parameter, only uses this (as far as I can tell) to convert values in the arrays before returning them. To use this with MetPy you need to attach proper units. I would do this using a small helper function:

from metpy.units import units

def metpy_getvar(file, name, units_str):
    return getvar(file, name, units=units_str) * units(units_str)

z = metpy_getvar(ncfile, "z", units="m")
ua = metpy_getvar(ncfile, "ua", units="m s-1")
va = metpy_getvar(ncfile, "va", units="m s-1")

That should eliminate the complaints about missing units.

EDIT: Fix name collision in hastily written function.