%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
Normally you can use the Table.read
method to read in FITS tables, but this shows how to do it directly in astropy.io.fits
# Read in Point Source Catalog
hdulist = fits.open('data/gll_psc_v08.fit')
psc = hdulist[1].data
# Extract Galactic Coordinates
l = hdulist[1].data['GLON']
b = hdulist[1].data['GLAT']
# Coordinates from 0 to 360, wrap to -180 to 180 to match image
l[l > 180.] -= 360.
# Plot the image
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, aspect='equal')
ax.scatter(l, b)
ax.set_xlim(180., -180.)
ax.set_ylim(-90., 90.)
ax.set_xlabel('Galactic Longitude')
ax.set_ylabel('Galactic Latitude')
# Read in Background Model
hdulist = fits.open('data/gll_iem_v02_P6_V11_DIFFUSE.fit')
bg = hdulist[0].data[0, :, :]
# Read in Point Source Catalog
hdulist = fits.open('data/gll_psc_v08.fit')
psc = hdulist[1].data
# Extract Galactic Coordinates
l = hdulist[1].data['GLON']
b = hdulist[1].data['GLAT']
# Coordinates from 0 to 360, wrap to -180 to 180 to match image
l[l > 180.] -= 360.
# Plot the image
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.imshow(bg ** 0.5, extent=[-180., 180., -90., 90.], cmap=plt.cm.gist_heat,
origin='lower', vmin=0, vmax=2e-3)
ax.scatter(l, b, s=10, edgecolor='none', facecolor='blue', alpha=0.5)
ax.set_xlim(180., -180.)
ax.set_ylim(-90., 90.)
ax.set_xlabel('Galactic Longitude')
ax.set_ylabel('Galactic Latitude')