astropy.coordinatesを用いた座標の取り扱い、その３ ― SkyCoodクラスの使用法

その１でも若干説明したようにSkyCoodクラスは、柔軟で使いやすいユーザーインターフェースを提供しており、天文座標の扱いや座標間の変換などを比較的簡単に取り扱うことができます。SkyCoodは、ICRSクラスやFK5クラス等の、特定の座標系を扱うクラスに比べてより多くの機能を備えており、ユーザーにやさしい設計となっていることから、「高水準クラス」と呼ばれています。一方で、ICRSクラスやFK5クラスなどは、どちらかというと機械よりの設計がなされており、使いづらい点はあるものの計算量の大きい作業などを効率よくこなすことができるため、「低水準クラス」と呼ばれています。

高水準クラスである SkyCoordと、ICRSクラスやFK5クラス等の低水準クラスの主な違いは次のとおりです。

• SkyCoodオブジェクトには、観測時間や観測場所といった細かな観測の属性を保存することができ、このような属性は例え座標系の変換を行っても保持されます。一方で低水準クラスを用いて、異なるクラス間（例えば、ICRSクラスからFK5クラス）で変換を行うと、その作業に必要な情報は変換先のクラスに受け継がれますがそれ以外の情報は失われてしまいます。したがって、逆方向に再度変換しても完全にもとの状態には戻りません。SkyCoodクラスの場合は再び座標変換をすると完全に元の状態に戻ります。
• SkyCoordクラスは、低水準クラスに比べると、座標値のフォーマットが柔軟で、さまざまな書式で座標値を定義することができます。
• SkyCoordには、沢山のメソッドが付随しており、典型的な解析を簡単に取り扱うことができます。
• SkyCoordオブジェクトは、現時点では（ver. 1.0.1）、astropy.coordinates.frame_transform_graph で定義された座標系のみを取り扱うことができます。

SkyCoordオブジェクトの作成

SkyCoordクラスでは、さまざまなフォーマットでSkyCoordオブジェクトを定義することができます。最低限必要なのは座標値と座標系ですが、場合によっては座標系の指定すら省略することも可能です。

SkyCoodオブジェクトの基本的なフォーマットは以下のとおりです。大文字で、COORDやFRAMEと示しているパラメータについては、以下の＊＊＊の項目で説明しています。[unit=UNIT] のように[...]内のパラメータはオプションです。

SkyCoord(COORD, [FRAME], keyword_args ...)
SkyCoord(LON, LAT, [frame=FRAME], [unit=UNIT], keyword_args ...)
SkyCoord([FRAME], =LON, =LAT, keyword_args ...)

SkyCoordオブジェクトの標準的な定義の例を以下に示します。ここに示す例は、全て球面座標の場合です。（球面座標は、全ての座標系においてデフォルトの座標となっています。）　SkyCoordでは、必要に応じて直交座標系や円柱座標系なども用いることができます。他の座標系を用いる場合は、以下の＊＊＊のセクションで説明します。SkyCoorオブジェクトを作成するには、まず必要なモジュールをインポートします。

>>> from astropy.coordinates import SkyCoord  # SkyCoodをインポートします。
>>> from astropy.coordinates import ICRS, Galactic, FK4, FK5  # その他低水準クラスもインポートします。
>>> from astropy.coordinates import Angle, Latitude, Longitude  # Angleクラスをインポートします。
>>> import astropy.units as u # astropy.unitをインポートします。
>>> import numpy as np # numpyをインポートします。

まず、赤道座標系（RA、Dec）でのSkyCoodオブジェクトの定義から見て行きましょう。

>>> SkyCoord(10, 20, unit='deg')  # デフォルトの座標系はICRSです。
<SkyCoord (ICRS): (ra, dec) in deg
    (10.0, 20.0)>

>>> SkyCoord([1, 2, 3], [-30, 45, 8], frame='icrs', unit='deg')
<SkyCoord (ICRS): (ra, dec) in deg
    [(1.0, -30.0), (2.0, 45.0), (3.0, 8.0)]>

最初の例では、座標系は明示的には定義していません。SkyCoodオブジェクトでは、デフォルトの座標系がICRSなので、座標系を明示的に定義しない場合は自動的にICRSが採用されます。しかし、基本的には使用する座標系を明示的に定義するほうが間違いが少ないでしょう。ICRS（国際天文基準座標系）は、元期J2000.0の赤道座標系と近似的にほぼ同じ座標値を与えます。

