カテゴリー

当サイトはアフィリエイトプログラムによる収益を得ています〈景品表示法に基づく表記です)

astropy

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

2021年4月8日

こんな人にオススメ


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

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

運営者のメガネです。YouTubeTwitterInstagram、自己紹介はこちら、お問い合わせはこちらから。

運営者メガネ

Pro天パで詳しく学ぶ

Pro天パにてリライト

本記事の内容はプログラミング専用のサイト「Pro天パ」の以下の記事でリライトを兼ねて解説している。より詳しく学びたい人はこちらを読んでほしい。

単位換算が大変

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

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

ガジェット

2023/11/11

【デスクツアー2022下半期】モノは少なく、でも効率的に Desk Updating #0

今回はガジェットブロガーなのにデスク環境を構築していない執筆者の ...

ライフハック

2023/9/16

【Audible vs YouTube Premium】耳で聴く音声学習コンテンツを比較

ワイヤレスイヤホンが普及し耳で学習することへのハードルが格段に下 ...

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

2023/9/18

【SENNHEISER MOMENTUM True Wireless 3レビュー】全てが整ったイヤホン

今回は高音質・高機能なSENNHEISERのフラグシップ完全ワイヤレスイヤホン「SENNH ...

ライフハック

2023/3/11

【YouTube Premiumとは】メリットしかないから全員入れ

今回はYouTube Premiumを実際に使ってみてどうなのか、どんなメリット/デメリット ...

マウス

2023/11/4

【Logicool MX ERGOレビュー】疲れない作業効率重視トラックボールマウス

こんな人におすすめ トラックボールマウスの王道Logicool MX ERGOが気になるけどऩ ...

ベストバイ

2023/10/29

【ベストバイ2022】今年買って良かったモノのトップ10

2022年ベストバイ この1年を振り返って執筆者は何を買ったのか。ガジェッ& ...

スマホ

2023/1/15

【楽天モバイル×povo2.0の併用】月1,000円の保険付きデュアルSIM運用

こんな人におすすめ 楽天モバイルとpovo2.0のデュアルSIM運用って実際のとこ ...

マウス

2023/9/16

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

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

macOSアプリケーション

2022/9/30

【Chrome拡張機能】便利で効率的に作業できるおすすめの拡張機能を18個紹介する

こんな人におすすめ Chromeの拡張機能を入れたいけど、調べても同じような ...

macOSアプリケーション

2023/5/3

【Automator活用術】Macで生産性を上げる作業の自動化術

今回はMacに標準でインストールされているアプリ「Automator」を使ってできる ...

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

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

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

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

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

  • この記事を書いた人

メガネ

Webエンジニア駆け出し。独学のPythonで天文学系の大学院を修了。常時金欠のガジェット好きでM2 Pro MacBook Pro(31万円) x Galaxy Z Fold5(25万円)使いの狂人。自己紹介と半生→変わって楽しいの繰り返しレビュー依頼など→お問い合わせ運営者情報、TwitterX@m_ten_pa、 YouTube@megatenpa、 Threads@megatenpa、 Instagram@megatenpa

-astropy
-,