カテゴリー

astropy

【astropy&単位】astropyで変数に単位を付与

2021年4月8日

こんな人にオススメ


数値計算する際に「m」と「mm」とかの換算が面倒。後で見返したときにどこでどういう変換をしているのかがわからない!自動で単位を計算してくれるようないい方法ってないの?

ということで、今回は執筆者が普段から使用しているastropyを用いて数値に単位を付与する方法について紹介していく。これを知ってから単に対しての心配事が減った。

スポンサーリンク
スポンサーリンク

運営者のメガネです。TwitterInstagramも運営してます。

自己紹介はこちらから、お問い合わせはこちら。

運営者メガネ

単位換算が大変

執筆者は理系出身で天文学系の専攻をしていた。なので、単位換算をする場面が多々あった。例えば

  • 1 Jy(ジャンスキー)と10-26 W m-2 Hz-1
  • 1 erg(エルグ)と10-7 J
  • 1''(1秒角)と4.85×10-6 rad

などなど。これらの単位換算を手計算で行うのは面倒だし計算間違いの危険性もある。例えば二つ目のergからJの関係を10-7ではなく107にしてしまうなど。

今回はこれを解消するために執筆者が使っている「astropy.unis」について紹介する。astropy.unitsでは値に単位を付与することができる。

python環境は以下。

  • Python 3.9.2
  • numpy 1.20.1
  • astropy 4.2

astropy.unitsがない時

astropy.unitsがない時は前述の通り、自分で計算しないといけない。例えば以下のように。

# 以下の値はerg単位で書いている
val_erg = 1.2345
# 以下の値はJ単位で書いている
val_joule = val_erg * 1E-7

毎回毎回これを行うのは面倒だし単位換算を行ったかどうか忘れてしまう。さらに、複雑な計算が混ざるものだとさらに面倒。例えば角度の単位換算。以下では一旦arcsecdegに変換した後にradに変換している。

# 以下の値は秒角単位で書いている
val_arcsec = 1.2345
# 以下の値はradの単位で書いている
val_rad = (val_arcsec / 3600) * (np.pi / 180)

astropy.unitsがある時

単位付与

一方でastropy.unitsがある時は上記の心配はいらない。astropy.unitsは以下の2種類のパーツに分かれており、これらのパーツを組み合わせることで単位を付与することができる。

  • value: 値の部分
  • unit: 単位の部分

以下で実際に試してみる。単位の付与はvalueunitの掛け算で行うことができる。なお、単位を付与したい値は全て以下のものを使用する。

val = 1.2345

それでは単位を付与していこう。出力は「#」で示す。

from astropy import units as u

# erg単位
erg = u.erg
val_erg = val * erg
print(val_erg)
# 1.2345 erg
print(type(val_erg))
# <class 'astropy.units.quantity.Quantity'>

# arcsec単位
arcsec = u.arcsec
val_arcsec = val * arcsec
print(val_arcsec)
# 1.2345 arcsec
print(type(val_arcsec))
# <class 'astropy.units.quantity.Quantity'>

このように単位を付与することができる。値を取り出したい時、単位を取り出したい時はそれぞれ.value.unitとすればいい。辞書型などではvaluessが末尾に着くが、astropy.unitsの場合はsがつかないのでそこは注意。

print(f"value: {val_erg.value}, unit: {val_erg.unit}")
# value: 1.2345, unit: erg
print(f"value: {type(val_erg.value)}, unit: {type(val_erg.unit)}")
# value: <class 'numpy.float64'>, unit: <class 'astropy.units.core.Unit'>

print(f"value: {val_arcsec.value}, unit: {val_arcsec.unit}")
# value: 1.2345, unit: arcsec
print(f"value: {type(val_arcsec.value)}, unit: {type(val_arcsec.unit)}")
# value: <class 'numpy.float64'>, unit: <class 'astropy.units.core.Unit'>

単位換算(1つの文字列の単位への変換)

本題の単位換算についてだが、大きく分けると3種類の方法がある。

  • .to('J').to('rad')など、直接単位を指定
  • .to(u.J).to(u.rad)など、astropy.units.core.Unitの形式で単位を指定
  • .to(u.Unit('rad'))など、単位を文字列で指定

それぞれの方法で実際に単位を換算してみる。

lst_joule = [val_erg.to(u.J), val_erg.to('J'), val_erg.to(u.Unit('J'))]
for i in lst_joule:
    print(f"ans: {i}, unit: {i.unit}")
# ans: 1.2345e-07 J, unit: J
# ans: 1.2345e-07 J, unit: J
# ans: 1.2345e-07 J, unit: J

lst_rad = [val_arcsec.to(u.rad), val_arcsec.to('rad'), val_arcsec.to(u.Unit('rad'))]
for i in lst_rad:
    print(f"ans: {i}, unit: {i.unit}")
# ans: 5.985024893297221e-06 rad, unit: rad
# ans: 5.985024893297221e-06 rad, unit: rad
# ans: 5.985024893297221e-06 rad, unit: rad

このように他に換算を簡単に行うことができる。

単位換算(2つ以上の文字列の単位への換算)

一方でJyからW m-2 Hz-1への変換のように、1つの文字列への単位変換ではない場合は単位換算の3種類の2つ目のastropy.units.core.Unitの形式で換算するということができない。今回の例ではJyを使用する。

# Jy単位
jansky = u.jansky
val_jansky = val * jansky
print(val_jansky)
# 1.2345 Jy
print(type(val_jansky))
# <class 'astropy.units.quantity.Quantity'>

JyをSI単位系に変換する。

lst_jansky_si = [
    val_jansky.to('W m-2 Hz-1'),
    val_jansky.to(u.Unit('W m-2 Hz-1')),
]
for i in lst_jansky_si:
    print(f"ans: {i}, unit: {i.unit}")
# ans: 1.2345e-26 W / (Hz m2), unit: W / (Hz m2)
# ans: 1.2345e-26 W / (Hz m2), unit: W / (Hz m2)

簡単に変換することができる。では、複数文字列の単位はどのように書いたら良いのか。上記の例では単位間をスペースで区切り、単位の累乗部分はそのまま並べて書いた。以下のコードでどの書き方ならエラーが出ないのかを検証した。

unit_tpl = (
    'W m-2 Hz-1',
    'Wm-2Hz-1',

    'W m^-2 Hz^-1',
    'W m^{-2} Hz^{-1}',

    'W (m2 Hz)-1',
    'W (m2 Hz)^{-1}',

    'W/m2/Hz',
    'W/m2/Hz1',

    'W /(m2 Hz)',
    'W /(m2 Hz1)',

    'W /m2 /Hz',
    'W /m2 /Hz1',

    'W/m2/Hz',
    'W/m2/Hz1',
)

# 単位として読み込むことができる文字列を入れる
ok_lst = []
# 単位として読み込むことができない文字列を入れる
ng_lst = []
for i in unit_tpl:
    try:
        ans = val_jansky.to(i)
        ok_lst.append(i)
    except ValueError:
        ng_lst.append(i)

ここで、このコードを動かすと以下の警告が出る。

WARNING: UnitsWarning: 'W/m2/Hz' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
WARNING: UnitsWarning: 'W/m2/Hz1' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
WARNING: UnitsWarning: 'W /m2 /Hz' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
WARNING: UnitsWarning: 'W /m2 /Hz1' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]

ざっと解説すると、スラッシュ/が複数あるとFITS推奨ではないということ。FITS形式については「天文学辞典(日本天文学会)」の「FITS形式」の項目に以下のようにある。

天体データ形式の一種。天体画像データを主として、その他の多くの天文データ(分光データ、高エネルギーイベントデータ、カタログ表データなど)を格納するために現在最も一般的に用いられる形式。FITSはflexible image transport systemの頭文字からきた名称。

要するにファイル形式ってこと。「〇〇.fits」形式でファイルが保存される。また、astropy.units.format.genericを含むastropyの単位については他の記事で書く予定なのでお待ちを。

単位として読み込むことができる文字列を入れたok_lstは以下。すなわち使える書き方。

print(*ok_lst, sep='\\n')
# W m-2 Hz-1
# W m^-2 Hz^-1
# W/m2/Hz
# W/m2/Hz1
# W /(m2 Hz)
# W /(m2 Hz1)
# W /m2 /Hz
# W /m2 /Hz1
# W/m2/Hz
# W/m2/Hz1

逆に単位として読み込むことができない文字列を入れたng_lstは以下。すなわちエラーを吐かれる文字列。

print(*ng_lst, sep='\\n')
# Wm-2Hz-1
# W m^{-2} Hz^{-1}
# W (m2 Hz)-1
# W (m2 Hz)^{-1}

これらの検証により以下のことがわかる。

  • スラッシュ/を複数個つけることで負の符号の単位を表すのは推奨されない
  • Wm-2Hz-1: 単位同士をくっつけるとエラー
  • W m^{-2} Hz^{-1}: 累乗を^{数字}の形式で表すのはエラー
  • W (m2 Hz)-1: 負の符号の単位をまとめて-1乗するのはエラー
  • W (m2 Hz)^{-1}: ^{-1}で負の符号をまとめてもエラー

執筆者的に一番簡単と思うのはW m-2 Hz-1のように負の符号を単位の横にそのまま並べて、単位全体で分数を使わずに並べる方法。シンプルでいい気がする。

いろいろと単位換算してみる

最後にここではastropy.unitsを使った単位換算を色々と試してみる。

単位同士の掛け算や割り算

いきなり単位換算ではないが、単純に単位を持った変数同士の掛け算や割り算はどうなるのか。単純に計算される、ただそれだけ。ここでは以下の2種類の計算を行う。

  • 速さ(velocity)と時間(time)から距離(distance)を求める
  • 距離(distance)と時間(time)から速さ(velocity)を求める
# 速さ: 20 m/s
velocity = 20 * (u.m / u.s)
# 時間: 5 sec
time = 5 * u.s
# 距離: 100 m
distance = velocity * time
print(distance)
# 100.0 m

# 距離: 10 m
distance = 10 * u.m
# 時間: 4 sec
time = 4 * u.s
# 速さ: 2.5 m/s
velocity = distance / time
print(velocity)
# 2.5 m / s

累乗があっても問題ない。面積を求めてみる。

# 長さ: 100 m
length = 100 * u.m
# 面積
area = length ** 2
print(area)
# 10000.0 m2

単位換算しまくる

ここでは色々単位換算をしてみる。まずは上述の速度の単位をm/sからcm/sに変換。

# 距離: 10 m
distance = 10 * u.m
# 時間: 4 sec
time = 4 * u.s
# 速さ: 2.5 m/s
velocity = distance / time
print(velocity)
print(velocity.to('cm/s'))
# 2.5 m / s
# 250.0 cm / s

次は地球と太陽の間の平均距離を表す「天文単位」。これをkmに変換してみる。

# 1天文単位
au = 1 * u.au
au_km = au.to('km')
print(f"{au_km:.3e}")
# 1.496e+08 km

最後に波長を周波数に変換する。この場合は長さと周波数の単位なので直接換算ができない。そこで.toの引数であるequivalenciesを光の関係であると宣言してあげると解決する。光速$c$、波長$\lambda$、周波数$f$とすると以下の関係が成り立つ。

$$c=f\lambda$$

# 波長を周波数に
wvlngth = 500 * u.nm
hz = (wvlngth).to(u.Hz, equivalencies=u.spectral())
print(f"{hz:.3e}")
# 5.996e+14 Hz

単位換算サイコー

今回はastropyの単位付与と変換について紹介した。執筆者自身この単位付与を使用してから計算ミスに対して不安が軽減されたし、いちいち計算を追わなくても結果に既に単位が書いているので余計なエネルギーを使うこともなくなった。皆様もこの機会に変数に単位を付与してみてはいかがだろう。

関連コンテンツ

スポンサーリンク

Amazonのお買い物で損したない人へ

1回のチャージ金額通常会員プライム会員
¥90,000〜2.0%2.5%
¥40,000〜1.5%2.0%
¥20,000〜1.0%1.5%
¥5,000〜0.5%1.0%

Amazonギフト券にチャージすることでお得にお買い物できる。通常のAmazon会員なら最大2.0%、プライム会員なら2.5%還元なのでバカにならない。

ゲットしたポイントは通常のAmazonでのお買い物に使えるからお得だ。一度チャージしてしまえば、好きなタイミングでお買いものできる。

なお、有効期限は10年だから安心だ。いつでも気軽にAmazonでお買い物できる。

Amazonチャージはここから出来るで

もっとお得なAmazon Prime会員はこちらから

30日間無料登録

執筆者も便利に使わせてもらってる

スポンサーリンク

  • この記事を書いた人

メガネ

独学でpythonを学び天文学系の大学院を修了。 ガジェット好きでMac×Android使い。色んなスマホやイヤホンを購入したいけどお金がなさすぎて困窮中。 元々、人見知りで根暗だったけど、人生楽しもうと思って良い方向に狂ったために今も人生めちゃくちゃ楽しい。 pythonとガジェットをメインにブログを書いていますので、興味を持たれましたらちょこちょこ訪問してくだされば幸いです🥰。 自己紹介→変わって楽しいの繰り返し

-astropy
-, ,