Using Astropy Units

Let’s start with some preliminary imports. To execute a cell in a Jupyter notebook, press Shift + Enter. If you need to restart the notebook, please execute this cell again.

If using Google Colab, click Run anyway when prompted. (If prompted again, select Restart runtime when the installation finishes.)

%matplotlib inline

import sys

if 'google.colab' in str(get_ipython()):
    if 'plasmapy' not in sys.modules:
        !pip install astropy matplotlib numpy

import astropy.units as u
from astropy.visualization import quantity_support
from astropy import constants
import matplotlib.pyplot as plt
import numpy as np

Motivation

In scientific computing, we often represent physical quantities as numbers.

distance_in_miles = 50
time_in_hours = 2
velocity_in_mph = distance_in_miles / time_in_hours
print(velocity_in_mph)

25.0

Representing a physical quantity as a number has risks. We might unknowingly perform operations with different units, like time_in_seconds + time_in_hours. We might even accidentally perform operations with physically incompatible units, like length + time, without catching our mistake. Unit conversion errors can be costly mistakes, such as the loss of spacecraft like the Mars Climate Orbiter.

We can avoid these problems by using a units package. This notebook introduces astropy.units with an emphasis on the functionality needed to work with plasmapy.particles and plasmapy.formulary. We typically import astropy.units subpackage as u.

Unit essentials

We can create a physical quantity by multiplying or dividing a number or array with a unit.

distance = 60 * u.km
print(distance)

60.0 km

This operation creates a Quantity object: a number, sequence, or array that has been assigned a physical unit.

type(distance)

astropy.units.quantity.Quantity

We can create an object by using the Quantity class itself.

time = u.Quantity(120, u.min)

We can create Quantity objects with compound units.

88 * u.imperial.mile / u.hour

$88 \; \mathrm{\frac{mi}{h}}$

We can even create Quantity objects that are explicitly dimensionless.

3 * u.dimensionless_unscaled

$3 \; \mathrm{}$

We can also create a Quantity based off of a NumPy array or a list.

np.array([2.5, 3.2, 1.1]) * u.kg

$[2.5,~3.2,~1.1] \; \mathrm{kg}$

[2, 3, 4] * u.m / u.s

$[2,~3,~4] \; \mathrm{\frac{m}{s}}$

Unit operations

Operations between Quantity objects handle unit conversions automatically. We can add Quantity objects together as long as their units have the same physical type.

1 * u.m + 25 * u.cm

$1.25 \; \mathrm{m}$

Units get handled automatically during operations like multiplication, division, and exponentiation.

velocity = distance / time
print(velocity)

0.5 km / min

area = distance**2
print(area)

3600.0 km2

Attempting an operation between physically incompatible units gives us an error, which we can use to find bugs in our code.

3 * u.m + 3 * u.s

---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/quantity_helper/helpers.py:77, in get_converters_and_unit(f, unit1, unit2)
     76 try:
---> 77     converters[changeable] = get_converter(unit2, unit1)
     78 except UnitsError:

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/quantity_helper/helpers.py:42, in get_converter(from_unit, to_unit)
     39 """Like Unit.get_converter, except returns None if no scaling is needed,
     40 i.e., if the inferred scale is unity.
     41 """
---> 42 converter = from_unit.get_converter(to_unit)
     43 return None if converter is unit_scale_converter else converter

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/core.py:1156, in UnitBase.get_converter(self, other, equivalencies)
   1154                 return lambda v: b(converter(v))
-> 1156 raise exc

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/core.py:1139, in UnitBase.get_converter(self, other, equivalencies)
   1138 try:
-> 1139     return self._apply_equivalencies(
   1140         self, other, self._normalize_equivalencies(equivalencies)
   1141     )
   1142 except UnitsError as exc:
   1143     # Last hope: maybe other knows how to do it?
   1144     # We assume the equivalencies have the unit itself as first item.
   1145     # TODO: maybe better for other to have a `_back_converter` method?

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/core.py:1090, in UnitBase._apply_equivalencies(self, unit, other, equivalencies)
   1088 other_str = get_err_str(other)
-> 1090 raise UnitConversionError(f"{unit_str} and {other_str} are not convertible")

UnitConversionError: 's' (time) and 'm' (length) are not convertible

During handling of the above exception, another exception occurred:

UnitConversionError                       Traceback (most recent call last)
Cell In[13], line 1
----> 1 3 * u.m + 3 * u.s

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/quantity.py:696, in Quantity.__array_ufunc__(self, function, method, *inputs, **kwargs)
    694     return NotImplemented
    695 else:
