天文データ解析入門 その8 (aplpy で RGB合成図を作成する)

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

TypeError: WCSAxes.__init__() got multiple values for argument 'wcs'

pip install matplotlib==3.5.2

などとしてダウングレードしてください。
[以上追記]

この記事では、aplpy を用いたRGB画像の作り方について記述します。

今回も例として、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
を用います。

まずは必要なものを import します。

import aplpy

これだけです。
fits の場所と名前を変数にしておきます。

data_fits_R = "~/your/fits/dir/W43_24um.fits"
data_fits_G = "~/your/fits/dir/GLM_03000+0000_mosaic_I4_reg.fits"
data_fits_B = "~/your/fits/dir/GLM_03000+0000_mosaic_I1_reg.fits"
fitss = [data_fits_R, data_fits_G, data_fits_B]

それぞれカラースケールを決めます。そして保存する名前も定義します。

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'

以下で png ファイルを生成します (少し時間がかかります)。

aplpy.make_rgb_image(fitss, save_png_name, vmin_r=color_min_R, vmax_r=color_max_R,vmin_g=color_min_G, vmax_g=color_max_G, vmin_b=color_min_B, vmax_b=color_max_B, stretch_r="log", stretch_g="log", stretch_b="log")

もし何かエラー (大抵 PIL) が出た場合は、エラーメッセージで google 検索すると解決できる可能性が高いです。どうしてもダメな場合は ipython などで実行するといいかもしれません。

座標をつけて plot したい場合は以下のようにします。

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(16, 12))
f = aplpy.FITSFigure(data_fits_R, slices=[0], figure=fig, convention='wells')
f.show_rgb(save_png_name)
f.ticks.set_color('w')
f.save('RGB_aplpy.pdf', dpi=300)

欠損部分について

24 µmのデータは強すぎるところがサチっていて NaN が入っています。
この NaN 部分をある値 (大抵 2000 Jy/str くらい) に置き換えたい場合は、

from astropy.io import fits
import numpy as np

hdu = fits.open("W43_24um.fits")[0]
data = hdu.data
data[data!=data] = 2000.0 # data!=dataは NaN の場所を示す

fits.PrimaryHDU(data, hdu.header).writeto("W43_24um_nan2000.fits", overwrite=True)

として、こちらを使用してください。

以上です。

リンク
目次