The library Astropy started as an effort to merge several astronomy related libraries (pyfits/asciitable/pywcs...) to read fits/ascii tables and convert astronomical coordinates. It evolved to become much more, including units/quantities cosmological calculation
In [1]:
# uncomment this line if you are using python 2
# from __future__ import print_function, division
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = 10, 10
matplotlib.rcParams['font.size'] = 20
matplotlib.rcParams['axes.labelsize'] = 'large'
In [53]:
from astropy.io import fits as pyfits
filename = "../data/data.fits"
Get quick information on the mutli-extension FITS file
In [52]:
pyfits.info(filename)
Read the 'image' extension into a numpy.array, together with the FITS header
In [55]:
image, header = pyfits.getdata(filename, 'image', header=True)
print(type(image), type(header))
The image contains the actual data
In [4]:
print(image[100:103,100:103]) # image contains the actual data
The header describes the image in terms of astrometry/unit/size
Here are the first 10 lines of header.
In [57]:
header[:10]
Out[57]:
Warning: axes in FITS files and NumPy arrays are inverted
In [51]:
print((header['naxis1'], header['naxis2'])) # this correpond to the size of the image
print(image.shape)
In [9]:
from astropy import wcs as pywcs
wcs = pywcs.WCS(header) # You can use the header to create a wcs (World Coordinate System) object.
In [10]:
ra, dec = wcs.wcs_pix2world(0.,0.,0) # sky coordinate of the center of the first pixel (0,0)
print(ra, dec)
x, y = wcs.wcs_world2pix(ra,dec,0) # back to pixel coordinate
print(x, y)
Raw/pixel index plot of the image with imshow
cmap: matplotlib colormap usedvmin, vmax: lower and upper boundaries for the colorbarorigin: where to display the [0, 0] value ('lower' = in the lower left corner)
In [11]:
plt.imshow(image, cmap='gray', origin='lower', interpolation='None', vmin=-2e-2, vmax=5e-2)
Out[11]:
Let's add information to the plot with interactive plotting (only available in notebooks)
In [12]:
%matplotlib notebook
plt.ion()
In [13]:
fig = plt.figure()
ax_wcs = fig.add_subplot(1, 1, 1, projection=wcs)
ax_wcs.imshow(image, cmap='gray', origin='lower', interpolation='None', vmin=-2e-2, vmax=5e-2)
Out[13]:
Change the frame labels
In [14]:
ax_wcs.coords['ra'].set_axislabel(r'$\alpha_\mathrm{J2000}$')
ax_wcs.coords['dec'].set_axislabel(r'$\delta_\mathrm{J2000}$')
Change tick style and add a coordinate grid
In [15]:
ax_wcs.coords['ra'].set_ticks(color='red')
ax_wcs.coords['dec'].set_ticks(color='red')
ax_wcs.coords.grid(color='red', linestyle='solid', alpha=0.7)
Add some overlay in a different frame (here galactic)
In [16]:
overlay = ax_wcs.get_coords_overlay('galactic')
overlay.grid(color='green', linestyle='solid', lw=1.0, alpha=0.5)
Add ticks and labels
In [17]:
overlay['l'].set_ticks(color='green')
overlay['l'].set_axislabel('Galactic Longitude')
overlay['b'].set_ticks(color='green')
overlay['b'].set_axislabel('Galactic Latitude')
This package can read and write in most of the ascii format that are used in astronomy http://cxc.harvard.edu/contrib/asciitable/
In [18]:
from astropy.io import ascii as asciitable
catalog = asciitable.read('../data/sources.txt') # ASCII w/wo header - CSV - IPAC - CdS - Daophot - LaTex
print(catalog)
Transform the (RA, DEC) position of the sources to pixel coordinate
In [19]:
(x,y) = wcs.wcs_world2pix(catalog['ra'], catalog['dec'], 0)
Overplot the sources from the catalog
In [20]:
ax_wcs.scatter(x, y, color='r', alpha=0.8)
Out[20]:
In [21]:
asciitable.write(catalog,'sources.tex', Writer=asciitable.Latex)
Quick check
In [22]:
!head -5 sources.tex
In [23]:
import astropy.io.votable as vo_table
votable = vo_table.from_table(catalog)
vo_table.writeto(votable, 'sources.xml')
In [24]:
!head -10 sources.xml
Astropy defines quantities as number with unit
In [25]:
from astropy import units as u
wavelength = 1.2 * u.millimeter
diameter = 30 * u.m
airy = 1.22 * wavelength/diameter
print(wavelength)
print(airy)
print(airy.decompose()) # unit change is made
Conversions are made easy
In [26]:
distance = 1.0 * u.parsec
print( distance.to(u.km) )
In the same way, astropy define usefull constants with their units
In [27]:
from astropy import constants as const
print(const.c)
print(const.c.to('km/s'))
In [28]:
frequency = 857 * u.GHz
wavelength = const.c / frequency
wavelength.to(u.micron)
Out[28]:
We can verify this works both ways
In [29]:
frequency.to(u.micron, equivalencies=u.spectral())
Out[29]:
Astropy allow for calculating the commonly used quantities as function of redshift, line distances, ages and loockback times
Cosmological constraints from different surveys are given
In [30]:
from astropy.cosmology import Planck15 as cosmo
In [31]:
print(cosmo)
We can build a redshift range and draw various cosmological distances
In [32]:
z = np.logspace(-2, 1, 100)
lumdist = cosmo.luminosity_distance(z)
angdist = cosmo.angular_diameter_distance(z)
Cosmological quantities from Astropy always come with units
In [33]:
print(lumdist[0:5])
Display the distances with a log-linear plot
In [39]:
fig, ax = plt.subplots()
ax.semilogy(z, lumdist, label='luminosity distance')
ax.semilogy(z, angdist, label='angular diameter distance')
Out[39]:
Add labels using internal units info
In [40]:
ax.set_xlabel('redshift')
ax.set_ylabel('distance [%s]' % lumdist.unit)
Out[40]:
Add a legend without frame around
In [42]:
ax.legend(frameon=False)
Out[42]:
In [78]:
import astropy.visualization as viz
Create routine to visualy compare plots
In [115]:
def plot_compare_stretch(stretchedimg, imglabel):
fig, ax = plt.subplots(1, 2, figsize=(10, 6))
ax[0].imshow(image, cmap='gray', origin='lower', interpolation=None)
ax[0].axis('off')
ax[0].set_title('Original image', fontsize=16)
ax[1].imshow(stretchedimg, cmap='gray', origin='lower', interpolation=None)
ax[1].set_title(imglabel, fontsize=16)
ax[1].axis('off')
The scale_image method to mimic DS9 scales: linear, log, sqrt, asinh..
In [119]:
sclimg = viz.scale_image(image, scale='log', min_cut=-2e-2, max_cut=5e-2)
plot_compare_stretch(sclimg, 'Log scaled image')
The logarithmic scaling can also be fine tuned..
In [118]:
a = 10**4
stretch = viz.LogStretch(a=a)
logimage = stretch(image.copy())
plot_compare_stretch(logimage, 'log-stretched image a=%d' % a)
astropy.modeling to provide a framework for representing 1-D and 2-D models and performing model fitting with parameter constraints.astropy.coordinates to work with various sky coordinates and anglesastropy.convolution to provide convolution functions and kernels that offers improvements compared to scipy routinesastropy.visualization to optimize the color stretch in your images
In [43]:
from IPython.core.display import HTML
HTML(open('../styles/notebook.css', 'r').read())
Out[43]: