如何用元胞法绘制锋生剖面

Question

我想做一个锋生的横截面我的代码是

domain_1 = os.path.abspath(filenames_in)
info = os.path.join(domain_1, filename)
data = xr.open_dataset(info)
data = data.metpy.parse_cf().squeeze()
print(data)
p1 = data['level'].values[:]* units.hPa
start = (56, -9)
end = (56.5, -6)
cross = cross_section(data, start, end).set_coords(('latitude', 'longitude'))
r,le,t,w =xr.broadcast(cross['r'],cross['level'],cross['t'],cross['w'])
# print('r',r.shape)
# print(le)
# print('w',w.shape)
# print(w)
theta=mpcalc.potential_temperature(level,t)
cross['Potential_temperature'] = xr.DataArray(theta, coords=t.coords, dims=t.dims)
ptemp=cross['Potential_temperature']
cross['u'] = cross['u']
u = cross['u']
cross['v'] = cross['v']
v = cross['v']
# Set subset slice for the geographic extent of data to limit download
lon_slice = slice(-40, 20)
lat_slice = slice(70, 40)
# Grab lat/lon values (GFS will be 1D)
lats = data.latitude.sel(latitude=lat_slice).values
lons = data.longitude.sel(longitude=lon_slice).values
# Compute dx and dy spacing for use in vorticity calculation
dx, dy = mpcalc.lat_lon_grid_deltas(lons,lats)
front = mpcalc.frontogenesis(ptemp, u, v, dx[None,:, :], dy[None,:, :], x_dim=-1, y_dim=-2)
tmpk = data.t.metpy.sel(
latitude=lat_slice, longitude=lon_slice).metpy.unit_array.squeeze()
uwnd = data['u'].metpy.sel(
latitude=lat_slice, longitude=lon_slice).metpy.unit_array.squeeze()
vwnd = data['v'].metpy.sel(
latitude=lat_slice, longitude=lon_slice).metpy.unit_array.squeeze()
potent_temp=mpcalc.potential_temperature(p1[:,None,None], tmpk)
fronto = mpcalc.frontogenesis(potent_temp, uwnd, vwnd, dx[None, :, :], dy[None, :, :], x_dim=-1, y_dim=-2)
# f=cross_section(fronto, start, end).set_coords(('latitude', 'longitude'))

总是有一些bug。如何解决?我不知道我是需要先计算锋面生，然后使用交叉函数还是先使用交叉函数，然后计算锋面生。我试了两种方法，但我不能得到横截面。

ValueError: operands could not be broadcast together with shapes (0,120,241) (19,100) 

如果我先计算锋面形成，然后使用交叉函数

lon_slice = slice(-40, 20) 
lat_slice = slice(70, 40) 
subset = data.metpy.sel(latitude=lat_slice, longitude=lon_slice)
subset['potential_temperature'] = mpcalc.potential_temperature(
subset['level'],
subset['t'])
subset['frontogenesis'] = mpcalc.frontogensis(
subset['potential_temperature'],
subset['u'],
subset['v'])
start = (56, -9) end = (56.5, -6)
cross = cross_section(subset, start, end).set_coords(('latitude', 'longitude'))

我总是遇到这个错误

ValueError: Data missing required coordinate information. Verify that your data have been parsed by MetPy with proper x and y dimension coordinates and added crs coordinate of the correct projection for each variable.

Answer 1

考虑到MetPy的frontogensis计算在2D水平网格上操作，您确实需要在取横截面之前计算frontogensis。类似地，子集应该在计算之前完成，因为MetPy还不支持延迟计算。此外，如果让MetPy为您处理网格增量和单元数组，而不是手动处理它们，您通常会遇到更少的错误，如下所示:

domain_1 = os.path.abspath(filenames_in)
info = os.path.join(domain_1, filename)
data = xr.open_dataset(info)
data = data.metpy.parse_cf().squeeze()
print(data)

lon_slice = slice(-40, 20)
lat_slice = slice(70, 40)
subset = data.metpy.sel(latitude=lat_slice, longitude=lon_slice)

subset['potential_temperature'] = mpcalc.potential_temperature(
subset['level'],
subset['t']
)
subset['frontogenesis'] = mpcalc.frontogensis(
subset['potential_temperature'],
subset['u'],
subset['v']
)

start = (56, -9)
end = (56.5, -6)
cross = cross_section(subset, start, end).set_coords(('latitude', 'longitude'))

(这假设你的数据有适当的"单位"属性和水平和垂直对齐的网格。如果不是这种情况，则需要在打开后进行一些额外的数据清理。)