编辑:我开始怀疑下面出现的问题是由于元数据，因为即使在纠正了关于单位mpcalc.geostrophic_wind(z)提出的问题之后，仍然发出关于坐标和顺序的警告。也许函数无法从文件中识别坐标?也许这是因为WRF输出数据不兼容cf ?
我想使用MetPy函数mpcalc.geostrophic_wind从WRF-ARW数据计算地转风和地转风。
我的尝试导致一堆错误，我不知道我做错了什么。有人能告诉我如何修改我的代码，以摆脱这些错误吗?
这是我到目前为止的尝试:
#
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
#
地转风的计算抛出了几个警告:
>>> # 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. '
谁能告诉我为什么我得到这些警告?
然后试图计算非地转风分量会导致一系列错误:
>>> # 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'
如有任何帮助，不胜感激。
(顺便说一下，我查看了https://github.com/Unidata/python-training/blob/master/pages/gallery/Ageostrophic_Wind_Example.ipynb上的脚本，并没有发现它有帮助，因为我不确定需要对WRF数据执行顶部附近的哪些数据操作。)
wrfpython的getvar函数，虽然它接受单位作为参数，但只使用它(据我所知)在返回数组之前转换数组中的值。要在MetPy中使用此功能，您需要附加适当的单元。我将使用一个小的辅助函数:
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")
这样应该可以消除关于缺少部件的投诉。
编辑:修复匆忙编写的函数名称冲突。
原始WRF-ARW数据集和通过wrf-python提取的变量所呈现的数据不具有与MetPy关于单位属性、坐标变量和网格投影(来自CF公约)的假设良好交互的元数据。相反，我建议使用xwrf，这是一个最近发布的软件包，用于以更适合cf - conventions的方式处理WRF数据。对于xwrf，您的示例如下:
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
我已经取得了一些进展:更新的脚本和生成的情节包含在下面。问题的一部分是我需要将dx, dy和lat传递到函数mepy .calc中。Geostrophic_wind，因为它们似乎没有自动从numpy数组中读取。
仍然(至少)有两个问题:
我传递了x_dim=-2和y_dim=-1来设置[X,Y]顺序。(这里的文档https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.geostrophic_wind.html说[…]的默认值是x_dim = -1和y_dim=-2。Y,X]顺序，但没有说明如何设置x_dim和y_dim为[…]X,Y]顺序，所以我只是猜测。)然而，我仍然得到' ' UserWarning:未找到水平尺寸数。默认为(…)， Y, X) order. "
其次，正如你在图中看到的，海岸线上的地转风成分发生了一些奇怪的事情。
300 mb地转风u分量
下面是我当前的脚本:
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')