2015年2月22日日曜日

astropy.coordinatesを用いた座標の取り扱い、その1 ― パッケージの概要

astropy.coordinatesを用いた天文学上の座標の扱いについてまとめて行きます。ここでは、python(x,y) 2.7.9.0とそれに付随してくるastropyを使用することを前提にしています。

天文学用のPythonパッケージであるastropyには、座標を扱うためのパッケージである coordinatesが含まれています。coordinatesを用いると、様々な座標の相互変換を簡単に行うことができます。まず最初は、coordinatesパッケージを使っておおよそどのようなことができるのか概要を見ていくことにします。

coordinatesを使うときに基本となるのが、SkyCoordというオブジェクトです。SkyCoordオブジェクトは、様々な座標表記に柔軟に対応しています。では、具体例を見て行きましょう。

まず、必要なパッケージを読み込みます。以下の例では、単位を扱うパッケージ unitsも使用するので、SkyCoordクラスと合わせて importしておきます。

>>> from astropy import units as u
>>> from astropy.coordinates import SkyCoord

必要なパッケージを読み込んだら、SkyCoordオブジェクトを定義してみます。SkyCoordオブジェクトは同じ内容を様々な形式で定義することができます。以下の定義は全て同じ内容です。

>>> c = SkyCoord(ra=10.625*u.degree, dec=41.2*u.degree, frame='icrs')
>>> c = SkyCoord(10.625, 41.2, frame='icrs', unit='deg')
>>> c = SkyCoord('00h42m30s', '+41d12m00s', frame='icrs')
>>> c = SkyCoord('00h42.5m', '+41d12m', frame='icrs')
>>> c = SkyCoord('00 42 30 +41 12 00', frame='icrs', unit=(u.hourangle, u.deg))
>>> c = SkyCoord('00:42.5 +41:12', frame='icrs', unit=(u.hourangle, u.deg))

各定義を行った後で、実際に c の内容を表示し、全て同じ内容であることを確認しましょう。

>>> c
<SkyCoord (ICRS): (ra, dec) in deg (10.625, 41.2)>

上記のように、座標値は特に変数名を指定しなくても記述することができます。変数名を明示して値を記述するときは、赤経、赤緯は、ra、dec、銀経銀緯は、l、b で指定します。座標系の種類にデフォルト値は存在せず必ず値を指定しなくてはいけません(以下で見るように、異なる座標間の変換や比較を行わない場合は省略しても構いません)。座標系の種類も、変数を使わず指定できますが、変数を明示して指定する場合は、frame という変数名を用います。

数値を与えるときは単位を指定する必要があります。例えば、10.5*u.degree、+41d12m00s の様な感じです。単位の指定は astropy.unitsパッケージを用いて行います。
(astropy.unitsに関しては、ここ参照。)

赤道座標系(赤経、赤緯)を表すSkyCoordオブジェクトは、LongitudeオブジェクトとLatitudeオブジェクトから構成されており、個別に情報を取り出すことができます。LongitudeオブジェクトとLatitudeオブジェクは、角度を扱う Angleクラスの一種なので、下記の用例のように角度表示の方法を様々に変更することが可能です。

>>> c = SkyCoord(ra=10.68458*u.degree, dec=41.26917*u.degree)
>>> c.ra  
<Longitude 10.68458 deg>
>>> c.ra.hour  
0.7123053333333335
>>> c.ra.hms  
hms_tuple(h=0.0, m=42.0, s=44.299200000000525)
>>> c.dec  
<Latitude 41.26917 deg>
>>> c.dec.degree  
41.26917
>>> c.dec.radian  
0.7202828960652683

SkyCoordオブジェクトの値は、文字列として取り出すこともできます。文字列への変換は、to_sting() メソッドを用いて行います。

>>> c = SkyCoord(ra=10.68458*u.degree, dec=41.26917*u.degree)
>>> c.to_string('decimal')
'10.6846 41.2692'
>>> c.to_string('dms')
'10d41m04.488s 41d16m09.012s'
>>> c.to_string('hmsdms')
'00h42m44.2992s +41d16m09.012s'

to_string() メソッドでは、上記のフォーマット以外にも、LaTeX形式など様々なフォーマットを利用することができます。

>>> c.ra.to_string(decimal=True)
'10.6846'
>>> c.dec.to_string(format='latex')
'$41^\\circ16{}^\\prime09.012{}^{\\prime\\prime}$'
>>> msg = 'My coordinates are: ra="{0}"" dec="{1}"'
>>> msg.format(c.ra.to_string(sep=':'), c.dec.to_string(sep=':'))
'My coordinates are: ra="10:41:04.488"" dec="41:16:09.012"'

ここまでに示した例では、座標の変換や比較は行う必要がなかったので frameの指定は行いませんでしたが、基本的には SkyCoordオブジェクトを定義するときには frameの指定を行う必要があります。

>>> c_icrs = SkyCoord(ra=10.68458*u.degree, dec=41.26917*u.degree, frame='icrs')

frameの指定を含めてSkyCoordオブジェクトを定義すると、座標値を他の座標における値に変換することができるようになります。座標値の変換にはいくつかの方法がありますが、一番応用が利く方法は、transform_to メソッドを用いる方法です。

>>> from astropy.coordinates import FK5
>>> c_icrs.galactic  
<SkyCoord (Galactic): (l, b) in deg
    (121.174241811, -21.5728855724)>
>>> c_fk5 = c_icrs.transform_to('fk5')  # c_icrs.fk5 does the same thing
>>> c_fk5  
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
    (10.6845915393, 41.2691714591)>
>>> c_fk5.transform_to(FK5(equinox='J1975'))  # precess to a different equinox  
<SkyCoord (FK5: equinox=J1975.000): (ra, dec) in deg
    (10.3420913461, 41.1323211229)>

transform_to メソッドを用いると、天体座標を AltAz座標系に変換することもできるので、例えば観測計画をたてるとき等にも SkyCoord オブジェクトを活用することができます。実際に SkyCoordオブジェクトを用いた観測計画の検討方法については機会をあらためて解説したいと思います。

SkyCoordは、Pythonの配列に対応しています。SkyCoordを用いて天体カタログなどの大量のデータを処理するときには、ループを用いて処理するよりも配列を用いて処理する方が遥かに高速です。

>>> SkyCoord(ra=[10, 11]*u.degree, dec=[41, -5]*u.degree)
<SkyCoord (ICRS): (ra, dec) in deg
    [(10.0, 41.0), (11.0, -5.0)]>

ここまでの例では球面座標を使用してきました。SkyCoordオブジェクトではデフォルトが球面座標なので特に明示的に定義は行いませんでした。SkyCoordオブジェクトでは、球面座標以外にも直交座標と円柱座標を使用することができます。球面座標以外の座標を用いる場合は、representation という変数に使用する座標の種類を明示的に定義する必要があります。

>>> c = SkyCoord(x=1, y=2, z=3, unit='kpc', frame='icrs', representation='cartesian')
>>> c
<SkyCoord (ICRS): (x, y, z) in kpc
    (1.0, 2.0, 3.0)>
>>> c.x, c.y, c.z
(<Quantity 1.0 kpc>, <Quantity 2.0 kpc>, <Quantity 3.0 kpc>)

>>> c.representation = 'cylindrical'
>>> c  
<SkyCoord (ICRS): (rho, phi, z) in (kpc, deg, kpc)
    (2.2360679775, 63.4349488229, 3.0)>
>>> c.phi
<Angle 63.434948... deg>
>>> c.phi.to(u.radian)
<Angle 1.107148... rad>

>>> c.representation = 'spherical'
>>> c  
<SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, kpc)
    (63.4349488229, 53.3007747995, 3.74165738677)>

>>> c.representation = 'unitspherical'
>>> c  
<SkyCoord (ICRS): (ra, dec) in deg
    (63.4349488229, 53.3007747995)>

SkyCoordオブジェクトを用いると、2つの天体の間の角距離の計算や、天体カタログどうしのクロスチェックなどを簡単に行うことができます。天体カタログどうしのクロスチェックについては別の機会に解説したいと思います。

>>> c1 = SkyCoord(ra=10*u.degree, dec=9*u.degree, frame='icrs')
>>> c2 = SkyCoord(ra=11*u.degree, dec=10*u.degree, frame='fk5')
>>> c1.separation(c2)  # c1とc2の間の各距離を求めています。  
<Angle 1.4045335865905868 deg>

Skycoordオブジェクトでは、天体までの距離も指定することができます。距離原点は座標系にもよりますが、基本的に地球からの距離を与えます。距離を指定すると空間における位置が決まるので、天体の位置をデカルト座標で表示することが可能となります。

>>> from astropy.coordinates import Distance
>>> c = SkyCoord(ra=10.68458*u.degree, dec=41.26917*u.degree, distance=770*u.kpc)
>>> c.cartesian.x  
<Quantity 568.7128654235232 kpc>
>>> c.cartesian.y  
<Quantity 107.3008974042025 kpc>
>>> c.cartesian.z  
<Quantity 507.88994291875713 kpc>

天体までの距離を与えると、2天体間の空間的な距離を計算するなど、SkyCoordオブジェクトを用いて、さらに様々な計算ができるようになります。

>>> c1 = SkyCoord(ra=10*u.degree, dec=9*u.degree, distance=10*u.pc, frame='icrs')
>>> c2 = SkyCoord(ra=11*u.degree, dec=10*u.degree, distance=11.5*u.pc, frame='icrs')
>>> c1.separation_3d(c2)  
<Distance 1.5228602415117989 pc>

インターネットへ接続できる環境があれば、from_nameというメソッドを使うことによって座標情報をネットから(おそらくSIMBAD)から取得することができます。

>>> SkyCoord.from_name("M42")  
<SkyCoord (ICRS): (ra, dec) in deg
    (83.82208, -5.39111)>

from_nameは便利ですが正確な位置情報の提供を必ずしも保証するものではないので信用できる位置情報を用いいる必要がある場合は主導で位置情報を与えるほうが無難です。