4
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

天文データ解析入門 その25 (aplpy.FITSFigure 設定ほぼ全部盛り)

Last updated at Posted at 2023-04-17

[以下追記] 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に張り付けられず、また、半透明部分がバグる。

qiita_aplpy_all.png

以上です。

リンク
目次

4
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?