Units and Quantities

Objectives

  • Use units
  • Create functions that accept quantities as arguments
  • Create new units

Basics

How do we define a Quantity and which parts does it have?

In [1]:
from astropy import units as u
In [2]:
# Define a quantity length
length = 26.2 * u.meter
# print it
print(length) # length is a quantity
26.2 m
In [3]:
# Type of quantity
type(length)
Out[3]:
astropy.units.quantity.Quantity
In [4]:
# Type of unit
type(u.meter)
Out[4]:
astropy.units.core.IrreducibleUnit
In [5]:
# Quantity
length
Out[5]:
$26.2 \; \mathrm{m}$
In [6]:
# value
length.value
Out[6]:
26.2
In [7]:
# unit
length.unit
Out[7]:
$\mathrm{m}$
In [8]:
# information
length.info
Out[8]:
dtype = float64
unit = m
class = Quantity
n_bad = 0

Quantities can be converted to other units systems or factors by using to()

In [9]:
# Convert it to: km, lyr
print(length.to(u.km))
print(length.to(u.lightyear))
0.0262 km
2.76934218514e-15 lyr

We can do arithmetic operations when the quantities have the compatible units:

In [10]:
# arithmetic with distances
distance_start = 10 * u.mm
distance_end = 23 * u.km
length = distance_end - distance_start
print(length)
22.99999 km

Quantities can also be combined, for example to measure speed

In [11]:
# calculate a speed
time = 15 * u.minute
speed = length / time
print(speed)
1.53333266667 km / min
In [12]:
# decompose it
print(speed.decompose())
print(speed.si)
25.5555444444 m / s
25.5555444444 m / s

Challenges

  1. Convert the speed in imperial units (miles/hour) using:
    ```from astropy.units import imperial```
  2. Calculate whether a pint is more than half litre
    You can compare quantities as comparing variables.
    Something strange? Check what deffinition of pint astropy is using.
  3. Does units work with areas? calculate the area of a rectangle of 3 km of side and 5 meter of width. Show them in m^2 and convert them to yards^2
In [13]:
#1
from astropy.units import imperial
print(speed.to(imperial.mile/u.hour))
57.166124831 mi / h
In [14]:
#2
imperial.pint > 0.5 * u.l
 # A liquid pint in US is 473 ml; in UK is 568 ml
Out[14]:
False
In [15]:
#3
rectangle_area = 3 * u.km * 5 * u.m
print(rectangle_area)
print(rectangle_area.decompose())
print(rectangle_area.to(imperial.yard ** 2))
15.0 km m
15000.0 m2
17939.8506945 yd2

Composed units

Many units are compositions of others, for example, one could create new combinationes for ease of use:

In [16]:
# create a composite unit
cms = u.cm / u.s
speed.to(cms)
Out[16]:
$2555.5544 \; \mathrm{\frac{cm}{s}}$
In [17]:
# and in the imperial system
mph = imperial.mile / u.hour
speed.to(mph)
Out[17]:
$57.166125 \; \mathrm{\frac{mi}{h}}$

and others are already a composition:

In [18]:
# what can be converted from s-1?
(u.s ** -1).compose()
Out[18]:
[Unit("Hz"), Unit("Bq"), Unit("3.7e+10 Ci")]
In [19]:
# or Jules?
(u.joule).compose()
Out[19]:
[Unit("J"), Unit("1e+07 erg"), Unit("4.58742e+17 Ry"), Unit("6.24151e+18 eV")]
In [20]:
# Unity of R
(13.605692 * u.eV).to(u.Ry)
Out[20]:
$0.99999996 \; \mathrm{R_{\infty}}$

Sometime we get no units quantitites

In [21]:
# no units
nounits = 20. * u.cm / (1. * u.m)
nounits
Out[21]:
$20 \; \mathrm{\frac{cm}{m}}$

What happen if we add a number to this?

