Hands On Astropy Tutorial - Solutions to the exercises¶

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

# distance sun - earth
distance_sun = 1 * u.au

# speed of light
speed = const.c

# time
time_sun = distance_sun / speed
print("Light travels from sun to earth in {:.2f}".format(time_sun.to('min')))

# distance to GC
distance_gc = 8 * u.kpc

# time
time_gc = distance_gc / speed
print("Light travels from the GC to earth in {:.2f}".format(time_gc.to('yr')))

import astropy.units as u

# define new unit baro-meter
bm = u.def_unit('baro-meter', 25 * u.cm)
u.add_enabled_units(bm)

# height of the empire state building
height = 381 * u.m

# convert to baro-meters
height_bm = height.to("baro-meter")
print("The height of the empire state building in baro-meters is {}".format(height_bm))

print("The numerical optimization gives: {[0]:.3e}".format(frequency_max))
print("Wien's displacement law gives   : {:.3e}".format(frequency_max_wien))

from astropy.coordinates import SkyCoord

# Coordinate of PKS2155-305
# coord_pks = SkyCoord.from_name('pks2155-304')
coord_pks = SkyCoord("21h58m52.7s", "-30d13m18s", frame="icrs")

# Coordinate of the Crab Nebula
# coord_crab = SkyCoord.from_name('crab')
coord_crab = SkyCoord(83.63 * u.deg,  22.01 * u.deg, frame='icrs')

separation_crab = coord_pks.separation(coord_crab)
print('The angular distance between PKS2155-204 and the Crab Nebula is {:.2f}'.format(separation_crab))

# Coordinate of the Galactic Center
coord_gc = SkyCoord(0, 0, unit='deg', frame='galactic')
separation_gc = coord_pks.separation(coord_gc)
print('The angular distance between PKS2155-204 and the Galactic Center is {:.2f}'.format(separation_gc))

from astropy.coordinates import EarthLocation, AltAz
from astropy.time import Time
import numpy as np
from astropy.visualization import quantity_support

# Define coordinate
coord_crab = SkyCoord(83.63 * u.deg,  22.01 * u.deg, frame='icrs')

# Set location
annecy = EarthLocation(lat=45.8058 * u.deg, lon=6.5726 * u.deg)

# Define array of times (30 minute intervals)
time_intervals = np.linspace(0 * u.day, 1 *u.day, 48)

time_start = Time('2019-04-12 00:00:00')
times = Time(time_start + time_intervals)

# Get AltAz coordinates
altaz = AltAz(obstime=times, location=annecy)
crab_altaz = coord_crab.transform_to(altaz)

# Get altitude over horizon
altitudes = crab_altaz.alt.to('deg')

# Get times when altitude > 0
above_horizon = altitudes > 0 * u.deg

with quantity_support():
    plt.plot_date(times.plot_date, altitudes, fmt='-')
    
    # orient date labels at a slant
    plt.gcf().autofmt_xdate() 
    
    alt_min, alt_max = -25 * u.deg, 75 * u.deg
    plt.fill_between(times.plot_date,
                     alt_min, alt_max,
                     where=above_horizon,
                     facecolor='green',
                     alpha=0.5)
    plt.ylim(alt_min, alt_max)

# Get example table
from astropy.table import Table

table = Table()
table['Source_Name'] = ['Crab', 'Sag A*', 'Cas A']
table['GLON'] = [184.55754381, 0, 111.74169477] * u.deg
table['GLAT'] = [-5.78427369, 0, -2.13544151] * u.deg

# Create SkyCoord object holding all 3 sky coordinates
coords = SkyCoord(table['GLON'], table['GLAT'], frame='galactic')

# Add new columns
table['RA'] = coords.icrs.ra
table['DEC'] = coords.icrs.dec

table

# Read table
fermi_2fhl = Table.read('data/fermi_2fhl_catalog.fits', hdu=1)

# Sort by brightness
fermi_2fhl.sort('Flux50')

# With -1 we get the last element
print('The brightest source in the catalog is: {} \n'.format(fermi_2fhl[-1]['ASSOC']))



# Defie SkyCoord for all objects in the table
coords = SkyCoord(fermi_2fhl['GLON'], fermi_2fhl['GLAT'], frame='galactic')

# Get Crab position
coord_crab = coords[fermi_2fhl["ASSOC"] == "Crab"]


# Find rows where separation < 1 deg
separation = coord_crab.separation(coords)

# Print result
max_separation = 10 * u.deg
fermi_2fhl[separation < max_separation]

source_classes = set(fermi_2fhl['CLASS'])

print('Source classes: {}'.format(source_classes))

from astropy.coordinates import Angle

plt.figure(figsize=(12, 5))
ax = plt.subplot(projection='aitoff')

glon = Angle(fermi_2fhl['GLON'].quantity)
glon = glon.wrap_at('180d')

glat = Angle(fermi_2fhl['GLAT'].quantity)

