Xarray Unit Troubles In 1.0

[jsillin]

Hi all,
I have finally been migrating over to 1.0 as it seems to have a fix for a bug in the wetbulb code. I think I've caught (or at least understand) the changing of some input names that might break scripts written in 0.12, but recently came across this error that has something to do with xarray and units which I'm vaguely aware may have changed going into 1.0 but I'm not sure what exactly this problem is or how to fix it.

Seems to happen when I'm shading CAPE/CIN

                for i in range(s,r):
                    levs_abv_ground = []
                    soundlon = sound_lons[i]
                    sound_temps = temp.interp(lat=soundlat,lon=soundlon)-273.15
                    sound_rh = rh.interp(lat=soundlat,lon=soundlon)
                    sound_pres = temp.lev
                    spres = sfc_pressure.interp(lat=soundlat,lon=soundlon)
                    
                    sound_dp = mpcalc.dewpoint_from_relative_humidity(sound_temps.data*units.degC,sound_rh.data*units.percent)
                    sound_wb = mpcalc.wet_bulb_temperature(sound_pres,sound_temps.data*units.degC,sound_dp)

                    #Only want data above the ground 
                    abv_sfc_temp = mbt.mask_below_terrain(spres,sound_temps,sound_pres)[0]
                    abv_sfc_dewp = mbt.mask_below_terrain(spres,sound_dp,sound_pres)[0]
                    abv_sfc_wetb = mbt.mask_below_terrain(spres,sound_wb,sound_pres)[0]
                    pres_abv_ground = mbt.mask_below_terrain(spres,sound_temps,sound_pres)[1]

                    #sound_wb = sound_wb*units.degC
                    skew = SkewT(fig=fig,rect=(0.75-(0.15*i),plot_elev,.15,.1))

                    parcel_prof = mpcalc.parcel_profile(sound_pres,sound_temps[0].data*units.degC,sound_dp[0])
                    cape = mpcalc.cape_cin(sound_pres,sound_temps.data*units.degC,sound_dp,parcel_prof)
                    capeout = int(cape[0].m)
                    cinout = int(cape[1].m)

                    #skew.ax.axvspan(-30, -10, color='cyan', alpha=0.4)

                    skew.plot(pres_abv_ground,abv_sfc_wetb,'b',linewidth=2)
                    skew.plot(pres_abv_ground,abv_sfc_dewp,'g',linewidth=3)
                    skew.plot(pres_abv_ground,abv_sfc_temp,'r',linewidth=3)

                    if capeout >100:
                        # Shade areas of CAPE and CIN
                        skew.shade_cin(pres_abv_ground, abv_sfc_temp.data*units.degC, parcel_prof)
                        skew.shade_cape(pres_abv_ground, abv_sfc_temp.data*units.degC, parcel_prof)
                        skew.plot(pres_abv_ground,parcel_prof,color='fuchsia',linewidth=1)
                        lcl_pressure, lcl_temperature = mpcalc.lcl(pres_abv_ground[0], abv_sfc_temp.data[0]*units.degC, abv_sfc_dewp[0])
                        skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='forestgreen')
                        skew.ax.axvline(-30, color='teal', linestyle='--', linewidth=1)
                        skew.ax.axvline(-10, color='teal', linestyle='--', linewidth=1)
                    skew.plot(1000,0, 'ko',markerfacecolor='k')

                    skew.ax.axvline(0, color='purple', linestyle='--', linewidth=3)
                    skew.ax.set_ylim((1000,300))
                    skew.ax.axis('off')

works for a while then

Traceback (most recent call last):
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\matplotlib\axis.py", line 1523, in convert_units
    ret = self.converter.convert(x, self.units, self)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\pint\matplotlib.py", line 36, in convert
    return [self._convert_value(v, unit, axis) for v in value]
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\pint\matplotlib.py", line 36, in <listcomp>
    return [self._convert_value(v, unit, axis) for v in value]
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\pint\matplotlib.py", line 46, in _convert_value
    return self._reg.Quantity(value, axis.get_units()).to(unit).magnitude
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\pint\quantity.py", line 203, in __new__
    raise TypeError(f"Quantity cannot wrap upcast type {type(value)}")
TypeError: Quantity cannot wrap upcast type <class 'xarray.core.dataarray.DataArray'>

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "get_soundingmaps.py", line 351, in <module>
    smap.plot_soundings(fig,ax1,prs_temps,prs_relh,sfc_pressure,centerlat,centerlon,domainsize,model,cape=cape_bool,wetbulb=wetbulb_bool)
  File "C:\Users\Jack\Documents\GitHub\plotting\dev\soundingmaps.py", line 147, in plot_soundings
    skew.shade_cin(pres_abv_ground, abv_sfc_temp.data*units.degC, parcel_prof)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\plots\skewt.py", line 740, in shade_cin
    return self.shade_area(pressure[idx], t_parcel[idx], t[idx], which='negative',
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\plots\skewt.py", line 669, in shade_area
    return self.ax.fill_betweenx(*arrs, **fill_args)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\matplotlib\__init__.py", line 1447, in inner
    return func(ax, *map(sanitize_sequence, args), **kwargs)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\matplotlib\axes\_axes.py", line 5315, in fill_betweenx
    return self._fill_between_x_or_y(
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\matplotlib\axes\_axes.py", line 5205, in _fill_between_x_or_y
    getattr(self, f"convert_{dep_dir}units")(dep1))
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\matplotlib\artist.py", line 175, in convert_xunits
    return ax.xaxis.convert_units(x)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\matplotlib\axis.py", line 1525, in convert_units
    raise munits.ConversionError('Failed to convert value(s) to axis '
matplotlib.units.ConversionError: Failed to convert value(s) to axis units: <xarray.DataArray (lev: 15)>
<Quantity([295.21072998 293.02791359 290.80367221 288.53606742 283.8622598
 281.02013484 278.29668786 275.29577887 271.9578587  268.20502548
 263.93479128 259.01299424 253.26804924 246.49220559 238.45676011], 'kelvin')>
Coordinates:
    time       datetime64[ns] 2021-03-10
  * lev        (lev) float64 975.0 950.0 925.0 900.0 ... 500.0 450.0 400.0 350.0
    metpy_crs  object Projection: latitude_longitude

[jthielen]

This looks to be a current sticking point in MetPy (I've opened #1752 to have an issue to keep track of it). Based on some quick exploration, I can get basic SkewT stuff working with just xarray, but not cape/cin shading. Mixing xarray and Pint Quantities reproduces the TypeError/ConversionError you found above, and it looks like you have some manual creation of quantities rather than just using source xarray objects, so I believe that is what going on in your example.

Unfortunately, this likely means either waiting for fixes in MetPy (or maybe even upstream in matplotlib if we can't fix the mixed object issues here) or casting all your xarray objects in the plotting calls to Pint quantities (i.e., data.metpy.unit_array).

[jsillin]

huh I tried hitting each of the arguments into shade_cin and shade_cape with .data.metpy.unit_array and got a whole assortment of errors. Will try to have another go at it tomorrow. You're saying I just need to make sure I'm only passing pint arrays to those functions? That seems doable. Trying to get these soundingmaps up and running before the tornadoes beat me to it haha

[jthielen]

I believe that should do it. Also, not sure if that initial . was a typo, but if not, note that is sound_field.metpy.unit_array (the unit_array attribute on the metpy accessor) rather than sound_field.data.metpy.unit_array (which will give errors since data attribute is a wrapped array (either ndarray or pint quantity)).

Good luck on getting the sounding maps all up and running, and thanks for continuing to bring up your MetPy issues here so that we can work on fixing them for everyone!

[jsillin]

This makes sense! Inexplicably something just worked this morning and I haven't had any issues with this all day. This is the code that works (for now)

                    levs_abv_ground = []
                    soundlon = sound_lons[i]
                    sound_temps = temp.interp(lat=soundlat,lon=soundlon)-273.15
                    sound_rh = rh.interp(lat=soundlat,lon=soundlon)
                    sound_pres = temp.lev
                    spres = sfc_pressure.interp(lat=soundlat,lon=soundlon)
                    sound_dp = mpcalc.dewpoint_from_relative_humidity(sound_temps.data*units.degC,sound_rh.data*units.percent)
                    sound_wb = mpcalc.wet_bulb_temperature(sound_pres,sound_temps.data*units.degC,sound_dp)

                    #Only want data above the ground
                    abv_sfc_temp = spt.mask_below_terrain(spres,sound_temps,sound_pres)[0]
                    abv_sfc_dewp = spt.mask_below_terrain(spres,sound_dp,sound_pres)[0]
                    abv_sfc_wetb = spt.mask_below_terrain(spres,sound_wb,sound_pres)[0]
                    pres_abv_ground = spt.mask_below_terrain(spres,sound_temps,sound_pres)[1]

                    #sound_wb = sound_wb*units.degC
                    skew = SkewT(fig=fig,rect=(0.75-(0.15*i),plot_elev,.15,.1))

                    parcel_prof = mpcalc.parcel_profile(pres_abv_ground,abv_sfc_temp[0].data*units.degC,abv_sfc_dewp[0])
                    cape = mpcalc.cape_cin(pres_abv_ground,abv_sfc_temp.data*units.degC,abv_sfc_dewp,parcel_prof)
                    capeout = int(cape[0].m)
                    cinout = int(cape[1].m)

                    #skew.ax.axvspan(-30, -10, color='cyan', alpha=0.4)

                    skew.plot(pres_abv_ground,abv_sfc_wetb,'b',linewidth=2)
                    skew.plot(pres_abv_ground,abv_sfc_dewp,'g',linewidth=3)
                    skew.plot(pres_abv_ground,abv_sfc_temp,'r',linewidth=3)

                    if capeout >100:
                        # Shade areas of CAPE and CIN
                        print(pres_abv_ground)
                        print(abv_sfc_temp.data*units.degC)
                        print(parcel_prof)
                        skew.shade_cin(pres_abv_ground, abv_sfc_temp.data*units.degC, parcel_prof)
                        skew.shade_cape(pres_abv_ground, abv_sfc_temp.data*units.degC, parcel_prof)
                        skew.plot(pres_abv_ground,parcel_prof,color='fuchsia',linewidth=1)
                        lcl_pressure, lcl_temperature = mpcalc.lcl(pres_abv_ground[0], abv_sfc_temp.data[0]*units.degC, abv_sfc_dewp[0])
                        skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='forestgreen')
                        skew.ax.axvline(-30, color='teal', linestyle='--', linewidth=1)
                        skew.ax.axvline(-10, color='teal', linestyle='--', linewidth=1)
                    skew.plot(975,0, 'ko',markerfacecolor='k')

                    skew.ax.axvline(0, color='purple', linestyle='--', linewidth=3)
                    skew.ax.set_ylim((1000,300))
                    skew.ax.axis('off')

[jsillin]

Also, thanks so much for all the help troubleshooting everything!