数据处理有很多方法，可以自己用C语言、Fortran写程序，也可以用各天文台提供的软件进行处理（IRAF、Ftool、AIPS、CLASS、 Starlink等）。其实用什么处理都可以，但是程序或软件好不好用，还是看相应的帮助文档和程序包是不是足够多。对于C语言来说，要处理天文中的数 据，得有相应的库，而且编译起来也麻烦。最近我逐渐发现，python是个不错的选择。
       几年前我考虑过用python进行数据处理，但是当时的理解是python执行的时候依赖于很多程序包，通用性可能有问题，而且当时读FITS文件也比较麻烦。但是最近再看，通用性的担心可能是多余的，现在几个常用程序包已经做得很好了，一般数据处理就只用到计算的包、画图的包，直接在python程序最前面加几句话就可以了。最关键的是，使用python的人很多，相应的社区很活跃，所以碰到问题的时候比较容易获取帮助。所以从现在开始尝试用python处理数据。
       下面是一个读入fits文件，画积分强度图，再把某个星表里的天体画到图上的python程序。
======================================================================
#!/usr/bin/env python

import pyfits
import math
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import mpl_toolkits
from matplotlib.patches import Ellipse

hdulist = pyfits.open('xxx.fits')
image = hdulist[0].data
nx = hdulist[0].header['naxis1']
ny = hdulist[0].header['naxis2']
nz = hdulist[0].header['naxis3']
crvalx = hdulist[0].header['crval1']
cdeltax = hdulist[0].header['cdelt1']
crpixx = hdulist[0].header['crpix1']
crvaly = hdulist[0].header['crval2']
cdeltay = hdulist[0].header['cdelt2']
crpixy = hdulist[0].header['crpix2']
crvalz = hdulist[0].header['crval3']
cdeltaz = hdulist[0].header['cdelt3']
crpixz = hdulist[0].header['crpix3']

x = np.arange(-crpixx*cdeltax+crvalx,(nx-1-crpixx)*cdeltax+crvalx,cdeltax)
y = np.arange(-crpixy*cdeltay+crvaly,(ny-1-crpixy)*cdeltay+crvaly,cdeltay)

# image[z,y,x]!!!! The first dimension is velocity channel!
Z = image[0,:,:]
Z = Z*0.0
for i in range(0,nz):
  Z = Z+image[i,:,:]

# set negative value to zero
Z = (Z+abs(Z))/2.0

ax = plt.subplot(111)
im = plt.imshow(Z, cmap=cm.gist_yarg,
                origin='lower', aspect='equal',
                extent=[max(x),min(x),min(y),max(y)])
locs, labels = plt.xticks()
xtv = locs
for i in range(0,len(locs)):
    xtv[i] = int((locs[i]-int(locs[i]/15.0)*15.0)/15.0*60.0)
x_ticks=[str(xtv[0]),str(xtv[1]),str(xtv[2]),str(xtv[3]),
          str(xtv[4]),'04h '+str(xtv[5])+' m' ,'04h '+str(xtv[6])+' m']
plt.gca().set_xticklabels(x_ticks)

ra,dec,major,minor,angle = np.loadtxt('catalog.txt',unpack=True,usecols=[16,17,18,19,20])
xr = ra
yd = dec
angle=90.0-angle
for i in range(len(ra)):
    yd[i] = dec[i]
    xr[i] = (ra[i]-crvalx)*math.cos(dec[i]*3.1415926/180.0)+crvalx

ells = [Ellipse(xy=[xr[i],yd[i]],width=major[i]*2,height=minor[i]*2,angle=angle[i])
        for i in xrange(len(ra))]

for e in ells:
    ax.add_artist(e)
    e.set_alpha(1.0)
    e.set_clip_box(ax.bbox)
    e.set_edgecolor([0.0,1.0,0.0])
    e.set_fill(False)
    e.set_linewidth(0.5)

plt.colorbar(im,orientation='vertical')
plt.show()
==========================================================

转载本文请联系原作者获取授权，同时请注明本文来自钱磊科学网博客。
链接地址：https://blog.sciencenet.cn/blog-117333-573394.html

上一篇：CLASS探索（一）
下一篇：吃亏与幸福感