In [22]:
# arithmetic with no units
nounits + 3
Out[22]:
$3.2 \; \mathrm{}$
In [23]:
# final value of a no unit quantity
nounits.decompose() # It's a unitless quantity
Out[23]:
$0.2 \; \mathrm{}$

Equivalencies

Some conversions are not done by a conversion factor as between miles and kilometers, for example converting between wavelength and frequency.

In [24]:
# converting spectral quantities
(656.281 * u.nm).to(u.Hz) # Fails because they are not compatible
---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
<ipython-input-24-1d7f56a8e756> in <module>()
      1 # converting spectral quantities
----> 2 (656.281 * u.nm).to(u.Hz) # Fails because they are not compatible

/home/dvd/.conda/envs/swc/lib/python2.7/site-packages/astropy/units/quantity.pyc in to(self, unit, equivalencies)
    631         unit = Unit(unit)
    632         new_val = self.unit.to(unit, self.view(np.ndarray),
--> 633                                equivalencies=equivalencies)
    634         return self._new_view(new_val, unit)
    635 

/home/dvd/.conda/envs/swc/lib/python2.7/site-packages/astropy/units/core.pyc in to(self, other, value, equivalencies)
    966             If units are inconsistent
    967         """
--> 968         return self._get_converter(other, equivalencies=equivalencies)(value)
    969 
    970     def in_units(self, other, value=1.0, equivalencies=[]):

/home/dvd/.conda/envs/swc/lib/python2.7/site-packages/astropy/units/core.pyc in _get_converter(self, other, equivalencies)
    867         except UnitsError:
    868             return self._apply_equivalencies(
--> 869                 self, other, self._normalize_equivalencies(equivalencies))
    870         return lambda val: scale * _condition_arg(val)
    871 

/home/dvd/.conda/envs/swc/lib/python2.7/site-packages/astropy/units/core.pyc in _apply_equivalencies(self, unit, other, equivalencies)
    858         raise UnitConversionError(
    859             "{0} and {1} are not convertible".format(
--> 860                 unit_str, other_str))
    861 
    862     def _get_converter(self, other, equivalencies=[]):

UnitConversionError: 'nm' (length) and 'Hz' (frequency) are not convertible
In [25]:
# but doing it right
(656.281 * u.nm).to(u.Hz, equivalencies=u.spectral())
Out[25]:
$4.5680502 \times 10^{14} \; \mathrm{Hz}$

Other built-in equivalencies are:

  • parallax()
  • Doppler (dopplr_radio, doppler_optical, doppler_relativistic)
  • spectral flux density
  • brigthness temperature
  • temperature energy
  • and you can build your own
In [26]:
# finding the equivalencies
u.Hz.find_equivalent_units()
Out[26]:
  Primary name | Unit definition | Aliases     
[
  Bq           | 1 / s           | becquerel    ,
  Ci           | 2.7027e-11 / s  | curie        ,
  Hz           | 1 / s           | Hertz, hertz ,
]
In [27]:
# but also using other systems
u.Hz.find_equivalent_units(equivalencies=u.spectral())
Out[27]:
  Primary name | Unit definition        | Aliases              
[
  AU           | 1.49598e+11 m          | au, astronomical_unit ,
  Angstrom     | 1e-10 m                | AA, angstrom          ,
  Bq           | 1 / s                  | becquerel             ,
  Ci           | 2.7027e-11 / s         | curie                 ,
  Hz           | 1 / s                  | Hertz, hertz          ,
  J            | kg m2 / s2             | Joule, joule          ,
  Ry           | 2.17987e-18 kg m2 / s2 | rydberg               ,
  cm           | 0.01 m                 | centimeter            ,
  eV           | 1.60218e-19 kg m2 / s2 | electronvolt          ,
  erg          | 1e-07 kg m2 / s2       |                       ,
  k            | 100 / m                | Kayser, kayser        ,
  lyr          | 9.46073e+15 m          | lightyear             ,
  m            | irreducible            | meter                 ,
  micron       | 1e-06 m                |                       ,
  pc           | 3.08568e+16 m          | parsec                ,
  solRad       | 6.95508e+08 m          | R_sun, Rsun           ,
]

Printing the quantities

In [28]:
# Printing values with different formats
print("{0.value:0.03f} {0.unit:FITS}".format(speed))
print("{0.value:0.03f} {0.unit:latex_inline}".format(speed))
1.533 km min-1
1.533 $\mathrm{km\,min^{-1}}$

Arrays

Quantities can also be applied to arrays

In [29]:
# different ways of defining a quantity for a single value
length = 44 * u.m
time = u.Quantity(23, u.s)
speed = length / time
speed
Out[29]:
$1.9130435 \; \mathrm{\frac{m}{s}}$
In [30]:
# now with lists
length_list = [1, 2, 3] * u.m

# and arrays
import numpy as np
time_array = np.array([1, 2, 3]) * u.s

# and its arithmetics
length_list / time_array
Out[30]:
$[1,~1,~1] \; \mathrm{\frac{m}{s}}$
In [31]:
# angles are smart!
angle = u.Quantity(np.arange(180), u.deg)
print(angle[[0, -1]])
print(np.sin(angle[[0, -1]]))
[   0.  179.] deg
[ 0.          0.01745241]

Plotting quantities

To work nicely with matplotlib we need to do as follows:

In [32]:
# allowing for plotting
from astropy.visualization import quantity_support
quantity_support()

# loading matplotlib
%matplotlib inline
from matplotlib import pyplot as plt
In [33]:
# Ploting the previous array
plt.plot(angle, np.sin(angle))
Out[33]:
[<matplotlib.lines.Line2D at 0x7f5ec58923d0>]

Creating functions with quantities as units

We want to have functions that contain the information of the untis, and with them we can be sure that we will be always have the right result.

In [34]:
# Create a function for the Kinetic energy
@u.quantity_input(mass=u.kg, speed=u.m/u.s)
def kinetic(mass, speed):
    return (mass * speed ** 2 / 2.).to(u.joule)
In [35]:
# run with and without units
kinetic(5, 10) # Fails! it doesn't specify the units.
kinetic(5 * u.kg, 100 * cms)
Out[35]:
$2.5 \; \mathrm{J}$

Challenges

  1. Create a function that calculates potential energy where *g* defaults to Earth value, but could be used for different planets. Test it for any of the *g* values for any other planet.
In [36]:
#4
@u.quantity_input(mass=u.kg, height=u.m, g=u.m / u.s ** 2)
def potential(mass, height, g=9.8 * u.m / u.s **2):
    return (0.5 * mass * g * height).to(u.joule)
In [37]:
# run it for some values
potential(5 * u.kg, 30 *u.cm)
Out[37]:
$7.35 \; \mathrm{J}$
In [38]:
# on Mars:
potential(5 * u.kg, 1 * u.m, g=3.75 * u.m/u.s**2)
Out[38]:
$9.375 \; \mathrm{J}$

Create your own units

Some times we want to create our own units:

In [39]:
# Create units for a laugh scale
titter = u.def_unit('titter')
chuckle = u.def_unit('chuckle', 5 * titter)
laugh = u.def_unit('laugh', 4 * chuckle)
guffaw = u.def_unit('guffaw', 3 * laugh)
rofl = u.def_unit('rofl', 4 * guffaw)
death_by_laughing = u.def_unit('death_by_laughing', 10 * rofl)
print((1 * rofl).to(titter))
240.0 titter

Challenges

  1. Convert the area calculated before `rectangle_area` in Hectare (1 hectare = 100 ares; 1 are = 100 m2).
In [40]:
#5
ares = u.def_unit('ares', (10 * u.m)**2)
hectar = u.def_unit('hectares', 100 * ares)
print(rectangle_area.to(hectar))
1.5 hectares