--> 696     raise e

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/quantity.py:641, in Quantity.__array_ufunc__(self, function, method, *inputs, **kwargs)
    636 # Determine required conversion functions -- to bring the unit of the
    637 # input to that expected (e.g., radian for np.sin), or to get
    638 # consistent units between two inputs (e.g., in np.add) --
    639 # and the unit of the result (or tuple of units for nout > 1).
    640 try:
--> 641     converters, unit = converters_and_unit(function, method, *inputs)
    643     out = kwargs.get("out", None)
    644     # Avoid loop back by turning any Quantity output into array views.

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/quantity_helper/converters.py:181, in converters_and_unit(function, method, *args)
    178 units = [getattr(arg, "unit", None) for arg in args]
    180 # Determine possible conversion functions, and the result unit.
--> 181 converters, result_unit = ufunc_helper(function, *units)
    183 if any(converter is False for converter in converters):
    184     # for multi-argument ufuncs with a quantity and a non-quantity,
    185     # the quantity normally needs to be dimensionless, *except*
   (...)
    188     # can just have the unit of the quantity
    189     # (this allows, e.g., `q > 0.` independent of unit)
    190     try:
    191         # Don't fold this loop in the test above: this rare case
    192         # should not make the common case slower.

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/quantity_helper/helpers.py:79, in get_converters_and_unit(f, unit1, unit2)
     77     converters[changeable] = get_converter(unit2, unit1)
     78 except UnitsError:
---> 79     raise UnitConversionError(
     80         f"Can only apply '{f.__name__}' function to quantities "
     81         "with compatible dimensions"
     82     )
     84 return converters, unit1

UnitConversionError: Can only apply 'add' function to quantities with compatible dimensions

Quantity arrays behave very similarly to NumPy arrays. In fact, Quantity is a subclass of numpy.ndarray.

balmer_series = [656, 486, 434, 410] * u.nm
Hα = balmer_series[0]
print(Hα)

656.0 nm

np.max(balmer_series)

$656 \; \mathrm{nm}$

⚠️ Most frequently encountered NumPy and SciPy functions can be used with Quantity objects. However, there are some functions that cause Quantity objects to lose their units!

Unit conversions

The .to method allows us to convert a Quantity to different units of the same physical type. This method accepts strings that represent a unit (including compound units) or a unit object.

velocity.to("m/s")

$8.3333333 \; \mathrm{\frac{m}{s}}$

velocity.to(u.m / u.s)

$8.3333333 \; \mathrm{\frac{m}{s}}$

The si and cgs attributes convert the Quantity to SI or CGS units, respectively.

velocity.si

$8.3333333 \; \mathrm{\frac{m}{s}}$

velocity.cgs

$833.33333 \; \mathrm{\frac{cm}{s}}$

Detaching units and values

The value attribute of a Quantity provides the number or array without the unit.

time.value

120.0

The unit attribute of a Quantity provides the unit without the value.

time.unit

$\mathrm{min}$

Equivalencies

Plasma scientists often use the electron-volt (eV) as a unit of temperature. This is a shortcut for describing the thermal energy per particle, or more accurately the temperature multiplied by the Boltzmann constant, \(k_B\). Because an electron-volt is a unit of energy rather than temperature, we cannot directly convert electron-volts to kelvin.

u.eV.to("K")

---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
Cell In[22], line 1
----> 1 u.eV.to("K")

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/core.py:1227, in UnitBase.to(self, other, value, equivalencies)
   1225     return UNITY
   1226 else:
-> 1227     return self.get_converter(Unit(other), equivalencies)(value)

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/core.py:1156, in UnitBase.get_converter(self, other, equivalencies)
   1153             else:
   1154                 return lambda v: b(converter(v))
-> 1156 raise exc

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/core.py:1139, in UnitBase.get_converter(self, other, equivalencies)
   1137 # if that doesn't work, maybe we can do it with equivalencies?
   1138 try:
-> 1139     return self._apply_equivalencies(
   1140         self, other, self._normalize_equivalencies(equivalencies)
   1141     )
   1142 except UnitsError as exc:
   1143     # Last hope: maybe other knows how to do it?
   1144     # We assume the equivalencies have the unit itself as first item.
   1145     # TODO: maybe better for other to have a `_back_converter` method?
   1146     if hasattr(other, "equivalencies"):

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/units/core.py:1090, in UnitBase._apply_equivalencies(self, unit, other, equivalencies)
   1087 unit_str = get_err_str(unit)
   1088 other_str = get_err_str(other)
-> 1090 raise UnitConversionError(f"{unit_str} and {other_str} are not convertible")

