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

[以下追記] matplotlib のバージョンが3.7.0等だと、aplpy でエラーが起きるようです。

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 を用います

そして、以前の記事 で regrid した 3 色分の fits
GLM_03000+0000_mosaic_I1_reg.fits
GLM_03000+0000_mosaic_I4_Reg.fits
W43_24um.fits
を用い、前回の記事 で作成した RGB 合成図の png ファイルも用います。

これらに、国立天文台の FUGIN プロジェクトで得られた野辺山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 します。

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にするため、以下のような関数を定義します。

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 を定義・作成します。

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 します。

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])

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

通常の 2D イメージの上に半透明のコントアを重ねる

透過度 (alpha) を指定して、for文でレベルごとに引いていきます。

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)

RGB イメージの上にコントアを重ねる

前回の記事 で作成した RGB 合成図の png ファイルを用います。

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])

以上です。

リンク
目次