論文やカタログなどでよく使われる標準的なフォーマットにであれば、座標値を文字列としてそのままSkyCoordオブジェクトに読み込むことができます。座標系は、FK4と大文字で書くこともできますし、小文字で与えることもできます（小文字の場合は'fk4'とします）。以下具体例を示します。

>>> coords = ["1:12:43.2 +1:12:43", "1 12 43.2 +1 12 43"]
>>> sc = SkyCoord(coords, frame=FK4, unit=(u.hourangle, u.deg), obstime="J1992.21")
>>> sc = SkyCoord(coords, frame=FK4(obstime="J1992.21"), unit=(u.hourangle, u.deg))
>>> sc = SkyCoord(coords, frame='fk4', unit='hourangle,deg', obstime="J1992.21")

>>> sc = SkyCoord("1h12m43.2s", "+1d12m43s", frame=Galactic)  # 単位は文字列から認識
>>> sc = SkyCoord("1h12m43.2s +1d12m43s", frame=Galactic)  # 単位は文字列から認識
>>> sc = SkyCoord(l="1h12m43.2s", b="+1d12m43s", frame='galactic')
>>> sc = SkyCoord("1h12.72m +1d12.71m", frame='galactic')

他にも、様々なフォーマットで ra と dec の値を書くことができます。

>>> sc = SkyCoord("15h17+89d15")
>>> sc = SkyCoord("275d11m15.6954s+17d59m59.876s")
>>> sc = SkyCoord("8 00 -5 00.6", unit=(u.hour, u.deg))
>>> sc = SkyCoord("J080000.00-050036.00", unit=(u.hour, u.deg))
>>> sc = SkyCoord("J1874221.31+122328.03", unit=u.deg)

前回登場したLongitudeオブジェクト、Latitudeオブジェクト、Angleオブジェクトを使ってSkyCoordオブジェクトを定義することもできます。astropyのunitを使って、普通のNumpy arrayなどからSkyCoordオブジェクトを定義することもできます。

>>> ra = Longitude([1, 2, 3], unit=u.deg)  # Angleオブジェクトも同様に使用可能
>>> dec = np.array([4.5, 5.2, 6.3]) * u.deg  # AstropyのUnitを使用してNumpy配列で定義
>>> sc = SkyCoord(ra, dec, frame='icrs')
>>> sc = SkyCoord(ra=ra, dec=dec, frame=ICRS, obstime='2001-01-02T12:34:56')

低水準クラスのオブジェクトを使ってSkyCoorオブジェクトを定義することもできます。

>>> c = FK4(1 * u.deg, 2 * u.deg)
>>> sc = SkyCoord(c, obstime='J2010.11', equinox='B1965')  # FK4のデフォルト設定を上書き

FK4などの低水準クラスは、それぞれに各種のデフォルト値を持っています。例えば、FK4の場合、equinox='B1950.0 や obstime=equinox がデフォルト値です。これらのデフォルト値は、低水準クラスをSkyCoordオブジェクトとして取り込む際に、別の値に書き換えても適切に処理されます。ただし、低水準クラスで明示的に値を与えている場合、例えば上の例で、c = FK4(1 * u.deg, 2 * u.deg, equinox='B1960') とした場合、異なる equinox でSkyCoordオブジェクトを定義しようとするとエラーの原因となります。

SkyCoordオブジェクトの書式

SkyCoordオブジェクトが扱うことのできるほとんどの座標系ではデフォルのと座標は球面座標です。球面座標を用いる場合、SkyCoordオブジェクトの書式は以下の形式になります。

SkyCoord(COORD, [FRAME | frame=FRAME], [unit=UNIT], keyword_args ...)
SkyCoord(LON, LAT, [DISTANCE], [FRAME | frame=FRAME], [unit=UNIT], keyword_args ...)
SkyCoord([FRAME | frame=FRAME], =LON, =LAT, [unit=UNIT],
         keyword_args ...)

上に示した書式の中で、大文字で書かれた要素（例えば、FRAME）はユーザーが与える情報です。括弧の中の要素はオプションなので必要に応じて与えます。（球面座表以外の場合の書式は＊＊＊で説明します。）

