I'm attempting to follow this training example to calculate QG omega on NCEP/NCAR data but I'm getting hung up on mpcalc.advection().
It appears as if my dx and dy variables are a different shape, but I'm directly following the routine in an online example that supposedly works.
import numpy as np
import xarray as xr
import metpy.calc as mc
import metpy.constants as mpconstants
from metpy.units import units
# CONSTANTS
# ---------
sigma = 2.0e-6 * units('m^2 Pa^-2 s^-2')
f0 = 1e-4 * units('s^-1')
Rd = mpconstants.Rd
path = './'
uf = 'uwnd.2018.nc'
vf = 'vwnd.2018.nc'
af = 'air.2018.nc'
ads = xr.open_dataset(path+af).metpy.parse_cf()
uds = xr.open_dataset(path+uf).metpy.parse_cf()
vds = xr.open_dataset(path+vf).metpy.parse_cf()
a700 = ads['air'].metpy.sel(
level=700 * units.hPa,
time='2018-01-04T12')
u700 = uds['uwnd'].metpy.sel(
level=700 * units.hPa,
time='2018-01-04T12')
v700 = vds['vwnd'].metpy.sel(
level=700 * units.hPa,
time='2018-01-04T12')
lats = ads['lat'].metpy.unit_array
lons = ads['lon'].metpy.unit_array
X, Y = np.meshgrid(lons,lats)
dx, dy = mc.lat_lon_grid_deltas(lons,lats)
avort = mc.absolute_vorticity(u700, v700,
dx, dy, lats[:,None])
print('Array shape:', avort.shape)
print('DX shape:', dx.shape)
print('DY shape:', dy.shape)
print('U700 shape:', u700.shape)
print('V700 shape:', v700.shape)
vortadv = mc.advection(avort, (u700,v700), (dx,dy)).to_base_units()
Here's the error message, it also looks like I may have a unit issue?
Found lat/lon values, assuming latitude_longitude for projection grid_mapping variable
Found lat/lon values, assuming latitude_longitude for projection grid_mapping variable
Found lat/lon values, assuming latitude_longitude for projection grid_mapping variable
/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/metpy/calc/basic.py:1033: UserWarning: Input over 1.5707963267948966 radians. Ensure proper units are given.
'Ensure proper units are given.'.format(max_radians))
/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/pint/quantity.py:888: RuntimeWarning: invalid value encountered in true_divide
magnitude = magnitude_op(new_self._magnitude, other._magnitude)
Array shape: (73, 144)
DX shape: (73, 143)
DY shape: (72, 144)
U700 shape: (73, 144)
V700 shape: (73, 144)
Traceback (most recent call last):
File "metpy.decomp.py", line 56, in <module>
vortadv = mc.advection(avort, (u700,v700), (dx,dy)).to_base_units()
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/metpy/xarray.py", line 570, in wrapper
return func(*args, **kwargs)
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/metpy/calc/kinematics.py", line 61, in wrapper
ret = func(*args, **kwargs)
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/metpy/calc/kinematics.py", line 320, in advection
wind = _stack(wind)
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/metpy/calc/kinematics.py", line 24, in _stack
return concatenate([a[np.newaxis] if iterable(a) else a for a in arrs], axis=0)
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/metpy/calc/kinematics.py", line 24, in <listcomp>
return concatenate([a[np.newaxis] if iterable(a) else a for a in arrs], axis=0)
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/dataarray.py", line 642, in __getitem__
return self.isel(indexers=self._item_key_to_dict(key))
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/dataarray.py", line 1040, in isel
indexers, drop=drop, missing_dims=missing_dims
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/dataset.py", line 2014, in _isel_fancy
name, var, self.indexes[name], var_indexers
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/indexes.py", line 106, in isel_variable_and_index
new_variable = variable.isel(indexers)
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/variable.py", line 1118, in isel
return self[key]
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/variable.py", line 766, in __getitem__
dims, indexer, new_order = self._broadcast_indexes(key)
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/variable.py", line 612, in _broadcast_indexes
return self._broadcast_indexes_outer(key)
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/variable.py", line 688, in _broadcast_indexes_outer
return dims, OuterIndexer(tuple(new_key)), None
File "/work1/jsa/anconda3/envs/earth/lib/python3.7/site-packages/xarray/core/indexing.py", line 410, in __init__
f"invalid indexer array, does not have integer dtype: {k!r}"
TypeError: invalid indexer array, does not have integer dtype: array(None, dtype=object)
Thanks in advance for any help!
So the problem is that in MetPy 1.0 the signature for advection changed, to be easier to use, to advection(scalar, u, v). Also, since you're working with Xarray data with MetPy 1.0, it can handle all of the coordinate stuff for you--you don't even need to deal with dx and dy manually:
subset = dict(level=700 * units.hPa, time='2018-01-04T12')
a700 = ads['air'].metpy.sel(**subset)
u700 = uds['uwnd'].metpy.sel(**subset)
v700 = vds['vwnd'].metpy.sel(**subset)
avort = mc.absolute_vorticity(u700, v700)
vortadv = mc.advection(avort, u700, v700).to_base_units()
We need to update that training example, but right now I'd recommend looking at the gallery examples 500 hPa Geopotential Heights, Absolute Vorticity, and Winds and 500 hPa Vorticity Advection.