First day: Dealing with spectra
Second day: Dealing with tables and imaging
Reading an ascii file
In [8]:
%matplotlib inline
In [9]:
# https://github.com/astropy/specutils/raw/master/specutils/io/tests/files/multispec_equispec.11.dat
In [10]:
%%bash
curl -O https://raw.githubusercontent.com/astropy/specutils/master/specutils/io/tests/files/multispec_equispec.11.dat
In [11]:
with open('multispec_equispec.11.dat','r') as fh: # fh is short for file handle
# file consists of two columns, e.g.:
# 14740.266391838 0.8220932
# 14743.8622868028 -1.856567
# declare empty list
all_lines = []
for line in fh:
# each line will split on a whitespace
all_lines.append(line.split())
It's more convenient to work with arrays. Also, we want the values to be floats intead of strings as they are now:
In [12]:
all_lines[:2]
Out[12]:
In [13]:
float_lines = [list(map(float,x)) for x in all_lines]
In [14]:
print("first two float lines: ",float_lines[:2])
print("length(float_lines): ",len(float_lines))
What if we want the array to have dimensions [2,1024] instead of [1024,2]?
In [15]:
float_lines_inverted = list(zip(*float_lines))
In [16]:
float_lines[0][:10]
Out[16]:
Numpy arrays are much more convenient to work with and are generally faster. As long as you have a list of numbers (not a list of strings), they are easy to use:
In [17]:
import numpy as np
float_lines_array = np.array(float_lines)
In [18]:
float_lines_array.shape
Out[18]:
For example, transposing an array is much easier with numpy:
In [19]:
float_lines_array.T.shape
Out[19]:
In [20]:
xaxis, yaxis = float_lines_array.T
With nested lists, you need to index each layer separately, whereas with numpy arrays you can index them together:
In [21]:
float_lines_array[:5,1]
Out[21]:
In [22]:
float_lines[:5]
Out[22]:
In [23]:
# difficult to access the second column:
list(zip(*float_lines[:5]))[1]
Out[23]:
Arrays can be manipulated like any other number, and arithmetic operations will be applied to each element:
In [24]:
5 * float_lines_array[:5,1]
Out[24]:
In [25]:
import pylab as pl
In [26]:
pl.plot(xaxis, yaxis)
pl.savefig("my_first_spectrum_plot.pdf")
In [27]:
import numpy as np
In [28]:
arr = np.loadtxt('multispec_equispec.11.dat')
arr
Out[28]:
In [29]:
arr = np.genfromtxt('multispec_equispec.11.dat')
arr
Out[29]:
In [30]:
arr = np.genfromtxt('multispec_equispec.11.dat', delimiter=" ", comments="#",
skip_header=0, skip_footer=0)
arr
Out[30]:
In [31]:
from astropy.table import Table
from astropy.io import ascii
In [32]:
tbl = Table.read('multispec_equispec.11.dat', format='ascii.no_header', delimiter=' ')
tbl
Out[32]:
In [33]:
import pandas as pd
ptbl = pd.read_csv('multispec_equispec.11.dat', delim_whitespace=True, header=None)
ptbl
Out[33]:
In [34]:
%timeit pd.read_csv('multispec_equispec.11.dat', delim_whitespace=True, header=None)
In [35]:
%timeit Table.read('multispec_equispec.11.dat', format='ascii.no_header', delimiter=' ')
Read this text file:
https://raw.githubusercontent.com/astropy/specutils/master/specutils/io/tests/files/AAO_11.txt
which contains both a header and two data columns. Plot it, and save as a .png and as a .pdf
Write your own function for file reading. Based on the original example, write a function that reads a 2-column (or n-column) space-separated text file into a numpy array. Compare its execution time to that of pandas.read_csv and astropy.io.table.Table.read.
URL for notebook from this session: goo.gl/EIbNDg
In [29]:
%%bash
curl -O https://raw.githubusercontent.com/astropy/specutils/master/specutils/io/tests/files/AAO_11.txt
In [45]:
!head -n 180 AAO_11.txt
How do we exclude the first N lines?
In [44]:
with open('AAO_11.txt','r') as fh:
for ii,line in enumerate(fh):
if 'END' in line:
last_header_line_number = ii
aao_data = np.genfromtxt('AAO_11.txt', skip_header=last_header_line_number+1)
print(last_header_line_number)
In [48]:
float('nan') # works
float('1.234') # works
float('blah') # fails
In [46]:
with open('AAO_11.txt','r') as fh:
for ii,line in enumerate(fh):
try:
float(line.split()[0])
break
except (ValueError, IndexError) as ex:
print(ex)
continue
last_header_line_number = ii
print(last_header_line_number)
In [37]:
ll = ['a','b','c']
for ii,x in enumerate(ll):
print(ii, x)
In [35]:
aao_data
Out[35]:
In [53]:
# inf vs nan
np.inf, -np.inf, np.nan
(np.isfinite([np.inf, -np.inf, np.nan, 0]),
np.isnan([np.inf, -np.inf, np.nan, 0]),
np.nan == np.nan)
Out[53]:
In [55]:
None, bool(None), None == np.nan
Out[55]:
In [59]:
def f():
pass
f() is None # preferred
f() == None # can result in problems in some situations
Out[59]:
In [40]:
%%bash
curl -O https://raw.githubusercontent.com/astropy/specutils/master/specutils/io/tests/files/gbt_1d.fits
FITs files can be read....
In [41]:
from astropy.io import fits
from astropy.wcs import WCS
In [42]:
fh = fits.open('gbt_1d.fits')
fh
Out[42]:
In [43]:
data_hdu = fh[0]
In [44]:
data_hdu.header
Out[44]:
In [45]:
data_hdu.data.shape
Out[45]:
In [46]:
pl.plot(data_hdu.data)
Out[46]:
Generate the x axis:
In [47]:
hdr = data_hdu.header
In [48]:
xarr = (np.arange(hdr['NAXIS1']) - hdr['CRPIX1'] + 1) * hdr['CDELT1'] + hdr['CRVAL1']
In [49]:
pl.plot(xarr, data_hdu.data)
Out[49]:
In [50]:
from astropy import units as u
xarr_u = xarr*u.Unit(hdr['CUNIT1'])
In [51]:
xarr_u
Out[51]:
In [52]:
xarr_u.unit.to_string(format='latex')
Out[52]:
In [53]:
pl.plot(xarr_u, data_hdu.data)
pl.xlabel(xarr_u.unit.to_string(format='latex'), fontsize=40)
pl.ylabel(u.Unit(hdr['BUNIT']).to_string(format='latex'), fontsize=40)
Out[53]:
In [54]:
pl.plot(xarr_u, data_hdu.data)
pl.xlabel(xarr_u.unit.to_string(format='latex_inline'), fontsize=40)
pl.ylabel(u.Unit(hdr['BUNIT']).to_string(format='latex'), fontsize=40)
Out[54]:
In [55]:
mywcs = WCS(hdr)
mywcs
Out[55]:
In [56]:
xarr_again = mywcs.wcs_pix2world(np.arange(hdr['NAXIS1']), 0) * u.Unit(mywcs.wcs.cunit[0])
In [57]:
xarr_u
Out[57]:
In [58]:
xarr_again
Out[58]:
In [59]:
np.isclose(xarr_again, xarr_u, atol=0)
Out[59]:
In [60]:
from specutils.io import fits
In [61]:
spec = fits.read_fits_spectrum1d('gbt_1d.fits')
In [62]:
import specutils
import numpy
import astropy
specutils.__version__, astropy.__version__, numpy.__version__
Out[62]:
In [63]:
spec.velocity
Out[63]:
In [64]:
pl.plot(spec.velocity, spec.flux)
pl.xlabel(spec.velocity.unit.to_string(format='latex_inline'), fontsize=40)
pl.ylabel(spec.flux.unit.to_string(format='latex'), fontsize=40)
Out[64]:
In [65]:
# you'll need to ``pip install pyspeckit`` to get access to this
import pyspeckit
In [66]:
sp = pyspeckit.Spectrum('gbt_1d.fits')
In [67]:
sp.plotter()
np.convolvescipy.ndimage.convolveastropy.convolutionpyspeckit.smooth
In [1]:
import specutils
from specutils.io import fits
In [60]:
spec = fits.read_fits_spectrum1d('gbt_1d.fits')
In [61]:
spec.flux
Out[61]:
In [62]:
spec.velocity
Out[62]:
In [6]:
from astropy import units as u
In [63]:
new_xarr = np.linspace(-50, 50, 1000)*u.km/u.s
In [66]:
new_xarr
Out[66]:
In [68]:
np.interp?
In [23]:
interpolated_data = np.interp(new_xarr[::-1], spec.velocity[::-1], spec.flux[::-1])
In [24]:
%matplotlib inline
import pylab as pl
In [27]:
pl.plot(spec.velocity, spec.flux)
pl.plot(new_xarr[::-1], interpolated_data)
Out[27]:
In [28]:
pl.plot(spec.velocity, spec.flux)
pl.plot(new_xarr[::-1], interpolated_data)
pl.xlim(-50, 50)
Out[28]:
In [ ]: