Python for shellscripters

解析フリークのみなさま、

High energy astro系の解析はNASAのHEASoftが主流だと思います。HEASoftに含まれるXSpecや各種ftoolsを使う上で、多くの人がshell scriptを書いていると思うのですが、個人的にshell scriptはクソだと思うので、即刻これを捨てて便利なpythonに移行しましょー！、というのが本ポストの趣旨です。

Pythonはちょいちょい使うものの、いざshell scriptを完全に置き換えたい、というとき具体的にどうすればいいんだー、と悩んだ個人的経験から、解析で頻出しそうなshell script的用法をまとめたいと思います。


 1. Shell commandや外部スクリプトの実行


os.systemを使います。

import os
os.system('cp aho1.txt aho2.txt')

()内はただのstringなのでpythonで使ってる変数も使えます。

import os
infile = 'aho1.txt'
outfile = 'aho2.txt'
os.system('cat '+infile+' > '+outfile)

これより複雑なことをしたければsubprocess.Popenを使います。

2. コマンドを実行して標準出力を返り値として得る


.communicate()を使います。

import subprocess
PIPE=subprocess.PIPE
eventfile = subprocess.Popen(['ls', 'temp.evt'], stdin=None, stdout=PIPE, shell=False).communicate()[0]

Popenの引数１つ目は実行するコマンドとその引数からなるリスト。shell=Trueにすれば、

import subprocess
PIPE=subprocess.PIPE
eventfile = subprocess.Popen('ls　temp.evt', stdin=None, stdout=PIPE, shell=True).communicate()[0]

とベタ書きしてもok。また、shell=Trueではwild cardが使えるので、

eventfiles = subprocess.Popen('ls　*_xis?.evt', stdin=None, stdout=PIPE, shell=True).communicate()[0]

なおこのままだと使い勝手がわるいし、最後に改行('\n')まで付いてきちゃうので、

eventfiles = subprocess.Popen('ls　*_xis?.evt', stdin=None, stdout=PIPE, shell=True).communicate()[0].split()

としてリストにしておくと便利。こう来ると大抵

eventfiles = subprocess.Popen('ls　*_xis?.evt', stdin=None, stdout=PIPE, shell=True).communicate()[0].split()
for eventfile in eventfiles:

とつづく。。

3. コマンド、プログラムを実行し、その標準入力とやり取りする

stdin.writeを使います。いわゆるhere documentです。XSpecやxselectを使うとき便利。
入力する各行末尾の改行を忘れずに。

import subprocess
PIPE=subprocess.PIPE
xsp = subprocess.Popen('xpec',stdin=PIPE,stdout=None,shell=False)
xsp.stdin.write(
    #ここにコマンド
    )
xsp.wait()

途中でfor回したりもできて便利。

import subprocess
PIPE=subprocess.PIPE
xse = subprocess.Popen('xselect',stdin=PIPE,stdout=None,shell=False)
xse.stdin.write(
    '\n'
    'no\n'
    'read event eventlist.txt\n'
    '.\n'
    'set image det\n'
    'set xybin 1\n'
    )

for seg in segs:

    subdir = 'xis'+xis+'_seg'+seg

    for grade in grades:
        qdpfile = str(tobs)+'_'+subdir+'_g'+grade+'.qdp'
        pifile = str(tobs)+'_'+subdir+'_g'+grade+'.pi'

        xse.stdin.write(
            'filter column "grade='+graderanges[grade]+' segment='+seg+':'+seg+' status=65536:131071 196608:262143 327680:393215 458752:524287"\n'
            'extract spec\n'
            'plot spec\n'
            'we '+subdir+'/'+qdpfile+'\n'
            'quit\n'
            'save spec '+subdir+'/'+pifile+'\n'
            'no\n'
            'clear column\n'
            )
xse.stdin.write(
    'exit\n'
    'no\n'
    '\n'
    )
xse.wait()

あとは得られたデータをnumpy, scipyで料理してmatplotlibでお絵描き、って流れです。
そこはまたおいおい。とりあえずここまでー。

さわだ