UnitConversionError: 'eV' (energy/torque/work) and 'K' (temperature) are not convertible

To handle non-standard unit conversions, astropy.units allows the use of equivalencies. The conversion from eV to K can be done by using the temperature_energy() equivalency.

(1 * u.eV).to("K", equivalencies=u.temperature_energy())

$11604.518 \; \mathrm{K}$

Radians are treated dimensionlessly when the dimensionless_angles() equivalency is in effect. Note that this equivalency does not account for the multiplicative factor of \(2π\) that is used when converting between frequency and angular frequency.

(3.2 * u.rad / u.s).to("1 / s", equivalencies=u.dimensionless_angles())

$3.2 \; \mathrm{\frac{1}{s}}$

Physical constants

few can use astropy.constants to access the most commonly needed physical constants.

print(constants.c)

  Name   = Speed of light in vacuum
  Value  = 299792458.0
  Uncertainty  = 0.0
  Unit  = m / s
  Reference = CODATA 2018

A Constant behaves very similarly to a Quantity. For example, we can use the Boltzmann constant to mimic the behavior of u.temperature_energy().

thermal_energy_per_particle = 0.6 * u.keV
temperature = thermal_energy_per_particle / constants.k_B
print(temperature.to("MK"))

6.962710872930049 MK

Electromagnetic constants often need the unit system to be specified, or will result in an exception.

2 * constants.e

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[27], line 1
----> 1 2 * constants.e

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/astropy/constants/constant.py:49, in ConstantMeta.__new__.<locals>.wrap.<locals>.wrapper(self, *args, **kwargs)
     47 if not self.system and name_lower in self._has_incompatible_units:
     48     systems = sorted(x for x in instances if x)
---> 49     raise TypeError(
     50         f"Constant {self.abbrev!r} does not have physically compatible "
     51         "units across all systems of units and cannot be "
     52         "combined with other values without specifying a "
     53         f"system (eg. {self.abbrev}.{systems[0]})"
     54     )
     56 return meth(self, *args, **kwargs)

TypeError: Constant 'e' does not have physically compatible units across all systems of units and cannot be combined with other values without specifying a system (eg. e.emu)

2 * constants.e.si

$3.2043533 \times 10^{-19} \; \mathrm{C}$

Code within PlasmaPy generally uses SI units.

Plotting quantities

Astropy has built-in support for plotting Quantity objects. Let’s plot the number density of electrons in the solar wind using an empirical formula given by Kruparova et al. (2023), which has a range of validity from 13 to 50 solar radii.

radii = np.linspace(13, 50, num=50) * constants.R_sun

Next we can apply the formula to get the electron density:

\[n_e(R) = \left( 343466\ \mbox{cm}^{-3} \right) × \left( \frac{R}{R_☉} \right)^{-1.87}\]

n_e = 343_466 * u.cm**-3 * (radii / constants.R_sun) ** -1.87

Will make use make use of astropy.visualization.quantity_support. This is a context manager, which means that we use the with statement.

with quantity_support():
    plt.figure()
    plt.plot(radii.to(u.R_sun), n_e)
    plt.draw()

Optimizing unit operations (if time)

Astropy’s documentation includes performance tips for using astropy.units in computationally intensive situations. We can test it with %timeit, which runs a command repeatedly to see how long it takes.

Putting compound units in parentheses speeds things up by reducing the number of copies made by the operation.

%timeit 1.6 * u.barn * u.Mpc
%timeit 1.6 * (u.barn * u.Mpc)

25 μs ± 159 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
16.3 μs ± 123 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

We can assign a unit to a value using the Quantity class directly.

%timeit u.Quantity(1.6, u.barn * u.Mpc)

13.6 μs ± 29.5 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

What else can we do to save time from the above operation?

volume_unit = u.barn * u.Mpc  # pre-define the unit
%timeit u.Quantity(1.6, volume_unit)

3.78 μs ± 14 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Physical types (if time)

A physical type corresponds to physical quantities with dimensionally compatible units. Astropy has functionality that represents different physical types. These physical type objects can be accessed using either the physical_type attribute of a unit or get_physical_type().

(u.m**2 / u.s).physical_type

PhysicalType({'diffusivity', 'kinematic viscosity'})

u.get_physical_type("number density")

PhysicalType('number density')

These physical type objects can be used for dimensional analysis.

energy_density = (u.J * u.m**-3).physical_type
velocity = u.get_physical_type("velocity")
print(energy_density * velocity)

energy flux/irradiance

There are some pretty obscure physical types too!

u.get_physical_type(u.m * u.s)

PhysicalType('absement')