2015年4月13日月曜日

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

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

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

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

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

wrap_angleで、ラッピングする範囲を指定します。例えば、角度を0°から360°の範囲で表示したい場合、wrap_angleは360°、-180°から180°の範囲で表したい場合は wrap_angleは180°となります。(つまり数式で表すと、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

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