python WRF资料的MetPy地转风

mbyulnm0 于 2022-12-02 发布在 Python

关注(0)|答案(3)|浏览(586)

编辑：我开始怀疑下面出现的问题是由元数据引起的，因为即使在纠正了有关单位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'

任何帮助都将不胜感激。
(By顺便说一句，我查看了https://github.com/Unidata/python-training/blob/master/pages/gallery/Ageostrophic_Wind_Example.ipynb处的脚本，发现它没有帮助，因为我不确定需要对WRF数据执行顶部附近的哪个数据操作。）

python

来源：https://stackoverflow.com/questions/74615766/metpy-geostrophic-wind-for-wrf-data

3条答案

按热度按时间

vof42yt11#

wrfpython的getvar函数，虽然它把units作为参数，但在返回数组中的值之前，只使用这个函数（据我所知）来转换它们。要在MetPy中使用这个函数，你需要附加适当的units。我将使用一个小的helper函数来完成这个任务：

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

这样就消除了对丢失单元的抱怨。
编辑：修正匆忙编写的函数中的名称冲突。

赞(0）回复(0）举报 2022-12-02

w6lpcovy2#

我已经取得了一些进展：更新后的脚本和生成的图如下所示。2问题的一部分是我需要将dx、dy和lat传递给函数metpy.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说[... Y，X]顺序的默认值是x_dim = -1和y_dim=-2，但是没有说[... X，Y]顺序的x_dim和y_dim设置为什么，所以我只是猜测。）但是，我仍然收到'UserWarning：未找到水平维编号。默认为（...，Y，X）顺序。''
其次，正如你在图中所看到的，海岸线上的地转风分量有些奇怪。
u-component of geostrophic wind at 300 mb
以下是我目前的脚本：

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

赞(0）回复(0）举报 2022-12-02

b09cbbtk3#

由原始WRF-ARW数据集和通过wrf-python提取的变量表示的数据没有与MetPy关于单位属性、坐标变量和格网投影的假设良好交互的元数据相反，我建议使用xwrf，这是一个最近发布的包，用于以更符合CF约定的方式处理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

赞(0）回复(0）举报 2022-12-02

我来回答

python WRF资料的MetPy地转风

3条答案

相关问题

热门标签

最新问答