LON、LAT

LON、LATでは球面座標における経度、緯度を指定します。経度、緯度を与えるとき、次のようなオプションを選択することができます。

• 単独の値を与える場合は以下の方法がある。（１）Quantityオブジェクトを用いる、（２）数値を与え、astropy.unitで単位を指定する、（３）その２で説明したAngle文字列で指定する。
• 複数の値を与える場合は、Pythonリスト、Quantityオブジェクトの配列、Numpy配列などが使用できる。
• その２で説明した、Angleオブジェクト、Latitudeオブジェクト、Longitudeオブジェクトを使って指定することも可能。これらのオブジェクトは、スカラー（単独値）としても配列としても使用可能。

DISTANCE

オプションとして、天体までの距離をSkyCoordオブジェクトの中に書き込むことができます。距離は、次のような大きく２つの方法で与えることができます。

• 単独の値として与える場合は、（１）Quantityオブジェクト、またはDistanceオブジェクトとして指定する、（２）無次元量としての距離を数値で与える、（３）数値で距離を与え、astropy.unitで単位を指定する。
• Pythonリスト、Quantityオブジェクトの配列、Distanceオブジェクトの配列、Numpy配列などで与える。

COORD

SkyCoordオブジェクトの座標値の情報は、COORDという単一の情報として与えることもできます。球面座標系の場合は、経度緯度の値のペアを単位として、１つのペア、もしくは複数のペアをCOORDとして与えることができます。COORDの情報を与える場合に使える書式は以下のとおりです。

• 経度と緯度の値をスペースで区切って並べた文字列。経度と緯度の書式は、その２の「角度の与え方」で説明した方法で与える。
• 座標値の文字列からなるPythonリストもしくはNumpy配列。
• (経度、緯度)の形のPythonタプルからなるリスト。ただし、経度、緯度はスカラーとして与える（つまり、配列では与えない）。
• N x 2 型のNumpy配列もしくはQuantity配列。第１コラムに経度、第２コラムに緯度を与える。例えば、 [[270, -30], [355, +85]] * u.deg のような形。
• (経度、緯度、距離)の形のタプルからなるPythonリスト。
• N x 3 型のNumpy配列もしくはQuantity配列。第１コラムに経度、第２コラムに緯度、第３コラムに距離を与える。

上記のような書式に加えて、以下の様な、より柔軟な表記にも対応しています。

• FK4などの低水準クラスで定義したオブジェクト。例えば、FK4(1*u.deg, 2*u.deg, obstime='J2012.2')
• SkyCoordオブジェクト（つまり、他のSkyCoordオブジェクトの座標情報をコピーすることになります）。
• SphericalRepresentation, CylindricalRepresentation, CartesianRepresentationなどの、BaseRepresentationサブクラスオブジェクト。

FRAME

ここでは座標系を指定します。 astropyには次のような座標系が用意されています、ICRS、FK5、 FK4、FK4NoETerms、Galactic、AltAz。パラメータとしてこれらの座標系を指定するときは小文字で書きます。

FRAMEを指定しない場合は、一応ICRSとして扱われます。ただし、具体的にFRAMEを指定しないと、異なった座標間における比較など、一部の機能が使用出来ません。

unit=UNIT

次のように単位を記法を指定することができます。

• 角度を表すUnitオブジェクト
• 角度の単位を表す文字列（スペースは含まない）。
• UnitオブジェクトからなるPythonタプル、もしくは単位を表す文字列からなるPythonタプル。単位の並べ方は、LON、LATの順。例えば、 ('hourangle', 'degree')
• 2つの単位を表す文字列をコンマで区切った単一の文字列。例えば、'hourangle,degree'

もし単位を一つだけ与えた場合は、LON、LATに同じ単位が指定されます。

その他の変数

経度、緯度を表すとき、一部の座標系ではその経度、緯度に対してその座標系特有の呼称を用いることができます。

例えば、ICRS、FK5、FK4、FK4NoETermsといった座標系では、ra、decで経度、緯度を表すことができます。同様に、銀河系座標（Galactic）では、l、ｂで経度、緯度を表すことができます。

次のキーワードはどのような座標系でも用いることができます。

