カテゴリー

astropy

【astropy】astropyを使って天文学の座標変換

2021年4月25日

Galaxy

こんな人にオススメ


天文学で座標系を変換したいんだけど、自分で計算するのは面倒だし、本に書いてある式を実装するのも面倒。何かパッとすぐに変換できる方法ってないの?

ということで、今回は究極にニッチな内容である「astropyを用いて天文座標系を変換する」ということについて解説する。というのも執筆者は天文学系の研究を行なっており、時折ある座標系から他の座標系への変換が必要な時があった。その時にいちいちサイトで座標を打って計算してもらってコピペして、というのが面倒すぎる。
そこで今回のようにpythonを用いて座標変換を行うことで自動的にすぐに値を取り出そうということだ。初めて知った時は衝撃的だった。

執筆者の人生的なものについては以下参照。

【M天パ(めがてんぱ)自己紹介】変わって楽しいの繰り返し

自己紹介 急にキツネの戦い ...

続きを見る

衝撃を受けたastropy関連の記事については以下参照。

【astropy&単位】astropyで変数に単位を付与

こんな人にオススメ 数値計 ...

続きを見る

【astropy&単位】astropyの単位出力書式

こんな人にオススメastropyで数Í ...

続きを見る

python環境は以下。

  • Python 3.9.2
  • astropy 4.2
  • numpy 1.20.1
  • pandas 1.2.3

今回は赤道座標を黄道座標、銀河座標に変換する。

運営者のメガネです。YouTubeTwitterInstagramも運営してます。

自己紹介はこちらから、お問い合わせはこちら。

運営者メガネ

座標として読み込む

基本的にはastropyのレファレンスに書いてある。

が、個人的にすごく読みにくい。したがって、サイトに書いてあるコードをここでも紹介することとする。

1つの座標を読み込む

座標を天文座標として読み込ませるために入力する方式は以下の6種類。

  1. (数値)*u.degreeで角度の単位を持つ値として入力
  2. (数値)と引数でunit='deg'で角度を指定して入力
  3. '00h42m30s'といったようにhmsで入力
  4. '00h42.5m'といったようにhmで入力
  5. '00 42 30 +41 12 00'といったふうにhmsの間をスペースで区切ってかつ、引数unitでそれぞれの角度を指定して入力
  6. '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種類もあるのだが、今回使用するのはbarycentricmeaneclipticgalacticの2つ。前者が黄道座標で後者が銀河座標に相当する。が、eclipticと名のつく項目が8項目もあり、さらにgalacと名のつく項目が4項目もある。以下にその「astropy.coordinates Package」で示されている各項目の説明を引用する。

  • ecliptic
    • barycentricmeanecliptic: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.
  • galac
    • galactic: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]

銀経・銀緯のままでも銀経だけ、銀緯だけとバラしても変換することが可能。

変換した座標系をpandasDataFrame

最後は黄道座標、銀河座標に変換した座標をpandasDataFrameに入れて整理する方法。ループを2重にし、それぞれDECRAで回しながら座標を変換して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.DataFrameDataFrameとして扱うことで見やすくする。最大表示列数や最大表示行数を設定することで、列を省略したり途中で改行されたりしないようにしている。

# 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座標ずつ手打ちで上記サイトで計算してメモしての繰り返しだった。しかし、自動化することで一瞬で作業が終わり他のことに時間を割り振ることができた。

今回の内容はかなりニッチな部類になるが、このように意外にもこんなこともできるんだというようなコードがあるかもしれない。そのコードを使用して作業を短縮することで一気に気持ちが楽になると思われる。

ぜひ自動化していただきたい。

スイッチボット

2022/11/28

【SwitchBotロックレビュー】これからのスタンダードになりうるスマートロック

こんな人にオススメ SwitchBotからスマートロック「SwitchBotロック」が発売された ...

生活に役立つ

2022/11/28

【メガネ厳選】クソ便利に使っているサービスやアイテム達

このページでは執筆者「メガネ」が実際に使って便利だと感じているサ ...

マウス

2022/9/11

【Logicool MX ERGO vs MX Master 3】ERGOをメインにした決定的な理由

こんな疑問・お悩みを持っている人におすすめ 執筆者はLogicoolのハイエンӠ ...

完全ワイヤレスイヤホン(TWS)

2022/11/21

【ながら聴きイヤホン比較】SONY LinkBuds、ambie、BoCoはどれがおすすめ?

こんな人におすすめ 耳を塞がない開放型のイヤホンに完全ワイヤレスӟ ...

macOSアプリケーション

2022/10/15

【M1 Mac】MacBook Proに入れている便利でニッチなアプリを21個紹介する

こんな人におすすめ MacBookを購入してLINEとか必要最低限のアプリは入れた。 ...

完全ワイヤレスイヤホン(TWS)

2022/10/23

【SENNHEISER MOMENTUM True Wireless 3レビュー】高レベルでバランス型の高音質イヤホン

こんな人におすすめ SENNHEISER MOMENTUM True Wireless 3って実際のところどうなの? 評判は良い ...

完全ワイヤレスイヤホン(TWS)

2022/11/21

【SONY WF-1000XM4レビュー】神とゴミのハーフ&ハーフ

こんな人におすすめ SONYのフラグシップモデル「SONY WF-1000XM4」ってどれくらい性 ...

完全ワイヤレスイヤホン(TWS)

2022/8/19

【Nothing ear (1)レビュー】ライトな完成度、アップデートに期待

こんな人にオススメ 完全ワイヤレスイヤホン(TWS)でスケルトンボディ ...

Pythonを学びたいけど独学できる時間なんてない人へのすゝめ

執筆者は大学の研究室・大学院にて独学でPythonを習得した。

でも社会人になったら独学で行うには時間も体力もなくて大変だ。

時間がない社会人だからこそプロの教えを乞うのが効率的。

ここでは色んなタイプに合ったプログラミングスクールの紹介をする。

  • この記事を書いた人

メガネ

ベンチャー企業のWebエンジニア駆け出し。独学のPythonで天文学系の大学院を修了→新卒を1.5年で辞める→転職→今に至る。
常時金欠のガジェット好きでM1 MacBook Pro x Galaxy S22 Ultraの狂人。
人見知りで根暗だったけど、人生楽しもうと思って良い方向に狂う→人生が楽しい

ガジェットのレビューとPythonコードを記事にしています。ぜひ楽しんでください🦊
自己紹介と半生→変わって楽しいの繰り返し

-astropy
-, , ,