编辑:我开始怀疑下面出现的问题是元数据造成的,因为即使在纠正了有关单元mpcalc.geostrophic_wind(z)提出的问题之后,mpcalc.geostrophic_wind(Z)仍然会发出关于坐标和排序的警告。也许函数无法从文件中识别坐标?这可能是因为WRF输出数据不符合CF?。
我想用MetPy函数mpcalc.geostrophic_wind从WRF资料中计算地转风和地转风。
我的尝试导致了许多错误,我不知道我做错了什么。有人能告诉我如何修改我的代码以消除这些错误吗?
以下是我迄今为止的尝试:
#
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'
任何帮助都将不胜感激。
(顺便说一句,我查看了Example.ipynb的脚本,发现它没有帮助,因为我不知道需要对WRF数据做哪些顶部的数据操作。)
发布于 2022-11-29 16:17:33
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")
这样就可以消除对失踪部队的投诉。
编辑:修正名称冲突在仓促编写的功能。
发布于 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)令
第二,正如你在图中所看到的,海岸线上的地转风分量发生了一些奇怪的事情。
下面是我当前的脚本:
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')
发布于 2022-12-01 21:15:16
原始WRF数据集和通过wrf-python
提取的变量提供的数据没有元数据,元数据与MetPy关于单元属性、坐标变量和网格投影(来自CF约定)的假设很好地交互。相反,我建议使用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
https://stackoverflow.com/questions/74615766
复制相似问题