📏🏋️☕️ Using Units¶

This page shows how to use astropy.units and astropy.constants to attach physical units to numbers and arrays. When doing calculations, numbers are just numbers, and Python will happily calculate whatever we tell it to, whether that calculation makes conceptual sense or not. To help remedy this problem, the community-developed astropy package provides a snazzy toolkit for attaching units to quantities, handling conversions between units, and making our code complain if we try to do something nonesensical.

What's the problem?¶

Let's say we have a calculation where we need to add two lengths together: $2\rm{km} + 3\rm{cm}$. In Python, if we're not paying attention, we might accidentally write this as...

In [1]:

a = 2
b = 3
a + b

a = 2 b = 3 a + b

Out[1]:

5

...and get an answer that doesn't make any sense whatsoever. We have so much else to think about when coding, so sometimes it'd be nice if we didn't have to expend too much mental energy on converting units between variables.

How can we do unit calculations more carefully?¶

The solution to this problem is to let variables have units associated with them and let Python do any necessary unit conversions for us. If we import some unit tools from astropy into the shortcut u, we can do this.

In [2]:

import astropy.units as u

import astropy.units as u

In [3]:

a = 2*u.km
b = 3*u.cm
a+b

a = 2*u.km b = 3*u.cm a+b

Out[3]:

$2.00003 \; \mathrm{km}$

Hooray! The unit conversion happened automatically and gave us a sensible answer. The result has a unit attached to it so we know what it is, and we can convert it to any other distance unit we want!

In [4]:

c = a+b

c = a+b

In [5]:

c.to('m')

c.to('m')

Out[5]:

$2000.03 \; \mathrm{m}$

In [6]:

c.to('au')

c.to('au')

Out[6]:

$1.3369375 \times 10^{-8} \; \mathrm{AU}$

In [7]:

c.to('lightyear')

c.to('lightyear')

Out[7]:

$2.1140334 \times 10^{-13} \; \mathrm{lyr}$

Let's try a slightly more complicated calculation, figuring out the gravitational acceleration near Earth's surface.

In [8]:

# define the inputs to our calculation
G = 6.7e-11*u.m**3/u.kg/u.s**2
M = 6.0e27*u.g
R = 6.4e3*u.km

# do the calculation
g = G*M/R**2

# define the inputs to our calculation G = 6.7e-11*u.m**3/u.kg/u.s**2 M = 6.0e27*u.g R = 6.4e3*u.km # do the calculation g = G*M/R**2

In [9]:

g

g

Out[9]:

$9.8144531 \times 10^{9} \; \mathrm{\frac{m^{3}\,g}{kg\,km^{2}\,s^{2}}}$

Ug! That seems to have the 9.8 part right, but the units are nasty. Let's ask Python to please convert that into more familiar and useful units.

In [10]:

g.to('m/s**2')

g.to('m/s**2')

Out[10]:

$9.8144531 \; \mathrm{\frac{m}{s^{2}}}$

If we don't know what the final units of our result may be, we can ask simply to "decompose" the units into something simpler. In this example, it gives us the same $\rm{m/s^2}$ we were hoping for.

In [11]:

g.decompose()

g.decompose()

Out[11]:

$9.8144531 \; \mathrm{\frac{m}{s^{2}}}$

Discuss!¶

What happens if we try to do an impossible unit conversion using astropy.units (like in the code line below)? Is that behavior helpful or annoying?

In [12]:

g.to('kg')

g.to('kg')

---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
Cell In[12], line 1
----> 1 g.to('kg')

File ~/opt/anaconda3/envs/2024/lib/python3.12/site-packages/astropy/units/quantity.py:943, in Quantity.to(self, unit, equivalencies, copy)
    939 unit = Unit(unit)
    940 if copy:
    941     # Avoid using to_value to ensure that we make a copy. We also
    942     # don't want to slow down this method (esp. the scalar case).
--> 943     value = self._to_value(unit, equivalencies)
    944 else:
    945     # to_value only copies if necessary
    946     value = self.to_value(unit, equivalencies)

File ~/opt/anaconda3/envs/2024/lib/python3.12/site-packages/astropy/units/quantity.py:896, in Quantity._to_value(self, unit, equivalencies)
    893     equivalencies = self._equivalencies
    894 if not self.dtype.names or isinstance(self.unit, StructuredUnit):
    895     # Standard path, let unit to do work.
