2015年6月25日木曜日

astropy.coordinatesを用いた座標の取り扱い、その7 ― 座標変換

astropy.coordinatesは、座標変換を行うための様々な機能を装備しています。天文学で用いられる主な座標系は予めastropyの中に組み込まれていますが、ユーザーが定義した特殊な座標系についても、独自の座標をastropy内で定義することによって他の座標系との間の座標変換を取り扱うことができます。独自の座標系の定義については別の機会に取り扱うとして、ここではastropy.coordinatesの中でどのように座標変換を実行できるかを見てゆくことにします。astropy.coordinatesに予め組み込まれた座標系と座標変換は、ここで確認することができます。

まず、一番単純な座標変換を見てみましょう。

>>> import astropy.units as u
>>> from astropy.coordinates import SkyCoord
>>> gc = SkyCoord(l=0*u.degree, b=45*u.degree, frame='galactic')
>>> gc.fk5  
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
    (229.272514629, -1.12844288043)>

より一般的には、transform_to() で変換を表します。transform_to() 用いることで入力方法が柔軟になり、なおかつ、より細かいオプションの設定が可能となります。

>>> from astropy.coordinates import FK5
>>> gc.transform_to('fk5')  
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
    (229.272514629, -1.12844288043)>
>>> gc.transform_to(FK5)  
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
    (229.272514629, -1.12844288043)>
>>> gc.transform_to(FK5(equinox='J1980.0'))  
<SkyCoord (FK5: equinox=J1980.000): (ra, dec) in deg
    (229.014693505, -1.05560349378)>

transform_to() には、SkyCoordオブジェクトをそのまま指定することもできます。SkyCoordオブジェクトを指定すると、そのSkyCoordオブジェクトで指定されている座標系に変換されます。

>>> sc = SkyCoord(ra=1.0, dec=2.0, unit='deg', frame=FK5, equinox='J1980.0')
>>> gc.transform_to(sc)  
<SkyCoord (FK5: equinox=J1980.000): (ra, dec) in deg
    (229.014693505, -1.05560349378)>

また、FK5, FK4, FK4NoETermsなどの座標系においては、同じ座標系の中での変換(例えば、分点を変更する)を行うことができます。

>>> fk5c = FK5('02h31m49.09s', '+89d15m50.8s')
>>> fk5c.equinox
<Time object: scale='utc' format='jyear_str' value=J2000.000>
>>> fk5c  
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
    (37.9545416667, 89.2641111111)>
>>> fk5_2005 = FK5(equinox='J2005')  # 分点を変更します
>>> fk5c.transform_to(fk5_2005)  
<SkyCoord (FK5: equinox=J2005.000): (ra, dec) in deg
    (39.3931763878, 89.2858442155)>

分点の変更は、以下のようにTimeオブジェクトを用いて行うこともできます。

>>> from astropy.time import Time
>>> fk5c = FK5('02h31m49.09s', '+89d15m50.8s',
...            equinox=Time('J1970', scale='utc'))
>>> fk5_2000 = FK5(equinox=Time(2000, format='jyear', scale='utc'))
>>> fk5c.transform_to(fk5_2000)  
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
    (48.0231710002, 89.386724854)>

2015年6月18日木曜日

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

実例1では、エイトフ図の上にランダムに点をプロットしましたが、今度はSkyCoordオブジェクトを使って、もう少し実際の星の分布に近い図を作ってみます。ここでは、多変量の正規分布(multivariate normal distribution)を用いて、銀河系のバルジとディスクの星の分布を再現してみます。まず最初に、作業に必要なパッケージを読み込みます。

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

ではまず、numpy.random.multivariate_normalを用いて、プロットするためのデータを作成します。生成されるデータは、直交座標系を用いて表されます。正規分布の、Z成分の幅を変えることで、銀河系のディスクとバルジを表現します。

>>> disk = np.random.multivariate_normal(mean=[0,0,0], cov=np.diag([1,1,0.5]), size=5000)
>>> bulge = np.random.multivariate_normal(mean=[0,0,0], cov=np.diag([1,1,1]), size=500)
>>> galaxy = np.concatenate([disk, bulge])

次に、生成したデータをSkycoordオブジェクトに変換し、直交座標系から赤道座標系に座標を変換します。

>>> c_gal = SkyCoord(galaxy, representation='cartesian', frame='galactic')
>>> c_gal_icrs = c_gal.icrs

実例1でやったのと同じように、角度の単位をラジアンに変換し、赤経を $-\pi$ から $\pi$ で表すように変換を行います。

>>> ra_rad = c_gal_icrs.ra.wrap_at(180 * u.deg).radian
>>> dec_rad = c_gal_icrs.dec.radian

最後に実例1と同じ設定で、matplotlibを用いてデータをプロットしてみます。

>>> import matplotlib.pyplot as plt
>>> plt.figure(figsize=(8,4.2))
>>> plt.subplot(111, projection="aitoff")
>>> plt.title("Aitoff projection of our random data")
>>> plt.grid(True)
>>> plt.plot(ra_rad, dec_rad, 'o', markersize=2, alpha=0.3)
>>> plt.subplots_adjust(top=0.95,bottom=0.0)
>>> plt.show()