distance：　座標系の原点からの距離を与えます。（Distanceクラスの設定値を初期化します。）
obstime：　観測時間を与えます。（Timeクラスの設定値を初期化します。）
equinox：　座標家の分点を与えます。（Timeクラスの設定値を初期化します。）

配列を使った操作

SkyCoordオブジェクトでは、座標値を配列として与えることができます。配列として座標値を与えておくと、単独の座標値を与えたSkyCoordオブジェクトを繰り返し処理するのに比べて、格段に高速に各種の作業を行うことができます。例えば、１０００個の座標値をICRSから銀河系座標（galactic）に変換する場合でその処理速度の違いを見てみましょう。

>>> ra = np.random.uniform(0, 360, size=1000) * u.deg
>>> dec = np.random.uniform(-90, 90, size=1000) * u.deg

>>> sc_list = [SkyCoord(r, d, frame='icrs') for r, d in zip(ra, dec)]
>>> timeit sc_gal_list = [c.galactic for c in sc_list]  
1 loops, best of 3: 7.66 s per loop

>>> sc = SkyCoord(ra, dec, frame='icrs')
>>> timeit sc_gal = sc.galactic  
100 loops, best of 3: 8.92 ms per loop

また、当然ながらPythonの配列と同様に選択やスライスなどの処理も行えます。

>>> north_mask = sc.dec > 0
>>> sc_north = sc[north_mask]
>>> len(sc_north)  
504
>>> sc[2:4]  
<SkyCoord (ICRS): (ra, dec) in deg
    [(304.304015..., 6.900282...),
     (322.560148..., 34.872244...)]>
>>> sc[2]  
<SkyCoord (ICRS): ra=304.304015... deg, dec=6.900282... deg>

SkyCoordの便利な機能

SkyCoordクラスには他にも多くの便利な機能（attribute、method）が備わっています。SkyCoordに備わった、attributeやmethodを使いこなすことで、より効率的に作業を進めることができます。

Attributeやmethodを使う際に知っておきたいのがtab補完（TAB-discovery）です。IPythonを使用している場合は、SkyCoordオブジェクトの名前を入力した後にピリオドをタイプしTABキーを押すと、どのようなattributeやmethodが使用可能が表示されます。

>>> sc = SkyCoord(1, 2, frame='icrs', unit='deg', obstime='2013-01-02 14:25:36')
>>> sc.  
sc.cartesian                           sc.match_to_catalog_3d
sc.data                                sc.match_to_catalog_sky
sc.dec                                 sc.name
sc.default_representation              sc.obstime
sc.distance                            sc.position_angle
sc.equinox                             sc.ra
sc.fk4                                 sc.realize_frame
sc.fk4noeterms                         sc.represent_as
sc.fk5                                 sc.representation
sc.frame                               sc.representation_component_names
sc.frame_attr_names                    sc.representation_component_units
sc.frame_specific_representation_info  sc.representation_info
sc.from_name                           sc.separation
sc.galactic                            sc.separation_3d
sc.get_frame_attr_names                sc.shape
sc.has_data                            sc.spherical
sc.icrs                                sc.time_attr_names
sc.is_frame_attr_default               sc.to_string
sc.is_transformable_to                 sc.transform_to
sc.isscalar

非常に多くの機能が表示されますが、attributeには比較的機能に準拠した名前が付けられているので、表示された文字列から機能は想像がつく場合が多いでしょう。簡単な例では、赤経、赤緯を表示するには、ra、dec を用います。

>>> sc.ra
<Longitude 1.0 deg>
>>> sc.dec
<Latitude 2.0 deg>

他にもいろいろと使えそうなattribute名が目につくはずです。例えば、座標系の名前、icrs, galactic, fk5 fk4, fk4noeterms が目につくと思います。これらのattributeを使用すると、元の座標値を該当する座標系における座標値に変換した後、その座標系における新たなSkyCoordオブジェクトが作成されます。

>>> sc_gal = sc.galactic
>>> sc_gal  
<SkyCoord (Galactic): (l, b) in deg
    (99.6378552814, -58.7096929334)>

この他にも、distance, equinox, obstime, shape、などは覚えておくと便利でしょう。

astropy.coordinatesを用いた座標の取り扱い、その２ ― 角度の扱い

