python の astroquery を用いた宇宙天文カタログの検索と銀河を用いた簡単なプロット方法

astroquery とは

宇宙天文のデータベースは山ほどあって調べるが大変な時代もあったが、今は astroquery からNED、ESAsky、Skyviewなど様々なデータベースへのアクセスが簡単になっている。ここではその使い方と簡単にその内容をプロットして確認する方法まで手短に紹介する。

インストール方法

anaconda3 系を入れていると、pip install astroquery だけで一式が揃う。python2系でも動く部分もあるが、データーベースに依存して動かない部分があるように見受けられるのでpython3系で使う方が良い。

ここで紹介するコードの一覧は、google Colab 上から閲覧できます。

簡単な例

ここでは例として、　An X-ray catalog and atlas of galaxies. Fabbiano+ 1992 を用いる。汎用性のため、テキストファイルへの保存とコラムの追加の方法なども例示的に含めた。

以下、ステップ毎に説明し、最後にコードを出力結果を示す。

1) astroquery.vizierを用いたデータ取得

from astroquery.vizier import Vizier 
v = Vizier(catalog="J/ApJS/80/531/gxfluxes",columns=['Name',"logLx","Bmag","dist","type"],row_limit=-1)
# http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=J/ApJS/80/531/gxfluxes
print(v.columns)
data = v.query_constraints()

で、カタログからデータを取得できる。

• カタログの名前は、"J/ApJS/80/531/gxfluxes"であり、この場合はVizierカタログから与えられた固有名で、カタログ毎に決まっていて、これは事前に取得しておく。
• row_limit = -1 がないと、デフォルトの50個しかダウンロードされない。Vizierオブジェクトを宣言した時点では、設定だけでまだダウンロードはされない。
• query_constraintsという関数を実行して初めてデータがダウンロードされる。返り値のdataはカタログのリストになる。

2) テキストファイルに書き出す

from astropy.io import ascii
ascii.write(data[0],"gxfluxes.txt",delimiter='|',overwrite=True)

でdataに保存された一つ目のカタログがテキストファイルに書き出される。

• delimiter を 有限の文字にしておくこと。天文系のデータは欠損が多いため、空白のdelimiterを使うと行数がずれるやすい。また文字に空白やカンマがある場合もあるので、無難な|を使う。
• dataはカタログリストであり、data[0]でカタログの一つ目のデータ群であり、astropy.table.table.Table という型になる。

astropy.table.table.Tableのまま解析してもよいが、汎用性のためにここでテキストファイルに書き出して、それをpandasで読み込む、というサンプルにしておく。

3) pandas で読み込む

import pandas as pd
df = pd.read_table('gxfluxes.txt', header=0,delimiter="|")

pandasのread_tableはデフォルトで一行目(header=0)をヘッダーとみなす。
この例は一行目がヘッダーなのでなくても良いが、忘れないように明示的に書いておいた。

4) コラムの追加、変更

df["gtype"] = df["type"]
df = df.replace({"gtype":hntype})

新しい行、"gtype" を追加して、その中身は "type" とする。その"gtype"の中身を、ハッブルタイプの文字で変更する。replaceの後に、上書きするコラム名をいれて、変換ルールを記した辞書を与えるだけでOK。

hntype={-8:"transition",-7:"dSph",-5:"E",-3:"S0-",-2:"S0o",-1:"S0+",0:"S0a",1:"Sa",\
2:"Sab",3:"Sb",4:"Sbc",5:"Sc",6:"Scd",7:"Sd",8:"Sdm",9:"Sm",10:"Ir"}

ここでは、ハッブルタイプと一対一対応するTtypeという整数値と文字を変換する例を示した。銀河のハッブルタイプに対応するTtypeについては例えば下記参照。

• Morphological type (numerical): t [http://leda.univ-lyon1.fr/leda/param/t.html]
  • The numerical code for the revised (de Vaucouleurs) morphological type is defined in RC2.

5) matplotlib の scatter でプロットする例

for i, sp in enumerate(set(df["gtype"])):
    c = scalarMap.to_rgba(i)
    df_species = df[df['gtype'] == sp]
    ax.scatter(data=df_species, x='logLx', y='Bmag', label=sp, color=c, s=5)

