こんな人にオススメ
数値計算する際に「m」と「mm」とかの換算が面倒。後で見返したときにどこでどういう変換をしているのかがわからない！自動で単位を計算してくれるようないい方法ってないの？
ということで、今回は執筆者が普段から使用している
astropyを用いて数値に単位を付与する方法について紹介していく。これを知ってから単に対しての心配事が減った。
単位換算が大変
執筆者は理系出身で天文学系の専攻をしていた。なので、単位換算をする場面が多々あった。例えば
- 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
毎回毎回これを行うのは面倒だし単位換算を行ったかどうか忘れてしまう。さらに、複雑な計算が混ざるものだとさらに面倒。例えば角度の単位換算。以下では一旦arcsecをdegに変換した後にradに変換している。
# 以下の値は秒角単位で書いている val_arcsec = 1.2345 # 以下の値はradの単位で書いている val_rad = (val_arcsec / 3600) * (np.pi / 180)
astropy.unitsがある時
単位付与
一方でastropy.unitsがある時は上記の心配はいらない。astropy.unitsは以下の2種類のパーツに分かれており、これらのパーツを組み合わせることで単位を付与することができる。
value: 値の部分unit: 単位の部分
以下で実際に試してみる。単位の付与はvalueとunitの掛け算で行うことができる。なお、単位を付与したい値は全て以下のものを使用する。
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とすればいい。辞書型などではvaluesとsが末尾に着くが、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の単位付与と変換について紹介した。執筆者自身この単位付与を使用してから計算ミスに対して不安が軽減されたし、いちいち計算を追わなくても結果に既に単位が書いているので余計なエネルギーを使うこともなくなった。皆様もこの機会に変数に単位を付与してみてはいかがだろう。