--> 896     return self.unit.to(
    897         unit, self.view(np.ndarray), equivalencies=equivalencies
    898     )
    900 else:
    901     # The .to() method of a simple unit cannot convert a structured
    902     # dtype, so we work around it, by recursing.
    903     # TODO: deprecate this?
    904     # Convert simple to Structured on initialization?
    905     result = np.empty_like(self.view(np.ndarray))

File ~/opt/anaconda3/envs/2024/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 ~/opt/anaconda3/envs/2024/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 ~/opt/anaconda3/envs/2024/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 ~/opt/anaconda3/envs/2024/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: 'm3 g / (kg km2 s2)' (acceleration) and 'kg' (mass) are not convertible

How do we include famous physical constants?¶

We use a lot of physical constants in astronomy, many with very strange units. Fortunately, astropy provides a database of many common astronomical constants that might be useful for our calculations. Let's repeat the above, using the astropy.constants module.

In [13]:

import astropy.constants as con

import astropy.constants as con

In [14]:

# define the inputs to our calculation
G = con.G
M = 1*u.M_earth
R = 1*u.R_earth

# do the calculation
g = G*M/R**2

# define the inputs to our calculation G = con.G M = 1*u.M_earth R = 1*u.R_earth # do the calculation g = G*M/R**2

In [15]:

g

g

Out[15]:

$6.6743 \times 10^{-11} \; \mathrm{\frac{m^{3}\,M_{\oplus}}{kg\,R_{\oplus}^{2}\,s^{2}}}$

In [16]:

g.decompose()

g.decompose()

Out[16]:

$9.7983981 \; \mathrm{\frac{m}{s^{2}}}$

That's more accurate and precise than our previous calculation, in which we provided each number just to 2 significant digits. In this particular example we used G from constants and Earth's mass and radius from units, but know that some quantities (like Earth's mass and radius) will appear in both units and constants.

In [17]:

print(con.G)

print(con.G)

  Name   = Gravitational constant
  Value  = 6.6743e-11
  Uncertainty  = 1.5e-15
  Unit  = m3 / (kg s2)
  Reference = CODATA 2018

Discuss!¶

Can we construct arrays of numbers that have units associated with them (for example, a 1D array of times, a 2D array of ADU per pixel, or a 3D array of density in some volume)? In what contexts might that be useful?

Ack! Can I get rid of the units?!?¶

We might find we want to do a calculation with units and constants but then convert our result back to a simple number or array. Some numpy, scipy, or plt functions don't play nicely with units, so we might want to drop the units from a quantity as some point. We can use the .to_value() method to return just the value of the quantity, in whatever units we specify.

In [18]:

g.to_value('m/s**2')

g.to_value('m/s**2')

Out[18]:

9.798398133669465

In [19]:

g.to_value(u.m/u.s**2)

g.to_value(u.m/u.s**2)

Out[19]:

9.798398133669465

We can call .to_value() without any arguments to return the value in its current unit, but be aware that current unit might be something confusing! For example, if we look at the undecomposed value for g, it will give us a value in strange units involved Earth radii and masses.

In [20]:

g.to_value()

g.to_value()

Out[20]:

6.6743e-11

In [21]:

g.decompose().to_value()

g.decompose().to_value()

Out[21]:

9.798398133669465

How can you apply this yourself?¶

Take what you've learned above and try out the following questions, using astropy units and/or constants in your calculations. Remember that in jupyter you can start typing a variable name and then hit <tab> to try to automatically complete the rest of its name.

1. If we drive a car continuously at a speed of 100 $\rm{km/hr}$ for 0.5 $\rm{fortnights}$, how many $\rm{lightseconds}$ have we traveled?
2. Calculate the energy $E$ of a photon with a wavelength $\lambda = 500 \rm{nm}$, using the equation $E = h\nu = hc/\lambda$. Print your result in units of Joules, ergs, electron volts, and kilwowatt hours.

The documentation for astropy.units and astropy.constants is extensive and helpful. You can do much more than the basics shown here; search for the documentation online to figure out what you need!