【Python载入数据】python极简教程
python是一种不一定简洁、但非常明了的语言,可读性强,虽然依赖的程序包有时候比较多。就天文数据的处理而言,python是非常合适的。之前已经写过一些python相关的博文(http://blog.sciencenet.cn/blog-117333-679285.html,http://blog.sciencenet.cn/blog-117333-679464.html,http://blog.sciencenet.cn/blog-117333-684226.html,blog.sciencenet.cn/blog-117333-694969.html,http://blog.sciencenet.cn/home.php?mod=space&uid=117333&do=blog&id=695303,http://blog.sciencenet.cn/blog-117333-573394.html,http://blog.sciencenet.cn/blog-117333-574518.html,http://blog.sciencenet.cn/blog-117333-574767.html)。和之前IDL的博文一样,这篇博文的目的也是起到总结的作用,以后不怎么会开新帖了。
python的各个版本可以在这里找到,http://www.python.org/download/releases/。更多python的程序包可以在一个非官方的python程序库网站http://www.lfd.uci.edu/~gohlke/pythonlibs/找到。python画图程序包matplotlib的相关内容可以参见http://matplotlib.org/。下面总结一些数据处理中经常用到的python命令。
注意:缩进在python里是有意义的,所以写代码的时候不能随便空格,要注意对齐。
1. 读取各种格式的文件
(1)FITS
>>> import pyfits
>>> hdulist = pyfits.open('xxx.fits') # 读头文件及数据,hdulist是一个结构体
>>> image = hdulist[0].data # 分离出数据
>>> nx = hdulist[0].header['naxis1'] # 读取头文件中相应的参数值
(2)文本文件
>>> import numpy as np
>>> ra,dec,major,minor,angle = np.loadtxt('catalog.txt',unpack=True,usecols=[16,17,18,19,20]) # unpack=True 表明对结果进行转置,使得ra,dec,major=...的格式是合法的
2. 定义数组
(1)定义一个未赋值数组
>>> import numpy as np
>>> a=np.array(10)
(2)定义一个所有元素为0的数组
>>> import numpy as np
>>> a=np.zeros(10)
>>> b=np.zeros((2,5))
(3)定义一个所有元素为1的数组
>>> import numpy as np
>>> a=np.ones(10)
>>> b=np.ones((2,5))
(4)定义一个指标数组
>>> import numpy as np
>>> a=np.arange(1,11,1)
>>> b=np.arange(1,11,1).reshape(2,5)
3. 计算
(1)傅里叶变换
>>> import numpy as np
>>> a=np.ones(100)
>>> b=np.fft.fftn(a)
>>> c=np.fft.ifftn(b).real
4. 画图
(1)积分强度图
# in intensity_map.py:
import pyfits
import numpy as np
import matplotlib.pyplot as plt
hdulist = pyfits.open('xxx.fits')
image = hdulist[0].data
nz = hdulist[0].header['naxis3']
for i in range(0,nz):
Z = Z+image[i,:,:]
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)])
plt.show()
(2)简单曲线图
# in simple_plot.py:
from pylab import *
t=arange(0.0, 2.0, 0.01)
s=sin(2*pi*t)
plot(t,s)
xlabel('time (s)')
ylabel('voltage (mV)')
savefig("test.eps")
show()
5. 查看帮助
>>> help(np.array) #查看np.array函数的定义和用法