参考来源： 这个神奇的英文教程, 据说是Innsbruck University的私货
注意
这篇是一篇搬运+翻译，以后可能也有陆陆续续的该教程的搬运+翻译。
据说来自Innsbruck University的私货。（黑人问号？？？），感觉挺不错的。
转载者或者引用者，请标注本站域名和该文作者。
该文作者：岐山凤鸣，真名熊楚原，域名ecohnoch.cn
Python天文包
这里我们介绍astropy库和附属的包astroquery
两个官方文档如下：
http://docs.astropy.org/en/stable/index.html#
https://astroquery.readthedocs.io/en/latest/
安装这两个包可以用conda：
conda install astropy
conda install -c astropy astroquery
有很多和astropy相关的包，你可能感兴趣，感兴趣就看看这个表吧：
http://www.astropy.org/affiliated/
FITS 文件
FITS文件是到目前为止天文学数据最喜欢被存储的格式。它来自于两种原料。
为了读取这两种文件，我们要导入两种不同的附属于astropy的包：
from astropy.io import fits
from astropy.table import Table
FITS 文件可以保存多维数据（2维，3维），任何给定的FITS文件都能包含多个图片或者表格，它们都被叫做扩展。每一个FITS扩展都包含一个头和数据。FITS头包含WCS（世界坐标系统）信息表示一个给定的坐标点在天空的什么地方。
不像Python，FITS习惯将索引值从1开始，并且很重视这一点。
读一个FITS文件
有几个很方便的函数可以让你很轻松读取一个FITS文件：
from astropy.io import fits
img1 = fits.getdata(filename) # 获取图片
head1 = fits.getheader(filename) # 获取头
这打开了一个numpy 数组（关于numpy的话，这是一个python 的科学计算库，有机会我再写博客解释），头是一个很像Python中字典数据结构的东西。
打开其他的扩展在FITS文件中：
img1 = fits.getdata(filename, 0) # 初级扩展
img2 = fits.getdata(filename, 1) # 二级扩展
img2 = fits.getdata(filename, ext = 1) # 等效
对FITS来说，用URL来访问也是可以的，可以和下载文件的函数协同使用：
from astropy.utils.data import download_file
from astropy.io import fits
image_file = download_file('http://data.astropy.org/tutorial/FITS-images/HorseHead.fits', cache = True) # 导入马头星云的图片
为了得到我们刚刚打开的这个文件的信息，我们可以用fits.info这个函数：
结果如下：
然后我们可以决定打开一个扩展，检查它的形状并且最终用matplotlib（是一款经常和numpy联合在一起使用的东西，可以想象用matlab画图的那种东西把）可视化
image_data = fits.getdata(image_file, ext = 0)
image_data.shape
结果：
（893， 891）
import matplotlib.pyplot as plt
%matplotlob inline
plt.figure()
plt.imshow(image_data, cmap = 'gist_heat', origin = 'lower')
plt.colorbar()
结果如下
用一种通常的方法来读FITS文件
在一种寻常的方法中，定义一个头数据单元(header data unit, hdu)可以保存所有的FITS文件中的goodis列表。
然后，得到我们只需要的那个。
首先，我们先打开列表：
hdulist = fits.open(image_file)
hdulist.info()
结果如下：
然后，我们可以提取出头和数据，从第一个扩展当中，我们可以用这两种方法：
- 通过特定识别扩展的编号
- 通过特定识别扩展的名字，如果这个名字定义了的话
header = hdulist['PRIMARY'].header
data = hudlist['PRIMARY'].data
FITS文件就是这种方法读取，从第一的元素（对天文数据经常是RA）到最后的元素的一个numpy数组，一定要二次确认你有你需要的元素。
我们就可以在这个点上，关闭这个文件了。
现在，让我们来探索这个头吧。为了显示出一个人类可以理解的东西，经常用repr这个函数来添加一些换行，在关键字之间：
print(repr(header[: 10])) # 头的开始部分
结果如下：
我们可以获取这个表中的关键字，值，和一些特定的关键字和注释：
print(header[:10].keys())
print(header[:10].values())
print(header['ORIGIN'])
print(header.comments['ORIGIN'])
为了提取天体测定数据，我们可以用wcs包（在astropy中有），这允许我们用图像来显示坐标后的天文数据。
from astropy.wcs import WCS
wcs = WCS(header)
print(wcs)
结果如下：
这里有一个可视化的显示结果，如果想要看更多细节，可以去这个页面
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection = wcs)
#ax = plt.subplot(projection = wcs)
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.imshow(data, cmap = 'gist_heat', origin = 'lower')
ra = ax.coords[0]
ra.set_major_formatter('hh:mm:ss')
dec = ax.coords[1]
dec.set_major_formatter('dd:mm:ss')
效果如下：
坐标转换
Astropy提供了一种方法解决坐标，并且自动处理变换：
from astropy.coordinates import SkyCoord
# 制作坐标
c1 = SkyCoord(ra, dec, frame = 'icrs', unit = 'deg')
c2 = SkyCoord(1, b, frame = 'galactic', unit = 'deg')
c3 = SkyCoord('00h12m30s', '+42d12m00s')
# 显示和变换
c1.ra, c1.dec, c1.ra.hour, c2.ra.hms, c3.dec.dms
c2.fk5, c1.galactic # 插入转换
c2.to_string('decimal'), c1.to_string('hmsdms')
举个栗子。让我们计算马头星云的中心坐标：
from astropy.coordinates import SkyCoord
c0 = SkyCoord('5h41m00s', '-2d27m00s', frame = 'icrs')
print(c0)
显示结果：
从点到坐标，反之亦然
这个wcs对象还包含从点到世界坐标的方法，并且反之亦然
# 从点到世界
ra, dec = w.all_pix2world(xpx, ypx, 0) # 可以是一个表
# 第三个参数表示你在开始
# 从0（py标准）或者1（FITS标准）
# 从世界到点
xpx, ypx = w.all_world2pix(ra, dec, 0)
center = wcs.all_world2pix(c0, ra, c0.dec, 0)
print(center)
显示结果：
[array(534.1235215073059), array(475.5504697035576)]
切割
我们偶尔有时候需要一张图像的一小部分，所以，我们可以提取一部分图像并且存储，这里可以用类Cutout2D
这个类允许你创造一个可以剪切的二维数组对象。如果一个WCS对象在输入，之后可以返回一个包含了一个来自源WCS复制品的对象，并且更新可以剪切的那个二维数组对象。
from astropy.nddata import Cutout2D
size = 400
cutout = Cutout2D(data, center, size, wcs = wcs)
print(cutout.bbox_original)
输出结果：
((276, 675), (334, 733))
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection = cutout.wcs)
#ax = plt.subplot(projection = wcs)
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.imshow(cutout.data, cmap = 'gist_heat', origin = 'lower')
ra = ax.coords[0]
ra.set_major_formatter('hh:mm:ss')
dec = ax.coords[1]
dec.set_major_formatter('dd:mm:ss')
结果如下：
保存一个新的FITS文件
为了保存一个新的FITS文件，我们需要创造一个头，然后创造扩展。最后打包所有的扩展到一个表里，然后将表写进文件。
# 制作一个初级头数据单元(header data unit, hdu)(一定要求)
primaryhdu = fits.PrimaryHDU(arr1) # 制作一个头
# 或者如果你有一个创造好的头
primaryhdu = fits.PrimaryHDU(arr1, header = head1)
# 如果你有其他的扩展
secondhdu = fits.ImageHDU(arr2)
# 制作一个新的HDU列表
hdulist1 = fits.HDUList([primaryHDU, secondhdu])
# 写入文件
hdulist1.writeto(filename, clobber = True)
这个clobber = True表明可以覆盖已有的文件，否则Python拒绝覆盖
cheader = cutout.wcs.to_header()
primaryhdu = fits.PrimaryHDU(cutout.data, cheader)
hdulist = fits.HDUList([primaryhdu])
hdulist.writeto('horse.fits', overwrite = True)
astropy中的表
当你在用FITS交互工具打开表时，astropy让这个过程很简单和方便。对于表的更多内容的帮助，你可以看这个页面
from astropy.table import Table
# 得到前一个表
t1 = Table.read(filename.fits)
# 得到第二个表
t2 = Table.read(filename.fits, hdu = 2)
这提供了一个很有灵活性的表来处理。它很方便的就可以访问不同类型的元素，并且读取和输出都有很广泛的格式，比如不一定就是FITS格式。让我们用前一个文件来打开它的扩展1吧。
hdulist = fits.open(image_file)
hdulist.info()
运行结果如下：
一次导入，一个表就能够在一个有趣的notebook中交互了。
from astropy.table import Table
t = Table.read(image_file, hdu = 1)
t[:10].show_in_notebook()
输出结果如下：
或者更多简单的输出：
```
print(t[:10])
```
格式可以改变：
```
t['ETA'].format = '4.1f'
print(t[:10])
```
一个表既可以是像字典那样，也可以是numpy数组那样，都可以通过键来访问值：
```
# 得到列的名字，行的数量
t1.colnames, len(t1)
# 返回特定的列
t1['name1'], t1['name1', 'name2']
# 返回特定的行
t1[0], t1[:3], t1[::-1]
# 直接搜寻也可以
inds = np.where(t1['name1'] > 5)
subtable = t1[inds] # 得到所有的列
```
举个栗子：
```
print(t[np,where(t['ETA_CORR'] > 0.8)])
```
# 制作一个表
制作一个表可以很简单用Numpy数组：
```
# 给了两个列的数组
t1 = Table([arr1, arr2], names = ['a', 'b'])
# 两个列有名字a和b
# 添加更多的列
col1 = Table.Column(name = 'c', data = arr3)
t1.add_column(col1)
# 添加更多的行
row = np.array([1, 2, 3])
t1.add_row(row)
```
```
import numpy as np
from astropy.table import Table
%matplotlob inline
import matplotlib.pyplot as plt
a = np.arrange(0, 10, 0.1)
b = a ** 2
t1 = Table([a, b], names = ('a', 'b'))
plt.plot(t1['a'], t1['b'])
```
为了在浏览器中显示表，你可以用这个方法：
```
t1.show_in_browser()
```
## 保存表
```
# 输出FITS
t1.write(filename.fits)
# 特定的文本格式
t1.write(filename.txt, format = 'ascii.tab')
# 输出到LaTeX
t1.write(filename.tex, format = 'ascii.latex')
```
```
t1.write('table.txt', format = 'ascii.tab', overwrite = True)
```
# 天文数据单元
astropy 提供了一种方法去调整数量，自动处理转换单元：
```
from astropy import units as u
# 定义单元的量
val1, val2 = 30.2 * u.cm, 2.2E4 * u.s
val3 = val1 / val2 # 单位cm/s
# 转换单元
val3km = val3.to(u.km/u.s)
# 简化单元
val4 = (10.3 * u.s / (3 * u.Hz)).decompose()
```
```
from astropy import units as u
val = 30.0 * u.cm
print(val.to(u.km))
# 插入
val1 = 10 * u.km
val2 = 100. * u.m
# 简化
print((val1/val2).decompose())
```
输出：
0.0003 km, 100.0
# 天文常量
astropy 还提供了天文常量：
```
from astropy import constants as c
# 一些常量
c.k_B, c.c, c.M_sun, c.L_sun
# 可以和单元一起用
energy = c.h * 30 * u.Ghz
# 可以插入单元
mass = (3.2E13 * u.kg).to(c.M_sun)
```
常量表可以在这个[网页得到](http://docs.astropy.org/en/stable/constants/)
```
from astropy import constants as c
print('solar mass: ', c.M_sun.value, c.M_sun.unit, '\n')
print(c, c)
```
输出结果：
大多数常量都可以用简单的用cgs方法来通过cgs单元插入
```
print(c.c.cgs)
```
结果为：
2.9979245800.0 cm/s
# 天文问题
这里有很多数据库对天文的问题，我们可以看看简单的几个例子，从SDSS问题中
为了访问SDSS，有一个包叫做astroquery.sdss, 我们将要导入这个包和坐标的包从astropy中，让我们看看特定的一些对象并且怎么从FITS文件中分析图像和光谱。我们要求这个对象只能从SDSS提供的光谱中选择。
```
from astroquery.sdss import SDSS
from astropy import coordinates as coords
pos = coords.SkyCoord('13h10m27.46s +18d26m17.4s', frame = 'icrs')
xid = SDSS.query_region(pos, spectro = True)
xid
```
运行结果：（图太大了，没有截全）
现在，我们可以用上面的命令得到表中的光谱和图像，我们可以得到一个对象的表从xid中，在这个例子中，只有一个对象。
```
sp = SDSS.get_spectra(matches = xid)
im = SDSS.get_images(matches = xid, band = 'r')
print(len(sp), len(im))
```
输出结果：
1, 1
我们也可以访问SDSS的模板库，举个栗子，我们得到qso模板用下面的命令：
```
template = SDSS.get_spectral_template('qso')
print(len(template))
```
输出：
1
让我们回到我们的图片，在这个例子中HDU标作为表的第一个元素，我们要探索里面的东西可以用下面的方法：
```
hdulist = im[0]
hdulist.info()
```
输出结果：
现在，我们来获取数据：
```
header = hdulist[0].header
data = hdulist[0].data # image in 1st extension
print (data.shape, data.dtype.name)
#data = hdulist['PRIMARY'].data
#print (data.shape, data.dtype.name)
import numpy as np
plt.imshow(np.sqrt(data+1.),origin='lower',
cmap='gist_heat',vmax=1.1,vmin=0.9)
plt.colorbar();
```
输出结果：
在这个例子中，我们还可以显示强度值的直方图：
```
fig,ax = plt.subplots()
ax.set_yscale('log')
ax.hist(data.ravel(),200)
ax.set_xlim([0,100]);
```
输出结果：
我们可能会感兴趣怎么将天文的图像显示，让我们考虑剪切一个目标银河的区域然后得到它的外形：
```
c0 = SkyCoord('13h10m27.46s','18d26m17.4s',frame='icrs')
wcs = WCS(header)
center = wcs.all_world2pix(c0.ra,c0.dec,0)
size=400
cutout = Cutout2D(data, center, size, wcs=wcs)
```
```
ax = plt.subplot(projection=cutout.wcs)
ra = ax.coords[0]
ra.set_major_formatter('hh:mm:ss')
dec = ax.coords[1]
dec.set_major_formatter('dd:mm:ss')
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.imshow(np.sqrt(cutout.data+1.), cmap='gist_heat',
origin='lower',vmax=1.1,vmin=0.9,aspect='auto');
a = np.sqrt(cutout.data+1.)
mina=np.min(a)
maxa=np.max(a)
levels = np.arange(mina,maxa,(maxa-mina)/20.)
labels = [item.get_text() for item in
ax.get_xticklabels()]
ax.contour(a,levels,color='cyan');
```
输出结果：
```
from astroquery.ukidss import Ukidss
import astropy.units as u
import astropy.coordinates as coord
image_ulrs = Ukidss.get_image_list(c0,frame_type='interleave',radius=5 * u.arcmin, waveband='K',programme_id='LAS')
```
# 重合两张天文图片
```
from astroquery.skyview import SkyView
survey = 'WISE 12'
sv = SkyView()
paths = sv.get_images(position='M 82',
survey=['WISE 12','GALEX Near UV'])
```
```
from astropy.wcs import WCS
wcs1 = WCS(paths[0][0].header)
wcs2 = WCS(paths[1][0].header)
```
```
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcs1)
ax.imshow(paths[0][0].data, origin='lower',
cmap='gist_heat_r')
ima2 = paths[1][0].data
levels = np.arange(np.nanmin(ima2),np.nanmax(ima2), 1.)
levels = np.nanmin(ima2)+[0.02,0.09,0.2]
ax.contour(ima2,levels, transform=ax.get_transform(wcs2),
colors='r')
plt.xlabel('RA')
plt.ylabel('Dec')
plt.show()
```
输出结果：
# 版本号
