Astropy is a the python package designed to contain the basic utilities needed for astronomy. It has been many years in the making and replicated some of the functionaily of other astronomy packages from IRAF and IDL. However, it still is missing many packages and utilities needed for say reduction, handling spectra, etc. As it is a community driven effort, its unclear how it will take for additional functionality to be included.
Data tables are convenient way of dealing with large data sets, particularly if you need to create, read, modifiy and/or write files (i.e. ASCII, FITS, HDF5, etc...). Both Pandas and Astropy have their own unique support of data tables and it is possible to convert back and forth betweeen them. The question arises why are there two different types of data table formats and what are the benifits of using one over the other.
Pandas strengths are within its statistical machinary available for data analysis. Whereas the Astropy table has a few astronomy specific advantages; 1) being able to work with astronomy formats (i.e. FITS, sextractor, CDS, VOTables, etc) and 2) being able to use units within the table. It may becomes necessary to switch back and forth between tables depending on the task at hand.
Here we will give you a very brief introduction to Astropy tables.
First you need to load the following module inorder to utilize an Astropy table.
In [1]:
from astropy.table import QTable
import astropy.units as u
import numpy as np
Here we will create a QTable which is a quantity table. The first column is the redshift of the source, the second is the star formation rate (SFR), the third is the stellar mass and the fourth is the rotational velocity. Since it is a QTable, we can define the units for each of the entries, if we like.
In [2]:
a = [0.10, 0.15, 0.2]
b = [10.0, 2.0, 100.0] * u.M_sun / u.yr
c = [1e10, 1e9, 1e11] * u.M_sun
d = [150., 100., 2000.] * u.km / u.s
When we create the table, we can assign the names to each column
In [3]:
t = QTable([a, b, c, d],
names=('redshift', 'sfr', 'stellar_mass', 'velocity'))
We can display the table by typing the variable into the notebook
In [4]:
t
Out[4]:
We can display the information about the table (i.e. data type, units, name)
In [5]:
t.info
Out[5]:
We can view a specific column by calling is name.
In [6]:
t['sfr']
Out[6]:
We can also group subsets of columns.
In [7]:
t['sfr','velocity']
Out[7]:
We can call a specific row, in this case the first row, which shows all the entries for the first object. Remember python indexing starts with zero!
In [8]:
t[0]
Out[8]:
We can write to a standard format if we like, binary FITS tables are often used for this. If could be an ASCII file (either sapace, tab or comma delimited) or whatever your favorite sharable file format is. There is a whole list of supported formats that can be read in, many astronomy specific.
In [9]:
t.write('my_table.fits', overwrite=True)
You can read it in just as easily too. As you can see, the units are stored as well, as long as you are using a QTable.
In [10]:
all_targets = QTable.read("my_table.fits")
In [11]:
all_targets
Out[11]:
We can add new entries to the table. In this case specific star formation rate (sSFR).
In [12]:
t['ssfr'] = t['sfr'] / t['stellar_mass']
In [13]:
t
Out[13]:
Maybe we are unhappy using years and we want to convert to Gyr. There are a whole list of available units here.
In [14]:
t['ssfr'].to(1./u.gigayear)
Out[14]:
If we want to estimate the H-alpha line lunminosity from the SFR using an emperical relation (such as one from Kennicutt et al. 1998).
In [15]:
t['lum'] = (t['sfr'] * u.erg / u.s )/(7.9e-42 * u.Msun / u.yr)
In [16]:
t
Out[16]:
Say we want to go from erg to photons, which are both "energy". You can see it is not quite as straight forward by the follwoing error.
In [17]:
(t['lum']).to(u.ph / u.s)
In order to convert to photons/s properly, we need to check the units. Our beginning and ending units match (i.e. energy/time), we just need to acount for the energy of the photon which is depends on the wavelength.
$$ E_{phot} = \frac{h c}{\lambda} $$(1) As an exercise for the reader. Modify the above table to estimate the number of photons/s from the H-alpha line luminosity the energy of the photon instead.
What is photometry? It is the measurement of light from a source, either represented as flux or magnitude. There are several ways to compute the flux from a source. In this tutorial we will go over aperture photometry, where you draw a circular region around the source and measure the flux inside.
Photoutils is a useful astropy package for performing photometry. This includes but not limited to aperture photometry, PSF photometry, star finding, source isophotes, background subtraction, etc...
Large parts of this tutorial are borrowed from Getting Started with Photutils documentation.
Here we load in a star field image provided with Photutils. We are interested in a small subsection of the image in which we want to do photometry.
In [18]:
import numpy as np
from photutils import datasets
hdu = datasets.load_star_image()
image = hdu.data[500:700, 500:700].astype(float)
head = hdu.header
The type function is useful for determining what data type you are dealing with, especially if it is an unfamilar one.
In [19]:
print(type(hdu))
print(type(image))
print(type(head))
Assuming your background is uniform (and it most certianly is not), here is a rough background subtraction.
This may not work for all astronomical fields (e.g. clusters, low surface brightness features, etc). There are more sophisticated ways to deal with the background subtraction, see the Background Estimation documentation.
In [20]:
image -= np.median(image)
We use DAOStarFinder to find all of the objects in the image.
The FWHM (Full Width Half Maximum) is a way of measureing the size of the source. Often for unresolved sources, like stars, the FWHM should be matched the resolution of the telescope (or the seeing).
The threshold is the value above the background in which the star finding algorithm will find sources. Often this should be done in a statistical sense in which a real detection is n-sigma above the background (or n standard deviations). We can measure the standard deviation of background.
Statistics aside: Assuming a Gaussian distribution for the noise, 1-sigma contains 68% of the distribution, 2-sigma contains 95% of the distribution and 3-sigma contains 99.7% of the distribution. Often 5-sigma is the golden standard in physics and astronomy for a "real" detection as it is contains 99.99994% background distribution.
Here we use 3-sigma (3 standard devations) above the background for a "real" detection.
In [21]:
from photutils import DAOStarFinder
from astropy.stats import mad_std
bkg_sigma = mad_std(image)
daofind = DAOStarFinder(fwhm=4., threshold=3.*bkg_sigma)
sources = daofind(image)
for col in sources.colnames: sources[col].info.format = '%.8g' # for consistent table output
print(sources[:10])
print(sources.info)
In [22]:
from photutils import aperture_photometry, CircularAperture
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))
apertures = CircularAperture(positions, r=4.)
phot_table = aperture_photometry(image, apertures)
for col in phot_table.colnames: phot_table[col].info.format = '%.8g' # for consistent table output
print(phot_table[:10])
print(phot_table.info)
In [23]:
import matplotlib.pyplot as plt
plt.imshow(image, cmap='gray_r', origin='lower')
apertures.plot(color='blue', lw=1.5, alpha=0.5)
In [24]:
print(head.keys())
Astroquery is a Astropy package for querying astronomical databases. There many, many, many catalogs and image services. Astroquery can be used to find catalogs and images. In this tutorial, we will show a few different databases and and ways to search for data.
NED is the NASA/IPAC Extragalacitc Database. It can be useful for finding measurements, properties and papers about your favorite galaxy (or galaxies). Like many databases you can search for galaxies via their name or coordinates.
First we need to import the modules that we will use for this tutorial.
In [25]:
from astroquery.ned import Ned
from astropy import coordinates
import astropy.units as u
In [26]:
result_table = Ned.query_object("m82")
print(result_table)
Print the columns returned from the search.
In [27]:
print(result_table.keys())
Print the redshift of the source
In [28]:
print(result_table['Redshift'])
Lets now search for objects around M82.
In [29]:
result_table = Ned.query_region("m82", radius=1 * u.arcmin)
print(result_table)
Gaia is a space observatory designed to observe the precise positions of stars in order to determine their distances and motions within the galaxy. Gaia star coordinates are the new golden standard for astrometry. Being able to querry their database to correct your catalogs and images to their coordinate system would be highly beneficial.
In [30]:
from astroquery.gaia import Gaia
In [31]:
coord = coordinates.SkyCoord(ra=280, dec=-60, unit=(u.degree, u.degree), frame='icrs')
radius = u.Quantity(1.0, u.arcmin)
j = Gaia.cone_search_async(coord, radius)
r = j.get_results()
r.pprint()
print(r.keys())
IRSA is the NASA/IPAC Infrared Science Archive. IRSA typically has a infrared mission focus, such as 2MASS, Herschel, Planck, Spitzer, WISE, etc...
In [32]:
from astroquery.irsa import Irsa
In [33]:
Irsa.list_catalogs()
Out[33]:
In [34]:
table = Irsa.query_region("m82", catalog="fp_psc", spatial="Cone",radius=10 * u.arcmin)
In [35]:
table
Out[35]:
In [36]:
table.keys()
Out[36]:
SDSS or the Sloan Digital Sky Survey is a massive survey of the sky which originally was optical imaging and spectroscopy. It has been operating for over 20 years and the well beyond the original planned survey as it has expanded to near-infrared observations of stars and IFUs (integral field units) for observing galaxies and eventually the whole sky LVM (Local Volume Mapper). It is an immense data set to draw from.
In [37]:
from astroquery.sdss import SDSS
Coordinates of your favorite object, in this case M82 in degrees.
In [38]:
ra, dec = 148.969687, 69.679383
In [39]:
co = coordinates.SkyCoord(ra=ra, dec=dec,unit=(u.deg, u.deg), frame='fk5')
In [40]:
xid = SDSS.query_region(co, radius=10 * u.arcmin)
# print the first 10 entries
print(xid[:10])
print(xid.keys())
In [41]:
print(xid['ra','dec'][:10]) # print the first 10 entries
Skyview is a virtual observatory that is connected to several surveys and imaging data sets spanning multiple wavelengths. It can be a handy way to get quick access to multi-wavelength imaging data for your favorite object.
In [42]:
from astroquery.skyview import SkyView
SkyView.list_surveys()
In [43]:
pflist = SkyView.get_images(position='M82', survey=['SDSSr'],radius=10 * u.arcmin)
In [44]:
ext = 0
pf = pflist[0] # first element of the list, might need a loop if multiple images
m82_image = pf[ext].data
Plot image
In [45]:
ax = plt.subplot()
ax.imshow(m82_image, cmap='gray_r', origin='lower', vmin=-10, vmax=20)
ax.set_xlabel('X (pixels)')
ax.set_ylabel('Y (pixels)')
Out[45]:
Plot the in RA and Declination coordinates space, instead of pixels. The WCS (or World Coordinates System) in a FITS file has a set of keywords that define the coordinate system.
Typical keywords are the following:
CRVAL1, CRVAL2 - RA and Declination (in degrees) coordinates of CRPIX1 and CRPIX2.
CRPIX1, CRPIX2 - Pixel positions of CRVAL1 and CRVAL2.
CDELT1, CDELT2 - pixel scale (or plate scale) in degrees/pixel (or CD matrix: CD1_1, CD1_2, CD2_1, CD2_2 which represent the pixel scale and rotation with respect to the position angle)
In [46]:
from astropy.wcs import WCS
Use the FITS header to define the WCS of the image.
In [47]:
head = pf[ext].header
wcs = WCS(head)
ax = plt.subplot(projection=wcs)
ax.imshow(m82_image, cmap='gray_r', origin='lower', vmin=-10, vmax=20)
#ax.grid(color='white', ls='solid')
ax.set_xlabel('Right Ascension (J2000)')
ax.set_ylabel('Declination (J2000)')
#ax.scatter(xid['ra'],xid['dec'],marker="o",s=50,transform=ax.get_transform('fk5'),edgecolor='b', facecolor='none')
http://learn.astropy.org/tutorials.html - list of all the Astropy tutorials
https://docs.astropy.org/en/stable/wcs/ - Astropy WCS implementation
https://docs.astropy.org/en/stable/table/ - Astropy data tables
https://docs.astropy.org/en/stable/units/ - Astropy units
https://astroquery.readthedocs.io - astroquery documentation
https://photutils.readthedocs.io - photutils documentation
In [ ]: