カラーバーを等級で表したい

Question

fig = plt.figure(figsize=(7.5, 7.5))
ax = fig.add_subplot(111)

X, Y = np.meshgrid(x[0:len(x)],y[0:len(y)]) #21 x 21
F = np.array(F).reshape(len(x)-1,-1) #20 x 20

m = ax.pcolormesh(X, Y, F, cmap='gray')
fig.colorbar(m, label='A')

plt.show()

用いるデータの一部を以下に示します。
列としては等級、フラックス、X, Yで並んでいます。
本来のデータは数十万行あります。

mag Flux X Y
16.78 0.0007045415735661775 329.399 42.232
16.806 0.0006878704054919076 329.023 42.052
13.317 0.017104351979956002 329.27 42.124
14.469 0.005919794420515654 328.675 41.868
15.511 0.0022672876636957917 329.196 41.957
16.706 0.0007542349013940496 329.374 42.217
16.051 0.0013788169862608637 329.501 42.265
16.962 0.0005958084044128775 329.141 42.167
15.238 0.0029154568968192433 328.801 41.978
16.209 0.0011920833578792421 329.246 42.019
14.663 0.004951156169049099 329.108 41.977
15.384 0.002548629026672342 329.15 42.032
15.111 0.0032772277608275092 328.878 41.979
15.681 0.0019386822097603555 329.284 42.092
16.151 0.0012574960369579535 329.339 42.157

上記のデータをnp.genfromtxtで読み込み、m = data[:,0], F_m = data[:,1]という風に読み込むと

m = np.array([16.78 , 16.806, 13.317, 14.469, 15.511, 16.706, 16.051, 16.962, 15.238, 16.209, 14.663, 15.384, 15.111, 15.681, 16.151])
F_m = np.array([0.00070454, 0.00068787, 0.01710435, 0.00591979, 0.00226729, 0.00075423, 0.00137882, 0.00059581, 0.00291546, 0.00119208, 0.00495116, 0.00254863, 0.00327723, 0.00193868, 0.0012575 ])
X = np.array([329.399, 329.023, 329.27 , 328.675, 329.196, 329.374, 329.501, 329.141, 328.801, 329.246, 329.108, 329.15 , 328.878, 329.284, 329.339])
Y = np.array([42.232, 42.052, 42.124, 41.868, 41.957, 42.217, 42.265, 42.167, 41.978, 42.019, 41.977, 42.032, 41.979, 42.092, 42.157])

となると思います。

少し天文学の知識が必要になる質問です。

上記のコードはデータなどかなり飛ばして抜粋したものですが、やりたいこととしては、
横Xと縦Yの値があり、その間に入る高さFの値を2次元マップで書きたいというものです。
Xの値は328330くらい、Yの値は4143くらい、Fの値は0.0, 0.001...くらいで全て小数(float)です。

ただ単にグラフを描くだけであればこれで終わりなのですが、厄介なことにこのFはフラックスで、描画自体はフラックスのままでいいのですが最終的にカラーバーでフラックスを等級に換算して見やすくしたいです。

なお、Fの中にはあるX, Yの中にデータがないということを意味する0も含まれており、さらに以下に示すフラックスと等級の換算の式で等級がたまたま0になったという場合もあります。

うまくして、フラックスFが元から0であった箇所は真っ黒で、かつ、等級を計算した結果0になったものも含めて等級の値が小さくなると白く、値が大きくなると黒に近づくというカラーバーにしたいです。
おそらくマップ自体はフラックスの値が大きい方が等級が小さくなる、すなわち明るく、フラックスでマップを描いた時に白で表され、逆もまた然りなので大丈夫かなと思っています。

フラックス→等級の換算式
m = -2.5 * log10(F_m/F_0)
m : F_mに対応する等級 /mag
F_m : フラックス /Jy(ジャンスキー)
F_0 : AB等級でのVEGAのフラックス(=3630 /Jy)

難題かもしれませんがよろしくお願いいたします。

Comments

[yamato_user]

難題ではないが、結局したいことは「フラックスFが元から0であった箇所は真っ黒で、かつ、等級を計算した結果0になったものも含めて等級の値が小さくなると白く、値が大きくなると黒に近づくというカラーバーにしたいです。」ということですか？

[mega_ne]

そうなります。

[yamato_user]

> 上記のコードはデータなどかなり飛ばして抜粋したものですが

ちゃんとデータを載せてください。グラフ描画用のデータ作成も回答者がやるのはめんどくさすぎます。

[mega_ne]

すみません。
修正いたしました。

[yamato_user]

DataFrameかNumpy配列で

[mega_ne]

すみません。
データの書き方を修正いたしました。

Answer 1

完全に仕様を理解しているか？？なのですが。

まず、「等級を計算した結果0になったものも含めて等級の値が小さくなると白く、値が大きくなると黒に近づく」とありますので、等級を計算した結果を mとすると、

c = ax.pcolormesh(X, Y, m, cmap='gray_r', vmin=0, vmax=20)

のようにcmapを gray_r を使い vmin と vmax パラメータにて範囲を指定して記述することで実現できます。

また、「フラックスFが元から0であった箇所は真っ黒で」の部分に関しては numpyにて フラックス→等級を計算すると Fの値が0の時は、-2.5*np.log10(0) にて結果が np.inf となり上記の vmax 以上の値となりますので特に何の処理も行わなくても自動的に黒にるかと思います。

以下簡単にサンプルを書きました（計算が合っているのかは不明ですが）ので御確認ください。

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(328,330,21)
y = np.linspace(41,43,21)
# F = np.random.choice([0, 0.001, 0.01, 0.1],20*20)
F = np.random.uniform(0,1,20*20)

X, Y = np.meshgrid(x,y)
F = np.array(F).reshape(len(x)-1,-1)
m = -2.5 * np.log10(F/3630)

fig = plt.figure(figsize=(7.5, 7.5))
ax = fig.add_subplot(111)
c = ax.pcolormesh(X, Y, m, cmap='gray_r', vmin=0, vmax=20)
fig.colorbar(c, label='A')

plt.show()

Comments

[magichan]

いつの間にかサンプルデータが追加されていたのですね・・・
私のサンプルでは、そのデータを使っておりません。スミマセン

[8524ba23]

当方環境(matplotlib=2.2.2)ではinfはなぜか白くなっちゃいますね。。。

[magichan]

コメントありがとうございました。
たしかにそのようになってますね・・。

[mega_ne]

データを示しておらず申し訳ございません。
こちらも白くなってしまいました...

Answer 2

「フラックスFが元から0であった箇所は真っ黒で、かつ、等級を計算した結果0になったものも含めて等級の値が小さくなると白く、値が大きくなると黒に近づくというカラーバーにしたいです。」

magichanさんの回答のとおりでよいと思うのですが、なぜか当方環境ではinfが白で描画されるので回避策の提案です。（matplotlibの仕様？）

幸い、元の０は変換後に特殊な値infになるので、変換後のinfを黒を表す大きな値に置き換えればよいかと思います。
下記コード中ではVMAX=20を指定していますが、適宜調整ください。

import numpy as np
import matplotlib.pyplot as plt

import matplotlib
print(matplotlib.__version__) # 2.2.2

x = np.linspace(328,330,3)
y = np.linspace(41,43,3)
F = [0, 3630, 363, 36.3]

X, Y = np.meshgrid(x,y)
F = np.array(F).reshape(2,-1)

m = -2.5 * np.log10(F/3630)
print(m)
"""
[[ inf -0. ]
 [ 2.5  5. ]]
"""

# 元から０(変換後のinf)は最大値(黒)に
VMAX = 20
m[m == np.inf] = VMAX

fig = plt.figure()
ax = fig.add_subplot(111)
m = ax.pcolormesh(X, Y, m, cmap='gray_r', vmin=0, vmax=VMAX)
fig.colorbar(m, label='A')

plt.show()

Comments

[magichan]

「なぜか当方環境ではinfが白で描画される」
マジか・・・・。確認しましたところ確かにそうですね。全く気づいておりませんでした。

[mega_ne]

こちらの方法でうまくいきました。ありがとうございます。
もっと勉強します。

Answer 3

コードを書くには２つの段階があります。

