Python で今日の月

Python で今日の月 (Skyfield版)

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

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

from skyfield.api import load, Topos
from skyfield.almanac import find_discrete, risings_and_settings, moon_phases, MOON_PHASES
from datetime import timedelta
from pytz import timezone
import numpy as np

Timescale オブジェクトを読み込み、 t0 に現在時刻を設定します。日時は、世界時で指定します。
また、タイムゾーンを日本に設定します。

ts = load.timescale()
t0 = ts.now()
tz = timezone('Asia/Tokyo')

天体暦のファイルを読み込み、太陽、月、地球のオブジェクト（ここではsun, moon, earth）を作成します。

eph = load('de421.bsp')
sun, moon, earth = eph['sun'], eph['moon'], eph['earth']

観測者の位置を指定する Topos オブジェクトを設定します。

osaka = Topos('34.6914 N', '135.4917 E')

t1 に現在時刻 t0 の1日後の値を設定します。
find_discrete(t0, t1, risings_and_settings(eph, moon, osaka)) で、t0, t1 の間に起こる月の出・月の入りの時刻を計算します。
結果は、t、updown に配列形式で代入されます。
t には月の出または月の入りの時刻が、updown には月の出の場合 True、月の入りの場合 False が、順番に配列形式で代入されます。

t1 = ts.utc(t0.utc_datetime() + timedelta(days=1))
t, updown = find_discrete(t0, t1, risings_and_settings(eph, moon, osaka))

t、updown に代入された結果を、順番に表示します。

for yi, ti in zip(updown, t):
    print('次の月の出:' if yi else '次の月の入り:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')

表示結果

次の月の入り: 2020-01-28 20:23:33 JST
次の月の出: 2020-01-29 09:36:35 JST

次回の新月・上弦・満月・下弦の日時は次のようにして求めます。
t2 に現在時刻 t0 から30日後の値を設定します。
find_discrete(t0, t2, moon_phases(eph)) で、t0, t2 の間に起こる新月・満月等の時刻を計算します。
結果は、t、y に配列形式で代入されます。
t には新月等になる時刻が、y には位相に応じて0～3の値（新月：0、上弦：1、満月：2、下弦：3）が、順番に配列形式で代入されます。

t2 = ts.utc(t0.utc_datetime() + timedelta(days=30))
t, y = find_discrete(t0, t2, moon_phases(eph))

t、y に代入された結果を、順番に表示します。
MOON_PHASES は、各位相の名前の配列になっています。

for yi, ti in zip(y, t):
    print(MOON_PHASES[yi], ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')

表示結果

First Quarter 2020-02-02 10:41:40 JST
Full Moon 2020-02-09 16:33:16 JST
Last Quarter 2020-02-16 07:17:12 JST
New Moon 2020-02-24 00:32:00 JST

次に月齢を計算します。そのために、前回の新月の日時を求めます。
t3 に現在時刻 t0 から30日前の値を設定します。
find_discrete(t3, t0, moon_phases(eph)) で、t3, t0 の間に起こる新月・満月等の時刻を計算します。
結果は、t、y に配列形式で代入されます。
t には新月等になる時刻が、y には位相に応じて0～3の値（新月：0、上弦：1、満月：2、下弦：3）が、順番に配列形式で代入されます。

t3 = ts.utc(t0.utc_datetime() - timedelta(days=30))
t, y = find_discrete(t3, t0, moon_phases(eph))

配列 y の中で、値が 0（新月）になっている項目の順番（インデックス）を求めます。
np.where(y == 0)[0] で、インデックスの配列が得られます。
2回新月が含まれる可能性があるので、[-1]で、配列の中の後ろのインデックスの値を index に代入します。
これより前回の新月の日時を t_newmoon に代入します。

index = (np.where(y == 0)[0])[-1]
t_newmoon = t[index]

前回の新月の日時と、月齢を表示します。
月齢の計算は、前回の新月から現在までの時間を引き算して求めます。

print('前回の新月:', t_newmoon.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
print('月齢:{0:.1f}'.format(t0 - t_newmoon))

表示結果

前回の新月: 2020-01-25 06:41:59 JST
月齢:3.4

月の満ち欠けを表示するため、太陽と月の離角を計算します。離角が0度の時が新月、180度の時が満月です。
ここでは、太陽と月の黄経の差を計算します。

_, slon, _ = earth.at(t0).observe(sun).apparent().ecliptic_latlon()
_, mlon, _ = earth.at(t0).observe(moon).apparent().ecliptic_latlon()
moon_elong = (mlon.degrees - slon.degrees) % 360

print('離角:{0:.1f}'.format(moon_elong))

表示結果

離角:37.7

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

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

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

プロットする領域を作り、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()

表示結果

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

from skyfield.api import load, Topos
from skyfield.almanac import find_discrete, risings_and_settings, moon_phases, MOON_PHASES
from datetime import timedelta
from pytz import timezone
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 初期時刻設定
ts = load.timescale()
t0 = ts.now()
tz = timezone('Asia/Tokyo')

# 太陽・月・地球
eph = load('de421.bsp')
sun, moon, earth = eph['sun'], eph['moon'], eph['earth']

# 観測地設定
osaka = Topos('34.6914 N', '135.4917 E')

# 月の出・月の入り計算
t1 = ts.utc(t0.utc_datetime() + timedelta(days=1))
t, updown = find_discrete(t0, t1, risings_and_settings(eph, moon, osaka))

for yi, ti in zip(updown, t):
    print('次の月の出:' if yi else '次の月の入り:', ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')

# 新月・満月等計算
t2 = ts.utc(t0.utc_datetime() + timedelta(days=30))
t, y = find_discrete(t0, t2, moon_phases(eph))

for yi, ti in zip(y, t):
    print(MOON_PHASES[yi], ti.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')

# 前回の新月・月齢計算
t3 = ts.utc(t0.utc_datetime() - timedelta(days=30))
t, y = find_discrete(t3, t0, moon_phases(eph))

index = (np.where(y == 0)[0])[-1]
t_newmoon = t[index]
    
print('前回の新月:', t_newmoon.astimezone(tz).strftime('%Y-%m-%d %H:%M:%S'), 'JST')
print('月齢:{0:.1f}'.format(t0 - t_newmoon))

# 月と太陽の離角を計算
_, slon, _ = earth.at(t0).observe(sun).apparent().ecliptic_latlon()
_, mlon, _ = earth.at(t0).observe(moon).apparent().ecliptic_latlon()
moon_elong = (mlon.degrees - slon.degrees) % 360

print('離角:{0:.1f}'.format(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 で今日の月 (Skyfield版)

　江越　航のホームページ