首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >WRF资料的MetPy地转风

WRF资料的MetPy地转风
EN

Stack Overflow用户
提问于 2022-11-29 14:40:53
回答 3查看 103关注 0票数 0

编辑:我开始怀疑下面出现的问题是元数据造成的,因为即使在纠正了有关单元mpcalc.geostrophic_wind(z)提出的问题之后,mpcalc.geostrophic_wind(Z)仍然会发出关于坐标和排序的警告。也许函数无法从文件中识别坐标?这可能是因为WRF输出数据不符合CF?

我想用MetPy函数mpcalc.geostrophic_wind从WRF资料中计算地转风和地转风。

我的尝试导致了许多错误,我不知道我做错了什么。有人能告诉我如何修改我的代码以消除这些错误吗?

以下是我迄今为止的尝试:

代码语言:javascript
运行
复制
#
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
#

地转风的计算提出了几个警告:

代码语言:javascript
运行
复制
>>> # 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. '

有人能告诉我为什么我会收到这些警告吗?

然后试图计算一个非地转风分量会产生一系列错误:

代码语言:javascript
运行
复制
>>> # 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'

任何帮助都将不胜感激。

(顺便说一句,我查看了Example.ipynb的脚本,发现它没有帮助,因为我不知道需要对WRF数据做哪些顶部的数据操作。)

EN

回答 3

Stack Overflow用户

发布于 2022-11-29 16:17:33

wrfpython的getvar函数虽然以单元作为参数,但在返回它们之前只使用它(据我所知)转换数组中的值。要将它与MetPy一起使用,您需要附加适当的单元。我会使用一个小助手函数来完成这个任务:

代码语言:javascript
运行
复制
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")

这样就可以消除对失踪部队的投诉。

编辑:修正名称冲突在仓促编写的功能。

票数 2
EN

Stack Overflow用户

发布于 2022-11-30 18:16:05

我已经取得了一些进展:下面包含了一个更新的脚本和相应的情节。部分问题是,我需要将dx、dy和lat传递到函数metpy.calc.geostrophic_wind中,因为它们似乎不是自动从numpy数组中读取的。

仍然(至少)有两个问题:

为了设置X,Y顺序,我已经通过了x_dim=-2和Y=-1。(这里的文档wind.html表示,对于...Y,X order,缺省值是x_dim = -1,y_dim=-2,但是没有说明如何为...X,Y顺序设置x_dim和y_dim,所以我猜对了。)但是,我仍然得到UserWarning:未找到水平维数。违抗(.,Y,X)令

第二,正如你在图中所看到的,海岸线上的地转风分量发生了一些奇怪的事情。

300 mb地转风的u分量

下面是我当前的脚本:

代码语言:javascript
运行
复制
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
EN

Stack Overflow用户

发布于 2022-12-01 21:15:16

原始WRF数据集和通过wrf-python提取的变量提供的数据没有元数据,元数据与MetPy关于单元属性、坐标变量和网格投影(来自CF约定)的假设很好地交互。相反,我建议使用xwrf,这是最近发布的一个包,用于以一种对CF更友好的方式处理WRF数据。使用xwrf,您的示例如下所示:

代码语言:javascript
运行
复制
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
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/74615766

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档