APLpyに頼らず，Astropyのみを使ったFITSのプロットを極めたい

# 更新履歴
- (2023/01/04) 

# Astropyとは

[Astropy](https://www.astropy.org)は，天文でよく用いられるFITSファイルを読み書きするためのパッケージ．FITSを画像として出力する関数が充実している．


# なぜ極めたいのか

[APLpy](https://aplpy.github.io)を用いてよくプロットを作成していたが，どうにも，macOSの更新やARMアーキテクチャに更新が追いついていないようなので，APLpyに頼りきりも良くないなと思ったことが主な要因．同じようなことを[これまでにも記事で書いた](https://qiita.com/R1ngNebula/items/a2b7c6b1f239fe62bb74)．

正確にはAPLpyが悪いわけではなく，APLpyが依存している別のパッケージの更新が間に合っていない様子．[PyPI](https://pypi.org/project/aplpy/) を見る限り，2022年にもAPLpyはリリースされている．そんなこんなでこれまではPyenvでギリギリAPLpyが使用できていた3.6.15をインストールし，APLpyを使用していた．


しかし先日，M2のMacBook Airを購入し，環境構築をしていたところ，3.6.15でAPLpyが入れられないことに気がつき，「もうこれはAPLpyを使っている場合じゃねえ！」と考えた．

そこで，2022年12月19日現在でも，最新のPython 3.11.1にちゃんと対応しているなど，更新され続けているAstropyで全てを完結させることにした．

# 使用したデータ

野辺山45 m鏡で近傍銀河を観測した[COMING](http://jvo.nao.ac.jp/portal/nobeyama/coming.do)プロジェクトのアーカイブデータを使用する．


# 手っ取り早くプロットの仕方を知りたい人向け

さっさと方法だけ知りたい人向けにコピペ欄を用意．
これでも論文とかに使えるように，
- 座標軸
- カラーバー
- カラーバーのラベル
- ビームサイズ (空間分解能)
- グリッド
など，最低限必要な情報を全てプロットしてくれる．

とりあえず実行してみて，おかしいところは適宜修正することをおすすめする．


```Python
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from astropy.io import fits
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.axes as maxes
from radio_beam import Beam
import astropy.units as u


## Figureのsize ## 適宜調整
width = 10
height = 10

fits_name = 'fits/NGC1087_12CO_RADEC.mom0.fits'
f = fits.open(fits_name)

try:
    from astropy.wcs import WCS
    from astropy.visualization.wcsaxes import WCSAxes
    wcs = WCS(f[0].header)
    fig = plt.figure(figsize=(width, height))
    plt.rcParams['font.size'] = 20 ## 適宜調整
    ax = WCSAxes(fig, [0.16, 0.1, 0.8, 0.8], wcs=wcs) ## 適宜調整
    fig.add_axes(ax)
except ImportError:
    ax = plt.subplot(111)

image = ax.imshow(f[0].data, cmap=cm.jet, vmin=400, vmax=900) ## 適宜調整


## colorbar
divider = make_axes_locatable(ax)
colorbar_ax = divider.append_axes("top", "3%", pad="1%", axes_class=maxes.Axes)
cbar = fig.colorbar(image, orientation='horizontal', cax=colorbar_ax)
colorbar_ax.xaxis.set_ticks_position('top')


## color barの単位 ## 適宜調整
plt.title('$\mathrm{km\ s^{-1}}$', x=0.5, y=3)

## 軸の名前 ## 適宜調整
ax.set_xlabel('RA (J2000)')
ax.set_ylabel('Dec (J2000)')

##====beam sizeを楕円として描画==========
## beam sizeがFITSに入っていないとエラーが出るはず
header = fits.getheader(fits_name)
beam_size = Beam.from_fits_header(header)
my_beam = Beam(beam_size.major.value*3600*u.arcsec,
               beam_size.minor.value*3600*u.arcsec,
               beam_size.pa.value*u.deg)
ycen_pix, xcen_pix = 3, 3 ## 適宜調整
pixscale = f[0].header["CDELT1"]*3600 * u.arcsec
ellipse_artist = my_beam.ellipse_to_plot(xcen_pix, ycen_pix, pixscale)
_ = ax.add_artist(ellipse_artist)
##===================================

##====グリッド線の追加 ## 適宜調整======
ax.grid(color='black', ls='dotted')
##==================================

fig.savefig('fig/test.png')

```

## 実行結果

![test.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/567602/ac1b9e66-40d4-2d17-3462-0ecd628e2baf.png)




それぞれの解説は気が向いたら...笑

# もう一つ別のパターン

いろいろ参考にしてもう一つ別のパターンを作ってみた．

```python
import matplotlib.pyplot as plt
from astropy.wcs import WCS
from astropy.io import fits
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.axes as maxes
from radio_beam import Beam
import astropy.units as u

fits_name = 'fits/NGC1087_12CO_RADEC.mom0.fits'

hdu = fits.open(fits_name)[0]
wcs = WCS(hdu.header)

fig = plt.figure()
ax = plt.subplot(projection=wcs, label='overlays')
image = ax.imshow(hdu.data, vmin=400, vmax=900, origin='lower')

ax.grid(color='black', ls='dotted')
ax.coords[0].set_axislabel('Right Ascension (J2000)')
ax.coords[1].set_axislabel('Declination (J2000)')
fig.add_axes(ax)

#colorbar
divider = make_axes_locatable(ax)
colorbar_ax = divider.append_axes("top", "3%", pad="1%", axes_class=maxes.Axes)
cbar = fig.colorbar(image, orientation='horizontal', cax=colorbar_ax)
colorbar_ax.xaxis.set_ticks_position('top')
plt.title('$\mathrm{km\ s^{-1}}$', x=0.5, y=3)

header = fits.getheader(fits_name)
beam_size = Beam.from_fits_header(header)
my_beam = Beam(beam_size.major.value*3600*u.arcsec,
               beam_size.minor.value*3600*u.arcsec,
               beam_size.pa.value*u.deg)
ycen_pix, xcen_pix = 3, 3
pixscale = hdu.header["CDELT1"]*3600*u.arcsec
ellipse_artist = my_beam.ellipse_to_plot(xcen_pix, ycen_pix, pixscale)
_ = ax.add_artist(ellipse_artist)
plt.savefig(fname='fig/test.png')
```
## 実行結果
![test2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/567602/94f65dc0-fff9-3cdd-197a-843b00a9f8b6.png)




# 参考

https://radio-beam.readthedocs.io/en/latest/plotting_beams.html

https://pyregion.readthedocs.io/en/latest/getting_started.html#draw-regions-with-matplotlib

※pyregionは使用していないがプロットの方法を参考にした．

https://docs.astropy.org/en/stable/visualization/wcsaxes/index.html

https://astropy-regions.readthedocs.io/en/stable/plotting.html

https://astropy-regions.readthedocs.io/en/stable/plot_reg.py