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

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

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

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

>>> 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>

天体までの距離を与えると、２天体間の空間的な距離を計算するなど、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は便利ですが正確な位置情報の提供を必ずしも保証するものではないので信用できる位置情報を用いいる必要がある場合は主導で位置情報を与えるほうが無難です。

【第５回】 Pythonの基本的な型： コンテナ、その３ （文字列の基本）

さて前回までリストについて見てきましたが（第３回、第４回）今回は同じくコンテナ型の代表格である「文字列」について見てゆきましょう。

Pythonで文字列を表す方法は以下のように四通りあります。

>>> s = 'Hello, how are you?'
>>> s = "Hi, what’s up"
>>> s = '''Hello, 
>>>           how are you''' 
>>> s = """Hi,
>>>        what’s up?"""

基本的には、「’」か「”」で文字列を囲うだけです。一重の場合は改行ができませんが「’’’」や「”””」のように三重に囲うと、文字列の途中で改行できます。

「’」と「”」は、文字列中に出てくる文字によって使い分けます。例えば、以下の様に、文中に「’」が出てくる場合は、「’」で文字列指定を行うとエラーの原因となります。

>>> 'Hello, what's up you?'
  File "", line 1
    'Hello, what's up you?'
                 ^
SyntaxError: invalid syntax

この場合、「”」を使うことでエラーを解消することができます。

>>> "Hello, what's up you?"
"Hello, what's up you?"

「’’’」や「”””」で囲むと、入力する文字列の途中に改行を入れることができます。改行は「/n」で表されます。（ちなみに、tab は「\t」で表されます。）

>>> '''Hello,
>>> How are you'''
'Hello,\nHow are you'

文字列はリストと同様にコンテナ型です。コンテナ型は第３回で説明したようにオブジェクトの集まりで、オブジェクトを操作することができます。文字列においてリストの各要素に相当するものが各文字ということになります。したがって、リストと全く同様、インデクス指定やスライス指定を行うことができます。以下例を見てゆきましょう。

インデクス指定：

>>> a = "hello"
>>> a[0]
'h'
>>> a[1]
'e'
>>> a[-1]
'o'

スライス指定：

>>> a = "hello, world!"
>>> a[3:6]
'lo,'
>>> a[2:10:2]
'lo o'
>>> a[::3]
'hl r!'

a[3:6]は、左から４番目から５番目の要素（文字列の場合は文字）を指定しています。a[2:10:2]は、２番めの文字から始まって、一つ置きに９文字目まで文字を取り出します。a[ : :3]は、最初から最後まで２つ置きに文字を取り出します。スライス構文の定義については第３回を参照してください。

複雑な記号などを扱いたいときは、クォーテーションの前に「u」をつけることで、ユニコードで文字を指定することもできます。

>>> u'Hello\u0020World !'
u'Hello World !'

\u0020 は、ユニコード値 0x0020 に相当するユニコードキャラクターを表します。0x0020 はスペースを表すので、Pythonはスペースに表示を置き換えてくれます。（ユニコードを使った文字の扱いに関するさらに詳しい情報はここ。）

次回は文字列の操作についてまとめます。

【第４回】 Pythonの基本的な型： コンテナ、その２ （リストの操作）

前回（第３回）では、コンテナのその１として、リスト、インデクス、スライスについてまとめました。今回は、コンテナのその２として、リストの要素のさらに詳しい操作方法を見てゆきましょう。

Pythonのリストの大変重要な性質として、「内容を変更できる」という点が挙げられます。（英語では mutable と言います。コンテナの中には後に紹介するように変更できないものもあります。）実際にどのように内容を変更するのか、例を見てみましょう。

>>> l = [1,2,3,4,5]
>>> l[0] = 28
>>> l
[28, 2, 3, 4, 5]
>>> l[2:4] = [3, 8]
>>> l
[28, 2, 3, 8, 5]

まず、一行目の入力で、l にリスト [1, 2, 3, 4, 5] を代入します。l[0] は前回説明したとおり、リスト l の一つ目の要素を表します。二つ目の入力で、l[0] = 28 とすると、一つ目の要素が 28 に置き換えられます。実際 l を表示してみると、一つ目の要素が置き換わっているのがわかります。

置き換えは、複数の要素に対して同時に行うこともできます。l [2:4} のスライスで表されるのは、前回学んだ通り、[3, 4] ですが、l[2:4] = [3, 8] とすることで、[3, 4] を [3, 8] に置き換えることができます。実際 l を表示してみると、3番目と4番目の要素が置き換わっているのが確認できます。

注：　同じ型を持つ数値データの集まりを扱うには、 numpy というモジュールが提供する array型を使うことがしばしばより効率的です。numpyについては後ほど詳しく取り上げます。

上記のような方法以外にもPythonにはリストを操作するためのさまざまな関数がデフォルトで組み込まれています。幾つか例を挙げてみましょう。まず、要素を追加するには次のように append を用います。

