こんな人にオススメ
天文学で座標系を変換したいんだけど、自分で計算するのは面倒だし、本に書いてある式を実装するのも面倒。何かパッとすぐに変換できる方法ってないの？
ということで、今回は究極にニッチな内容である「astropyを用いて天文座標系を変換する」ということについて解説する。というのも執筆者は天文学系の研究を行なっており、時折ある座標系から他の座標系への変換が必要な時があった。その時にいちいちサイトで座標を打って計算してもらってコピペして、というのが面倒すぎる。
そこで今回のようにpythonを用いて座標変換を行うことで自動的にすぐに値を取り出そうということだ。初めて知った時は衝撃的だった。
執筆者の人生的なものについては以下参照。
-
-
【M天パ（めがてんぱ）自己紹介】変わって楽しいの繰り返し
続きを見る
衝撃を受けたastropy関連の記事については以下参照。
-
-
【astropy&単位】astropyで変数に単位を付与
続きを見る
-
-
【astropy&単位】astropyの単位出力書式
続きを見る
python環境は以下。
- Python 3.9.2
- astropy 4.2
- numpy 1.20.1
- pandas 1.2.3
今回は赤道座標を黄道座標、銀河座標に変換する。
座標として読み込む
基本的にはastropyのレファレンスに書いてある。
が、個人的にすごく読みにくい。したがって、サイトに書いてあるコードをここでも紹介することとする。
1つの座標を読み込む
座標を天文座標として読み込ませるために入力する方式は以下の6種類。
（数値）*u.degreeで角度の単位を持つ値として入力（数値）と引数でunit='deg'で角度を指定して入力'00h42m30s'といったようにhmsで入力'00h42.5m'といったようにhmで入力'00 42 30 +41 12 00'といったふうにhmsの間をスペースで区切ってかつ、引数unitでそれぞれの角度を指定して入力'00:42.5 +41:12'といったふうにh:mで入力かつ、引数unitでそれぞれの角度を指定して入力
どれでも結果は同じ。以下に例を示す。
import numpy as np
import pandas as pd
from astropy.coordinates import SkyCoord
from astropy import units as u
lst = [
SkyCoord(ra=10.625 * u.degree, dec=41.2 * u.degree, frame='icrs'),
SkyCoord(10.625, 41.2, frame='icrs', unit='deg'),
SkyCoord('00h42m30s', '+41d12m00s', frame='icrs'),
SkyCoord('00h42.5m', '+41d12m'),
SkyCoord('00 42 30 +41 12 00', unit=(u.hourangle, u.deg)),
SkyCoord('00:42.5 +41:12', unit=(u.hourangle, u.deg)),
]
for c in lst:
print(c)
# <SkyCoord (ICRS): (ra, dec) in deg
# (10.625, 41.2)>
# <SkyCoord (ICRS): (ra, dec) in deg
# (10.625, 41.2)>
# <SkyCoord (ICRS): (ra, dec) in deg
# (10.625, 41.2)>
# <SkyCoord (ICRS): (ra, dec) in deg
# (10.625, 41.2)>
# <SkyCoord (ICRS): (ra, dec) in deg
# (10.625, 41.2)>
# <SkyCoord (ICRS): (ra, dec) in deg
# (10.625, 41.2)>
全て(10.625, 41.2)として読み込まれていることがわかる。値の取り出しについては後述。
配列として読み込む
先ほどは単発でデータを読み込ませたが、配列の中に入れて複数個同時に読み込むということも可能。
# 配列でも大丈夫
c = SkyCoord(ra=[10, 11, 12, 13] * u.degree, dec=[41, -5, 42, 0] * u.degree)
print(c)
# <SkyCoord (ICRS): (ra, dec) in deg
# [(10., 41.), (11., -5.), (12., 42.), (13., 0.)]>
c = SkyCoord(
ra=['00h42m30s', '10h42m30s', '20h42m30s'],
dec=['+41d12m00s', '+51d12m00s', '+61d12m00s']
)
print(c)
# <SkyCoord (ICRS): (ra, dec) in deg
# [( 10.625, 41.2), (160.625, 51.2), (310.625, 61.2)]>
c = SkyCoord(
['00:42.5 +41:12', '10:42.5 +51:12', '20:42.5 +61:12'],
unit=(u.hourangle, u.deg)
)
print(c)
# <SkyCoord (ICRS): (ra, dec) in deg
# [( 10.625, 41.2), (160.625, 51.2), (310.625, 61.2)]>
座標変換
前章で座標を読み込むという部分を行なった。ここでは本題の座標変換について解説する。
変換可能な座標系
まずは変換可能な座標系について。さっき使用したのはframe='ircs'のicrsいうもの。これは大雑把には赤道座標系。その他にも複数種類のframeが用意されている。確認方法は読み込みを失敗すること。以下のグラフでは例えばicrsを大文字でICRSと書くことでエラー発生させている。
# frameを確認
print(SkyCoord('00h42m30s', '+41d12m00s', frame='ICRS'))
# raise ValueError('Coordinate frame name "{}" is not a known '
# ValueError: Coordinate frame name "ICRS" is not a known coordinate frame (['altaz', 'barycentricmeanecliptic', 'barycentrictrueecliptic', 'cirs', 'custombarycentricecliptic', 'fk4', 'fk4noeterms', 'fk5', 'galactic', 'galacticlsr', 'galactocentric', 'gcrs', 'geocentricmeanecliptic', 'geocentrictrueecliptic', 'hcrs', 'heliocentriceclipticiau76', 'heliocentricmeanecliptic', 'heliocentrictrueecliptic', 'icrs', 'itrs', 'lsr', 'lsrd', 'lsrk', 'precessedgeocentric', 'supergalactic', 'teme', 'tete'])
27種類もあるのだが、今回使用するのはbarycentricmeaneclipticとgalacticの2つ。前者が黄道座標で後者が銀河座標に相当する。が、eclipticと名のつく項目が8項目もあり、さらにgalacと名のつく項目が4項目もある。以下にその「astropy.coordinates Package」で示されている各項目の説明を引用する。
eclipticbarycentricmeanecliptic：Barycentric mean ecliptic coordinates.barycentrictrueecliptic：Barycentric true ecliptic coordinates.custombarycentricecliptic：Barycentric ecliptic coordinates with custom obliquity.geocentricmeanecliptic：Geocentric mean ecliptic coordinates.geocentrictrueecliptic：Geocentric true ecliptic coordinates.heliocentriceclipticiau76：Heliocentric mean (IAU 1976) ecliptic coordinates.heliocentricmeanecliptic：Heliocentric mean ecliptic coordinates.heliocentrictrueecliptic：Heliocentric true ecliptic coordinates.
galacgalactic：A coordinate or frame in the Galactic coordinate system.galacticlsr：A coordinate or frame in the Local Standard of Rest (LSR), axis-aligned to the Galactic frame.galactocentric：A coordinate or frame in the Galactocentric system.supergalactic：Supergalactic Coordinates (see Lahav et al. 2000, https://ui.adsabs.harvard.edu/abs/2000MNRAS.312..166L, and references therein).
galacticに関してはシンプルにgalacticでいったが、eclipticに関してはエラーになるものもあった。例えば、heliocentrictrueeclipticは以下のエラーとなる。
raise UnitsError(_NEED_ORIGIN_HINT.format(from_coo.__class__.__name__)) astropy.coordinates.errors.UnitsError: The input ICRS coordinates do not have length units. This probably means you created coordinates with lat/lon but no distance. Heliocentric<->ICRS transforms cannot function in this case because there is an origin shift.
今回はエラーの出なかったbarycentricmeaneclipticを使用してコードを作成した。barycentricmeaneclipticは以下のような説明がされている。
Barycentric mean ecliptic coordinates. These origin of the coordinates are the barycenter of the solar system, with the x axis pointing in the direction of the mean (not true) equinox as at the time specified by the equinox attribute (as seen from Earth), and the xy-plane in the plane of the ecliptic for that date.
座標の原点は太陽系の重心である、重心平均黄道座標系のようだ。まあそれぞれの状況で最適なものを使用すればいいだろう。
座標変換
実際に座標変換をするための下準備として以下のコードを作成した。
def trans(ra, dec, ra_unit='deg', dec_unit='deg'):
"""入力した赤道座標の赤経・赤緯の座標とその単位をSkyCoordの方に変換する
Parameters
----------
ra : float or str
赤経の値。deg単位でもhms表記でも大丈夫。
dec : float or str
赤緯の値。deg単位でもhms表記でも大丈夫
ra_unit : str or astropy.units.core.Unit, optional
赤経の単位, by default 'deg'
dec_unit : str or astropy.units.core.Unit, optional
赤緯の単位, by default 'deg'
Returns
-------
astropy.coordinates.sky_coordinate.SkyCoord
SkyCoord用に変換された座標系
"""
# RA, DECから任意の座標へ変換する関数
result = SkyCoord(ra, dec, frame='icrs', unit=(ra_unit, dec_unit))
return result
引数として入力した赤経・赤緯の値を座標系として読み込む、上で行なった作業のdef版だ。実際に使用してみる。
RA = [0, 30, 90, 153, 360] # RAのdeg表記 DEC = [-90, -30, 0, 60, 90] # DECのdeg表記 # 赤道座標で読み込み equatorial = trans(RA[0], DEC[0]) print(equatorial) # <SkyCoord (ICRS): (ra, dec) in deg # (0., -90.)> print(type(equatorial)) # <class 'astropy.coordinates.sky_coordinate.SkyCoord'>
これを銀河座標に変換する。.galacticとすることで変換完了。黄道座標も同じように計算することができる。
# 銀河座標に変換 galactic = equatorial.galactic print(galactic) # <SkyCoord (Galactic): (l, b) in deg # (302.93192526, -27.12825241)> print(type(galactic)) # <class 'astropy.coordinates.sky_coordinate.SkyCoord'>
銀経・銀緯の抽出は.l、.bで、黄経・黄緯の抽出は.lon、.latで行う。
# 銀経を計算 print(galactic.l) # 302d55m54.9309s print(type(galactic.l)) # <class 'astropy.coordinates.angles.Longitude'> # 銀緯を計算 print(galactic.b) # -27d07m41.7087s print(type(galactic.b)) # <class 'astropy.coordinates.angles.Latitude'>
degの角度表示したい場合は.degにする。
# 銀経を角度表示 print(galactic.l.deg) # 302.9319252554167 print(type(galactic.l.deg)) # <class 'numpy.float64'> # 銀緯を計算 print(galactic.b.deg) # -27.12825241496801 print(type(galactic.b.deg)) # <class 'numpy.float64'>
銀経・銀緯をバラさずにそのまま角度表示するには.to_string()を使用する。この時、引数として以下の3種類が用意されている。
dms: 経度、緯度ともにdms表記をするhmsdms: 経度はhms、緯度はdms表記をするdecimal: 「°」のままスペースを使って区切るだけ
全て文字列として出力される。
# 銀経・銀緯のセットのまま値を取り出す
print(galactic.to_string('dms'))
# 302d55m54.9309s -27d07m41.7087s
print(type(galactic.to_string('dms')))
# <class 'str'>
print(galactic.to_string('hmsdms'))
# 20h11m43.6621s -27d07m41.7087s
print(type(galactic.to_string('hmsdms')))
# <class 'str'>
print(galactic.to_string('decimal'))
# 302.932 -27.1283
print(type(galactic.to_string('decimal')))
# <class 'str'>
座表変換の確認
変換した座標の計算が果たして合っているのかということを確認するために執筆者は「NASA/IPAC EXTRAGALACTIC DATABASE Coordinate Transformation & Galactic Extinction Calculator」というサイトを使用している。このサイトに座標を入れることで座標を変換してくれる。
このサイトで計算する際に「Calculate」を押すのだが、新規タブで開くなどして計算元を残しておかないと、変換後のページの「Back to NED Home」で戻ると大元のサイトに戻ってしまうので注意。
座標変換を一気に行う
最後に座標変換を一気に行う方法について紹介する。基本的には配列にしておけば大丈夫だ。
シンプルに配列の座標を変換
まずはシンプルに配列の座標を変換する方法。これは初めに行なった、配列のまま変換するもの。
# 一気に座標変換
equatorial = trans(RA, DEC)
print(equatorial)
# <SkyCoord (ICRS): (ra, dec) in deg
# [( 0., -90.), ( 30., -30.), ( 90., 0.), (153., 60.), ( 0., 90.)]>
# 銀河座標も大丈夫
galactic = equatorial.galactic
print(galactic.to_string('dms'))
# ['302d55m54.9309s -27d07m41.7087s', '227d46m29.4597s -74d41m23.7897s', '206d59m20.8719s -11d25m28.1675s', '151d12m28.9474s 47d25m56.4996s', '122d55m54.9309s 27d07m41.7087s']
print(galactic.to_string('hmsdms'))
# ['20h11m43.6621s -27d07m41.7087s', '15h11m05.964s -74d41m23.7897s', '13h47m57.3915s -11d25m28.1675s', '10h04m49.9298s +47d25m56.4996s', '08h11m43.6621s +27d07m41.7087s']
print(galactic.to_string('decimal'))
# ['302.932 -27.1283', '227.775 -74.6899', '206.989 -11.4245', '151.208 47.4324', '122.932 27.1283']
# 銀経だけ、銀緯だけも大丈夫
print(galactic.l.deg)
# [302.93192526 227.77484991 206.98913108 151.20804095 122.93192526]
print(galactic.b.deg)
# [-27.12825241 -74.68994158 -11.42449097 47.43236101 27.12825241]
銀経・銀緯のままでも銀経だけ、銀緯だけとバラしても変換することが可能。
変換した座標系をpandasのDataFrameに
最後は黄道座標、銀河座標に変換した座標をpandasのDataFrameに入れて整理する方法。ループを2重にし、それぞれDECとRAで回しながら座標を変換して1DECごとにlistの最後尾に変換後の座標を格納する。
coord = {'Ecliptic': [], 'Galactic': []} # 座標変換した後の結果入れる用
for dec in DEC:
coord['Ecliptic'].append([])
coord['Galactic'].append([])
for ra in RA:
# 黄道座標系での黄経
ra_e = trans(ra, dec).barycentricmeanecliptic.lon.deg
# 黄道座標系での黄緯
dec_e = trans(ra, dec).barycentricmeanecliptic.lat.deg
coord['Ecliptic'][-1].append(f"({ra_e:.3f}, {dec_e:.3f})")
# 銀河座標系での銀経
ra_g = trans(ra, dec).galactic.l.deg
# 銀河座標系での銀緯
dec_g = trans(ra, dec).galactic.b.deg
coord['Galactic'][-1].append(f"({ra_g:.3f}, {dec_g:.3f})")
print(coord)
# {'Ecliptic': [['(270.000, -66.561)', '(270.000, -66.561)', '(270.000, -66.561)', '(270.000, -66.561)', '(270.000, -66.561)'], ['(347.066, -27.306)', '(14.817, -39.123)', '(90.000, -53.439)', '(168.155, -37.962)', '(347.066, -27.306)'], ['(0.000, -0.000)', '(27.911, -11.472)', '(90.000, -23.439)', '(154.945, -10.404)', '(0.000, -0.000)'], ['(34.566, 52.614)', '(52.963, 44.037)', '(90.000, 36.561)', '(128.868, 44.771)', '(34.566, 52.614)'], ['(90.000, 66.561)', '(90.000, 66.561)', '(90.000, 66.561)', '(90.000, 66.561)', '(90.000, 66.561)']], 'Galactic': [['(302.932, -27.128)', '(302.932, -27.128)', '(302.932, -27.128)', '(302.932, -27.128)', '(302.932, -27.128)'], ['(15.640, -78.354)', '(227.775, -74.690)', '(235.858, -23.549)', '(266.360, 21.325)', '(15.640, -78.354)'], ['(96.337, -60.189)', '(157.005, -58.262)', '(206.989, -11.424)', '(241.572, 43.092)', '(96.337, -60.189)'], ['(116.538, -2.232)', '(131.410, -1.738)', '(153.616, 17.209)', '(151.208, 47.432)', '(116.538, -2.232)'], ['(122.932, 27.128)', '(122.932, 27.128)', '(122.932, 27.128)', '(122.932, 27.128)', '(122.932, 27.128)']]}
これをpd.DataFrameでDataFrameとして扱うことで見やすくする。最大表示列数や最大表示行数を設定することで、列を省略したり途中で改行されたりしないようにしている。
# pandas.DataFrameのFrame自体の最大表示列数（Noneは無制限）
pd.options.display.max_columns = None
# pandas.DataFrameの出力時の最大表示幅（Noneは無制限）
pd.options.display.width = None
df_dct = {}
for X in coord:
df_dct[X] = pd.DataFrame(np.array(coord[X]), columns=RA, index=DEC)
for name, val in df_dct.items():
print(f"{name}\\n{val}\\n")
# Ecliptic
# 0 30 90 153 360
# -90 (269.982, -66.560) (269.982, -66.560) (269.982, -66.560) (269.982, -66.560) (269.982, -66.560)
# -30 (347.059, -27.308) (14.813, -39.126) (90.006, -53.440) (168.154, -37.958) (347.059, -27.308)
# 0 (359.995, -0.000) (27.908, -11.473) (90.002, -23.440) (154.944, -10.403) (359.995, -0.000)
# 60 (34.566, 52.618) (52.965, 44.040) (90.003, 36.561) (128.871, 44.769) (34.566, 52.618)
# 90 (90.010, 66.562) (90.010, 66.562) (90.010, 66.562) (90.010, 66.562) (90.010, 66.562)
# Galactic
# 0 30 90 153 360
# -90 (302.932, -27.128) (302.932, -27.128) (302.932, -27.128) (302.932, -27.128) (302.932, -27.128)
# -30 (15.640, -78.354) (227.775, -74.690) (235.858, -23.549) (266.360, 21.325) (15.640, -78.354)
# 0 (96.337, -60.189) (157.005, -58.262) (206.989, -11.424) (241.572, 43.092) (96.337, -60.189)
# 60 (116.538, -2.232) (131.410, -1.738) (153.616, 17.209) (151.208, 47.432) (116.538, -2.232)
# 90 (122.932, 27.128) (122.932, 27.128) (122.932, 27.128) (122.932, 27.128) (122.932, 27.128)
上にあげたCoordinate Transformation & Galactic Extinction Calculatorで実際に計算していただきたい。ほぼほぼ同じ値になっていると思う。あとは日付などを調整すればいいだろう。
座標変化を一括で行うと時間が生まれる
今回の記事の内容を知る前は、1座標ずつ手打ちで上記サイトで計算してメモしての繰り返しだった。しかし、自動化することで一瞬で作業が終わり他のことに時間を割り振ることができた。
今回の内容はかなりニッチな部類になるが、このように意外にもこんなこともできるんだというようなコードがあるかもしれない。そのコードを使用して作業を短縮することで一気に気持ちが楽になると思われる。
ぜひ自動化していただきたい。
