本記事では、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 はファイル名が長すぎるとエラーになるので注意です。)
以上です。
リンク
目次