matplotlibのscatterという関数を使ってプロットする例を示した。pythonのsetという関数は重複を除いたリストを生成する関数で、今のようにtypeが整数値で何種類かあるか具体的にわからない時に重宝する関数。

6) seborn の sns.pairplot でプロットする例

sns.pairplot(df)
sns.pairplot(df,hue='gtype',markers='.')

sns.pairplotは df の中に含まれる数値データの相関を一気にプロットしてくれる便利な関数。文字列はプロットできない。

ただし、カットして色分けする場合には、hueというオプションにコラム名を与えるとカットをかけることができる。

コード

コードは、google Colab 上からも閲覧できます。

qiita_astroquery_gxfluxes.py

#!/usr/bin/env python

from astroquery.vizier import Vizier 
v = Vizier(catalog="J/ApJS/80/531/gxfluxes",columns=['Name',"logLx","Bmag","dist","type"],row_limit=-1)
# http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=J/ApJS/80/531/gxfluxes
print(v.columns)
data = v.query_constraints()

# dump text
from astropy.io import ascii
ascii.write(data[0],"gxfluxes.txt",delimiter='|',overwrite=True)
# type(data[0])                                                                                          
# Out[7]: astropy.table.table.Table

import pandas as pd
df = pd.read_table('gxfluxes.txt', header=0,delimiter="|")

hntype={-8:"transition",-7:"dSph",-5:"E",-3:"S0-",-2:"S0o",-1:"S0+",0:"S0a",1:"Sa",\
2:"Sab",3:"Sb",4:"Sbc",5:"Sc",6:"Scd",7:"Sd",8:"Sdm",9:"Sm",10:"Ir"}
#http://leda.univ-lyon1.fr/leda/param/t.html

df["gtype"] = df["type"]
df = df.replace({"gtype":hntype})

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.rcParams['font.family'] = 'serif'
import matplotlib.colors as colors

fig = plt.figure(figsize=(10,7.))
ax = fig.add_subplot(1, 1, 1)

num=len(set(df["gtype"]))
usercmap = plt.get_cmap('jet')
cNorm  = colors.Normalize(vmin=0, vmax=num)
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=usercmap)

for i, sp in enumerate(set(df["gtype"])):
    c = scalarMap.to_rgba(i)
    df_species = df[df['gtype'] == sp]
    ax.scatter(data=df_species, x='logLx', y='Bmag', label=sp, color=c, s=5)

plt.xlabel("logLx")
plt.ylabel("Bmag")
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=8, mode="expand", borderaxespad=0.,fontsize=12)
plt.savefig("gxfluxes_scat.png")
plt.show()


# plot pair plot. Note that the column of "string" is omitted. 
import seaborn as sns
sns.pairplot(df)
plt.savefig("gxfluxes_pairplot1.png")

sns.pairplot(df,hue='gtype',markers='.')
plt.savefig("gxfluxes_pairplot2.png")

出力結果

matplotlibで生成されたscatter plot

sns.pairplotのオプションなしで生成された相関図

sns.pairplotのhueオプションで生成されたカテゴリ分けされた相関図

別のカタログの情報と結合してプロットする方法

OpenNGCのデータとの結合方法の例

ここではOpenNGCのデータと、上記で取得したデータを結合する例を示す。
再現性のため、下記のサンプルで用いるテキストファイルも置いておく。OpenNGCについては、
pythonでOpenNGCデータベースを用いて銀河の可視光データをプロットする方法を参照。
ここでは、上で http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=J/ApJS/80/531/gxfluxes にあるAn X-ray catalog and atlas of galaxies. (Fabbiano+, 1992) から抽出した gxfluxes.txt を用いる。

• http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/gxfluxes.txt
• http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/NGC.csv

中身はどちらもNameというコラムがあり、天体名が入っている。

>head -2 gxfluxes.txt                                     
Name|logLx|Bmag|dist|type
I0342|40.37|7.86|6.300|6

