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

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

例として、国立天文台の [FUGIN](https://nro-fugin.github.io/) プロジェクトで得られた野辺山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 で開いて編集してください。

```python
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](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve.html#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)
![スクリーンショット 2021-06-18 9.36.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/d57b6850-a6a9-31ad-5e43-3552d9b9e7f8.png)

# 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](https://www.atnf.csiro.au/computing/software/miriad/doc/convol.html) を参照してください。
(miriad はファイル名が長すぎるとエラーになるので注意です。)

以上です。

リンク
[目次](https://qiita.com/Shinji_Fujita/items/ecb869e5b6a9fa468483)