Astronomy I & II sessions

[Astropy Basics]
- Constants and units
- Celestial Coordinates
- Time and dates
- FITS files
  - FITS Tables (and other common formats)
  - Spectra
  - Images & WCS
[Astroquery]
- Query to CVS/Vizier catalogues
[Photutils]
- Source detection (DAO vs SExtractor)
- Modelling: Measuring the PSF FHWM with a moffat profile
- Background modelling and aperture Photometry
[Pyephem/Astroplan]
- Earth, Time and Fixed Bodies
- Sun rising, setting and dark time
- Night planning: objects visibility and distance to moon

Other interesting packages not covered in this workshop

Astroscrappy: Cosmic ray detection / rejection / masking
Astropy-Cosmology: Cosmology utility functions
Astropy-ccdproc: Data reduction of CCD images
Sewpy: Sextractor wrapper
Naima:Derivation of non-thermal particle distributions through MCMC spectral fitting
Reproject: Astronomical image reprojection
gammapy: Gamma-ray astronomy
sncosmo: Supernova cosmology



In [ ]:

    
%matplotlib inline
import numpy as np
import math
import matplotlib.pyplot as plt
import seaborn
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 14
plt.rcParams['lines.linewidth'] = 2
plt.rcParams['xtick.labelsize'] = 13
plt.rcParams['ytick.labelsize'] = 13
plt.rcParams['axes.titlesize'] = 14
plt.rcParams['legend.fontsize'] = 13

Astropy Basics

In this section, we will explain the basics of what can be done with Astropy, such as working with internal units, opening FITS files, tables, spectra and WCS.

Constants and Units

Astropy provides a large amount of astronomical constants...

but warning: The use of units can slow down the processing of a large data set.



In [ ]:

    
from astropy import constants as const
from astropy import units as u

By default astropy constants uses S.I. units...



In [ ]:

    
print(const.c)

It can be transformed to any units...



In [ ]:

    
const.c.to('km/s')



In [ ]:

    
const.c.to('pc/yr')

You can also define your own constant using astropy Units



In [ ]:

    
my_emission_line_flux = 12.32 * u.erg / u.cm ** 2 / u.s
my_emission_line_flux

Here we can compute earth's orbit speed using astropy constants...



In [ ]:

    
speed_of_earth = const.au * 2 * math.pi / u.yr
speed_of_earth.to('km/s')

Exercise Astropy 1: Working with astronomy constants

Compute (approximately) the speed of the earth's orbit using the gravitational force between the two.



In [ ]:

    
# %load -r 2-3 solutions/10_Astronomy.py

Celestial Coordinates

The simplest coordinate we can define is a single point in the sky, by default in the ICRS frame.



In [ ]:

    
from astropy.coordinates import SkyCoord
c = SkyCoord(ra=10.625*u.degree, dec=41.2*u.degree, frame='icrs')

Definition

It can be defined in almost any format used in astronomy (and there are many, as usual...) all representing the same location.



In [ ]:

    
c = SkyCoord(10.625, 41.2, frame='icrs', unit='deg')
c = SkyCoord('00h42m30s', '+41d12m00s', frame='icrs')
c = SkyCoord('00h42.5m', '+41d12m')
c = SkyCoord('00 42 30 +41 12 00', unit=(u.hourangle, u.deg))
c = SkyCoord('00:42.5 +41:12', unit=(u.hourangle, u.deg))
c

Astropy also has a significantly large list of sources than can be retrieved by its name:



In [ ]:

    
a_big_blue_star = SkyCoord.from_name('rigel')
print (a_big_blue_star.ra, a_big_blue_star.dec)

Transformation

We can easily convert to other coordinate systems, like the galactic...



In [ ]:

    
c.galactic

Or even get what is the closest constellation to the object, very useful for astronomers as you know...



In [ ]:

    
c.get_constellation()

Distances

Coordinates allow also to define distances:



In [ ]:

    
c = SkyCoord(ra=10.68458*u.degree, dec=41.26917*u.degree, distance=770*u.kpc)
print (c.cartesian.x, c.cartesian.y, c.cartesian.z)

If we define one or more coordinates we can compute the distance between the two objects:



In [ ]:

    
c1 = SkyCoord(ra=10*u.degree, dec=9*u.degree, distance=10*u.pc, frame='icrs')
c2 = SkyCoord(ra=11*u.degree, dec=10*u.degree, distance=11.5*u.pc, frame='icrs')

print ("Angular Separation: %s" % c1.separation(c2))
print ("Distance between objects: %s" % c1.separation_3d(c2))

Catalogue of sources

A catalogue of positions can also be created using numpy arrays:



In [ ]:

    
ras = np.array([0-.7, 21.5, 120.9]) * u.deg  
decs = np.array([4.5, -5.2, 6.3]) * u.deg   
catalogue = SkyCoord(ras, decs, frame='icrs')
catalogue.galactic

