(搬运+翻译)天文Python包

Apr 14, 2017


参考来源: 这个神奇的英文教程, 据说是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开始,并且很重视这一点。

image

读一个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这个函数:

fits.info(image_file)

结果如下:

image

然后我们可以决定打开一个扩展,检查它的形状并且最终用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()

结果如下

image

用一种通常的方法来读FITS文件

在一种寻常的方法中,定义一个头数据单元(header data unit, hdu)可以保存所有的FITS文件中的goodis列表。

然后,得到我们只需要的那个。

首先,我们先打开列表:

hdulist = fits.open(image_file)
hdulist.info()

结果如下:

image

然后,我们可以提取出头和数据,从第一个扩展当中,我们可以用这两种方法:

  • 通过特定识别扩展的编号
  • 通过特定识别扩展的名字,如果这个名字定义了的话
header = hdulist['PRIMARY'].header
data = hudlist['PRIMARY'].data

FITS文件就是这种方法读取,从第一的元素(对天文数据经常是RA)到最后的元素的一个numpy数组,一定要二次确认你有你需要的元素。

我们就可以在这个点上,关闭这个文件了。

hdulist.close()

现在,让我们来探索这个头吧。为了显示出一个人类可以理解的东西,经常用repr这个函数来添加一些换行,在关键字之间:

print(repr(header[: 10]))  # 头的开始部分

结果如下:

image

我们可以获取这个表中的关键字,值,和一些特定的关键字和注释:

print(header[:10].keys())
print(header[:10].values())
print(header['ORIGIN'])
print(header.comments['ORIGIN'])

image

为了提取天体测定数据,我们可以用wcs包(在astropy中有),这允许我们用图像来显示坐标后的天文数据。

from astropy.wcs import WCS
wcs = WCS(header)
print(wcs)

结果如下:

image

这里有一个可视化的显示结果,如果想要看更多细节,可以去这个页面

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')

效果如下:

image

坐标转换

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)

显示结果:

image

从点到坐标,反之亦然

这个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')

结果如下:

image

保存一个新的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()

运行结果如下:

image

一次导入,一个表就能够在一个有趣的notebook中交互了。

from astropy.table import Table
t = Table.read(image_file, hdu = 1)
t[:10].show_in_notebook()

输出结果如下:

<Table length = 10>

image

或者更多简单的输出:

print(t[:10])

image

格式可以改变:

t['ETA'].format = '4.1f'
print(t[:10])

image

一个表既可以是像字典那样,也可以是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)])

image

制作一个表

制作一个表可以很简单用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'])

image

为了在浏览器中显示表,你可以用这个方法:

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)

常量表可以在这个网页得到

from astropy import constants as c

print('solar mass: ', c.M_sun.value, c.M_sun.unit, '\n')
print(c, c)

输出结果:

image

大多数常量都可以用简单的用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

运行结果:(图太大了,没有截全)

image

现在,我们可以用上面的命令得到表中的光谱和图像,我们可以得到一个对象的表从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()

输出结果:

image

现在,我们来获取数据:

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();

输出结果:

image

在这个例子中,我们还可以显示强度值的直方图:

fig,ax = plt.subplots()
ax.set_yscale('log')
ax.hist(data.ravel(),200)
ax.set_xlim([0,100]);

输出结果:

image

我们可能会感兴趣怎么将天文的图像显示,让我们考虑剪切一个目标银河的区域然后得到它的外形:

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');

输出结果:

image

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()

输出结果:

image

版本号

image