astropy.coordinatesの中で、座標系を表すオブジェクトは、Angle というクラスによって構成されています。Angle クラスはastropy.coordinatesの中で特に重要ですが、astropyの他のパッケージでも角度を扱うときに必ず出てくるクラスなのでその機能を知っておくことは重要です。

Angleオブジェクトの作成

Angleオブジェクトは様々な形で定義することができます。インプットとして与える角度は、単独の値はもちろんのこと、Python配列、タプル、文字列など様々な形をとることができます。Angleオブジェクトの定義の例を以下に示します。

>>> import numpy as np
>>> from astropy import units as u
>>> from astropy.coordinates import Angle

>>> Angle('10.2345d')              # d（°）を単位として文字列で定義
<Angle 10.2345 deg>
>>> Angle(['10.2345d', '-20d'])    # 文字列の配列で定義
<Angle [ 10.2345,-20.    ] deg>
>>> Angle('1:2:30.43 degrees')     # d（°）を60進数で定義
<Angle 1.0417861111111113 deg>
>>> Angle('1 2 0 hours')           # hoursを60進数で定義 
<Angle 1.0333333333333334 hourangle>
>>> Angle(np.arange(1., 8.), unit=u.deg)  # Numpy array from 1..7 in degrees
<Angle [ 1., 2., 3., 4., 5., 6., 7.] deg>
>>> Angle(u'1°2′3″')               # Unicodeの記号をそのまま使うこともできる
<Angle 1.0341666666666667 deg>
>>> Angle('1d2m3.4s')              # Degree, arcmin, arcsec.  
<Angle 1.0342777777777779 deg>
>>> Angle('-1h2m3s')               # Hour, minute, second  
<Angle -1.0341666666666667 hourangle>
>>> Angle((-1, 2, 3), unit=u.deg)  # (degree, arcmin, arcsec)  
<Angle -1.0341666666666667 deg>
>>> Angle(10.2345 * u.deg)         # From a Quantity object in degrees  
<Angle 10.2345 deg>
>>> Angle(Angle(10.2345 * u.deg))  # AngleオブジェクトでAngleオブジェクトを定義
<Angle 10.2345 deg>

角度の与え方

Angleオブジェクトは、様々な方法で角度の値を与えることができます。角度の値は文字列でも浮動小数点数のどちらでも構いません。

>>> a = Angle(1, u.radian)
>>> a
<Angle 1.0 rad>
>>> a.radian
1.0
>>> a.degree  
57.29577951308232
>>> a.hour  
3.8197186342054885
>>> a.hms  
hms_tuple(h=3.0, m=49.0, s=10.987083139758766)
>>> a.dms  
dms_tuple(d=57.0, m=17.0, s=44.806247096362313)
>>> a.signed_dms  
signed_dms_tuple(sign=1.0, d=57.0, m=17.0, s=44.806247096362313)
>>> (-a).dms  
dms_tuple(d=-57.0, m=-17.0, s=-44.806247096362313)
>>> (-a).signed_dms  
signed_dms_tuple(sign=-1.0, d=57.0, m=17.0, s=44.806247096362313)
>>> a.arcminute  
3437.7467707849396
>>> a.to_string()
u'1rad'
>>> a.to_string(unit=u.degree)
u'57d17m44.8062s'
>>> a.to_string(unit=u.degree, sep=':')
u'57:17:44.8062'
>>> a.to_string(unit=u.degree, sep=('deg', 'm', 's'))
u'57deg17m44.8062s'
>>> a.to_string(unit=u.hour)
u'3h49m10.9871s'
>>> a.to_string(unit=u.hour, decimal=True)
u'3.81972'

Angleオブジェクトの使用例

Angleオブジェクトは、直接Python内の各種演算に対応しています。

>>> a = Angle(1.0, u.radian)
>>> a + 0.5 * u.radian + 2 * a
<Angle 3.5 rad>
>>> np.sin(a / 2)  
<Quantity 0.479425538604203>
>>> a == a
array(True, dtype=bool)
>>> a == (a + a)
array(False, dtype=bool)

もちろんAngleオブジェクトは天体の座標を指定する際にも使われます。

>>> from astropy.coordinates import ICRS
>>> ICRS(Angle(1, u.radian), Angle(0.5, u.radian)) 
<ICRS Coordinate: (ra, dec) in deg (57.2957795131, 28.6478897565)>