Time and Date

The astropy.time package provides functionality for manipulating times and dates used in astronomy, such as UTC or MJD.

Definition



In [ ]:

    
from astropy.time import Time
times = ['2017-09-13T00:00:00', '2017-09-15T11:20:15.123456789',]
t1 = Time(times)
t1

Default format is ISOT and scale UTC, but it can be set to others.



In [ ]:

    
times = [58009, 58011.47239726]
t2 = Time(times, format='mjd', scale='tai')
t2

Format conversion



In [ ]:

    
print ("To julian date: %s" % t1[0].jd)
print ("To modified julian date: %s" % t1[0].mjd)
print ("To FITS: %s" % t1[0].fits)
print ("To GPS: %s" % t1[0].gps)
print ("To Bessel Epoch Year: %s" % t1[0].byear_str)
print ("To Julian Epoch Year: %s" % t1[0].jyear_str)

Timezones

When a Time object is constructed from a timezone-aware datetime, no timezone information is saved in the Time object. However, Time objects can be converted to timezone-aware datetime objects:



In [ ]:

    
from datetime import datetime
from astropy.time import Time, TimezoneInfo
import astropy.units as u
utc_plus_one_hour = TimezoneInfo(utc_offset=1*u.hour)
dt_aware = datetime(2000, 1, 1, 0, 0, 0, tzinfo=utc_plus_one_hour)
t = Time(dt_aware)  # Loses timezone info, converts to UTC
print(t)            # will return UTC



In [ ]:

    
print(t.to_datetime(timezone=utc_plus_one_hour)) # to timezone-aware datetime

FITS files

Tables in FITS

We will read a FITS table containing a catalogue, in this case a custom collection of Gaia stars created with CosmoHub.

With two instructions I can open the fits file and preview the content of it. For this file we find a list of two units, a primary HDU and the binary table:



In [ ]:

    
from astropy.io import fits
gaia_hdulist = fits.open('../resources/cosmohub_catalogue.fits')

gaia_hdulist.info()

To access the primary HDU, we can directly call it by its name, and preview the header like this:



In [ ]:

    
gaia_hdulist['PRIMARY'].header

As the second extension has no name, we can access to it from its index:



In [ ]:

    
gaia_hdulist[1].header

The data is contained in the Binary table and can be accessed very similarly to a numpy/pandas table:



In [ ]:

    
plt.scatter(gaia_hdulist[1].data['ra'], gaia_hdulist[1].data['dec'])
plt.xlabel('Right Ascension (deg)')
plt.ylabel('Declination (deg)')

Tables not in FITS

Even FITS is widely used in astronomy, there are a several other widely used formats for storing table and catalogue data. The module astropy.io.ascii provides read capability for most of them. Find the list of supported formats in astropy's documentation: http://docs.astropy.org/en/stable/io/ascii/index.html#supported-formats



In [ ]:

    
from astropy.io import ascii
data = ascii.read("../resources/sources.dat")  
print(data)

The read method tries to identify the file format automatically, but it can be specified in the format input parameter:



In [ ]:

    
data = ascii.read("../resources/sources.csv", format='csv')  
print(data)

Any catalogue can then be exported (in this case to screen) to any format:



In [ ]:

    
import sys
ascii.write(data, sys.stdout, format='latex')

Spectra data

Let's read a fits file containing spectra for a QSO observed with SDSS.

First we want to open the fits and inspect what it's in there...



In [ ]:

    
sdss_qso_hdulist = fits.open('../resources/sdss_qso_spec-0501-52235-0313.fits')

sdss_qso_hdulist.info()

The coadd seems to have the coadded spectra from several observations.

Let's now inspect what columns we get in the spectra:



In [ ]:

    
sdss_qso_hdulist['COADD'].columns

We can now have a look at the spectra data itself, using a scatter plot.



In [ ]:

    
plt.scatter(10**sdss_qso_hdulist['COADD'].data['loglam'], sdss_qso_hdulist['COADD'].data['flux'], s=2)
plt.xlabel('Wavelengths (Angstroms)')
plt.ylabel(r'f$\lambda$ (erg/s/cm2/A)')

The previous spectra seems to have some bad measurements, but we can make use of the OR mask included to discard those measurements.

To better visualize the spectra file we will apply a gaussian filtering...



In [ ]:

    
from scipy.ndimage.filters import gaussian_filter
from scipy.interpolate import interp1d

y_values_masked = np.ma.masked_where(sdss_qso_hdulist['COADD'].data['or_mask'], 
                                     sdss_qso_hdulist['COADD'].data['flux'])
x_values_masked = np.ma.masked_where(sdss_qso_hdulist['COADD'].data['or_mask'], 
                                     sdss_qso_hdulist['COADD'].data['loglam'])

