カテゴリー

astropy

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

2021年4月8日

こんな人にオススメ


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

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

運営者のメガネとです。YouTubeTwitterInstagramも運営中。

自己紹介はこちらから、お問い合わせはこちらからお願いいたします。

運営者メガネ

単位換算が大変

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

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

スイッチボット

2022/11/28

【SwitchBotロックレビュー】これからのスタンダードになりうるスマートロック

こんな人にオススメ SwitchBotからスマートロック「SwitchBotロック」が発売された ...

生活に役立つ

2022/11/28

【メガネ厳選】クソ便利に使っているサービスやアイテム達

このページでは執筆者「メガネ」が実際に使って便利だと感じているサ ...

マウス

2022/9/11

【Logicool MX ERGO vs MX Master 3】ERGOをメインにした決定的な理由

こんな疑問・お悩みを持っている人におすすめ 執筆者はLogicoolのハイエンӠ ...

完全ワイヤレスイヤホン(TWS)

2022/11/21

【ながら聴きイヤホン比較】SONY LinkBuds、ambie、BoCoはどれがおすすめ?

こんな人におすすめ 耳を塞がない開放型のイヤホンに完全ワイヤレスӟ ...

macOSアプリケーション

2022/10/15

【M1 Mac】MacBook Proに入れている便利でニッチなアプリを21個紹介する

こんな人におすすめ MacBookを購入してLINEとか必要最低限のアプリは入れた。 ...

完全ワイヤレスイヤホン(TWS)

2022/10/23

【SENNHEISER MOMENTUM True Wireless 3レビュー】高レベルでバランス型の高音質イヤホン

こんな人におすすめ SENNHEISER MOMENTUM True Wireless 3って実際のところどうなの? 評判は良い ...

完全ワイヤレスイヤホン(TWS)

2022/11/21

【SONY WF-1000XM4レビュー】神とゴミのハーフ&ハーフ

こんな人におすすめ SONYのフラグシップモデル「SONY WF-1000XM4」ってどれくらい性 ...

完全ワイヤレスイヤホン(TWS)

2022/8/19

【Nothing ear (1)レビュー】ライトな完成度、アップデートに期待

こんな人にオススメ 完全ワイヤレスイヤホン(TWS)でスケルトンボディ ...

Pythonを学びたいけど独学できる時間なんてない人へのすゝめ

執筆者は大学の研究室・大学院にて独学でPythonを習得した。

でも社会人になったら独学で行うには時間も体力もなくて大変だ。

時間がない社会人だからこそプロの教えを乞うのが効率的。

ここでは色んなタイプに合ったプログラミングスクールの紹介をする。

  • この記事を書いた人

メガネ

ベンチャー企業のWebエンジニア駆け出し。独学のPythonで天文学系の大学院を修了→新卒を1.5年で辞める→転職→今に至る。
常時金欠のガジェット好きでM1 MacBook Pro x Galaxy S22 Ultraの狂人。
人見知りで根暗だったけど、人生楽しもうと思って良い方向に狂う→人生が楽しい

ガジェットのレビューとPythonコードを記事にしています。ぜひ楽しんでください🦊
自己紹介と半生→変わって楽しいの繰り返し

-astropy
-, ,