python怎么算积分_科学网—python数据处理笔记(一)积分强度图 - 钱磊的博文

python数据处理笔记(一)积分强度图

已有 12600 次阅读

2012-5-21 14:38

|个人分类:知识|系统分类:科研笔记|

Python, 数据处理

数据处理有很多方法，可以自己用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()

==========================================================

转载本文请联系原作者获取授权，同时请注明本文来自钱磊科学网博客。

链接地址：http://blog.sciencenet.cn/blog-117333-573394.html

上一篇：CLASS探索（一）

下一篇：吃亏与幸福感