plt.scatter(10**x_values_masked, y_values_masked, s=2, label='masked')
plt.plot(10**sdss_qso_hdulist['COADD'].data['loglam'], gaussian_filter(y_values_masked, sigma=16), 
         color='orange', linewidth=3, label='masked and filtered')
plt.xlabel('Wavelengths (Angstroms)')
plt.ylabel(r'f$\lambda$ (erg/s/cm2/A)')
plt.legend()

Another information that is included in this spectra file is the emission lines measured by SDSS. We can inspect the columns of that extension:



In [ ]:

    
sdss_qso_hdulist['SPZLINE'].data.columns

Exercise Astropy 2: Working with spectra

Display the emission lines available in the SPZLINE extension over the QSO spectra



In [ ]:

    
# %load -r 7-18 solutions/10_Astronomy.py

FITS images and WCS



In [ ]:

    
hst_hdulist = fits.open('../resources/hst_656nmos.fits')
hst_hdulist.info()



In [ ]:

    
plt.imshow(hst_hdulist['PRIMARY'].data)
plt.xlabel('X pixels')
plt.ylabel('Y pixels')
plt.colorbar()



In [ ]:

    
from astropy.visualization import ZScaleInterval

norm = ZScaleInterval()
vmin, vmax = norm.get_limits(hst_hdulist['PRIMARY'].data)

plt.imshow(hst_hdulist[0].data, vmin=vmin, vmax=vmax, interpolation='none', origin='lower')
plt.xlabel('X pixels')
plt.ylabel('Y pixels')
plt.colorbar()

The image header contains the World Coordinate System (WCS) information, stored in a set of keywords (CD, CRVAL, CRPIX and optionally some distortion parameters). The WCS provides the projection of the image in the sky, allowing to work with pixels and sky coordinates.



In [ ]:

    
print ("WCS projection type:")
print (hst_hdulist['PRIMARY'].header['CTYPE1'])
print (hst_hdulist['PRIMARY'].header['CTYPE2'])

print ("WCS reference values:")
print (hst_hdulist['PRIMARY'].header['CRVAL1'])
print (hst_hdulist['PRIMARY'].header['CRVAL2'])

print ("WCS reference pixel:")
print (hst_hdulist['PRIMARY'].header['CRPIX1'])
print (hst_hdulist['PRIMARY'].header['CRPIX2'])

print ("WCS pixel to sky matrix:")
print (hst_hdulist['PRIMARY'].header['CD1_1'])
print (hst_hdulist['PRIMARY'].header['CD1_2'])
print (hst_hdulist['PRIMARY'].header['CD2_1'])
print (hst_hdulist['PRIMARY'].header['CD2_2'])

We can load the WCS of the image directly like this:



In [ ]:

    
from astropy import wcs

hst_image_wcs = wcs.WCS(hst_hdulist['PRIMARY'].header)
hst_image_wcs.printwcs()

Once loaded the WCS, we can retrieve the corners of the image footprint:



In [ ]:

    
hst_image_wcs.calc_footprint()

We can infer its pixel scale from the CD matrix:



In [ ]:

    
hst_pixelscale = np.mean(wcs.utils.proj_plane_pixel_scales(hst_image_wcs) * u.degree).to('arcsec')
hst_pixelscale

It is also useful to know the coordinates of a specific pixel in the image:



In [ ]:

    
# Origin of the pixel coordinates convention:
# Set 0 when first pixel is 0 (c/python-like)
# Set 1 when first pixel is 1 (fortran-like)
origin = 0  

# convert the pixels
lon, lat = hst_image_wcs.all_pix2world(20, 30, origin)
print (lon, lat)

In the same way, sky coordinates can be transformed to pixel positions in the image.



In [ ]:

    
x, y = hst_image_wcs.all_world2pix(lon, lat, origin)
print (x, y)

Note that the function we used is called all_XXXXXX. This is the method to use all distortion information (such as SIP, TPV,...). To use only the WCS without the distortion, use the equivalent method wcs_XXXXXXX.

Exercise Astropy 3: Plot sources on image

Use tha gaia catalogue loaded previously and plot the stars over the HST image.

TIP: a list of coordinates can be passed directly to the WCS function.



In [ ]:

    
# %load -r 22-38 solutions/10_Astronomy.py

Using WCS axes

X and Y in the plot correspond to the X and Y of the image. WCS axes allows you to plot sky coordinates without remapping the image pixels:



In [ ]:

    
ax = plt.subplot(projection=hst_image_wcs)

ax.imshow(hst_hdulist[0].data, vmin=vmin, vmax=vmax, origin='lower')

overlay = ax.get_coords_overlay('fk5')
overlay.grid(color='white', ls='dotted')
overlay[0].set_axislabel('Right Ascension (J2000)')
overlay[1].set_axislabel('Declination (J2000)')