Python で今日の月

Python で今日の月 (PyEphem版)

Python で月の出入り時刻を計算したり、満ち欠けの様子を表示することができます。

必要モジュールをインポートします。

import ephem
import datetime

観測者のインスタンス（ここでは osaka）を作成します。
作成した観測者 osaka に、現在日時（世界時）を設定ます。

osaka = ephem.Observer()
osaka.lat = '34.6914'
osaka.lon = '135.4917'
osaka.date = datetime.datetime.utcnow()

月のインスタンス（ここでは moon）を作成します。

moon = ephem.Moon()

osaka.next_rising(moon) とすれば、次の月の出の時刻が計算されます。
計算される時刻は世界時のため、日本標準時で表示するには、ephem.localtime メソッドを使用します。

print('次の月の出　 :', ephem.localtime(osaka.next_rising(moon)))
print('次の月の入り :', ephem.localtime(osaka.next_setting(moon)))

表示結果

次の月の出　 : 2019-09-09 15:29:13
次の月の入り : 2019-09-09 00:46:42

ephem.next_full_moon(osaka.date) とすれば、次の満月等の日時が計算されます。
月齢の計算は、前回の新月から現在までの時間を引き算して求めます。

print('次の新月 :', ephem.localtime(ephem.next_new_moon(osaka.date)))
print('次の上弦 :', ephem.localtime(ephem.next_first_quarter_moon(osaka.date)))
print('次の満月 :', ephem.localtime(ephem.next_full_moon(osaka.date)))
print('次の下弦 :', ephem.localtime(ephem.next_last_quarter_moon(osaka.date)))
print('現在の月齢 :', osaka.date - ephem.previous_new_moon(osaka.date))

表示結果

次の新月 : 2019-09-29 03:26:20.000006
次の上弦 : 2019-10-06 01:47:04.000006
次の満月 : 2019-09-14 13:32:44.000005
次の下弦 : 2019-09-22 11:40:52.000006
現在の月齢 : 9.182557198269933

月の満ち欠けの様子を表示するのは、python のグラフ描画ライブラリ Matplotlib を利用します。
ここでは三次元で描写した球面を平面に投射して、満ち欠けを表現します。

必要モジュールをインポートします。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

太陽と月の離角を計算します。離角が0度の時が新月、180度の時が満月です。
（正確には新月や満月の時でも、月と太陽はぴったり重ならないため、離角も0度や180度にはなりません。）
太陽と月の離角は、moon.elong メソッドで求められます。
rad2deg メソッドで、ラジアンから度の単位に変更しておきます。

moon.compute(osaka)
moon_elong = np.rad2deg(moon.elong)

プロットする領域を作り、3次元で表示するよう設定します。

fig = plt.figure(figsize=(5,5))
ax = fig.gca(projection='3d')

x, y, z軸の表示範囲を設定します。

ax.set_xlim([-1., 1.])
ax.set_ylim([-1., 1.])
ax.set_zlim([-1., 1.])

x, y, z軸や目盛は表示しないようにします。

for a in [ax.xaxis, ax.yaxis, ax.zaxis]:
    a.set_ticklabels([])
    a._axinfo['grid']['linewidth'] = 0
    a._axinfo['tick']['linewidth'] = 0

基準で表示される背景の x, y, z面も、何も表示しないよう、透明にします。

for a in [ax.w_xaxis, ax.w_yaxis, ax.w_zaxis]:
    a.line.set_linewidth(0)
    a.set_pane_color((0., 0., 0., 0.))

グラフ背面を灰色にします。

ax.set_facecolor('lightgray')

メッシュ状の球面 (u, v) を用意します。
u は接線方向の値として 0～2π を50分割、v は動経方向の値として 0～π を25分割して入れます。
各 (u, v) に対して、3次元空間の (x, y, z) 値を計算します。

u, v = np.mgrid[0:2*np.pi:50j, 0:np.pi:25j]
x = np.cos(u) * np.sin(v)
y = np.sin(u) * np.sin(v)
z = np.cos(v)

メッシュの球面に貼りつける色を用意します。
半分だけ黄色（R,G,B）=(1,1,0)、残りは黒にします。

colors = np.zeros((50, 25, 3))
for i in range(0, 25):
    for j in range(0, 25):
        colors[i][j][0] = 1
        colors[i][j][1] = 1
        colors[i][j][2] = 0

球面をプロットします。立体の影は表示しないようにします。

ax.plot_surface(x, y, z, facecolors = colors, shade = False)

グラフを見る方向を設定します。
z軸方向から見た迎角は0度、方位角は太陽と月の離角に応じた値にします。

ax.view_init(elev = 0, azim = moon_elong - 90)

実際に表示します。

plt.show()

表示結果

全体は次のようになります。

import ephem
import datetime
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 観測地設定
osaka = ephem.Observer()
osaka.lat = '34.6914'
osaka.lon = '135.4917'
osaka.date = datetime.datetime.utcnow()

# 月
moon = ephem.Moon()

# 次回　月の出・月の入り等を表示
print('次の月の出　 :', ephem.localtime(osaka.next_rising(moon)))
print('次の月の入り :', ephem.localtime(osaka.next_setting(moon)))

print('次の新月 :', ephem.localtime(ephem.next_new_moon(osaka.date)))
print('次の上弦 :', ephem.localtime(ephem.next_first_quarter_moon(osaka.date)))
print('次の満月 :', ephem.localtime(ephem.next_full_moon(osaka.date)))
print('次の下弦 :', ephem.localtime(ephem.next_last_quarter_moon(osaka.date)))
print('現在の月齢 :', osaka.date - ephem.previous_new_moon(osaka.date))

# 月と太陽の離角を計算
moon.compute(osaka)
moon_elong = np.rad2deg(moon.elong)

# 描画領域を準備
fig = plt.figure(figsize=(5,5))
ax = fig.gca(projection='3d')

# x, y, z軸の範囲設定
ax.set_xlim([-1., 1.])
ax.set_ylim([-1., 1.])
ax.set_zlim([-1., 1.])

# x, y, z軸や目盛を非表示に
for a in [ax.xaxis, ax.yaxis, ax.zaxis]:
    a.set_ticklabels([])
    a._axinfo['grid']['linewidth'] = 0
    a._axinfo['tick']['linewidth'] = 0

# 背景の x, y, z面を非表示に
for a in [ax.w_xaxis, ax.w_yaxis, ax.w_zaxis]:
    a.line.set_linewidth(0)
    a.set_pane_color((0., 0., 0., 0.))

# 背面を灰色に
ax.set_facecolor('lightgray')

# メッシュ状の球面 (u, v) を準備し、(x, y, z) 値を計算
u, v = np.mgrid[0:2*np.pi:50j, 0:np.pi:25j] # u:接線方向　v:動経方向
x = np.cos(u) * np.sin(v)
y = np.sin(u) * np.sin(v)
z = np.cos(v)

# メッシュの球面に貼りつける色を準備（半分だけ黄色に）
colors = np.zeros((50, 25, 3))
for i in range(0, 25):
    for j in range(0, 25):
        colors[i][j][0] = 1
        colors[i][j][1] = 1
        colors[i][j][2] = 0

# 球面をプロット
ax.plot_surface(x, y, z, facecolors = colors, shade = False)
# グラフを見る方向を設定
ax.view_init(elev = 0, azim = moon_elong - 90)

plt.show()

江越 航のホームページ

Python で今日の月 (PyEphem版)

　江越　航のホームページ