天文データ解析入門 その16 (fitsのsmoothing (convolution))

本記事では、fits のsmoothing (convolution) の方法について紹介します。異なる分解能のデータの比較の際に必要となります。天文データ解析入門 その7 (fitsのregrid) も合わせて参照ください。

例として、国立天文台の FUGIN プロジェクトで得られた野辺山45m電波望遠鏡の CO 輝線のアーカイブデータを使用します。データは
http://jvo.nao.ac.jp/portal/nobeyama/
から
FGN_03100+0000_2x2_12CO_v1.00_cube.fits
をダウンロードします (重いです)。

casa

ターミナルで casa を起動します。パスが通っている場合は、casa と打つだけです。通していない場合は、

your/casa/dir/casa-X.X.X-XXX/bin/casa

で起動します。
まず、fitsを読み込みます。

importfits(fitsimage="FGN_03100+0000_2x2_12CO_v1.00_cube.fits", imagename="FGN_03100+0000_2x2_12CO_v1.00_cube.im/", overwrite=True)

これで casa 形式 (image) になりました。最後に true と表示されていればおそらく成功です。
この fits には、header に 分解能の長軸, 短軸, 傾き, そして静止周波数が書かれています。

BMAJ    = 0.005611111111111111                                                  
BMIN    = 0.005611111111111111                                                  
BPA     =                  0.0 
RESTFRQ =       115271202000.0 / Rest Frequency (Hz)

書かれていない fits の場合は、以下のように、astropy.io.fits で開いて編集してください。

from astropy.io import fits

hdu = fits.open("~/your/fits/dir/xxx.fits")[0]
h = hdu.header
h["BMAJ"] = xx.x
h["BMIN"] = xx.x
h["BPA"] = xx.x
h["RESTFRQ"] = xx.x

fits.PrimaryHDU(hdu.data, h).writeto("~/your/fits/dir/xxx_header.fits", overwrite=True)

空間方向

ある大きさの gaussian をかけるとき

それでは、30''の gaussian で smoothing する場合を考えます。元々分解能が20''ですので、合成して約 36''(sqrt(20^2+30^2)) の分解能になります。
casa の imsmooth を実行します。(時間かかります。)

imsmooth(imagename="FGN_03100+0000_2x2_12CO_v1.00_cube.im/", major="30arcsec", minor="30arcsec", pa="0deg", targetres=False, outfile="FGN_03100+0000_2x2_12CO_v1.00_cube.smooth1.im/") 

それでは、30''の gaussian で smoothing する場合を考えます。元々分解能が20''ですので、合成して約 36''(sqrt(20^2+30^2)) の分解能になります。
casa の imsmooth を実行します。(kernel はデフォルトで gaussian になっています。計算には時間がかかります。)

imsmooth(imagename="FGN_03100+0000_2x2_12CO_v1.00_cube.im/", major="30arcsec", minor="30arcsec", pa="0deg", targetres=False, outfile="FGN_03100+0000_2x2_12CO_v1.00_cube.smooth1.im/") 

ある大きさの gaussian の分解能になるようにかけるとき

今度は、30''の分解能 に なるように smoothing する場合を考えます。元々分解能が20''ですので、約 22''(sqrt(30^2-20^2)) の gaussian を convolution すると30''の分解能になります。ただし、casa の場合、それを自分で計算する必要はなく、option の targetres を True にすることで major と minor で指定した分解能になるように header から内部的に計算してくれます。
以下のように casa の imsmooth を実行します。(時間かかります。)

imsmooth(imagename="FGN_03100+0000_2x2_12CO_v1.00_cube.im/", major="30arcsec", minor="30arcsec", pa="0deg", targetres=True, outfile="FGN_03100+0000_2x2_12CO_v1.00_cube.smooth2.im/") 

速度方向

残念ながら、casa の task の中に速度方向の smoothing ができるものは見つけられませんでした。ただし、imrebin という rebin (まとめあげ) する task はありました。これを速度方向だけに適用します。factor で３軸分のまとめあげサイズ (channel 数) を指定します。(時間かかります。)

imrebin(imagename="FGN_03100+0000_2x2_12CO_v1.00_cube.im/", outfile="FGN_03100+0000_2x2_12CO_v1.00_cube.rebin.im/", factor=[1,1,3])

どうしても速度方向に smoothing をしたい場合は、scipy.signal.convolve などを使用すると良いかと思います。

最後に、image を fits に戻して終わりです。

exportfits(imagename="FGN_03100+0000_2x2_12CO_v1.00_cube.smooth1.im/", fitsimage="FGN_03100+0000_2x2_12CO_v1.00_cube.smooth1.fits", velocity=True, dropstokes=True, overwrite=True) 

分解能は落としましたが、ノイズが軽減されているのがわかります。
(左: original, 右: 30''smoothing)

miriad

fits in=FGN_03100+0000_2x2_12CO_v1.00_cube.fits out=FGN_03100+0000_2x2_12CO_v1.00_cube.mir op=xyin

convol map=FGN_03100+0000_2x2_12CO_v1.00_cube.mir fwhm=30 out=FGN_03100+0000_2x2_12CO_v1.00_cube.mirconv.mir

fits in=FGN_03100+0000_2x2_12CO_v1.00_cube.mirconv.mir  out=FGN_03100+0000_2x2_12CO_v1.00_cube.mirconv.fits op=xyout

上と同様の fits が出てくると思います。
casa と同様、final という option を入れることで tagetres を指定することができます。詳細は CONVOL を参照してください。
(miriad はファイル名が長すぎるとエラーになるので注意です。)

以上です。

リンク
目次