>head -2 NGC.csv
Name;Type;RA;Dec;Const;MajAx;MinAx;PosAng;B-Mag;V-Mag;J-Mag;H-Mag;K-Mag;SurfBr;Hubble;Cstar U-Mag;Cstar B-Mag;Cstar V-Mag;M;NGC;IC;Cstar Names;Identifiers;Common names;NED notes;OpenNGC notes
IC0001;**;00:08:27.05;+27:43:03.6;Peg;;;;;;;;;;;;;;;;;;;;;

これをNameを軸にして結合するには、

df1 = pd.read_table('NGC.csv', header=0,delimiter=";")
df2 = pd.read_table('gxfluxes.txt', header=0,delimiter="|")
df2['Name'] = df2['Name'].str.replace('N', 'NGC').str.replace('I', 'IC')
df12 = pd.merge(df1, df2, on='Name',indicator="indicator")

とするだけでOKである。df2のNameの方は、NをNGCにreplaceして、IをICにreplaceして、df1に合うように調整している。その後で、mergeをする。onというオプションにどのコラムを軸に結合するかを明記うする。indicatorはなくても良いが、和集合か積集合か確認するために使われる。

サンプルコード

上からNGC.csvとgxfluxes.txtをダウンロードして下記を実行すると、２つをNameを軸にマージできる。参考までに、その結果として生成可能になるX線と可視、近赤のプロットの例も示した。

qiita_merge.py

#!/usr/bin/env python

import pandas as pd
df1 = pd.read_table('NGC.csv', header=0,delimiter=";")
df2 = pd.read_table('gxfluxes.txt', header=0,delimiter="|")
df2['Name'] = df2['Name'].str.replace('N', 'NGC').str.replace('I', 'IC')
df12 = pd.merge(df1, df2, on='Name',indicator="indicator")

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.rcParams['font.family'] = 'serif'
import matplotlib.colors as colors

fig = plt.figure(figsize=(14,7.))
ax = fig.add_subplot(1, 2, 1)
plt.ylabel("logLx")
plt.xlabel("B-mag")
ax.scatter(data=df12[df12["type"]<0],  x='B-Mag', y='logLx', label="E+S0", s=4)
ax.scatter(data=df12[df12["type"]>=0], x='B-Mag', y='logLx', label="S", s=4)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left', borderaxespad=0.,fontsize=8)

ax = fig.add_subplot(1, 2, 2)
plt.ylabel("logLx")
plt.xlabel("K-mag")
ax.scatter(data=df12[df12["type"]<0],  x='K-Mag', y='logLx', label="E+S0", s=4)
ax.scatter(data=df12[df12["type"]>=0], x='K-Mag', y='logLx', label="S", s=4)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left', borderaxespad=0.,fontsize=8)

plt.savefig("fig_BKX.png")
plt.show()

生成される図

２つのデータベースを結合したことで生成できるX線とBバンド、Kバンドの相関図を例に示した。ハッブルタイプに正負で場合分けしている。

統計情報の確認

Nameの書き方や、ちょっとした記法によってはうまく働かない場合がある。結合後には、ちゃんと意図したように結合できているか、info()で確認しよう。

>df12.info()                                                                                                            
<class 'pandas.core.frame.DataFrame'>
Int64Index: 444 entries, 0 to 443
Data columns (total 31 columns):
Name             444 non-null object
Type             444 non-null object
RA               444 non-null object
Dec              444 non-null object
Const            444 non-null object
MajAx            441 non-null float64
MinAx            436 non-null float64
PosAng           436 non-null float64
B-Mag            437 non-null float64
V-Mag            276 non-null float64
J-Mag            427 non-null float64
H-Mag            427 non-null float64
K-Mag            428 non-null float64
SurfBr           436 non-null float64
Hubble           436 non-null object
Cstar U-Mag      0 non-null float64
Cstar B-Mag      0 non-null float64
Cstar V-Mag      0 non-null float64
M                45 non-null float64
NGC              21 non-null object
IC               15 non-null object
Cstar Names      0 non-null object
Identifiers      441 non-null object
Common names     31 non-null object
NED notes        86 non-null object
OpenNGC notes    9 non-null object
logLx            443 non-null float64
Bmag             444 non-null float64
dist             444 non-null float64
type             444 non-null int64
indicator        444 non-null category
dtypes: category(1), float64(16), int64(1), object(13)
memory usage: 108.1+ KB

天体の見かけの等級を絶対等級に変換する方法

天体の見かけの明るさと真の明るさは距離がわかってないと変換できない。
距離が分かっている場合は簡単な演算で計算できる。( https://physmemo.shakunage.net/phys/magnitude/magnitude.html
など参考。) 注意点としては、一般的に天体の距離を計測するのは難しいため、距離には誤差が必ず含まれるが、データベースには誤差までは含まれてないことが多い。天体の物理量の絶対値で議論する方が理論との比較はしやすいが、天文学においては距離を介することで不定性が増すことは意識して解析したい。(理論屋さんは観測値ではなくて絶対値を知りたいのであるが、観測屋は観測値の方を明記する理由の一つである。)

absB = df12["B-Mag"] + 5 - 5 * np.log10(df12["dist"] * 1e6)
df12["absB"] = absB
bv = df12["B-Mag"] - df12["V-Mag"]
df12["bv"] = bv

変換した後は、コラムを追加しておこう。コラムの追加は、Dataframe["追加したいコラム名"]=配列、で追加できる。

サンプルコード

#!/usr/bin/env python

import numpy as np
import pandas as pd
df1 = pd.read_table('NGC.csv', header=0,delimiter=";")
df2 = pd.read_table('gxfluxes.txt', header=0,delimiter="|")
df2['Name'] = df2['Name'].str.replace('N', 'NGC').str.replace('I', 'IC')
df12 = pd.merge(df1, df2, on='Name',indicator="indicator")

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.rcParams['font.family'] = 'serif'
import matplotlib.colors as colors

absB = df12["B-Mag"] + 5 - 5 * np.log10(df12["dist"] * 1e6)
df12["absB"] = absB
bv = df12["B-Mag"] - df12["V-Mag"]
df12["bv"] = bv

df12 = df12.loc[:,["bv","absB","type","dist","B-Mag"]]
df12 = df12.dropna(how='any')

fig = plt.figure(figsize=(12,8.))
ax = fig.add_subplot(2, 2, 1)
plt.ylabel("relative B-Mag")
plt.xlabel("B-V")
ax.scatter(data=df12[df12["type"]<0],  x='bv', y='B-Mag', label="E+S0", s=4)
ax.scatter(data=df12[df12["type"]>=0], x='bv', y='B-Mag', label="S", s=4)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left', borderaxespad=0.,fontsize=8, ncol = 8)

ax = fig.add_subplot(2, 2, 3)
plt.ylabel("absolute B-Mag")
plt.xlabel("B-V")
ax.scatter(data=df12[df12["type"]<0],  x='bv', y='absB', label="E+S0", s=4)
ax.scatter(data=df12[df12["type"]>=0], x='bv', y='absB', label="S", s=4)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left', borderaxespad=0.,fontsize=8, ncol = 8)

ax = fig.add_subplot(2, 2, 2)
plt.ylabel("relative B-Mag")
plt.xlabel("distance [Mpc]")
ax.scatter(data=df12[df12["type"]<0],  x='dist', y='B-Mag', label="E+S0", s=4)
ax.scatter(data=df12[df12["type"]>=0], x='dist', y='B-Mag', label="S", s=4)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left', borderaxespad=0.,fontsize=8, ncol = 8)

ax = fig.add_subplot(2, 2, 4)
plt.ylabel("absolute B-Mag")
plt.xlabel("distance [Mpc]")
ax.scatter(data=df12[df12["type"]<0],  x='dist', y='absB', label="E+S0", s=4)
ax.scatter(data=df12[df12["type"]>=0], x='dist', y='absB', label="S", s=4)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left', borderaxespad=0.,fontsize=8, ncol = 8)

plt.savefig("fig_CM.png")
plt.show()

生成されるCM図

天文用語でCM図とは、フラックスの比と光度の相関図のことである。(color-Magnitude図という意味。たぶん)
このような距離補正を行なったときは、それが正しく変換できているか、その前後のパラメータをプロットして確認しよう。