角度値のラッピングと範囲内かどうかの判定

角度を扱う上で便利な機能を２つ紹介しておきます。１つは角度をラッピングするための wrap_at() 、もう１つは角度が特定の範囲内に収まっているかどうかを判断する is_within_bounds() です。

wrap_at() は次のような形式で使用します。wrap_at(wrap_angle, inplace=False) 

wrap_angleで、ラッピングする範囲を指定します。例えば、角度を０°から３６０°の範囲で表示したい場合、wrap_angleは３６０°、－１８０°から１８０°の範囲で表したい場合は wrap_angleは１８０°となります。（つまり数式で表すと、wrap_angle - 360d <= angle < wrap_angle ということになります。）　inplaceでは、結果を別のAngleオブジェクトとして返すか入力として与えた角度を上書きするかを指定します。デフォルトでは inplace=Falseとなっており、結果を新たなAngleオブジェクトとしてして返します。入力した角度を上書きしたい場合は、inplace=True とします。以下、具体的な使用例を見てみましょう。

>>> from astropy.coordinates import Angle
>>> import astropy.units as u
>>> a = Angle([-20.0, 150.0, 350.0] * u.deg)

>>> a.wrap_at(360 * u.deg).degree  # 0°から360°の範囲でラッピング。inplaceはデフォルト値。
array([ 340.,  150.,  350.])

>>> a.wrap_at('180d', inplace=True)  # -180°から180°でラッピング。inplace=Trueなのでaが上書きされる。
>>> a.degree  # aをdegree単位で表示。
array([ -20.,  150.,  -10.])

is_within_bounds()は、インプットとして当てる角度が指定した範囲内に入っているかどうかを判定します。is_within_dounds()は、次の形式で使用します。is_within_bounds(lower=None, upper=None)

与えたlowerとupperの値に対して、与えたangleが、lower <= angle < upper を満たしているかどうかを判定します。もし、lower の値を指定しない場合（もしくは、None を指定した場合）は、upperのみが考慮されます。これは upper に対しても同じです。以下、使用例を見てみましょう。

>>> from astropy.coordinates import Angle
>>> import astropy.units as u
>>> a = Angle([-20, 150, 350] * u.deg)
>>> a.is_within_bounds('0d', '360d')
False
>>> a.is_within_bounds(None, '360d')
True
>>> a.is_within_bounds(-30 * u.deg, None)
True

LongitudeオブジェクトとLatitudeオブジェクト

LongitudeオブジェクトとLatitudeオブジェクトはAngleオブジェクトと類似したオブジェクトですが、球面座標系における座標値を表すことに特化したオブジェクトです。Longitudeオブジェクトは、例えば、赤道座標系における赤経、銀河座標系における銀系、Alt-Az座標系における方位角を表す場合に用います。一方、Latitudeオブジェクトは、赤道座標系における赤緯、銀河座標系における銀緯、Alt-Az座標系における仰角を表す場合に用いられます。

Longitudeオブジェクト

Longitudeオブジェクトは、wrap_angle というパラメータを設定することができる点が普通のAngleオブジェクトと異なるところです。wrap_angle を指定することによって、あらゆる角度の値を次のような一定の範囲内の角度値で表現します。

wrap_angle - 360 * u.deg <= angle(s) < wrap_angle

デフォルトの wrap_angle は360°です。wrap_angle=180 * u.deg と設定すると、値の範囲が-180°から+180°に変更されます。Longitudeオブジェクトは、wrap_angleを変更すると新しい値はインプットのAngleオブジェクトに上書きされます。以下で使用例を見てみましょう。

>>> from astropy.coordinates import Longitude
>>> a = Longitude([-20, 150, 350, 360] * u.deg)
>>> a.degree
array([ 340., 150., 350.,   0.])
>>> a.wrap_angle = 180 * u.deg
>>> a.degree
array([ -20., 150., -10.,   0.])

Latitudeオブジェクト

Latitudeオブジェクトは、値が以下のような範囲に限定される点で通常のAngleオブジェクトと異なります。

-90.0 * u.deg <= angle(s) <= +90.0 * u.deg

この範囲外で値を定義しようとするとエラーが返ってきます。