metpy.calc.dewpoint_from_relative_humidity w/ GFS data : ValueError: operands could not be broadcast together with shapes (31,) (34,)

Question

我正在尝试根据 GFS 数据执行 metpy.plots.SkewT()。

当我尝试计算 Td（从 T、rh）时，压力水平不匹配。

是否有一些聪明的方法（w/ xarray？）来切片和切块以便它们排成一行？

以下代码绘制了倾斜 t 上的温度，但是当 metpy.calc.dewpoint_from_relative_humidity 未注释时，它会抱怨：ValueError: operands could not be broadcast together with shapes (31,) (34,)

也欢迎任何有关改进代码的提示....

import xarray
xds0 = xarray.open_dataset(f'https://thredds.ucar.edu/thredds/dodsC/'
    'grib/NCEP/GFS/Global_onedeg/Best')

subset0 = {'method' : 'nearest'}
subset0['lon'] = 238. ; subset0['lat'] = 38.

import metpy.plots
xda1 = xds0.metpy.parse_cf('Temperature_isobaric')
import datetime
vtname0 = xda1.metpy.time.name # eg : time1
subset0[vtname0] = datetime.datetime.utcnow()

T = xda1.metpy.sel(**subset0).values * metpy.units.units(xda1.units) #.squeeze()

vname0 = xda1.metpy.vertical.name # isobaric6 = 0..33
p = xda1[vname0].values * metpy.units.units(xda1[vname0].units)

skewt1 = metpy.plots.SkewT()
skewt1.plot(p, T, 'r')

xda1 = xds0.metpy.parse_cf('Relative_humidity_isobaric') # isobaric = 0..30
rh = xda1.metpy.sel(**subset0).values #* metpy.units.units(xda1.units) #.squeeze()
#Td = metpy.calc.dewpoint_from_relative_humidity(T, rh) # ValueError: operands could not be broadcast together with shapes (31,) (34,) 
#skewt1.plot(p, Td, 'g')

import matplotlib.pyplot as plt
plt.show()

Answer 1

这个过程可以通过使用 MetPy 的 XArray accessor to simplify selection based on generic dimensions like 'time' and 'vertical'. We can also get the values on common levels formed using np.intersect1d:

import datetime
import metpy
import numpy as np
import xarray

xds0 = xarray.open_dataset('https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_onedeg/Best')
temp_da = xds0.metpy.parse_cf('Temperature_isobaric')
rh_da = xds0.metpy.parse_cf('Relative_humidity_isobaric')

# Formulate subset
subset0 = {'method' : 'nearest',
           'lon': 238.,
           'lat': 38.,
           # Can use 'vertical' here if we pass to metpy's version of sel
           # Use the array of values for pressure that are present in both
           'vertical': np.intersect1d(temp_da.metpy.vertical, rh_da.metpy.vertical),
           'time': datetime.datetime.utcnow()}

temp = temp_da.metpy.sel(**subset0)
rh = rh_da.metpy.sel(**subset0)

# Need to rename the coordinate on RH to match that of temperature
rh = rh.rename({rh_da.metpy.vertical.name: temp_da.metpy.vertical.name})

# Get vertical coordinate (pressure) as a unitted array
p = temp.metpy.vertical.metpy.unit_array

# Calculate dewpoint
td = metpy.calc.dewpoint_from_relative_humidity(temp, rh)

# Plot on SkewT
skewt1 = metpy.plots.SkewT()

# For some reason matplotlib doesn't like the temp array with units inside
# a DataArray
skewt1.plot(p, temp.metpy.unit_array, 'r')
skewt1.plot(p, td, 'g')

上面的代码需要 MetPy 1.0 才能原生支持 XArray DataArrays。对于早期版本，对 dewpoint_from_relative_humidity 的调用应该改为：

td = metpy.calc.dewpoint_from_relative_humidity(temp.metpy.unit_array,
                                                rh.metpy.unit_array)

metpy.calc.dewpoint_from_relative_humidity w/ GFS data : ValueError: operands could not be broadcast together with shapes (31,) (34,)

metpy.calc.dewpoint_from_relative_humidity w/ GFS data : ValueError: operands could not be broadcast together with shapes (31,) (34,)

python-xarray

metpy