fitsをアタマから作りたい

はじめに

天文データでよく使われるfitsと仲良くなれたと思ったときに記録を残していきます。
これはfitsとはなんぞやというのを勉強するために1から作ってみました。

まず、fitsは天文データによく使われるファイル形式であり、座標、周波数といった情報を格納するヘッダー部と、観測から得られた情報を格納するデータ部の2から構成される。詳しくはwikipediaで...(逃げ)
https://ja.wikipedia.org/wiki/FITS

#fitsの作成
なんだかんだで悩んだfitsの作成。本当にこの方法でいいのか。。。いいと信じたい。
とりあえず[100, 100]の配列を持つfitsを作るスクリプトを示します。

make_fits_1.py

import numpy as np
from astropy.io import fits

x = 100      #x_axis pixel number
y = 100      #y_axis pixel number
z = 100      #z_axis pixel number

###配列の作成###
#3次元の時
data = np.zeros((z, y, x), dtype=np.float32)
#2次元の時
#data = np.zeros((y, x), dtype=np.float32)

hdu = fits.PrimaryHDU(data = data)

#作ったfitsの保存
hdu.writeto('make_fits.fits',overwrite=True)

3次元cubeのfitsを作るときは「配列の作成」の部分を変更します。
この状態ではまだヘッダーを入れていないのでds9などで表示しても味気ないのしか出てこない。

ヘッダーを入れていきましょう。
ヘッダーの追加はこれでできる。

ヘッダーの追加

hdu.header['keyword'] = 'something' #keywordの時
hdu.header['keyword'] = something   #座標といった数値情報を入れる時

お試しに、100x100 の配列をもち、
WCS: Galactic
center position 00:00:00.00 00:00:00.00 gal
beamsize 20.0 arcsec
gridsize 10.0 arcsec
のfitsを作ってみる。

make_fits_2.py

import numpy as np
from astropy.io import fits

x = 100      #x_axis pixel number
y = 100      #y_axis pixel number
z = 100      #z_axis pixel number

beamsize = 20.0 #arcsec

###配列の作成###
#3次元の時
data = np.zeros((z, y, x), dtype=np.float32)
#2次元の時
data = np.zeros((y, x), dtype=np.float32)

hdu = fits.PrimaryHDU(data = data)

#入れたいヘッダ一覧                                      
hdu.header['BITPIX']  =                  -64                                                 
hdu.header['NAXIS']   =                    2                                                 
hdu.header['NAXIS1']  =                  x                                                  
hdu.header['NAXIS2']  =                  y
hdu.header['BPA']     = 0.000000
hdu.header['BMAJ']    = beamsize/3600                                                  
hdu.header['BMIN']    = beamsize/3600
hdu.header['EPOCH']   =                 2000                                                                                                    
hdu.header['BUNIT']  = 'K       '                                                                                                    
hdu.header['CTYPE1']  = 'GLON-GLS'                                                            
hdu.header['CRVAL1']  =                  0.0                                                  
hdu.header['CDELT1']  =     -beamsize/2/3600                                                  
hdu.header['CRPIX1']  =                  x/2                                                  
hdu.header['CROTA1']  =                  0.0                                                  
hdu.header['CUNIT1']  = 'deg     '                                                            
hdu.header['CTYPE2']  = 'GLAT-GLS'                                                            
hdu.header['CRVAL2']  =                  0.0                                                  
hdu.header['CDELT2']  =      beamsize/2/3600                                                  
hdu.header['CRPIX2']  =                  y/2                                                  
hdu.header['CROTA2']  =                  0.0                                                  
hdu.header['CUNIT2']  = 'deg     '          

#作ったfitsの保存
hdu.writeto('make_fits_2.fits',overwrite=True)

これでちゃんと座標情報が入ったfitsが作れました。

#fitsを作る注意点(自分への戒め)
ヘッダーを間違えると間違ったデータになっちゃうので確認はしっかりしましょう。