Supporting the potential vorticity unit PVU

[marcowurth]

As I am not sure how would be the order in supporting this new unit I'll ask here.

I checked on the Default Pint units definition file and the PVU unit does not seems to exist yet in pint.

When then later supported it would be nice apart from being able to convert between the SI units and PVU to also have the calculation function metpy.calc.potential_vorticity_baroclinic() return PVU units. 1PVU equals 1.0 × 10-6 m2 s-1K kg-1 (see here)

[dopplershift]

I think this has been requested a few times. I think we could add another call to units.define() defining "PVU" (or should it be "pvu"?), just like did for gpm: https://github.com/Unidata/MetPy/blob/9b8d85dd9bf000178a6c362764a6dd9a306e491f/src/metpy/units.py#L59

You could test this out right now with our registry:

from metpy.units import units
units.define(...)

If it works like you want it to, I'd be happy to merge a PR that adds it.

[marcowurth]

Okay, now this way by adding

units.define('PVU = 1e-6 m^2 s^-1 K kg^-1')

unit conversion from and to SI units already work great! I have seen the lower case "pvu" as well in a few documents but the upper case "PVU" seems to be used mostly and as well by the highly influential paper of Hoskins et al. (1985, QJRMS).

Apart from that I think it makes most sense to have 1PVU defined on both hemispheres with the same sign as the above and then go on calculating and plotting the "Cyclonic PV" in PVU for both hemispheres (like here). I will look a bit into the metpy.calc.potential_vorticity_baroclinic() function so that it returns on both hemispheres the cyclonic baroclinic (=Ertel's) potential vorticity in PVUs. Maybe @kgoebber could approve this.

After this I will open my first ever pull request wohoo. Thanks a lot already for linking the great resources on learning to request PRs!

[jthielen]

I'll let @kgoebber chime in on those more technical details, but you should also be able to support the lower case "pvu" as an alias:

units.define('PVU = 1e-6 m^2 s^-1 K kg^-1 = _ = pvu')

(the _ is to skip defining a symbol for the unit)

[marcowurth]

Oh yes, this is a good idea to implement both lower and upper case versions! I'm not sure if the extra "= _" is needed tough. I tried the following without it and it seems to be working:

>>> from metpy.units import units
>>> units.define('PVU = 1e-6 m^2 s^-1 K kg^-1 = pvu')
>>> (1*units.pvu).to_base_units()
<Quantity(1e-06, 'kelvin * meter ** 2 / kilogram / second')>
>>> (1*units.PVU).to_base_units()
<Quantity(1e-06, 'kelvin * meter ** 2 / kilogram / second')>

[kgoebber]

I've had a brief chance to look and it appears that there are some instances where individuals multiply the PV in the Southern Hemisphere by -1 to get a "Cyclonic PV" in the literature. My suggestion would be that we add an option that would allow getting back PV with the Southern Hemisphere data multiplied by -1.

The option could be: return_cyclonic_PV and be a boolean variable, defaulting to False.

The question on implementation would be how to apply this correction. It would seem easiest to define a variable of all ones, based on the size of the latitude array that is input, or comes from the xarray dataset, and modify based on the those values if the option is set to True and then add that variable as a multiplication to the returned values (not changing the values by default). Something like,

hemisphere_pv_multiplier = np.ones(latitude.shape)
if return_cyclonic_PV:
    hemisphere_pv_multiplier[latitude < 0] *= -1

return hemisphere_pv_multiplier * (PV_equation)

Variable naming and implementation open to suggestions, but main points of adding the option, but it is off by default would be an acceptable path forward.

[DWesl]

Oh yes, this is a good idea to implement both lower and upper case versions! I'm not sure if the extra "= _" is needed tough. I tried the following without it and it seems to be working:

From the description, the relevant test would be:

>> from metpy.units import units
>> units.define("PVU = 1e-6 m^2 s^-1 K kg^-1 = pvu")
>> data = units.Quantity(1, "PVU")
>> str(data)
'1 PVU'
>> "{:~P}".format(data)
'1 pvu'

Including the underscore means f"{data:~P}" produces 1 PVU instead of 1 pvu.

If you want to get verbose, you could do

units.define("potential_vorticity_unit = 1e-6 m^2 s^-1 K kg^-1 = PVU = pvu")

which might help people who have never run into the unit before.