こんな人にオススメ
天文学で座標系を変換したいんだけど、自分で計算するのは面倒だし、本に書いてある式を実装するのも面倒。何かパッとすぐに変換できる方法ってないの?
ということで、今回は究極にニッチな内容である「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」で示されている各項目の説明を引用する。
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]
銀経・銀緯のままでも銀経だけ、銀緯だけとバラしても変換することが可能。
変換した座標系を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座標ずつ手打ちで上記サイトで計算してメモしての繰り返しだった。しかし、自動化することで一瞬で作業が終わり他のことに時間を割り振ることができた。
今回の内容はかなりニッチな部類になるが、このように意外にもこんなこともできるんだというようなコードがあるかもしれない。そのコードを使用して作業を短縮することで一気に気持ちが楽になると思われる。
ぜひ自動化していただきたい。