使用 mpcalc.advection( ) 计算标量场的平流问题

Question

我正在尝试按照 this training example 计算 NCEP/NCAR 数据的 QG 欧米茄，但我被 mpcalc.advection() 挂断了。

看起来好像我的 dx 和 dy 变量的形状不同，但我直接遵循了一个应该有效的在线示例中的例程。

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

这是错误消息，看起来我可能遇到了单位问题？

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)

在此先感谢您的帮助！

Answer 1

所以问题在于，在 MetPy 1.0 中，advection 的签名更改为 advection(scalar, u, v)，以便于使用。此外，由于您正在使用 MetPy 1.0 处理 Xarray 数据，它可以为您处理所有坐标内容——您甚至不需要手动处理 dx 和 dy：

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

我们需要更新该训练示例，但现在我建议您查看图库示例 500 hPa Geopotential Heights, Absolute Vorticity, and Winds and 500 hPa Vorticity Advection。

使用 mpcalc.advection( ) 计算标量场的平流问题

Issue using mpcalc.advection( ) to calculate advection of a scalar field

python

metpy