In this lecture we will introduce the astropy library and the affiliated package astroquery. The official documents of these packages are available at:
To install these packages with conda:
conda install astropy
conda install -c astropy astroquery
There are many more packages affiliated to astropy which can be of interest to you. You can find the list at:
FITS files are by far the favorite format to store and distribute astronomical data. They come in two flavors:
To read these two types of files we will need to import from astropy two different sub-libraries:
from astropy.io import fits
from astropy.table import Table
FITS files can store multi-dimensional data (commonly 2 or 3 dimensions). Any given FITS file can contain multiple images (or tables) called extensions. Every FITS extension contains a header and data. FITS headers can contain World Coordinate System (wcs) information that indicates where a given pixel is on the sky.
Unlike Python, the FITS convention is indexing starting at 1. Generally, astropy takes this into account.
Convenience functions make reading FITS images easy.
from astropy.io import fits
img1 = fits.getdata(filename) # Getting the image
head1 = fits.getheader(filename) # and the Header
This opens the image as a Numpy array, and the header as a “dictionary-like” object (i.e., you can access the individual header keywords through “head1[‘key’]” ).
To open other extensions in the fits file:
img1 = fits.getdata(filename, 0) # Primary Ext
img2 = fits.getdata(filename, 1) # Second Ext
img2 = fits.getdata(filename, ext=1) # Equivalent
It is possible to import a FITS file also using an URL. This is done with the download_file function.
In [19]:
from astropy.utils.data import download_file
from astropy.io import fits
image_file = download_file('http://data.astropy.org/tutorials/FITS-images/HorseHead.fits',
cache=True)
To have some information about the file we just open, we can use the fits.info method:
In [20]:
fits.info(image_file)
Then, we can decide to open an exstension, check its shape and eventually plot it with matplolib imshow.
In [21]:
image_data = fits.getdata(image_file, ext=0)
image_data.shape
Out[21]:
In [22]:
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure()
plt.imshow(image_data, cmap='gist_heat',origin='lower')
plt.colorbar();
In [23]:
hdulist = fits.open(image_file)
hdulist.info()
Next, we can extract header and data from the first extension. We can do this in two ways:
In [24]:
header = hdulist['PRIMARY'].header
data = hdulist['PRIMARY'].data
FITS files are read in such a way that the first axis (often the RA for astronomical images) is read in as the last axis in the numpy array. Be sure to double check that you have the axis you need.
We can, at this point, close the file.
In [25]:
hdulist.close()
Now, let's explore the header. To print in a human readable way, it's useful to use the repr function which adds line breaks in between the keywords:
In [26]:
print(repr(header[:10])) # Beginning of the header
We can access the list of keywords, values, a specific keyword or comment:
In [27]:
print (header[:10].keys())
print (header[:10].values())
print (header['ORIGIN'])
print (header.comments['ORIGIN'])
To extract the astrometry, we will use wcs package inside astropy. This will allow us to display the image with astronomical coordinates
In [28]:
from astropy.wcs import WCS
wcs = WCS(header)
print wcs
In the following, a plot with astrometry and some label customization. For more details, have a look at this page: https://github.com/astropy/astropy-api/blob/master/wcs_axes/wcs_api.md
In [29]:
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcs)
#ax = plt.subplot(projection=wcs)
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.imshow(data, cmap='gist_heat',origin='lower')
ra = ax.coords[0]
ra.set_major_formatter('hh:mm:ss')
dec = ax.coords[1]
dec.set_major_formatter('dd:mm:ss');
Astropy provides a way of dealing with coordinates, and automatically deal with conversions:
from astropy.coordinates import SkyCoord
# Making Coordinates:
c1 = SkyCoord(ra, dec, frame=‘icrs’, unit=‘deg’)
c2 = SkyCoord(l, b, frame=‘galactic’, unit=‘deg’)
c3 = SkyCoord(’00h12m30s’, ‘+42d12m00s’)
# Printing and Conversions:
c1.ra, c1.dec, c1.ra.hour, c2.ra.hms, c3.dec.dms
c2.fk5, c1.galactic # Converting Coordinates
c2.to_string(‘decimal’), c1.to_string(‘hmsdms’)
For instance, let's compute the coordinates of the center of the horse head:
In [30]:
from astropy.coordinates import SkyCoord
c0 = SkyCoord('5h41m00s','-2d27m00s',frame='icrs')
print c0
The wcs object contains functions that conversion from pixel to world coordinates and vice versa.
# From pixel => world:
ra, dec= w.all_pix2world(xpx, ypx, 0)# Can be lists
# The third parameter indicates if you’re starting
# from 0 (Python-standard) or 1 (FITS-standard)
# From world => pixel:
xpx, ypx= w.all_world2pix(ra, dec, 0)
In [31]:
center = wcs.all_world2pix(c0.ra,c0.dec,0)
print (center)
It is not infrequent that we need only a part of an image. So, we would like to extract this part and save another FITS file with correct astrometry. We can do this using the class Cutout2D.
This class allows one to create a cutout object from a 2D array. If a WCS object is input, then the returned object will also contain a copy of the original WCS, but updated for the cutout array.
In [32]:
from astropy.nddata import Cutout2D
size=400
cutout = Cutout2D(data, center, size, wcs=wcs)
print cutout.bbox_original
In [33]:
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8],
projection=cutout.wcs)
#ax = plt.subplot(projection=wcs)
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.imshow(cutout.data, cmap='gist_heat',origin='lower')
ra = ax.coords[0]
ra.set_major_formatter('hh:mm:ss')
dec = ax.coords[1]
dec.set_major_formatter('dd:mm:ss');
To save the new fits we have to create the header, then the extensions, finally pack all the extensions in a list and write the list to a file.
# Making a Primary HDU (required):
primaryhdu= fits.PrimaryHDU(arr1)# Makes a header
# or if you have a header that you’ve created:
primaryhdu= fits.PrimaryHDU(arr1, header=head1)
# If you have additional extensions:
secondhdu= fits.ImageHDU(arr2)
# Making a new HDU List:
hdulist1 = fits.HDUList([primaryhdu, secondhdu])
# Writing the file:
hdulist1.writeto(filename, clobber=True)
The clobber=True instruction is given to allow the rewriting of a file over an existing file. Otherwise, Python refuses to overwrite.
In [34]:
cheader = cutout.wcs.to_header()
primaryhdu = fits.PrimaryHDU(cutout.data, cheader)
hdulist = fits.HDUList([primaryhdu])
hdulist.writeto('horse.fits', overwrite=True)
While you can use the FITS interface to open tables, Astropy makes it very easy and convienientwith the astropy.table interface. For an extensive help on Tables, have a look to the documentation page:
http://docs.astropy.org/en/stable/table/
from astropy.table import Table
# Getting the first table
t1 = Table.read(filename.fits)
# Getting the second table
t2 = Table.read(filename.fits, hdu=2)
This provides a really flexible Table object that is a pleasure to deal with. It is easy to access different types of data, and read in and output to a wide variety of formats (not just FITS). Let's open the table in the extension 1 of the previous file:
In [35]:
hdulist = fits.open(image_file)
hdulist.info()
Once imported, a table can be shown with a fancy notebook interface:
In [36]:
from astropy.table import Table
t = Table.read(image_file, hdu=1)
t[:10].show_in_notebook()
Out[36]:
Or more simply printed:
In [37]:
print(t[:10])
The format can be fixed:
In [38]:
t['ETA'].format = '4.1f'
print(t[:10])
A table is both a dictionary-like and numpy array-like data type that can either be accessed by key (for columns) or index (for rows):
# Getting column names, number of rows:
t1.colnames, len(t1)
# Getting specific columns:
t1[‘name1’], t1[[‘name1’, ‘name2’]]
# Getting specific rows (all normal indexing works):
t1[0], t1[:3], t1[::-1]
# Where searching also works:
inds= np.where(t1[‘name1’] > 5)
subtable= t1[inds] # Gets all columns
For instance:
In [39]:
print t[np.where(t['ETA_CORR'] > 0.8)]
To make a table manually is easy with Numpy arrays:
# Given two columns (1D) arr1 and arr2:
t1 = Table([arr1, arr2], names=('a', 'b'))
# The columns are named “a” and “b”.
# Adding an additional column:
col1 = Table.Column(name='c', data=arr3)
t1.add_column(col1)
# Adding an additional row:
row = np.array([1, 2, 3])
t1.add_row(row)
In [40]:
import numpy as np
from astropy.table import Table
%matplotlib inline
import matplotlib.pyplot as plt
a = np.arange(0,10,0.1)
b = a**2
t1 = Table([a, b], names=('a', 'b'))
plt.plot(t1['a'],t1['b']);
To show the table in a browser in a nicely formatted manner, you can do:
t1.show_in_browser()
In [41]:
t1.write('table.txt',format='ascii.tab',overwrite=True)
Astropy provides a way to manipulate quantities, automatically taking care of unit conversions automatically.
from astropy import units as u
# Defining Quantities with units:
val1, val2 = 30.2 * u.cm, 2.2E4 * u.s
val3 = val1/val2 # Will be units cm / s
# Converting Units
val3km = val3.to(u.km/u.s)
# Simplifying Units
val4 = (10.3 * u.s/ (3 * u.Hz)).decompose()
In [42]:
from astropy import units as u
val = 30.0 * u.cm
print val.to(u.km)
# convert
val1 = 10 * u.km
val2 = 100. * u.m
# simplify
print (val1/val2).decompose()
Astropy also provides constants (with units).
from astropy import constants as c
# Some constants
c.k_B, c.c, c.M_sun, c.L_sun
# Can use with units
energy = c.h* 30 * u.Ghz
# Can convert units
mass = (3.2E13 * u.kg).to(c.M_sun)
The list of available constant is on: http://docs.astropy.org/en/stable/constants/
In [43]:
from astropy import constants as c
print 'solar mass: ', c.M_sun.value, c.M_sun.unit,'\n'
print (c.c)
Most constant can be converted in cgs units simply using the "cgs" method:
In [44]:
print c.c.cgs
There are lots of possible databases to query with astroquery. Let's see an example with the SDSS query.
To access the SDSS, there is a package called astroquery.sdss. We will import this and also the coordinate package from astropy. Let's look for a particular object and explore its FITS files for imaging and spectroscopy. We require the object to be selected only if the spectrum is available in SDSS.
In [45]:
from astroquery.sdss import SDSS
from astropy import coordinates as coords
pos = coords.SkyCoord('13h10m27.46s +18d26m17.4s',
frame='icrs')
xid = SDSS.query_region(pos, spectro=True)
xid
Out[45]:
Now, we can get the spectra and images for this list of objects using the following commands. We will obtain a list with as many objects a the list from xid. In this case, only one object.
In [46]:
sp = SDSS.get_spectra(matches=xid)
im = SDSS.get_images(matches=xid, band='r')
print len(sp), len(im)
We can also access the SDSS template library. For instance, we will get qso template with the command:
In [47]:
template = SDSS.get_spectral_template('qso')
print len(template)
Let's go back to our image. In this case the HDU list is the first element of the list. We can explore what is inside using the .info method:
In [48]:
hdulist = im[0]
hdulist.info()
Now, let's get the data.
In [49]:
header = hdulist[0].header
data = hdulist[0].data # image in 1st extension
print (data.shape, data.dtype.name)
#data = hdulist['PRIMARY'].data
#print (data.shape, data.dtype.name)
import numpy as np
plt.imshow(np.sqrt(data+1.),origin='lower',
cmap='gist_heat',vmax=1.1,vmin=0.9)
plt.colorbar();
In the case we want to display the histogram of intensity values:
In [50]:
# How to display an histogram of the intensity values
fig,ax = plt.subplots()
ax.set_yscale('log')
ax.hist(data.ravel(),200)
ax.set_xlim([0,100]);
We can be interested in displaying the image with astrometry. Let's consider a cutout around the target galaxy and overlap the contours.
In [51]:
c0 = SkyCoord('13h10m27.46s','18d26m17.4s',frame='icrs')
wcs = WCS(header)
center = wcs.all_world2pix(c0.ra,c0.dec,0)
size=400
cutout = Cutout2D(data, center, size, wcs=wcs)
In [52]:
ax = plt.subplot(projection=cutout.wcs)
ra = ax.coords[0]
ra.set_major_formatter('hh:mm:ss')
dec = ax.coords[1]
dec.set_major_formatter('dd:mm:ss')
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.imshow(np.sqrt(cutout.data+1.), cmap='gist_heat',
origin='lower',vmax=1.1,vmin=0.9,aspect='auto');
a = np.sqrt(cutout.data+1.)
mina=np.min(a)
maxa=np.max(a)
levels = np.arange(mina,maxa,(maxa-mina)/20.)
labels = [item.get_text() for item in
ax.get_xticklabels()]
ax.contour(a,levels,color='cyan');
In [53]:
from astroquery.ukidss import Ukidss
import astropy.units as u
import astropy.coordinates as coord
image_ulrs = Ukidss.get_image_list(c0,frame_type='interleave',radius=5 * u.arcmin, waveband='K',programme_id='LAS')
In [54]:
from astroquery.skyview import SkyView
survey = 'WISE 12'
sv = SkyView()
paths = sv.get_images(position='M 82',
survey=['WISE 12','GALEX Near UV'])
In [55]:
from astropy.wcs import WCS
wcs1 = WCS(paths[0][0].header)
wcs2 = WCS(paths[1][0].header)
In [56]:
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcs1)
ax.imshow(paths[0][0].data, origin='lower',
cmap='gist_heat_r')
ima2 = paths[1][0].data
levels = np.arange(np.nanmin(ima2),np.nanmax(ima2), 1.)
levels = np.nanmin(ima2)+[0.02,0.09,0.2]
ax.contour(ima2,levels, transform=ax.get_transform(wcs2),
colors='r')
plt.xlabel('RA')
plt.ylabel('Dec')
plt.show()
In [57]:
%reload_ext version_information
%version_information numpy, astropy
Out[57]: