天文データ解析入門 その9 (イメージにコントアを重ねる)

[以下追記] matplotlib のバージョンが3.7.0等だと、aplpy でエラーが起きるようです。
```python
TypeError: WCSAxes.__init__() got multiple values for argument 'wcs'
```

```
pip install matplotlib==3.5.2
```
などとしてダウングレードしてください。
[以上追記]

この記事では、イメージにコントアを重ねる方法について記述します。

今回も例として、_Spitzer_ 3.6 µm, 8.0 µm, 24 µm のアーカイブデータを使用します。
https://irsa.ipac.caltech.edu/data/SPITZER/GLIMPSE/images/I/1.2_mosaics_v3.5/
から
GLM_03000+0000_mosaic_I1.fits
GLM_03000+0000_mosaic_I4.fits 
をダウンロードします。
また、
https://irsa.ipac.caltech.edu/data/SPITZER/MIPSGAL/images/mosaics24/
から
MG0290n005_024.fits
MG0290p005_024.fits
MG0300n005_024.fits
MG0300p005_024.fits
MG0310n005_024.fits
MG0310p005_024.fits
の計６つの fits をダウンロードし、つなぎ合わせた W43_24um.fits を用います

そして、[以前の記事](https://qiita.com/Shinji_Fujita/items/bdd114823efb7c6cee42) で regrid した 3 色分の fits 
GLM_03000+0000_mosaic_I1_reg.fits
GLM_03000+0000_mosaic_I4_Reg.fits 
W43_24um.fits
を用い、[前回の記事](https://qiita.com/Shinji_Fujita/items/10c872a99bc0ecc7c8ed) で作成した RGB 合成図の png ファイルも用います。

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


# 通常の 2D イメージの上に通常のコントアを重ねる
例として、8.0 µm のデータの上に CO 輝線の強度を plot します。
まずは必要なものを import します。

```python
from astropy.io import fits
import numpy as np
from matplotlib import pyplot as plt
import aplpy
from astropy.wcs import WCS
```

CO の cube を 2D fitsにするため、以下のような関数を定義します。

```python
def v2ch(v, w): # v(km/s)をchに変える
    x_tempo, y_tempo, v_tempo   = w.wcs_pix2world(0, 0, 0, 0)
    x_ch, y_ch, v_ch   = w.wcs_world2pix(x_tempo, y_tempo, v*1000.0, 0)
    v_ch = int(round(float(v_ch), 0))
    return v_ch

def del_header_key(header, keys): # headerのkeyを消す
    import copy
    h = copy.deepcopy(header)
    for k in keys:
        try:
            del h[k]
        except:
            pass
    return h
def make_new_hdu_integ(hdu, v_start_wcs, v_end_wcs, w): # 積分強度のhduを作る
    data = hdu.data
    header = hdu.header
    start_ch, end_ch = v2ch(v_start_wcs, w), v2ch(v_end_wcs, w)
    new_data = np.nansum(data[start_ch:end_ch+1], axis=0)*header["CDELT3"]/1000.0
    header = del_header_key(header, ["CRVAL3", "CRPIX3", "CRVAL3", "CDELT3", "CUNIT3", "CTYPE3", "CROTA3", "NAXIS3", "PC1_3", "PC2_3", "PC3_3", "PC3_1", "PC3_2"])
    header["NAXIS"] = 2
    new_hdu = fits.PrimaryHDU(new_data, header)
    return new_hdu
def make_new_hdu_peak(hdu, v_start_wcs, v_end_wcs, w): # 積分強度のhduを作る
    data = hdu.data
    header = hdu.header
    start_ch, end_ch = v2ch(v_start_wcs, w), v2ch(v_end_wcs, w)
    new_data = np.nanmax(data[start_ch:end_ch+1], axis=0)
    header = del_header_key(header, ["CRVAL3", "CRPIX3", "CRVAL3", "CDELT3", "CUNIT3", "CTYPE3", "CROTA3", "NAXIS3", "PC1_3", "PC2_3", "PC3_3", "PC3_1", "PC3_2"])
    header["NAXIS"] = 2
    new_hdu = fits.PrimaryHDU(new_data, header)
    return new_hdu
```

以下のように 2D fits を定義・作成します。

```python
I4_hdu = fits.open("~/your/fits/dir/GLM_03000+0000_mosaic_I4_reg.fits")[0]
CO_hdu = fits.open("~/your/fits/dir/FGN_03100+0000_2x2_12CO_v1.00_cube.fits")[0]
w = WCS(CO_hdu.header)
peak_hdu = make_new_hdu_peak(CO_hdu, 75.0, 125.0, w)
```

例えば以下のように plot します。

```python
fig = plt.figure(figsize=(12, 12))
f = aplpy.FITSFigure(I4_hdu, figure=fig, slices=[0],convention='wells')
f.show_colorscale(vmin=10, vmax=1000, stretch='log', cmap="Greens", aspect="equal")
f.add_colorbar()
f.colorbar.show()
f.colorbar.set_axis_label_text("8.0 $\mu$m [Jy str$^{-1}$]")    
f.recenter(30.75, 0.0, width=0.5, height=0.5)
f.show_contour(peak_hdu, levels=[20, 25, 30, 35, 40, 45, 50], convention='wells', colors=["w"], slices=[0], linewidths=[1,2,2,2,2,2,2,2])
```
![ダウンロード (16).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/e9128b82-16a3-ac98-a5f8-0abc8014f3c1.png)

色や線幅をコントアレベルごとに指定できます。

# 通常の 2D イメージの上に半透明のコントアを重ねる
透過度 (alpha) を指定して、for文でレベルごとに引いていきます。

```python
fig = plt.figure(figsize=(12, 12))
f = aplpy.FITSFigure(I4_hdu, figure=fig, slices=[0],convention='wells')
f.show_colorscale(vmin=10, vmax=1000, stretch='log', cmap="Greys", aspect="equal")
f.add_colorbar()
f.colorbar.show()
f.colorbar.set_axis_label_text("8.0 $\mu$m [Jy str$^{-1}$]")    
f.recenter(30.75, 0.0, width=0.5, height=0.5)
levels = [20, 25, 30, 35, 40, 45, 50]
for l in levels:
    f.show_contour(peak_hdu, levels=[l, 99999.9], convention='wells', colors="r", slices=[0], linewidths=[1,2,2,2,2,2,2,2], filled=True, alpha=0.15)
```
![ダウンロード.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/a6f7fc2e-75db-6765-b692-500d756a6e49.png)

# RGB イメージの上にコントアを重ねる
[前回の記事](https://qiita.com/Shinji_Fujita/items/10c872a99bc0ecc7c8ed) で作成した RGB 合成図の png ファイルを用います。

```python
color_min_R = 30.0
color_max_R = 1500.0
color_min_G = 30.0
color_max_G = 1500.0
color_min_B = 10.0
color_max_B = 100.0
colorval = "%.1f_%.1f_%.1f_%.1f_%.1f_%.1f"%(color_min_R, color_max_R, color_min_G, color_max_G, color_min_B, color_max_B)
save_png_name = "RGB_%s"%(colorval)+'.png'

fig = plt.figure(figsize=(12, 12))
f = aplpy.FITSFigure(I4_hdu, slices=[0], figure=fig, convention='wells')
f.recenter(30.75, 0.0, width=0.5, height=0.5)
f.show_rgb(save_png_name)
f.ticks.set_color('w')
f.show_contour(peak_hdu, levels=[20, 25, 30, 35, 40, 45, 50], convention='wells', colors="w", slices=[0], linewidths=[1,2,2,2,2,2,2,2])
```
![ダウンロード (7).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/86192f25-5f0b-99ca-637c-6aa5cbb0dfa2.png)
以上です。

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