1

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

n0rthern.dancer
  • 11
  • 3
  • Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. – Community Nov 29 '22 at 15:33
  • I have attempted to clarify -- what I need is to know why I'm getting these errors and how to modify my script to successfully compute the geostrophic winds. – n0rthern.dancer Nov 29 '22 at 15:53
  • Does this answer your question? [Getting error message while ploting metpy SkewT](https://stackoverflow.com/questions/73077017/getting-error-message-while-ploting-metpy-skewt) – Michael Delgado Nov 29 '22 at 16:18
  • Thanks. Alas no, it does not fix the problem. Doing z = units.Quantity(z, 'm') resulted in the error `...raise TypeError(f"PlainQuantity cannot wrap upcast type {type(value)}") TypeError: PlainQuantity cannot wrap upcast type ' I think the problem may be that MetPy's calculation functions need to be working with the type DataArray not Dataset, as per this answer: https://stackoverflow.com/questions/68114074/errors-when-using-metpy-to-calculate-specific-humidity Now I need to calculate full model height without using getvar. – n0rthern.dancer Nov 29 '22 at 18:11

3 Answers3

2

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.

DopplerShift
  • 5,472
  • 1
  • 21
  • 20
  • Thanks for your suggestion. Alas, it did not work. Entering `z = metpy_getvar(ncfile, "z", units="m")' resulted in the error. ``....Traceback (most recent call last): File "", line 1, in File "", line 2, in metpy_getvar TypeError: 'str' object is not callable....'' Based on this answer https://stackoverflow.com/questions/68114074/errors-when-using-metpy-to-calculate-specific-humidity I think the problem may be that MetPy's calculation functions need the data type to be DataArray. To test this, I need to get the full model height without using getvar. – n0rthern.dancer Nov 29 '22 at 18:17
  • No, the problem is a name collision in the function I wrote above. I've edited to correct. MetPy will work fine with both `DataArray` and numpy arrays wrapped with Pint units. – DopplerShift Nov 29 '22 at 22:20
  • The corrected code you've written above does seem to successfully add the units. However, when I compute the geostrophic wind from z, I still get the same warnings, and the geostrophic wind component arrays are full of nans, so it's still not working. – n0rthern.dancer Nov 30 '22 at 09:31
  • I got some sensible-looking winds out with the following: #### PH = ncfile.metpy.parse_cf('PH') PHB = ncfile.metpy.parse_cf('PHB') height = ( PH + PHB ) / (9.81 * units.meter / (units.second * units.second)) geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(height) geo_wind_u = geo_wind_u * units.meter / units.second geo_wind_v = geo_wind_v * units.meter / units.second ### However, there are two problems: the height field is then not on the mass grid (while U and V are), and I don't know whether the wrf output is xy or yx ordering (and how to change to xy ordering if needed). – n0rthern.dancer Nov 30 '22 at 09:46
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
Jon Thielen
  • 479
  • 2
  • 7
  • Hi, thanks for your reply! I ran your script and I'm still getting the following warnings after doing geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(z): ``More than one latitude coordinate present for variable. warnings.warn('More than one ' + axis + ' coordinate present for variable' ValueError: Must provide dx/dy arguments or input DataArray with latitude/longitude coordinates.'' It seems like it's not identifying XLAT and XLONG as latitude and longitude coordinates. Should I pass in dx, dy, and lat explicitly? – n0rthern.dancer Dec 02 '22 at 09:24
  • @n0rthern.dancer Sorry about that! This is intended to work seamlessly between xwrf and metpy, so if you're still encountering those issues, I'd consider that a bug on xwrf's end. Would you be able to open an issue for xwrf on GitHub (https://github.com/xarray-contrib/xwrf/issues) with your file and list of installed versions of the libraries used? – Jon Thielen Dec 02 '22 at 20:58
  • Sure, I'll do it tomorrow morning! Thanks for your help! – n0rthern.dancer Dec 02 '22 at 23:33
  • I've opened an issue here: https://github.com/xarray-contrib/xwrf/issues/120 – n0rthern.dancer Dec 03 '22 at 15:43
0

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')
n0rthern.dancer
  • 11
  • 3
  • For `x_dim` and `y_dim`, the values are the index of that dimension within the overall shape. So for data that are e.g. (time, vertical, x, y), `x_dim` should be either 2 or -2 (depends if you count from the front or the end), `y_dim` similarly should be 3 or -1. Honestly, though, if `z` is a DataArray and has the right metadata, you should be able to just do `mpcalc.geostrophic_wind(z)`. I would also recommend calculating geostrophic wind *after* interpolating to a single level. – DopplerShift Nov 30 '22 at 20:26
  • Thanks for the explanation of x_dim and y_dim. Doing ``mpcalc.geostrophic_wind(z)'' is not working even though z is a DataArray, so I suspect the metadata must not be right. – n0rthern.dancer Dec 01 '22 at 09:00
  • Outputting z (and suppressing the array entries) yields `` Coordinates: XLONG (south_north, west_east) float32 -44.18 -44.14 ... 53.16 53.22 XLAT (south_north, west_east) float32 25.65 25.66 25.67 ... 67.91 67.89 XTIME float32 4.32e+03 Time datetime64[ns] 2016-10-04T12:00:00 Dimensions without coordinates: bottom_top, south_north, west_east'' Can you see what's wrong from this? I'm really stuck. – n0rthern.dancer Dec 01 '22 at 09:02
  • If I use getvar to extract the u-wind and then output the array, some additional metadata (not included for z) appears: ``Attributes: FieldType: 104 MemoryOrder: XYZ description: destaggered u-wind component units: m s-1 stagger: coordinates: XLONG XLAT XTIME projection: LambertConformal(stand_lon=-17.0, moad_cen_lat=55.000007629...'' Could the problem be that this metadata is missing from the z DataArray? – n0rthern.dancer Dec 01 '22 at 09:20
  • I added attributes to z using z.attrs[] and it's still not working, so apparently that isn't the problem either. I'm so stuck I'm considering writing my own function to calculate geostrophic wind. – n0rthern.dancer Dec 01 '22 at 09:43