[以下追記] matplotlib のバージョンが3.7.0等だと、aplpy でエラーが起きるようです。
TypeError: WCSAxes.__init__() got multiple values for argument 'wcs'
pip install matplotlib==3.5.2
などとしてダウングレードしてください。
[以上追記]
この記事では、(備忘録として、) aplpy を用いて画像の作る際の設定について記述します(随時更新する可能性あり)。
今回は例として、NGC628 の、JWST と ALMA (PHANGS) の公開データを使用します。
このサイトから JWST の中間赤外線のデータ (jw02107-o039_t018_miri_f770w_i2d.fits) をダウンロードします。
また、このサイトのリンクから ALMA のCO のデータ (ngc0628_12m+7m+tp_co21.fits) をダウンロードします。
まずは必要なものを import します。
from astropy.io import fits
from astropy.wcs import WCS
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.coordinates import ICRS, Galactic, FK4, FK5
import matplotlib.pyplot as plt #v3.5.2
import numpy as np
import aplpy #v2.0.3
CO の積分強度の hdu を作るための関数を適当に定義します。
def make_integ_hdu(hdu):
data = hdu.data
header = hdu.header
new_data = np.nansum(data[:], axis=0)*np.abs(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 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
必要なパラメータを適当に定義します。
map_center = SkyCoord("1:36:41.7863 +15:47:00.801", frame=FK5, unit=(u.hourangle, u.deg))
random_star_ra = map_center.ra.degree + (np.random.random(10)-0.5)*60.0/3600.0
random_star_dec = map_center.dec.degree + (np.random.random(10)-0.5)*60.0/3600.0
integ_hdu = make_integ_hdu(fits.open("./ngc0628_12m+7m+tp_co21.fits")[0])
hdu_color = fits.open("jw02107-o039_t018_miri_f770w_i2d.fits")[1]
以下が aplpy の設定ほぼ全部盛りです。適宜変更してご使用ください。
plt.rcParams['xtick.direction'] = 'out' # inで内側
plt.rcParams['ytick.direction'] = 'out'
#fig, ax = plt.subplots(figsize=(16, 16))
fig = plt.figure(figsize=(16, 16))
f = aplpy.FITSFigure(hdu_color, slices=[0], convention='wells', figure=fig, subplot=[0.1, 0.1, 0.8, 0.8])
f.show_colorscale(vmin=10, vmax=40, stretch='log', cmap="Purples", aspect="equal")
f.recenter(map_center.ra.degree, map_center.dec.degree, width=120.0/3600.0, height=120.0/3600.0)
f.ticks.set_xspacing(0.002) # degree
f.ticks.set_yspacing(0.002)
f.ticks.set_length(5, minor_factor=0.5) # points
f.ticks.set_color('k')
f.ticks.set_linewidth(1.5) # points
f.ticks.set_minor_frequency(2)
f.tick_labels.set_xformat('hh:mm:ss.ss') # FK5
f.tick_labels.set_yformat('dd:mm:ss.ss')
#f.tick_labels.set_xformat('dd.d') # GAL
#f.tick_labels.set_yformat('dd.d')
f.tick_labels.set_font(size=15, family='serif')
f.axis_labels.set_font(size=15, family='serif')
f.add_scalebar(np.rad2deg(np.arcsin(100.0/11e6))) # 距離 11e6 pc で 100 pc の bar
f.scalebar.set_label('100 pc')
f.scalebar.set_font_family('serif')
f.scalebar.set_font_size(15)
f.scalebar.set_color("k")
f.scalebar.set_linewidth(3)
#f.scalebar.set_corner("bottom right")
f.add_colorbar()
f.colorbar.set_axis_label_text(r"f770w (MJy str$^{-1}$)")
f.colorbar.set_axis_label_font(size=20, family='serif')
#f.colorbar.set_ticks([0,1,2,3,4,5])
f.colorbar.set_font(size=20, family='serif')
f.colorbar.set_location("right")
f.add_beam(major=0.3*u.arcsec, minor=0.3*u.arcsec, angle=0.0*u.degree) # header に書いてあれば引数の記述は必要なし
f.beam.set_color('white')
f.beam.set_hatch('+')
f.add_grid()
f.grid.set_color('lightgray')
f.grid.set_alpha(1)
levels = [5, 10, 16, 32]
f.show_contour(integ_hdu, levels=levels, linewidths=1.5, colors="orange", zorder=100) # もしくはファイル名
for l in levels: # 半透明コントア。999999 はとりあえず大きな数。
f.show_contour(integ_hdu, levels=[l, 999999], convention='wells', colors="orange", slices=[0], filled=True, alpha=0.2, zorder=101)
f.show_markers(random_star_ra, random_star_dec, marker="*", c="c", s=100, zorder=102)
f.show_circles(map_center.ra.degree, map_center.dec.degree, 10.0/3600.0, color="w", linestyle="--", linewidth=3, zorder=103)
f.show_rectangles(map_center.ra.degree, map_center.dec.degree, 20.0/3600.0, 40.0/3600.0, angle=0.0, color="w", linestyle=(0, (5, 3, 1, 3, 1, 3)), linewidth=3, zorder=103)
f.add_label(map_center.ra.degree, map_center.dec.degree, "neko", weight="bold", color="m", family="serif", size=24, zorder=104)
x_start, y_start, x_length, y_length = map_center.ra.degree + 20.0/3600.0, map_center.dec.degree + 20.0/3600.0, 20.0/3600.0, 20.0/3600.0
f.show_arrows(x_start, y_start, x_length, y_length, color="k", width=2, head_width=10, head_length=10, zorder=105)
f.show_lines([np.array([[map_center.ra.degree - 20.0/3600.0, map_center.ra.degree - 40.0/3600.0], [map_center.dec.degree - 20.0/3600.0, map_center.dec.degree - 40.0/3600.0]])], color='k', linewidth=2, alpha=0.9, linestyle="-.", zorder=106)
f.set_nan_color("lightgray")
f.set_title("NGC628", color="k", family="serif", size=24, y=1.00)
plt.tight_layout()
plt.rcParams['figure.facecolor'] = "w"# 背景を白にしたい時
plt.rcParams['figure.edgecolor'] = "w"
plt.show()
f.save("qiita_aplpy_all.png", dpi=100) # pdf や eps でも出力可能。ただし、eps はそのままではTeXに張り付けられず、また、半透明部分がバグる。
以上です。
リンク
目次