Give an overview about Astropy, with the focus on using it for Gamma-Ray Astronomy
Give illustrative examples
Make some advertisement for a great open source software package
Show what I've worked on during the summer (Better late than never...)
Two weeks ago Astropy v0.3 was released
git clone https://github.com/adonath/gamma_astropy_talk.git
ipython notebook gamma_astropy_talk.ipynb --pylab inline
ipython nbconvert gamma_astropy_talk.ipynb --to slides --post serve
Catalog pipeline
Where former PyFITS and PyWCS packages are used
Installed on the cluster, can be used by:
source ~hess/software/init.sh
activate-python
Further information on Open Source Astronomy software on HOST confluence page
from astropy import constants as const
print(const.c)
Name = Speed of light in vacuum Value = 299792458.0 Error = 0.0 Units = m / s Reference = CODATA 2010
print(const.c.to('km/s'))
299792.458 km / s
print(const.c.to('pc/yr'))
0.306601393788 pc / yr
from astropy import units as u
dist_sun = 8.5 * u.kpc
dist_sun
dist_sun.to(u.m)
(u.count / u.m ** 2 / u.s / u.deg ** 2 * (0.1 * u.deg) ** 2 / u.pixel).to(u.count / u.cm ** 2 / u.s / u.pixel)
from astropy.table import Table
fermi_2FGL = Table.read('Data/gll_psc_v08.fit')
print(fermi_2FGL)
Source_Name RAJ2000 DEJ2000 ... ASSOC2 Flags ------------------ -------- -------- ... -------------------------- ----- 2FGL J0000.9-0748 0.233711 -7.8155 ... 0 2FGL J0001.7-4159 0.438849 -41.9965 ... 0 2FGL J0002.7+6220 0.679812 62.3396 ... 0 2FGL J0004.2+2208 1.05569 22.1365 ... 0 2FGL J0004.7-4736 1.18021 -47.6116 ... 0 2FGL J0006.1+3821 1.52516 38.3502 ... 0 2FGL J0007.0+7303 1.77352 73.0545 ... 0 2FGL J0007.7+6825c 1.92505 68.4232 ... 544 2FGL J0007.8+4713 1.97432 47.2298 ... 0 2FGL J0008.7-2344 2.19605 -23.7363 ... 0 2FGL J0009.0+0632 2.26231 6.54235 ... 0 ... ... ... ... ... ... 2FGL J2353.3+6643c 358.341 66.7193 ... 40 2FGL J2353.5-3034 358.384 -30.5801 ... 0 2FGL J2354.2-6615 358.565 -66.2627 ... 0 2FGL J2356.0-5256 359.021 -52.9389 ... 2 2FGL J2356.1+4034 359.028 40.5814 ... 0 2FGL J2356.3+0432 359.091 4.541 ... 384 2FGL J2358.4-1811 359.613 -18.1964 ... 0 2FGL J2358.9+6325 359.745 63.4218 ... 152 2FGL J2359.0-3037 359.759 -30.6252 ... 0 2FGL J2359.4+6751c 359.86 67.8633 ... 41 2FGL J2359.6+6543c 359.907 65.7305 ... 32
from astropy.coordinates import ICRS, Galactic
# Get crab postion from online database
crab_position = ICRS.from_name('Crab')
# Convert to galactic coordinates
crab_position.transform_to(Galactic)
<Galactic l=184.55752 deg, b=-5.78434 deg>
# Print as string
print(crab_position.to_string())
5h34m31.9399s 22d00m52.2s
from astropy.wcs import wcs
from astropy.io import fits
# Read WCS information from fits header
header = fits.getheader('Data/galactic_center_fermi.fits')
# Set up WCS transformation
wcs_transformation = wcs.WCS(header)
# Transform pixel to world coordinates
wcs_transformation.wcs_pix2world(100, 100, 0)
[array(359.95), array(0.05)]
# Transform world to pixel coordinates
wcs_transformation.wcs_world2pix(359.95, 0.05, 0)
[array(100.00000000000011), array(100.0)]
import numpy as np
from astropy.modeling import models, fitting
# Generate fake data
y, x = np.indices((101, 101))
model = models.Gaussian2D(10, 50, 50, 10, 10)
data = model(x, y) + 0.5 * np.random.randn(101, 101)
# Fit model
fitter = fitting.NonLinearLSQFitter()
fitted_model = fitter(model, x, y, data)
# Plot data and model
figsize(16, 5)
subplot(1, 2, 1)
imshow(data, origin='lower')
title('Faked data')
subplot(1, 2, 2)
imshow(data - fitted_model(x, y), origin='lower')
title('Residuals')
show()
from astropy.io import fits
# Read fits file and print info
hdulist = fits.open('Data/galactic_center_fermi.fits')
hdulist.info()
Filename: Data/galactic_center_fermi.fits No. Name Type Cards Dimensions Format 0 PRIMARY PrimaryHDU 266 (200, 200) int32 1 GTI BinTableHDU 85 26215R x 2C [D, D]
# Get and show image data
data = hdulist['PRIMARY'].data
figsize(16, 5)
imshow(data, origin='lower')
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x9ad158c>
from astropy.convolution import Tophat2DKernel, convolve
from astropy.io import fits
# Read data
data = fits.getdata('Data/galactic_center_fermi.fits')
# Define tophat kernel
tophat_kernel = Tophat2DKernel(10)
# Filter data
filtered_data_top = convolve(data, tophat_kernel)
# Plot image and kernel
figsize(16, 5)
subplot(1, 2, 1)
imshow(filtered_data_top, origin='lower')
title("Filterd data")
subplot(1, 2, 2)
imshow(tophat_kernel, origin='lower', interpolation='none')
title("Kernel Image")
<matplotlib.text.Text at 0xa76036c>
from astropy.cosmology import WMAP9, Planck13
# Redshift
z = 4
# Calculate distance using WMAP9
print "WMAP: {0}".format(WMAP9.comoving_distance(z))
# Calculate distance using Planck13
print "Planck: {0}".format(Planck13.comoving_distance(z))
WMAP: 7344.44477641 Mpc Planck: 7339.03732136 Mpc
import aplpy
from astropy.io import fits
# Set up FITS figure, draw fermi diffuse model in the first energy bin
ff = aplpy.FITSFigure('Data/gll_iem_v02_P6_V11_DIFFUSE.fit', slices=[0])
ff.show_colorscale(vmin=0)
ff.show_arrows(-10, 30, 10, -30, color='white')
ff.add_label(0, 40, "Galactic Center", size=22, color="white")
ff.show_colorbar()
INFO: Auto-setting vmax to 1.326e-05 [aplpy.aplpy]
from astroquery import fermi
# Request fermi crab data
results = fermi.FermiLAT.query_object('M1', energyrange_MeV='1000, 100000', obsdates='2013-01-01 00:00:00, 2013-01-02 00:00:00')
print(result)
['http://fermi.gsfc.nasa.gov/FTP/fermi/data/lat/queries/L1312021943384C91964326_SC00.fits', 'http://fermi.gsfc.nasa.gov/FTP/fermi/data/lat/queries/L1312021943384C91964326_PH00.fits']
# Open fits file
from astropy.io import fits
crab_event_list = fits.open(results[1])
from gammapy.stats import significance_on_off
# Compute significance
significance_on_off(n_on=18, n_off=97, alpha=0.1, method='lima')
2.2421704424844875
from gammapy.spectrum import crab
for ref in crab.refs.keys():
emin, emax = 1, 10 # TeV
int_flux = crab.int_flux(emin, emax, ref)
print("{0:15s}: {1:.3g}".format(ref, float(int_flux)))
hess_pl : 2.07e-11 hegra : 1.71e-11 hess_ecpl : 2.24e-11 meyer : 2.04e-11
Seems to be not that much for 3 months, but code has to be also tested and documented, which takes usually much more time than implementing the code itself.
from astropy.convolution import Gaussian2DKernel, Tophat2DKernel, convolve
from astropy.io import fits
# Read data
data = fits.getdata('Data/galactic_center_fermi.fits')
# Define filter kernels
gauss_kernel = Gaussian2DKernel(5)
tophat_kernel = Tophat2DKernel(5)
# Filter data
filtered_data_gauss = convolve(data, gauss_kernel, boundary='extend')
filtered_data_top = convolve(data, tophat_kernel, boundary='extend')
# Plot data
figsize(16, 12)
subplot(1, 2, 1)
imshow(filtered_data_gauss, origin='lower')
title('Gauss smooting')
subplot(1, 2, 2)
imshow(filtered_data_top, origin='lower')
title('Tophat smoothing')
show()
from IPython.display import Math
# The King or beta function
Math("$K(x, \sigma, \gamma) = \\frac{1}{2 \pi \sigma^2} \\left(1 - \\frac{1}{\gamma}\\right)"
"\\cdot \\left[1 + \\frac{x^2}{2\gamma\sigma^2}\\right]^{-\gamma}$")
# Fermi PSF can be described by core and tail component
Math("PSF(x) = f_{core}K(x, \sigma_{core}, \gamma_{core}) + (1 - f_{core})K(x, \sigma_{tail}, \gamma_{tail})")
# Lets try to model the Fermi PSF with astropy
from numpy import pi
from astropy.modeling.models import Beta2D
from astropy.convolution import Model2DKernel
# Define core and tail parameters
# Values are taken from:
# http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_PSF.html
f_core = 0.05
sigma_core, sigma_tail = 0.5399, 1.063
gamma_core, gamma_tail = 2.651, 2.932
A_core = 1 / (2 * pi * sigma_core ** 2) * (1 - 1 / gamma_core)
A_tail = 1 / (2 * pi * sigma_tail ** 2) * (1 - 1 / gamma_tail)
psf_core = Beta2D(A_core, 0, 0, 2 * gamma_core * sigma_core ** 2, gamma_core)
psf_tail = Beta2D(A_tail, 0, 0, 2 * gamma_tail * sigma_tail ** 2, gamma_tail)
# Define PSF kernel
psf_kernel = f_core * Model2DKernel(psf_core, x_size=91) + (1 - f_core) * Model2DKernel(psf_tail, x_size=91)
# Plot PSF image and profile
figsize(16, 5)
subplot(1, 2, 1)
imshow(log(psf_kernel))
colorbar()
title('PSF image')
subplot(1, 2, 2)
plot(psf_kernel.array[psf_kernel.center[0], psf_kernel.center[0]:-1])
loglog()
ylim(1e-6, 1)
xlim(1, 5e1)
title('PSF radial profile')
show()
import numpy as np
from astropy.io import fits
from astropy.wcs import wcs
from astropy.convolution import Gaussian2DKernel, convolve
# Read significance, counts and background map
counts = fits.getdata('Data/galactic_center_fermi.fits')
header = fits.getheader('Data/galactic_center_fermi.fits')
# Set up WCS transformation
wcs_transformation = wcs.WCS(header)
# Smooth data with Gaussian kernel of width 0.1°
kernel = Gaussian2DKernel(header['CDELT2'] / 0.1)
smoothed_counts = convolve(counts, kernel)
# Find peak significance
y_max, x_max = np.unravel_index(np.argmax(smoothed_counts), smoothed_counts.shape)
# Transform pixel to Galactic coordinates
lon_max, lat_max = wcs_transformation.wcs_pix2world(x_max, y_max, 0)
print("GLON: {0}, GLAT: {1}".format(lon_max, lat_max))
# Sum up excess in a radius of 0.3°
y, x = np.indices(smoothed_counts.shape)
r2 = (x - x_max) ** 2 + (y - y_max) ** 2
mask = r2 < (0.3 / header['CDELT2']) ** 2
total_counts = counts[mask].sum()
print("Total counts in 0.3 deg: {0}".format(total_counts))
GLON: 359.95, GLAT: -0.05 Total counts in 0.3 deg: 249
import numpy as np
from astropy.io import fits
from astropy.table import Table
from astropy.wcs import WCS
# Read Fermi LAT catalog source positions
fermi_catalog = Table.read('Data/gll_psc_v08.fit')
lon = fermi_catalog['GLON']
lat = fermi_catalog['GLAT']
# Read ROSAT map
image, header = fits.getdata('Data/ROSAT.fits', header=True)
# Instantiate WCS object
wcs = WCS(header)
# Find pixel positions of LAT sources. Note we use ``0`` here for the last
# argument, since we want zero based indices (for Numpy), not the FITS
# pixel positions.
x, y = wcs.wcs_world2pix(lon, lat, 0)
# Find the nearest integer pixel
x = np.round(x).astype(int)
y = np.round(y).astype(int)
# Find the ROSAT values (note the reversed index order)
values = image[y, x]
# Print out the values
print(values)
[ 123.7635498 163.27642822 221.76609802 ..., 255.07995605 100.35219574 87.62506104]
Most examples were taken from the Astropy documentation: http://docs.astropy.org/
Astropy (v0.2) paper on ADS: http://adsabs.harvard.edu/abs/2013A%26A...558A..33A
Official Astropy Website: http://www.astropy.org/
Astropy on GitHub: https://github.com/astropy/astropy
Recently version 0.3 was released: http://docs.astropy.org/en/stable/whatsnew/0.3.html
More examples and tutorials can be found at:
https://astropy4mpik.readthedocs.org
http://scipy-lectures.github.io
http://python4astronomers.github.io