使用 mpcalc 计算锋生时出错
Error computing frontogenesis with mpcalc
我的目标是使用 mpcalc.frontogenesis 方法计算给定压力水平下的锋生。
我在NOMADS上读到了相关数据:
t7 = data['temp'].sel(lev=700.0,lat=lats,lon=lons)
u7 = data['u'].sel(lev=700.0,lat=lats,lon=lons).squeeze()*1.94384449
v7 = data['v'].sel(lev=700.0,lat=lats,lon=lons).squeeze()*1.94384449
h7 = data['gph'].sel(lev=700.0,lat=lats,lon=lons).squeeze()
计算出相关的 coordinate/grid 间距值:
x, y = t7.metpy.coordinates('x', 'y')
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)
lat, lon = xr.broadcast(y, x)
用单位标记每个数组:
t7u = t7.values * units.degC
u7u = u7.values * units.knots
v7u = v7.values * units.knots
然后我去计算FGEN:
h7_fgen = mpcalc.frontogenesis(t7u,u7u,v7u,dx,dy,dim_order='xy')
产生以下错误:
ValueError: operands could not be broadcast together with shapes (359,219) (358,220)
温度、u、v数组的形状分别为:
(220, 360)
(220, 360)
(220, 360)
因此,虽然可以理解 dx 和 dy 的长度为 n-1,但我不确定为什么这不起作用。
如果直接从 NOMADS GFS 源中提取数据,看起来您的数据应该是 yx 维度顺序(这也是 MetPy 的 lat_lon_grid_deltas
函数假定的顺序)。因此,我认为脚本中关于维度排序的东西搞砸了。解决这个问题,并额外使用 MetPy 的一些 xarray 集成功能来清理实现,我们有:
Pre MetPy v1.0
# Imports
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
# Open data and parse for CRS
data = xr.open_dataset('http://nomads.ncep.noaa.gov/dods/gfs_0p25_1hr/gfs20201230/gfs_0p25_1hr_12z').metpy.parse_cf()
# Subset out the desired data
# Note that this leaves in time, but that can also be selected out if only
# one time is needed and you want to have smaller data arrays to compute on
t7 = data['tmpprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310))
u7 = data['ugrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
v7 = data['vgrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
h7 = data['hgtprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
# Clean up the xarray metadata for MetPy's use, since provided data is missing
# the units attribute
t7.attrs['units'] = 'K'
u7.attrs['units'] = v7.attrs['units'] = 'm/s'
h7.attrs['units'] = 'm'
# Get grid deltas
dx, dy = mpcalc.grid_deltas_from_dataarray(t7)
# Compute frontogenesis
h7_fgen = mpcalc.frontogenesis(t7, u7, v7, dx, dy)
Post MetPy v1.0
在 MetPy v1.0 中,网格参数是自动从 DataArray 输入中提取的,因此,这可以简化为
# Imports
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
# Open data and parse for CRS
data = xr.open_dataset('http://nomads.ncep.noaa.gov/dods/gfs_0p25_1hr/gfs20201230/gfs_0p25_1hr_12z').metpy.parse_cf()
# Subset out the desired data
t7 = data['tmpprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310))
u7 = data['ugrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
v7 = data['vgrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
h7 = data['hgtprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
# Clean up the xarray metadata for MetPy's use, since provided data is missing
# the units attribute
t7.attrs['units'] = 'K'
u7.attrs['units'] = v7.attrs['units'] = 'm/s'
h7.attrs['units'] = 'm'
# Compute frontogenesis
h7_fgen = mpcalc.frontogenesis(t7, u7, v7)
我的目标是使用 mpcalc.frontogenesis 方法计算给定压力水平下的锋生。
我在NOMADS上读到了相关数据:
t7 = data['temp'].sel(lev=700.0,lat=lats,lon=lons)
u7 = data['u'].sel(lev=700.0,lat=lats,lon=lons).squeeze()*1.94384449
v7 = data['v'].sel(lev=700.0,lat=lats,lon=lons).squeeze()*1.94384449
h7 = data['gph'].sel(lev=700.0,lat=lats,lon=lons).squeeze()
计算出相关的 coordinate/grid 间距值:
x, y = t7.metpy.coordinates('x', 'y')
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)
lat, lon = xr.broadcast(y, x)
用单位标记每个数组:
t7u = t7.values * units.degC
u7u = u7.values * units.knots
v7u = v7.values * units.knots
然后我去计算FGEN:
h7_fgen = mpcalc.frontogenesis(t7u,u7u,v7u,dx,dy,dim_order='xy')
产生以下错误:
ValueError: operands could not be broadcast together with shapes (359,219) (358,220)
温度、u、v数组的形状分别为:
(220, 360)
(220, 360)
(220, 360)
因此,虽然可以理解 dx 和 dy 的长度为 n-1,但我不确定为什么这不起作用。
如果直接从 NOMADS GFS 源中提取数据,看起来您的数据应该是 yx 维度顺序(这也是 MetPy 的 lat_lon_grid_deltas
函数假定的顺序)。因此,我认为脚本中关于维度排序的东西搞砸了。解决这个问题,并额外使用 MetPy 的一些 xarray 集成功能来清理实现,我们有:
Pre MetPy v1.0
# Imports
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
# Open data and parse for CRS
data = xr.open_dataset('http://nomads.ncep.noaa.gov/dods/gfs_0p25_1hr/gfs20201230/gfs_0p25_1hr_12z').metpy.parse_cf()
# Subset out the desired data
# Note that this leaves in time, but that can also be selected out if only
# one time is needed and you want to have smaller data arrays to compute on
t7 = data['tmpprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310))
u7 = data['ugrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
v7 = data['vgrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
h7 = data['hgtprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
# Clean up the xarray metadata for MetPy's use, since provided data is missing
# the units attribute
t7.attrs['units'] = 'K'
u7.attrs['units'] = v7.attrs['units'] = 'm/s'
h7.attrs['units'] = 'm'
# Get grid deltas
dx, dy = mpcalc.grid_deltas_from_dataarray(t7)
# Compute frontogenesis
h7_fgen = mpcalc.frontogenesis(t7, u7, v7, dx, dy)
Post MetPy v1.0
在 MetPy v1.0 中,网格参数是自动从 DataArray 输入中提取的,因此,这可以简化为
# Imports
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
# Open data and parse for CRS
data = xr.open_dataset('http://nomads.ncep.noaa.gov/dods/gfs_0p25_1hr/gfs20201230/gfs_0p25_1hr_12z').metpy.parse_cf()
# Subset out the desired data
t7 = data['tmpprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310))
u7 = data['ugrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
v7 = data['vgrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
h7 = data['hgtprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
# Clean up the xarray metadata for MetPy's use, since provided data is missing
# the units attribute
t7.attrs['units'] = 'K'
u7.attrs['units'] = v7.attrs['units'] = 'm/s'
h7.attrs['units'] = 'm'
# Compute frontogenesis
h7_fgen = mpcalc.frontogenesis(t7, u7, v7)