with quantity_support():
    for source_class in sorted(source_classes):
        selection = fermi_2fhl['CLASS'] == source_class
        ax.scatter(glon[selection].rad, glat[selection].rad, label=source_class)
    ax.grid()
    plt.legend(bbox_to_anchor=(1.15, 1), loc=1)

from astropy.io import fits

filenames = ['data/casa_0.5-1.5keV.fits.gz', 'data/casa_1.5-3.0keV.fits.gz',
             'data/casa_4.0-6.0keV.fits.gz']

data_list = []

for filename in filenames:
    image_hdu = fits.open(filename)['PRIMARY']
    data_list.append(image_hdu.data)
    
data_summed = np.sum(data_list, axis=0)
plt.imshow(data_summed, origin='lower', cmap='afmhot')

from astropy.nddata.utils import Cutout2D
from astropy.wcs import WCS

image_hdu = fits.open('data/casa_0.5-1.5keV.fits.gz')['PRIMARY']
wcs = WCS(image_hdu.header) 
image_data = image_hdu.data

pos_casa = SkyCoord('23h23m27.94s', '+58d48m42.4s', frame='icrs')

data_cutout = Cutout2D(image_data, position=pos_casa,
                       wcs=wcs, size=0.05 * u.deg)

# Find index of maximum value (this will return the
# index on the flattened array)
idx_flat = data_cutout.data.argmax()

# Transform back to the original shape
idx = np.unravel_index(idx_flat, data_cutout.data.shape)

# Transform to SkyCoord
pos = SkyCoord.from_pixel(*idx, wcs=data_cutout.wcs)
print("The position of the brightest pixel in the image is: {}".format(pos))

from astropy.wcs import WCS
from astropy.io import fits 
from astropy.coordinates import SkyCoord
from astropy.visualization import quantity_support
    
# open image
image_hdu = fits.open('data/casa_1.5-3.0keV.fits.gz')['PRIMARY']
wcs = WCS(image_hdu.header)
image_data = image_hdu.data

# define position and integration radius
pos_casa = SkyCoord('23h23m27.94s', '+58d48m42.4s', frame='icrs')

# create an array of sky positions for each pixel
yp, xp = np.indices(image_data.shape)
sky_positions = SkyCoord.from_pixel(xp=xp, yp=yp, wcs=wcs)

# calculate separation image
separation = pos_casa.separation(sky_positions)

# initialize radius arrays
radii = np.linspace(0, 5, 50).reshape(-1, 1, 1) * u.arcmin
radii_min, radii_max = radii[:-1], radii[1:]

# mask out the annulus regions
mask = (separation > radii_min) & (separation < radii_max)

data_masked = mask * image_data.reshape(-1, 1024, 1024)
data_summed = data_masked.sum(axis=-1).sum(axis=-1)

radius_center = (radii_min[:, 0, 0] + radii_max[:, 0, 0]) / 2.

r_dr = radius_center * (radii_max - radii_min).squeeze()

with quantity_support():
    plt.plot(radius_center, data_summed / (2 * np.pi * r_dr))
    plt.xlabel('Radius (arcmin)')
    plt.ylabel('Profile Amplitude (A.U.)')

from astropy.io import fits
from astropy.wcs import WCS

# We define a figure size, that has an aspect ratio to accomodate three plots side by side
fig = plt.figure(figsize=(12, 3))
filenames = ['data/casa_0.5-1.5keV.fits.gz', 'data/casa_1.5-3.0keV.fits.gz', 'data/casa_4.0-6.0keV.fits.gz']

# For convenience we use a Python loop here, but the same can be achieved
# by copy and pasting the code for one energy band
for idx, filename in enumerate(filenames):
    image_hdu = fits.open(filename)['PRIMARY']
    wcs = WCS(image_hdu.header)
    ax = plt.subplot(1, 3, idx + 1, projection=wcs)
    ax.imshow(image_hdu.data, origin='lower', cmap='afmhot')
    
    # We extract the energy range from the filename, by setting it "by hand" is also fine
    energy_range = filename[10:20]
    ax.set_title('Energy range: {}'.format(energy_range))
    
    ax.set_xlabel('RA (deg)')
    ax.set_ylabel('DEC (deg)')
    

from astropy.visualization import AsinhStretch

# use asinh stretching to brighten up the image
stretch = AsinhStretch(0.1)

data_rgb = []

for filename in filenames:
    image_hdu = fits.open(filename)['PRIMARY']
    data = image_hdu.data
    data /= data.max()
    data = stretch(data)
    data_rgb.append(data)

data_rgb_stacked = np.stack(data_rgb, axis=2)

plt.figure(figsize=(10, 10))
ax = plt.subplot(projection=wcs)
ax.imshow(data_rgb_stacked, origin='lower')

ax.set_xlabel('RA (deg)')
ax.set_ylabel('DEC (deg)')