(A) 解こうとする問題の定式化
(B) 定式化した問題を特定の機能・文法を持つ言語・ライブラリーを前提として計算機の言葉（コード）に翻訳すること

当然ながら(A)を行うにはその問題の専門領域の知識が必要ですが、問題が天文学であろうと何であろうとその定式化の中には必ずと言ってよいほど「数学の知識」が含まれていると思います。プログラマーは定式化の中の数学的部分についてはそれほど難しいとは感じないことが多いです。何故ならあらゆる問題に取り組む際に共通の道具として使う数学の応用に慣れているからです。もちろんその数学の難易度のレベルは様々であり義務教育を終えた人なら必ず教わっているはずのことから数学の専門家にしか理解が困難な高度なものまであるとは思いますが、多くの問題は義務教育～高校・大学の数学の中で解決できることがら(α)であると思えます。

少なくとも本件については(α)の範囲内に収まった数学の知識と計算機の特徴の知識が必要ではないでしょうか。

• ある範囲の値[a, b]を別の値[f(a), f(b)]へ写像する手段
  これは(α)の範囲内の知識だと思います。天文学とは直接関係ありません。

• 計算機のディスプレイで白色光は輝度(フラックスと考えてよいと思います)の量で表す

• 計算機のディスプレイでは輝度は256段階であるのが普通
  0～255の整数で輝度を表したり、輝度の分解能を意識する必要がないよう0.0～1.0までの浮動小数点数で表現するなど少しバリエーションはあります。

といったことを踏まえていればコードは作れると思います。

フラックスFの値に基づき、等級mを表す画素値をどう計算するかですが、色を画素値[0.0, 1.0]で表現する(ただし0.0は完全な黒, 1.0は完全な白)としますと

・Fが0なら
=>画素値は0.0(黒)
・Fが0以外なら
=>等級を計算し等級の範囲に応じて画素値へ写像

とすればよいです。画素値[0.0, 1.0]への写像の計算式をどうするかですが、それは等級を輝度へどう変換したいかの質問者さん自身の意思により決まることです。それ(定式化)が決まらないとコードに表現することはできません。

ググりますと太陽の等級は-26.7, 月は-12.7らしいので、等級を輝度へ線形に写像するなら

最小の等級	最大の等級	等級から画素値への写像
6	太陽	pixel = (m - 6) / (26.7 + 6)
10	月	pixel = (m - 10) / (12.7 + 10)

といった具合になるでしょう。写像の方法には目的によっていろいろ考えることができ、同じ線形写像でも

• (1) 前述のように特定のFの範囲を定め最小値を画素値0.0とし最大値を画素値1.0へ写像する
• (2) 輝度の最小値は自分が望む値にすることもできましょう。
  最小輝度を0.5, 最大を1.0にしたってよいのです。要は結果としての画像がどう見えてほしいかです。
• (3) ある写像を行い、画素値が範囲外になったら[0.0, 1.0]の範囲に丸めるという考え方もある
• (4) (1)のようにあらかじめ範囲を定めるのではなく、対象画像のFの中で非0のものの最小値を画素値0.0とし最大値を画素値1.0へ写像するといった考え方。

等々のバリエーションが考えられます。いずれにせよ質問者さんの課題は「ある範囲の値を別の値へ写像する」という数学的な応用を本件に適用できないでいることのように思います。

---

数学は多くの科学（プログラミングもその一種といって差し支えないと思います）の基礎(道具)です。プログラミング初心者の方の場合、大抵は義務教育～高校数学の範囲に有用なテクニックが含まれているにもかかわらずそれをうまく適用できないでいるという状況が目立つ気がします。

もしそうであれば数学を(単なる知識として知っていても)具体的な場面で応用・適用することに慣れていないだけだと思いますので、問題を数学的な面から捉える(数学的アプローチをする)練習をするのがよいのではないかと思います。

Comments

[KSwordOfHaste]

回答中「輝度(フラックスと考えてよいと思います)」は質問者さんを混乱させかねない表現かも知れません。ここでいうフラックスとは「ディスプレイを見る人間の感覚に与えるフラックス」のことで星のフラックスのことを言っているのではありません。

[mega_ne]

そうですね...
数学を全然使えていないと思います。
定式化や応用法を考えます。