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

I'm trying to do metpy.plots.SkewT() from GFS data.

When I try to calculate Td (from T, rh), the pressure levels don't match up.

Is there some clever way (w/ xarray ?) to slice and dice so they line up ?

The following code plots the temperature on a skew-t, but when the metpy.calc.dewpoint_from_relative_humidity is uncommented, it complains : ValueError: operands could not be broadcast together with shapes (31,) (34,)

Any tips on making the code better are welcome as well ....

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()

Solution

This process can get a little cleaner by making use of MetPy's 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')

The code above requires MetPy 1.0 for its native support of XArray DataArrays. For earlier versions, the call to dewpoint_from_relative_humidity should instead be:

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