カテゴリー

当サイトはアフィリエイトプログラムによる収益を得ています〈景品表示法に基づく表記です)

astropy

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

2021年4月25日

Galaxy

こんな人にオススメ


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

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

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

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

続きを見る

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

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

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

ガジェット

2023/9/18

【デスクツアー2022下半期】モノは少なく、でも効率的に Desk Updating #0

今回はガジェットブロガーなのにデスク環境を構築していない執筆者の ...

ライフハック

2023/9/16

【Audible vs YouTube Premium】耳で聴く音声学習コンテンツを比較

ワイヤレスイヤホンが普及し耳で学習することへのハードルが格段に下 ...

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

2023/9/18

【SENNHEISER MOMENTUM True Wireless 3レビュー】全てが整ったイヤホン

今回は高音質・高機能なSENNHEISERのフラグシップ完全ワイヤレスイヤホン「SENNH ...

ライフハック

2023/3/11

【YouTube Premiumとは】メリットしかないから全員入れ

今回はYouTube Premiumを実際に使ってみてどうなのか、どんなメリット/デメリット ...

マウス

2023/8/17

【Logicool MX ERGOレビュー】疲れない作業効率重視トラックボールマウス

こんな人におすすめ トラックボールマウスの王道Logicool MX ERGOが気になるけどऩ ...

ベストバイ

2023/9/18

【ベストバイ2022】今年買って良かったモノのトップ10

2022年ベストバイ この1年を振り返って執筆者は何を買ったのか。ガジェッ& ...

スマホ

2023/1/15

【楽天モバイル×povo2.0の併用】月1,000円の保険付きデュアルSIM運用

こんな人におすすめ 楽天モバイルとpovo2.0のデュアルSIM運用って実際のとこ ...

マウス

2023/9/16

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

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

macOSアプリケーション

2022/9/30

【Chrome拡張機能】便利で効率的に作業できるおすすめの拡張機能を18個紹介する

こんな人におすすめ Chromeの拡張機能を入れたいけど、調べても同じような ...

macOSアプリケーション

2023/5/3

【Automator活用術】Macで生産性を上げる作業の自動化術

今回はMacに標準でインストールされているアプリ「Automator」を使ってできる ...

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

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

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

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

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

  • この記事を書いた人

メガネ

Webエンジニア駆け出し。独学のPythonで天文学系の大学院を修了。常時金欠のガジェット好きでM2 Pro MacBook Pro(30万円) x Galaxy S22 Ultra(17万円)使いの狂人。自己紹介と半生→変わって楽しいの繰り返しレビュー依頼など→お問い合わせ運営者情報、TwitterX@m_ten_pa、 YouTube@megatenpa、 Threads@megatenpa

-astropy
-,