One of the most powerful, central features of Astropy is the unit conversion library, and the Quantity
object that associates numerical values (and arrays) with units. Effort is made to integrate units throughout Astropy so it is important to understand their usage.
Since we may want to use a number of units in expressions, it is easiest and most concise to import the units module with:
from astropy import units as u
though note that this will conflict with any variable called u
.
Units can then be accessed with:
u.m
In the notebook this just displays a LaTeX representation of the unit:
u.Angstrom
They all have a docstring defining them:
print(u.m.__doc__)
meter: base unit of length in SI
An associated physical type:
print(u.m.physical_type)
length
print(u.Angstrom.physical_type)
length
And other various associated information that can be explored through the full documentation.
To create a quantity we just multiply a value by a unit. The value can be either a scalar or a Numpy array (or an object that can be converted to an array, like a Python list):
1.0 * u.m
[1.0, 2.0, 3.0] * u.m
<Quantity [ 1., 2., 3.] m>
Since u.m
and u.Angstrom
have the same physical type we can convert a quantity in meters to a quantity in Angstroms:
q = [1.0, 2.0, 3.0] * u.m
q.to(u.Angstrom)
<Quantity [ 1.00000000e+10, 2.00000000e+10, 3.00000000e+10] Angstrom>
Quantities can then be combined:
q1 = 3. * u.m
q1
q2 = 5. * u.cm / u.s / u.g**2
q2
q1 * q2
and converted to different units:
(q1 * q2).to(u.m**2 / u.kg**2 / u.s)
We can also incorporate physical constants from the astropy.constants
library. For example, ff we want to compute
the Gravitational force felt by a 100. kg space probe by the Sun, at a distance of 3.2au, we can do:
from astropy.constants import G
print repr(G)
<Constant name=u'Gravitational constant' value=6.67384e-11 error=8e-15 units='m3 / (kg s2)' reference=u'CODATA 2010'>
F = (G * 1. * u.M_sun * 100. * u.kg) / (3.2 * u.au)**2
F.unit.find_equivalent_units()
Primary name | Unit definition | Aliases [ N | kg m / s2 | Newton, newton , dyn | 1e-05 kg m / s2 | dyne , ]
F.to(u.N)
Units don't necessarily have to have the same physical type to be converted between each other. Astropy includes a number of unit equivalencies that can be enabled (with the possibility of defining your own). For example, without equivalencies the following fails:
(450. * u.nm).to(u.GHz)
--------------------------------------------------------------------------- UnitsError Traceback (most recent call last) <ipython-input-216-e4fbcb033257> in <module>() ----> 1 (450. * u.nm).to(u.GHz) /usr/local/lib/python2.7/dist-packages/astropy/units/quantity.pyc in to(self, unit, equivalencies) 575 unit = Unit(unit) 576 new_val = np.asarray( --> 577 self.unit.to(unit, self.value, equivalencies=equivalencies)) 578 return self._new_view(new_val, unit) 579 /usr/local/lib/python2.7/dist-packages/astropy/units/core.pyc in to(self, other, value, equivalencies) 910 If units are inconsistent 911 """ --> 912 return self.get_converter(other, equivalencies=equivalencies)(value) 913 914 def in_units(self, other, value=1.0, equivalencies=[]): /usr/local/lib/python2.7/dist-packages/astropy/units/core.pyc in get_converter(self, other, equivalencies) 846 except UnitsError: 847 return self._apply_equivalences( --> 848 self, other, self._normalize_equivalencies(equivalencies)) 849 return lambda val: scale * _condition_arg(val) 850 /usr/local/lib/python2.7/dist-packages/astropy/units/core.pyc in _apply_equivalences(self, unit, other, equivalencies) 809 raise UnitsError( 810 "{0} and {1} are not convertible".format( --> 811 unit_str, other_str)) 812 813 def get_converter(self, other, equivalencies=[]): UnitsError: 'nm' (length) and 'GHz' (frequency) are not convertible
But with the spectral equivalency:
(450. * u.nm).to(u.GHz, equivalencies=u.spectral())
(450. * u.eV).to(u.nm, equivalencies=u.spectral())
Most of the Numpy functions also understand Astropy quantities (with some notable exceptions you can find in the documentation). For example, by default np.sin
takes input in radians:
import numpy as np
np.sin(90)
0.89399666360055785
But if we pass in a quantity in degrees:
np.sin(90 * u.degree)
More examples and exercises with units in this notebook.
The astropy Table class provides an extension of NumPy structured arrays for storing and manipulating heterogeneous tables of data. A few notable features of this package are:
from astropy.table import Table
import numpy as np
There is great deal of flexibility in the way that a table can be initially constructed:
Read an existing table from a file or web URL
Add columns of data one by one
Add rows of data one by one
From an existing data structure in memory:
See the documentation section on Constructing a table for the gory details and plenty of examples.
The easiest (though maybe not always most useful in practice) way is to just create an empty table object and add columns to it by assigning arrays to column names, in a syntax reminiscent on dictionaries:
t = Table()
t['name'] = ['source 1', 'source 2', 'source 3', 'source 4']
t['flux'] = [1.2, 2.2, 3.1, 4.3]
Tables have numerous useful attributes for introspection. Just a couple examples are:
t.colnames
['name', 'flux']
t.dtype
dtype([('name', '|S8'), ('flux', '<f8')])
There are several options for viewing tables, but for a small table it will just print in full:
t
name | flux |
---|---|
source 1 | 1.2 |
source 2 | 2.2 |
source 3 | 3.1 |
source 4 | 4.3 |
t.show_in_browser()
<open file '<fdopen>', mode 'w+b' at 0x529b300>
We can access single rows of columns, single columns, or selection sof multiple rows/columns:
t['flux']
<Column name='flux' unit=None format=None description=None> array([ 1.2, 2.2, 3.1, 4.3])
t[1]
<Row 1 of table values=('source 2', 2.2) dtype=[('name', '|S8'), ('flux', '<f8')]>
t[1:3]
name | flux |
---|---|
source 2 | 2.2 |
source 3 | 3.1 |
Table columns can also have units associated with them for storing quantities to tables. This is useful, for example, when we read FITS tables--if a column has an associated TUNIT
keyword we can return dimensionful quantities from the table.
Masked tables are another powerful feature of the Astropy Table class:
t2 = Table([['x', 'y', 'z'],
[1.1, 2.2, 3.3]],
names=['name', 'value'],
masked=True)
t2
name | value |
---|---|
x | 1.1 |
y | 2.2 |
z | 3.3 |
t2['value'].mask = [False, True, False]
t2
name | value |
---|---|
x | 1.1 |
y | -- |
z | 3.3 |
Numerical operations on masked columns ignore the masked values:
t2['value'].mean()
2.2000000000000002
t2['value'].fill_value
1e+20
We can set our own fill value which can be used when outputting the table to a format that does not natively support masked values:
t2['value'].fill_value = -99
t2.filled()
name | value |
---|---|
x | 1.1 |
y | -99.0 |
z | 3.3 |
Tables support a number of high-level operations for combining tables, including merging, concatenating, and grouping. For example, recall our first table of flux values:
t
name | flux |
---|---|
source 1 | 1.2 |
source 2 | 2.2 |
source 3 | 3.1 |
source 4 | 4.3 |
Now say that we now got some additional flux values from a different reference for a different, but overlapping sample of sources:
t2 = Table()
t2['name'] = ['source 1', 'source 3', 'source 8']
t2['flux2'] = [1.4, 3.5, 8.6]
Now we can get a master table of flux measurements which are joined matching the values the name
column. This includes every row from each of the two tables, which is known as an outer
join:
from astropy.table import join
t3 = join(t, t2, keys=['name'], join_type='outer')
t3
name | flux | flux2 |
---|---|---|
source 1 | 1.2 | 1.4 |
source 2 | 2.2 | -- |
source 3 | 3.1 | 3.5 |
source 4 | 4.3 | -- |
source 8 | -- | 8.6 |
There are several other join options. For example we could choose to keep only rows where both tables had a valid measurement using an inner
join:
t3 = join(t, t2, keys=['name'], join_type='inner')
t3
name | flux | flux2 |
---|---|---|
source 1 | 1.2 | 1.4 |
source 3 | 3.1 | 3.5 |
Tables can be written to and read in from a variety of data formats:
t3.write('example.fits', overwrite=True)
Tables can be read back in using the Table.read
interface, and then written to another format (this is an easy way to do format conversion):
t4 = Table.read('example.fits')
t4
name | flux | flux2 |
---|---|---|
source 1 | 1.2 | 1.4 |
source 3 | 3.1 | 3.5 |
t4.write('example.vot', format='votable', overwrite=True)
!cat example.vot
<?xml version="1.0" encoding="utf-8"?> <!-- Produced with astropy.io.votable version 0.4.2 http://www.astropy.org/ --> <VOTABLE version="1.2" xmlns="http://www.ivoa.net/xml/VOTable/v1.2" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://www.ivoa.net/xml/VOTable/v1.2"> <RESOURCE type="results"> <TABLE> <FIELD ID="name" arraysize="8" datatype="char" name="name"/> <FIELD ID="flux" datatype="double" name="flux"/> <FIELD ID="flux2" datatype="double" name="flux2"/> <DATA> <TABLEDATA> <TR> <TD>source 1</TD> <TD>1.2</TD> <TD>1.3999999999999999</TD> </TR> <TR> <TD>source 3</TD> <TD>3.1000000000000001</TD> <TD>3.5</TD> </TR> </TABLEDATA> </DATA> </TABLE> </RESOURCE> </VOTABLE>
Tables can also be written to a variety of ASCII formats, including custom formats. Before trying to hack together a hand-written table reader/writer look into what is possible with the astropy.io.ascii
.
There is also an API for developing new table readers/writers for other formats, though the built-in formats cover much of the existing data out there.
t = Table()
t['name'] = ['source 1', 'source 2', 'source 3', 'source 4']
t['flux'] = [1.2, 2.2, 3.1, 4.3]
t
name | flux |
---|---|
source 1 | 1.2 |
source 2 | 2.2 |
source 3 | 3.1 |
source 4 | 4.3 |
t.write('mytable.tex', format='latex')
t.write(sys.stdout, format='latex')
\begin{table} \begin{tabular}{cc} name & flux \\ source 1 & 1.2 \\ source 2 & 2.2 \\ source 3 & 3.1 \\ source 4 & 4.3 \\ \end{tabular} \end{table}
text = """
\begin{table}
\begin{tabular}{cc}
name & flux \\
source 1 & 1.2 \\
source 2 & 2.2 \\
source 3 & 3.1 \\
source 4 & 4.3 \\
\end{tabular}
\end{table}
"""
Table.read('mytable.tex', format='latex')
name | flux |
---|---|
source 1 | 1.2 |
source 2 | 2.2 |
source 3 | 3.1 |
source 4 | 4.3 |
import sys
sys.stdout
<IPython.kernel.zmq.iostream.OutStream at 0x2cc3e10>
We briefly saw how to save and load binary tables from FITS files, but that was through a high-level interface that does not give us a great deal of access to the lower-level details of FITS files. We would also in many cases like to be able to manipulate the contents of FITS headers and the arrangement of HDUs.
And although there is a high-level generic interface for tabular data in the Table
class, the equivalent for gridded data like images is still under development (see: NDData
).
The full interface to FITS images that allows direct access to their image arrays and headers is via the astropy.io.fits
package. If you have ever used the venerable PyFITS library, then you have already used astropy.io.fits
, which is just an import of PyFITS into Astropy. By moving PyFITS into Astropy it is able to evolve new capabilities, such as the improved table interface, and improved support for units.
To open a FITS file import the fits
package and use the fits.open
function:
from astropy.io import fits
# The following FITS file lives in the lesson repository, relative
# to this file. However, if you get a "file not found" error you can download and open the file
# by uncommenting these lines (and commenting the last line)
# from astropy.utils import download_file
# filename = download_file('https://github.com/abostroem/2015-01-03-aas/raw/gh-pages/intermediate/astropy/data/gll_iem_v02_P6_V11_DIFFUSE.fit')
# hdulist = fits.open(filename)
hdulist = fits.open('data/gll_iem_v02_P6_V11_DIFFUSE.fit')
The returned object, hdulist
, (an instance of the HDUList class) behaves like a Python list, and each element
maps to a Header-Data Unit (HDU) in the FITS file. You can view more
information about the FITS file with:
hdulist.info()
Filename: data/gll_iem_v02_P6_V11_DIFFUSE.fit No. Name Type Cards Dimensions Format 0 PRIMARY PrimaryHDU 34 (720, 360, 30) float32 1 ENERGIES BinTableHDU 19 30R x 1C [D]
As we can see, this file contains two HDUs. The first contains an image cube, the second a table with one column of double-precision floating point values. To access the primary HDU, which contains the main data, you can then do:
primary_hdu = hdulist[0]
primary_hdu = hdulist['PRIMARY']
To access the header use the .header
attribute:
primary_hdu.header
SIMPLE = T / Written by IDL: Thu Jan 20 07:19:05 2011 BITPIX = -32 / NAXIS = 3 / number of data axes NAXIS1 = 720 / length of data axis 1 NAXIS2 = 360 / length of data axis 2 NAXIS3 = 30 / length of data axis 3 EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H FLUX = 8.42259635886 / CRVAL1 = 0. / Value of longitude in pixel CRPIX1 CDELT1 = 0.5 / Step size in longitude CRPIX1 = 360.5 / Pixel that has value CRVAL1 CTYPE1 = 'GLON-CAR' / The type of parameter 1 (Galactic longitude in CUNIT1 = 'deg ' / The unit of parameter 1 CRVAL2 = 0. / Value of latitude in pixel CRPIX2 CDELT2 = 0.5 / Step size in latitude CRPIX2 = 180.5 / Pixel that has value CRVAL2 CTYPE2 = 'GLAT-CAR' / The type of parameter 2 (Galactic latitude in C CUNIT2 = 'deg ' / The unit of parameter 2 CRVAL3 = 50. / Energy of pixel CRPIX3 CDELT3 = 0.113828620540137 / log10 of step size in energy (if it is logarith CRPIX3 = 1. / Pixel that has value CRVAL3 CTYPE3 = 'photon energy' / Axis 3 is the spectra CUNIT3 = 'MeV ' / The unit of axis 3 CHECKSUM= '3fdO3caL3caL3caL' / HDU checksum updated 2009-07-07T22:31:18 DATASUM = '2184619035' / data unit checksum updated 2009-07-07T22:31:18 DATE = '2009-07-07' / FILENAME= '$TEMPDIR/diffuse/gll_iem_v02.fit' /File name with version number TELESCOP= 'GLAST ' / INSTRUME= 'LAT ' / ORIGIN = 'LISOC ' /LAT team product delivered from the LISOC OBSERVER= 'MICHELSON' /Instrument PI HISTORY Scaled version of gll_iem_v02.fit for use with P6_V11_DIFFUSE
The Header
object works very similarly to the Python dict
type--the values for individual keywords can be looked up like so:
primary_hdu.header['OBSERVER']
'MICHELSON'
primary_hdu.header['DATE']
'2009-07-07'
The Header
interface has a wealth of other capabilities that can be read about in the full documentation.
The array data associated with an HDU can be accessed via the .data
attribute, which returns a view of the data as a Numpy array. In this example the data is an image cube for a series of wavelengths. To display the first one:
%matplotlib inline
from matplotlib.pyplot import imshow
imshow(primary_hdu.data[0], origin='lower')
<matplotlib.image.AxesImage at 0x66887d0>
from aplpy import FITSFigure
fig = FITSFigure('data/gll_iem_v02_P6_V11_DIFFUSE.fit', slices=[0])
fig.show_colorscale()
INFO:astropy:Auto-setting vmin to -1.114e-06 INFO:astropy:Auto-setting vmax to 1.337e-05
INFO: Auto-setting vmin to -1.114e-06 [aplpy.core] INFO: Auto-setting vmax to 1.337e-05 [aplpy.core]
We will use this for a brief demonstration of how to write a new FITS file. Say we wanted to save just this slice to a file. We can create a new primary HDU with this data like so:
new_hdu = fits.PrimaryHDU(data=primary_hdu.data[0], header=primary_hdu.header)
This copies over the header from the original file, but some of the keywords are no longer applicable, such as the WCS for the spectral axis, or the filename keyword:
del new_hdu.header['FILENAME']
We can use a wildcard expression to take care of the WCS:
del new_hdu.header['C*3']
Note, however, that we did not have to fuss around with keywords related to the basic structure of the file, such as BITPIX or the NAXIS keywords. It is not necessary to manipulate such basic structural keywords by hand--their consistency with the data is managed automatically (in the future there will be further enhancements to make WCS manipulation more transparent too):
new_hdu.header
SIMPLE = T / conforms to FITS standard BITPIX = -32 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 720 NAXIS2 = 360 COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H FLUX = 8.42259635886 / CRVAL1 = 0. / Value of longitude in pixel CRPIX1 CDELT1 = 0.5 / Step size in longitude CRPIX1 = 360.5 / Pixel that has value CRVAL1 CTYPE1 = 'GLON-CAR' / The type of parameter 1 (Galactic longitude in CUNIT1 = 'deg ' / The unit of parameter 1 CRVAL2 = 0. / Value of latitude in pixel CRPIX2 CDELT2 = 0.5 / Step size in latitude CRPIX2 = 180.5 / Pixel that has value CRVAL2 CTYPE2 = 'GLAT-CAR' / The type of parameter 2 (Galactic latitude in C CUNIT2 = 'deg ' / The unit of parameter 2 CHECKSUM= '3fdO3caL3caL3caL' / HDU checksum updated 2009-07-07T22:31:18 DATASUM = '2184619035' / data unit checksum updated 2009-07-07T22:31:18 DATE = '2009-07-07' / TELESCOP= 'GLAST ' / INSTRUME= 'LAT ' / ORIGIN = 'LISOC ' /LAT team product delivered from the LISOC OBSERVER= 'MICHELSON' /Instrument PI HISTORY Scaled version of gll_iem_v02.fit for use with P6_V11_DIFFUSE
Writing the file out with the checksum=True
option will also update the CHECKSUM
and DATASUM
keywords:
new_hdu.writeto('lat_background_model_slice.fits', checksum=True,
clobber=True)
WARNING: Overwriting existing file 'lat_background_model_slice.fits'. [astropy.io.fits.file] WARNING:astropy:Overwriting existing file 'lat_background_model_slice.fits'.
This is a good place for a brief mention of the coordinates
interface; specifically the SkyCoord
object, a high level interface for representation and manipulation of celestial coordinates. Astropy 0.4 only includes a few common coordinate systems (ICRS, FK4, FK5, and Galactic), but future versions will include more built-in coordinate systems, and users can already define their own systems.
Creating coordinate objects is straightforward. We use units--often u.deg
(a shorthand for u.degree
) to specify the units of the angles:
from astropy.coordinates import SkyCoord
from astropy import units as u
SkyCoord(ra=10.68458 * u.deg, dec=41.26917 * u.deg, frame='icrs')
<SkyCoord (ICRS): ra=10.68458 deg, dec=41.26917 deg>
In addition to passing angular Quantity objects, the SkyCoord class can parse coordinate strings (this can be used to pass initialize a SkyCoord
object from, say, a column of coordinates read from an ASCII table with the Table
interface):
c = SkyCoord('00h42m44.3s +41d16m9s', frame='icrs')
c
<SkyCoord (ICRS): ra=10.6845833333 deg, dec=41.2691666667 deg>
We can access the different components as Angle
objects and convert them to other units:
c.ra
c.ra.hour
0.7123055555555556
...and so on.
It is also possible to convert our SkyCoord
object to other coordinate frames:
c.galactic
<SkyCoord (Galactic): l=121.174305193 deg, b=-21.5728034839 deg>
c.transform_to('galactic')
<SkyCoord (Galactic): l=121.174305193 deg, b=-21.5728034839 deg>
Or we can create custom coordinate frames in order to adjust other parameters, such as the equinox:
from astropy.coordinates import FK5
FK5?
frame = FK5(equinox='J2010')
c.transform_to(frame)
<SkyCoord (FK5): equinox=J2010.000, ra=10.8218390257 deg, dec=41.3238612021 deg>