>>> l = [1, 2, 3, 4, 5]
>>> l.append(6)
>>> l
[1, 2, 3, 4, 5, 6]

l.append(6)とすることで、l の最後に新しい要素 6 が追加されます。リストの最後の要素削除するには pop() を用います。

>>> l.pop()
6
>>> l
[1, 2, 3, 4, 5]

複数の要素を一度に末尾に付け足すには extend() を用います。

>>> l.extend([6, 7])
>>> l
[1, 2, 3, 4, 5, 6, 7]

付け足した末尾2つの要素を取り除くには pop() を2回適用しても構いませんが、次のような方法もあります。

>>> l = l[:-2]
>>> l
[1, 2, 3, 4, 5]

l[:-2] は1番目の要素から最後から3番目までの要素を表すスライス表現です。つまり、l[:-2] は後ろ2つの要素が l から除かれたリストになるので、これを l に代入することで extend() で追加した2つの要素を削除することができます。

また、以下のように操作するとリストの要素を逆順に並べることができます。

>>> r = l[::-1]
>>> r
[5, 4, 3, 2, 1]

l[ : :-1] はスライスで step が -1 に指定されたものです。step を -1 にすることで最後尾の要素から逆順に一つづつ前に戻るという指示ができます。リストは + 演算子を使って結合することができます。

>>> r + l
[5, 4, 3, 2, 1, 1, 2, 3, 4, 5]

また、* 演算子を使うと、同じ内容を繰り返したリストを作ることができます。

>>> 2 * r
[5, 4, 3, 2, 1, 5, 4, 3, 2, 1]

sort() 関数を使うと、要素の数字を昇順に並べ直すことができます。先ほど逆順に並べた r をソートしてみましょう。

>>> r.sort()
>>> r
[1, 2, 3, 4, 5]

r.sort() のような書き方のプログラミングをオブジェクト指向プログラミングと呼びます。（オブジェクト指向プログラミングがどういうものかはここでは気にする必要はありません。）　今、r はリストなので、r に対しては、リストに対して有効なメソッドを適用することができます。iPythonシェルでは、r に対して有効な全メソッドを次のようにして表示することができます。

>>> l. # ここで tabキーをタイプ
 append count extend index insert pop remove reverse sort

メソッドを様々に駆使することでより効率的にリストの操作を行うことができます。ここで、全てのメソッドを紹介するわけにも行かないので、自分でいろいろ実験してみることをおすすめします。次回は、別のコンテナの例として、「文字列」についてまとめる予定です。

【第３回】 Pythonの基本的な型： コンテナ、その１ （リストの基本）

Pythonの型で重要なものの１つに「コンテナ」があります。コンテナというのは、オブジェクトというものの集合です。オブジェクトが何かというのは現時点では解りにくいかもしれませんが、例えば文字、数値などの場合もありますし、それ以外の場合もあります。

これからPythonを学習していく中で徐々にわかってくると思いますが、コンテナの扱いに慣れることがPythonを使いこなす上での肝となります。今回は、コンテナ型の初回として、コンテナ型の１つである「リスト」と、リストを操作する上で基本となる「インデクス」と「スライス」についてまとめてみます。

リストは順序つきのオブジェクトの集まりです。（順序付きというところが大切です。後に順序がつかないコンテナ型も登場します）。要素となるオブジェクトは異なる型を取ることができます。以下に例を示します。

>>> l = [1, 2, 3, 4, 5]
>>> type(l)
list

[1, 2, 3, 4, 5] がリストで、このリストを l に代入した上で type() で型を確認しています。リストは [ ] で括り、各要素はコンマで区切ります。

リストの中の要素を取り出したいときは、次のようにします。

>>> l[2]
3

l[2] をインデクス指定と呼びます。リストの中の要素には番号がふられており、番号を指定することで書く要素にアクセスすることができます。注意すべきなのは、インデクスの番号は 0, 1, 2, 3…、というように、0 から始まることです。上記の例の場合、l[2] なので、左から三番目の要素、つまり 3 が指定されます。

インデクスには負の数も使うことができます。以下の例を見てみましょう。

>>> l[-1]
5
>>> l[-2]
4

負号を付けると、右から何番目かを指定することができます。ただし、左からかぞえるときは、0, 1, 2, 3…、と 0 から数え始めましたが、右から数えるときは、-1, -2, -3…、と-1から始まることに注意してください。上記の例では、l[-1] は一番右、l[-2] は右から2番めの要素を表すので、5, 4 が指定されます。

インデクスを使うと、何番目から何番目の要素、というようにリストの一部分を指定できます。具体例で見てみましょう。

>>> l
[1, 2, 3, 4, 5]
>>> l[2:4]
[3, 4]