実例2

2015年6月16日火曜日

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

ここでは、SkyCoordオブジェクトを使って、ランダムに生成された座標値を、エイトフ図(Aitoff projection)上にプロットするという作業をやってみます。 

まず最初に必要なパッケージを読み込みます。ここでは、関連するastropyのパッケージ以外に、図を作成するためのパッケージである matplotlib、ランダムな数値データを生成するための numpy を読み込みます。

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

まず、座標値をランダムに生成し、プロットするためのデータを作成します。RA(赤経)は、0°から360°の間で値を生成し、ra_random という名前で保存します。赤緯(DEC)に関しては、-90°から+90°の間で値を生成し、dec_random という名前で保存します。

>>> ra_random = np.random.rand(100)*360.0 * u.degree
>>> dec_random = (np.random.rand(100)*180.0-90.0) * u.degree

次に ra_random、dec_random から、SkyCoordオブジェクトを生成します。

>>> c = SkyCoord(ra=ra_random, dec=dec_random, frame='icrs')

ここで使用する matplotlib という作図パッケージでは、単位として $-\pi$ から $\pi$ の間のラジアンが表記($0$ から $2 \pi$ ではない)が用いられるので、表記を変換する必要があります。この変換は、astropy.coordinate.Angle のラッピング機能を用いて行います。

>>> ra_rad = c.ra.wrap_at(180 * u.deg).radian
>>> dec_rad = c.dec.radian

最後に、matplotlibを用いて、エイトフ図法(Aitoff projection)でデータ点をプロットします。ここでは、タイトルを図の上に置き、座標を表すグリッドを加えます。データ点はサイズ2の塗りつぶした丸とします。(他にも細部の設定は、matplotlib のテキストを参照のこと。)

>>> import matplotlib.pyplot as plt
>>> plt.figure(figsize=(8,4.2))
>>> plt.subplot(111, projection="aitoff")
>>> plt.title("Aitoff projection of our random data")
>>> plt.grid(True)
>>> plt.plot(ra_rad, dec_rad, 'o', markersize=2, alpha=0.3)
>>> plt.subplots_adjust(top=0.95,bottom=0.0)
>>> plt.show()


出来上がり図

2015年6月10日水曜日

astropy.coordinatesを用いた座標の取り扱い、その4 ― SkyCoodクラスの使用法(上級編)

その4の最後で、SkyCoordオブジェクトの座標系を銀河座標系(galactic)へ変換する方法を紹介しました。一旦この変換を行うと、経度と緯度の座標値は、銀河座標系を表す際の慣習に従ってSkyCoordオブジェクト内で銀系と銀緯に対応した、、というラベルが付けられます。天文学の背景知識を持っているユーザーにとっては、(l, b) だけで何を意味しているか即座にわかると思いますが、SkyCoordオブジェクトには、(l, b) の意味も保存されており、状況に応じて適宜情報を取り出すことができます。

>>> sc_gal.representation_component_names
OrderedDict([(u'l', u'lon'), (u'b', u'lat'), (u'distance', u'distance')])

>>> sc_gal.representation_component_units
OrderedDict([(u'l', Unit("deg")), (u'b', Unit("deg"))])

>>> sc_gal.representation
<class 'astropy.coordinates.representation.SphericalRepresentation'>


上記のようなコマンドラインによって、(l, b) が経度と緯度を表していること、(l, b) の単位が度で有ること、直交座標系ではなく球面座標系であることが確認できます。また、上記で使用した representation_component_names という attribute を使うと、SkyCoodオブジェクトの中で使用する座標値の独自の呼称を定義することもできます。

別の重要な attribute としては、frame_attr_names があります。SkyCoordオブジェクトを作成する際、SkyCoordオブジェクトの中でどのようなパラメータが設定可能か確認したい場合があります。この attribute を使うと、簡単に未設定のパラメータを調べることができます。

>>> sc_fk4 = SkyCoord(1, 2, frame='fk4', unit='deg')
>>> sc_fk4.get_frame_attr_names()  
{u'equinox': 

明示的にパラメータを設定していない場合はデフォルト値が設定される場合もあります。上の例では、1行目でFK4座標系におけるSkyCoordオブジェクトを定義していますが、ユーザー側で完全に全パラメータを明示的に定義するには、equinoxobstime というパラメータに値を与える必要があることが分かります。

SkyCoordオブジェクトのパラメータ設定を行う際に若干ややこしのは、SkyCoordオブジェクトが、frameオブジェクトなど、いくつかの別のオブジェクトを土台として設計されていることです。このため、例えばSkyCoordオブジェクトの attribute で実行できる作業が、他の attribute に内包される別の attribute でも実行できてしまう場合があります。以下に例を示します。


>>> sc.frame
<ICRS Coordinate: (ra, dec) in deg
    (1.0, 2.0)>

>>> sc.has_data is sc.frame.has_data
True

>>> sc.frame.  
sc.frame.cartesian                           sc.frame.ra
sc.frame.data                                sc.frame.realize_frame
sc.frame.dec                                 sc.frame.represent_as
sc.frame.default_representation              sc.frame.representation
sc.frame.distance                            sc.frame.representation_component_names
sc.frame.frame_attr_names                    sc.frame.representation_component_units
sc.frame.frame_specific_representation_info  sc.frame.representation_info
sc.frame.get_frame_attr_names                sc.frame.separation
sc.frame.has_data                            sc.frame.separation_3d
sc.frame.is_frame_attr_default               sc.frame.shape
sc.frame.is_transformable_to                 sc.frame.spherical
sc.frame.isscalar                            sc.frame.time_attr_names
sc.frame.name                                sc.frame.transform_to
>>> sc.frame.name
'icrs'

この様に、has_data というSkyCoord自身の attribute でも frame という attribute以下に存在する、has_data でも全く同じ結果が得られます。このような性質は、知っておくと、利用方法によっては便利な場合もあります。

座標変換

座標変換については、実際に天文学のデータを扱う際に重要な項目でもあるので、場所を改めて詳しく紹介します。ここでは、SkyCoordクラスの概要の一部として、簡単な例をいくか紹介しておきます。

>>> from astropy.coordinates import FK5
>>> sc = SkyCoord(1, 2, frame='icrs', unit='deg')
>>> sc.galactic  
<SkyCoord (Galactic): (l, b) in deg
    (99.6378552814, -58.7096929334)>

>>> sc.transform_to('fk5')  # sc.fk5、sc.transform_to(FK5)としても同様の結果が得られます。  
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
        (1.00000655566, 2.00000243092)>

>>> sc.transform_to(FK5(equinox='J1975'))  # FK5において分点を変更しています。  
<SkyCoord (FK5: equinox=J1975.000): (ra, dec) in deg
        (0.679672818323, 1.86083014099)>

2つのSkyCoordオブジェクトがあった時、2つのオブジェクト間で座標系を一致させる操作は次のような方法で簡単に行えます。

>>> sc2 = SkyCoord(3, 4, frame='fk4', unit='deg', obstime='J1978.123', equinox='B1960.0')
>>> sc.transform_to(sc2) 
<SkyCoord (FK4: equinox=B1960.000, obstime=J1978.123): (ra, dec) in deg
    (0.48726331438, 1.77731617297)>

球面座標以外の座標の利用

ここまでの例では、ずっと球面座標を使ってきました。SkyCoordに組み込まれている座標系ではデフォルトが球面座標なので、球面座標の扱いを知っていれば、実際多くの天文学上の仕事を処理することができます。しかし、球面座標以外にも、直交座標や円柱座標の扱い方を知っておくと便利な場面は少なくありません。ここでは、SkyCoordオブジェクトで使用する座標をどうやって変更するのか簡単に紹介します。座標に関するより詳しい説明は、Using and Designing Coordinate Representationsを見て下さい。

初期化

SkyCoordオブジェクトの座標を扱う際に必要なことは、以下の実例を見ればほぼ十分だと思います。座標を初期化するのに必要なのは、representationというキーワードに、使用する座標を書き込むだけです。直交座標を使う場合の例は以下のとおりです。

>>> c = SkyCoord(x=1, y=2, z=3, unit='kpc', 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>)
他に、円柱座標を使ったり、各座標軸の単位を変更したりするには以下のようにします。
>>> SkyCoord(1, 2*u.deg, 3, representation='cylindrical')
<SkyCoord (ICRS): (rho, phi, z) in (, deg, )
    (1.0, 2.0, 3.0)>

>>> SkyCoord(rho=1*u.km, phi=2*u.deg, z=3*u.m, representation='cylindrical')
<SkyCoord (ICRS): (rho, phi, z) in (km, deg, m)
    (1.0, 2.0, 3.0)>

>>> SkyCoord(rho=1, phi=2, z=3, unit=(u.km, u.deg, u.m), representation='cylindrical')
<SkyCoord (ICRS): (rho, phi, z) in (km, deg, m)
    (1.0, 2.0, 3.0)>

>>> SkyCoord(1, 2, 3, unit=(None, u.deg, None), representation='cylindrical')
<SkyCoord (ICRS): (rho, phi, z) in (, deg, )
    (1.0, 2.0, 3.0)>

SkyCoordオブジェクトのフォーマットを一般的な形で書き下すと以下のようになります。

SkyCoord(COORD, [FRAME | frame=FRAME], [unit=UNIT], [representation=REPRESENTATION],
         keyword_args ...)
SkyCoord(COMP1, COMP2, [COMP3], [FRAME | frame=FRAME], [unit=UNIT],
         [representation=REPRESENTATION], keyword_args ...)
SkyCoord([FRAME | frame=FRAME], =COMP1, =COMP2,
         =COMP3, [representation=REPRESENTATION], [unit=UNIT],
         keyword_args ...)