本記事では、カタログの扱い方について紹介します。astroquery や pandas などを用います。numpy (loadtxt) だけでも一応大丈夫ですが、あまりメリットがないので pandas を使いましょう。
インストール済みでない方は
pip install astroquery
pip install pandas
を実行しましょう。
今回、例として、Spitzer 8.0 µm のアーカイブデータを使用します。
https://irsa.ipac.caltech.edu/data/SPITZER/GLIMPSE/images/I/1.2_mosaics_v3.5/GLON_30-53/
から GLM_03000+0000_mosaic_I4.fits をダウンロードしました。
astroquery
table の作成
よく用いられる Simbad についての例を紹介します。
まず必要なものを import します。
from astroquery.simbad import Simbad
from astropy import coordinates
import astropy.units as u
from matplotlib import pyplot as plt
from astropy.io import fits
import aplpy
W43 を取得します。
result_table = Simbad.query_object("W43")
print(result_table)
# MAIN_ID RA DEC ... COO_WAVELENGTH COO_BIBCODE
# "h:m:s" "d:m:s" ...
#--------------- -------- -------- ... -------------- -------------------
#SNR G030.8-00.0 18 47 32 -01 56.5 ... O 2013A&A...558A..53K
W40 - W49 までを取得したい場合は、wildcard で指定できます。
result_table = Simbad.query_object("W 4[0-9]", wildcard=True) # なぜかこっちはスペースがいる
print(result_table)
# MAIN_ID RA ... COO_WAVELENGTH COO_BIBCODE
# "h:m:s" ...
#------------------ ------------- ... -------------- -------------------
# LBN 028.77+03.43 18 31 26.5 ... R 1992ApJS...79..331W
# SNR G023.3-00.3 18 34 51.80 ... 2006ApJ...651..190L
#GRS G024.80 +00.10 18 38 14.8 ... 1984A&AS...58..197R
# SNR G030.8-00.0 18 47 32 ... O 2013A&A...558A..53K
# SNR G034.6-00.5 18 56 10.650 ... 2015ApJS..217....2M
# CTB 61 18 56.4 ... R 2012MNRAS.422.2429A
# IRAS 18549+0913 18 57 20.7234 ... O 2018yCat.1345....0G
# CTB 64 19 00 59 ...
# SNR G035.2-01.8 19 01.6 ... R 2012MNRAS.422.2429A
# W 49 19 10 20 ... 2003A&A...397..133R
次は、W43 の周囲 0.5 度の天体を取得します。
result_table = Simbad.query_region("W43", radius=0.5 * u.deg)
print(result_table)
# MAIN_ID RA ... COO_WAVELENGTH COO_BIBCODE
# "h:m:s" ...
#------------------------- ------------- ... -------------- -------------------
# SNR G030.8-00.0 18 47 32 ... O 2013A&A...558A..53K
# SDC G30.762-0.015 18 47 32.06 ... I 2009A&A...505..405P
#SSTOERC G030.7626-00.0243 18 47 33.97 ... M 2017ApJ...839..108S
#SSTOERC G030.7566-00.0121 18 47 30.71 ... M 2017ApJ...839..108S
# [HBM2005] G30.760-0.027 18 47 34.2 ... m 2005MNRAS.363..405H
#SSTOERC G030.7689-00.0238 18 47 34.57 ... M 2017ApJ...839..108S
# [CPA2006] N52 18 47 31 ... 2006ApJ...649..759C
#SSTOERC G030.7478-00.0185 18 47 31.11 ... M 2017ApJ...839..108S
# AGAL G030.763-00.031 18 47 35.2 ... s 2015ApJS..218....1M
# [MSL2003] W43-MM45 18 47 34.4 ... 2003ApJ...582..277M
# ... ... ... ... ...
# 2MASS J18454366-0209104 18 45 43.6263 ... O 2018yCat.1345....0G
# AGAL G031.254+00.057 18 48 10.2 ... s 2015ApJS..218....1M
#SSTOERC G030.4851-00.4363 18 48 31.64 ... M 2017ApJ...839..108S
#SSTOERC G030.3573-00.3154 18 47 51.78 ... M 2017ApJ...839..108S
#SSTOERC G031.2587-00.0357 18 48 30.75 ... M 2017ApJ...839..108S
#SSTOERC G031.2505+00.0732 18 48 06.59 ... M 2017ApJ...839..108S
# BGPS G030.315+00.212 18 45 54.5631 ... m 2010ApJS..188..123R
# BGPS G030.512-00.453 18 48 38.0883 ... m 2010ApJS..188..123R
#SSTOERC G030.2889-00.1880 18 47 17.05 ... M 2017ApJ...839..108S
#SSTOERC G030.2597-00.0367 18 46 41.50 ... M 2017ApJ...839..108S
#Length = 6684 rows
# もちろん座標指定でも取得できます。
result_table = Simbad.query_region(coordinates.SkyCoord(l=30.7593, b=-0.0186, frame='galactic', unit=u.deg), radius=0.5 * u.deg)
print(result_table)
# MAIN_ID RA ... COO_WAVELENGTH COO_BIBCODE
# "h:m:s" ...
#------------------------- ------------- ... -------------- -------------------
# SNR G030.8-00.0 18 47 32 ... O 2013A&A...558A..53K
# SDC G30.762-0.015 18 47 32.06 ... I 2009A&A...505..405P
#SSTOERC G030.7626-00.0243 18 47 33.97 ... M 2017ApJ...839..108S
#SSTOERC G030.7566-00.0121 18 47 30.71 ... M 2017ApJ...839..108S
# [HBM2005] G30.760-0.027 18 47 34.2 ... m 2005MNRAS.363..405H
#SSTOERC G030.7689-00.0238 18 47 34.57 ... M 2017ApJ...839..108S
# [CPA2006] N52 18 47 31 ... 2006ApJ...649..759C
#SSTOERC G030.7478-00.0185 18 47 31.11 ... M 2017ApJ...839..108S
# AGAL G030.763-00.031 18 47 35.2 ... s 2015ApJS..218....1M
# [MSL2003] W43-MM45 18 47 34.4 ... 2003ApJ...582..277M
# ... ... ... ... ...
# 2MASS J18454366-0209104 18 45 43.6263 ... O 2018yCat.1345....0G
# AGAL G031.254+00.057 18 48 10.2 ... s 2015ApJS..218....1M
#SSTOERC G030.4851-00.4363 18 48 31.64 ... M 2017ApJ...839..108S
#SSTOERC G030.3573-00.3154 18 47 51.78 ... M 2017ApJ...839..108S
#SSTOERC G031.2587-00.0357 18 48 30.75 ... M 2017ApJ...839..108S
#SSTOERC G031.2505+00.0732 18 48 06.59 ... M 2017ApJ...839..108S
# BGPS G030.315+00.212 18 45 54.5631 ... m 2010ApJS..188..123R
#SSTOERC G030.2889-00.1880 18 47 17.05 ... M 2017ApJ...839..108S
#SSTOERC G030.2597-00.0367 18 46 41.50 ... M 2017ApJ...839..108S
# BGPS G030.512-00.453 18 48 38.0883 ... m 2010ApJS..188..123R
#Length = 6684 rows
# 論文の bibcode からも取ってこれます。
Simbad.TIMEOUT = 300 # 天体数が多いと timeout するので、長くしておく
result_table = Simbad.query_bibobj('2017ApJ...839..108S') # Saral et al. 2017
# MAIN_ID RA ... COO_WAVELENGTH COO_BIBCODE
# "h:m:s" ...
#------------------------- ------------- ... -------------- -------------------
# IC 1848 02 51 06 ... O 2013A&A...558A..53K
# NAME G305 Complex 13 13 ... 2012MNRAS.419.1871D
# NAME RCW 106 GMC 16 21.0 ... R 2006MNRAS.367.1609B
#SSTOERC G029.0464+00.3218 18 43 11.73 ... M 2017ApJ...839..108S
#SSTOERC G030.5437+01.0894 18 43 12.06 ... M 2017ApJ...839..108S
#SSTOERC G030.8961+01.2696 18 43 12.20 ... M 2017ApJ...839..108S
#SSTOERC G030.5271+01.0799 18 43 12.25 ... M 2017ApJ...839..108S
# 2MASS J18431248-0213585 18 43 12.4780 ... O 2018yCat.1345....0G
#SSTOERC G029.2830+00.4386 18 43 12.74 ... M 2017ApJ...839..108S
#SSTOERC G029.0774+00.3319 18 43 12.96 ... M 2017ApJ...839..108S
# ... ... ... ... ...
#SSTOERC G049.2320-01.0191 19 25 31.43 ... M 2017ApJ...839..108S
#SSTOERC G049.2997-00.9834 19 25 31.60 ... M 2017ApJ...839..108S
# 2MASS J19253166+1417349 19 25 31.667 ... I 2003yCat.2246....0C
#SSTOERC G049.2291-01.0225 19 25 31.82 ... M 2017ApJ...839..108S
#SSTOERC G050.0326-00.5913 19 25 32.24 ... M 2017ApJ...839..108S
#SSTOERC G049.4660-00.9020 19 25 33.39 ... M 2017ApJ...839..108S
#SSTOERC G049.7610-00.7431 19 25 33.42 ... M 2017ApJ...839..108S
# NAME Cyg X 20 28 41 ... G 2013ApJS..209...34A
# NAME Sct-Cen ...
# NAME Galactic Bar ...
#Length = 15673 rows
SNR だけ欲しい場合など、こう言う書き方もあります。詳しくはドキュメントを参照。
result = Simbad.query_criteria('region(box, GAL, 0 +0, 3d 1d)', otype='SNR')
# もしくは
script = '(region(box, GAL, 0 +0, 0.5d 0.5d) | region(box, GAL, 43.3 -0.2, 0.25d 0.25d))'
result = Simbad.query_criteria(script, otype='SNR')
print(result)
# MAIN_ID RA ... COO_WAVELENGTH COO_BIBCODE
# "h:m:s" ...
#----------------- ------------ ... -------------- -------------------
# SNR G359.9-00.9 17 45.8 ...
# NAME Sgr A East 17 45 41 ... 2010ApJS..188..405A
# W 49b 19 11 09.000 ... 2015ApJS..217....2M
#SNR G000.13-00.12 17 46.4 ... 2013MNRAS.434.1339H
table の使い方
使いやすいように、coordinates.SkyCoord のオブジェクトにします。
ra, dec = result_table["RA"], result_table["DEC"]
c = coordinates.SkyCoord(ra=ra, dec=dec, frame="fk5", unit=(u.hourangle, u.deg))
# ValueError: Syntax error parsing angle ''
何かエラー出ました。おそらく座標が適切でないものが混じっているっぽいです。たぶん後ろの 2 行だと思います。
ra, dec = result_table["RA"][:-2], result_table["DEC"][:-2]
c = coordinates.SkyCoord(ra=ra, dec=dec, frame="fk5", unit=(u.hourangle, u.deg))
うまくいきました。
これを Spitzer 8.0 µm の上に plot してみます。
hdu = fits.open("~/your/fits/dir/GLM_03000+0000_mosaic_I4.fits")[0]
fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(vmin=10, vmax=1000, stretch='log', cmap="Greens", aspect="equal")
f.recenter(30.5, 0.0, width=1.0, height=1.0)
f.show_markers(c.galactic.l.deg, c.galactic.b.deg, s=2)
plt.show()
何かカタログの中の値 (result_table.columns で確認可能) で mask したい場合などは、
fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(vmin=10, vmax=1000, stretch='log', cmap="Greens", aspect="equal")
f.recenter(30.5, 0.0, width=1.0, height=1.0)
mask = result_table["COO_QUAL"][:-2]=="A"
f.show_markers(c.galactic.l.deg[mask], c.galactic.b.deg[mask], s=2)
plt.show()
COO_QUAL が A のものだけを plot しました。
自前の csv ファイルを読み取る
手元にある天体のリストのテキストファイルを python で読み取る場合について紹介します。
手っ取り早いのは、excel か何かに貼って、重複する column 名を変更して (pandas なら変更しなくても ok)、テキストエディタにコピペして textfile で保存し、それを読みこむ方法です。
例として、
https://ui.adsabs.harvard.edu/abs/2017yCat..18390108S/abstract
の mysos を取ってきて Saral_test.txt と言う名前で保存しました。
import pandas as pd
import numpy as np
# delimiter で区切りを指定します。タブ (\t) 区切りの場合は
Saral_df = pd.read_csv("Saral_test.txt", delimiter="\t")
numpy array にしたい場合などは、
ra, dec, lumi, mass = np.array(Saral_df["RAJ2000"]), np.array(Saral_df["DEJ2000"]), np.array(Saral_df["Lum"]), np.array(Saral_df["Mass"])
c = coordinates.SkyCoord(ra=ra, dec=dec, frame="fk5", unit=(u.deg, u.deg))
例えばluminosity で大きさを変えて、mass で色を変えて plot してみます。
fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(vmin=10, vmax=1000, stretch='log', cmap="Greys", aspect="equal")
mask = Saral_df["Reg"]=="W43"
f.show_markers(c.galactic.l.deg[mask], c.galactic.b.deg[mask], s=lumi[mask]*30, c=mass[mask], cmap="jet")
plt.show()
ちなみに、論文のテーブルがうまく選択できなくてコピペできない場合などは、Google Chrome や Adobe Acrobat Reader 等を使えば選択できる可能性があります。また、論文が古すぎて画像しかない場合は、文字起こしソフト (Google にもあった気がする) を使えば、手打ちするよりは早いと思います。
以上です。
リンク
目次