最初の入力行では確認のためリスト l の内容を表示しています。二つ目の入力行で範囲の指定をしています。この場合、3番目から4番目までの要素が指定されます。スライスの構文は、一般的には次のように書くことができます。

スライス構文： l[start:stop:step]

最初の start が範囲の最初の要素、stop が終わりの要素、step が増加量を表しています。step は省略すると 1 として扱われます。ここで覚えておきたいのは、部分を取り出す始めの要素が左から n番目のとき、start=n-1 となる一方で、終わりの要素が左から n 番目とすると、stop=n となることです。（この点、しばしば勘違いしてエラーの原因になりますのでしっかり覚えておきましょう。）

なぜこのような一見わかりにくそうな仕様になっているかというと、stop - start が取り出す要素の個数になるからです。ここでは解りにくいと思う方もいるかもしれませんが、Pythonでは要素の個数が分かりやすい方がプログラミングがやりやすいのでこのような仕様になっています。

以下、もう少しスライスの具体例を見てみましょう。

>>> l[3:]
[4, 5]
>>> l[:3]
[1, 2, 3]
>>> l[::2]
[1, 3, 5]

最初の入力 l[3:] は、start=3 で、stop と step が省略されています。この場合、範囲の最後は元のリストの最後までとなり、step は上述の通り 1 となります。

二つ目の入力 l[:3] は同様に、start=0, stop=3, step=1 と同義です。

三つ目の入力 l[::2] は、start と stop が省略されていて、step=2 です。この場合、元のリストの最初の要素から一つ置きに順次要素を取り出すことになります。 ここに上げた例以外にも、スライスは負の数を使った範囲指定もできますが、それはまた追々説明することにします。次回は、リストの操作の基本事項をまとめる予定です。

【第２回】 Pythonの基本的な型： 数値型

Pythonの変数にはいくつかの基本的な型があります。今回はそのうちの１つである「数値型」についてまとめます。数値型には以下の様なものがあります。

整数：

>>> 1 + 1
2
>>> a = 4

最後の行では、aに整数4を代入しています。


浮動小数点数：

>>> c = 2.1

小数点を入れて数値を書くと自動的に浮動小数点数となります。2.0の場合は、2.と省略して書くこともできます。上の例では、cに2.1という浮動小数点数を代入しています。


複素数：

>>> a = 1.5+0.5j
>>> a.real
1.5
>>> a.imag
0.5

複素数もPythonの組み込み型なので特別の配慮なしにいきなり使うことができます。虚数単位は数学では i ですが、Pythonでは j を用います。上の1行目で、a に複素数 1.5+0.5j を代入し、a.real で a の実部、a.imag で a の虚部を表示します。


ブール値：

>>> 3>4
False
>>> test=(3>4)
>>> test
False
>>> type(test)
bool

条件判断における真（true）か偽（false）をブール値と呼び、これも数値型の1つです。3>4 は偽なので、False と返ってきます。ブール値は文字列に代入することも可能です。上記の例では、test にFalse が代入されます。

type() を使うと型を調べることができます。上記の例では test にブール値が代入されているので、type(test) とすると bool と返ってきます。以下、type() の使用例をもう少し上げておきます。

>>> type(1)
int
>>> type(1.)
float
>>> type(1. + 0j)
complex
>>> a = 3
>>> type(a)
int

その他、数値型に関する基礎知識

数値型の種類は上記のとおりですが、さらに知っておいたほうが良い基礎知識を紹介します。

その（１）　基本的な演算子

一つ目は基本的な演算子（+, -, *, /, %）です。Python シェルは基本的な算術演算 +, -, *, /, % がそのまま使えるため以下の例のように電卓代わりに使うこともできます。

>>> 7 * 3.
21.0
>>> 2**10
1024
>>> 8 % 3
2

加減乗除は + - * / です。** と * を2つ続けるのは累乗を表します。上の例では$2^{10}$ということになります。% は剰余を表します。


その（２）　計算するときのコツ

3/2 を計算するとき、答えとして 1.5 を期待していると、以下の実行例のように、「オヤ？」という答えをPythonは返してきます。

>>> 3 / 2
1

これは、3 も 2 も整数なので答えも小数点以下が切り捨てられ整数で返ってくるためです。小数点以下を切り捨てないで完全に計算するためには以下のようにします。

>>> 3 / 2.
1.5

2の後にピリオドを付けることで期待した答えが得られます。3にピリオドをつけても同じ結果が得られますが、計算式中の1つの数値が浮動小数点型になれば他の数値も浮動小数点型として扱われるので省略してかまいません。

また、数値を a b c などの文字変数に代入した後に計算する場合は、float()  という関数を用いることで、以下の例のように整数を浮動小数に変換することができます。

>>> a = 3
>>> b = 2
>>> a / b
1
>>> a / float(b)
1.5

float() を単体で実行すると以下のようになります。

>>> float(1)
1.0

次回は、Pythonの重要な型の1つ「コンテナ」についてまとめます。