この1ヶ月間、日没の方位角と時刻を写真撮影して観察してきた。
日没の位置がかなりのスピードで淡路島の妙見山付近から南に移動してきて、ユニバーサルシティに近づくほどにゆっくりとなり、ついに極値（最も南）に達したのが冬至の数日前。
しかし、日没時刻は12月のはじめより少しずつ遅くなり始めて日が長くなっている。
11月16日から12月05日の日没の方向 （前回の記事で掲載した写真）
12月09日から12月19日の日没の方向
Google Gemini AIに日出・日没の時刻と方位角を計算しグラフ化するPythonスクリプトを作るようお願いしたら、数秒で「モノになる」スクリプトを作ってよこした。
ソースコードを読んで、改善箇所をいくつか指摘して完成したスクリプトと出力結果を掲載する。
日出・日没の時刻と方位角（11月15日から1月15日）
日出・日没の時刻と方位角を描画するPythonスクリプト
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
================================================================================
日の出・日没・日中時間の総合計算・グラフ表示スクリプト
================================================================================
概要:
指定された緯度・経度、および期間に基づき、pyephemライブラリを使用して
日の出・日没・日中時間・方位角を計算し、リスト出力およびグラフ表示する。
実行環境:
Python 3.x
必要なライブラリ: pyephem, numpy, matplotlib
作成者: Google Gemini3
================================================================================
"""
from typing import List, Dict
import ephem
from datetime import datetime, timedelta
import numpy as np
import matplotlib.pyplot as plt
import sys
from typing import List, Tuple, Dict
# =================================================================
# グローバル設定変数 (ここを編集してください)
# =================================================================
TARGET_LATITUDE_DEG = 34.67
# TARGET_LATITUDE_DEG = 0.0
TARGET_LONGITUDE_DEG = 135.5
# TARGET_LONGITUDE_DEG = 135.0
START_YEAR = 2025
START_MONTH = 11
START_DAY = 15
NUM_DAYS = 60
# =================================================================
def calculate_sun_data(
lat: float, lon: float, start_y: int, start_m: int, start_d: int, days_count: int
) -> List[Dict]:
"""
指定された条件に基づき、太陽のデータを計算して辞書のリストで返す。
"""
obs = ephem.Observer()
obs.lat = str(lat)
obs.lon = str(lon)
results = []
jst_offset = timedelta(hours=9)
# 計算開始の起点（JST 00:00 = UTC 前日 15:00）
start_date = datetime(start_y, start_m, start_d)
for i in range(days_count):
current_date = start_date + timedelta(days=i)
# その日の開始時刻（JST 00:00）をセット
obs.date = current_date - jst_offset
sun = ephem.Sun() # type: ignore[reportAttributeAccessIssue]
try:
# --- 日の出 ---
rise_time_ephem = obs.next_rising(sun)
obs.date = rise_time_ephem
sun.compute(obs)
rise_az = np.degrees(sun.az)
rise_jst = rise_time_ephem.datetime() + jst_offset
# --- 日没 ---
# 日没は日の出の後に来るように設定
set_time_ephem = obs.next_setting(sun)
obs.date = set_time_ephem
sun.compute(obs)
set_az = np.degrees(sun.az)
set_jst = set_time_ephem.datetime() + jst_offset
# --- 日中時間 (日没 - 日の出) ---
duration_td = set_time_ephem.datetime() - rise_time_ephem.datetime()
total_seconds = int(duration_td.total_seconds())
d_h, rem = divmod(total_seconds, 3600)
d_m, d_s = divmod(rem, 60)
duration_str = f"{d_h:02d}:{d_m:02d}:{d_s:02d}"
duration_hours = total_seconds / 3600.0
results.append({
"date": current_date.strftime('%Y/%m/%d'),
"rise_t": rise_jst.strftime('%H:%M:%S'),
"rise_az": rise_az,
"set_t": set_jst.strftime('%H:%M:%S'),
"set_az": set_az,
"duration_t": duration_str,
"duration_h": duration_hours,
"rise_h_num": rise_jst.hour + rise_jst.minute / 60.0 + rise_jst.second / 3600.0,
"set_h_num": set_jst.hour + set_jst.minute / 60.0 + set_jst.second / 3600.0
})
except (ephem.AlwaysUpError, ephem.NeverUpError):
results.append({
"date": current_date.strftime('%Y/%m/%d'),
"rise_t": "--:--:--", "rise_az": np.nan,
"set_t": "--:--:--", "set_az": np.nan,
"duration_t": "--:--:--", "duration_h": np.nan,
"rise_h_num": np.nan, "set_h_num": np.nan
})
return results
def plot_data(lat: float, lon: float, data: List[Dict]) -> None:
"""
4段のグラフを表示する。方位角はデータ範囲に合わせて動的にレンジを調整。
"""
dates = [d["date"][5:] for d in data]
rise_times = [d["rise_h_num"] for d in data]
set_times = [d["set_h_num"] for d in data]
durations = [d["duration_h"] for d in data]
rise_azs = [d["rise_az"] for d in data]
set_azs = [d["set_az"] for d in data]
fig, axes = plt.subplots(4, 1, figsize=(12, 14), sharex=True)
fig.suptitle(f"太陽データ: 緯度 {lat}°, 経度 {lon}°", fontsize=16)
# 共通プロット設定: 点のみ、線なし
plot_style = {'marker': '.', 'markersize': 4, 'linestyle': 'None'}
# 余白設定 (%)
margin = 0.05
# 1. 日の出時刻
axes[0].plot(rise_times, color='orange', **plot_style)
axes[0].set_ylabel("日の出時刻 (JST)")
axes[0].grid(True, linestyle='--', alpha=0.7)
# 2. 日没時刻
axes[1].plot(set_times, color='darkblue', **plot_style)
axes[1].set_ylabel("日没時刻 (JST)")
axes[1].grid(True, linestyle='--', alpha=0.7)
# 3. 日中時間
axes[2].plot(durations, color='green', **plot_style)
axes[2].set_ylabel("日中時間 (時間)")
axes[2].grid(True, linestyle='--', alpha=0.7)
# 4. 方位角 (左軸:日の出, 右軸:日没)
# 左軸 (日の出)
ax4_left = axes[3]
ax4_left.plot(rise_azs, color='orange', label='日の出方位角', **plot_style)
ax4_left.set_ylabel("日の出方位角 (度/左軸)", color='orange')
ax4_left.tick_params(axis='y', labelcolor='orange')
# 日の出方位角の動的レンジ設定 (NaNを除外して計算)
if not np.all(np.isnan(rise_azs)):
r_min, r_max = np.nanmin(rise_azs), np.nanmax(rise_azs)
ax4_left.set_ylim(r_min - abs(r_max-r_min)*margin, r_max + abs(r_max-r_min)*margin)
ax4_left.grid(True, linestyle='--', alpha=0.7)
# 右軸 (日没)
ax4_right = ax4_left.twinx()
ax4_right.plot(set_azs, color='darkblue', label='日没方位角', **plot_style)
ax4_right.set_ylabel("日没方位角 (度/右軸)", color='darkblue')
ax4_right.tick_params(axis='y', labelcolor='darkblue')
# 日没方位角の動的レンジ設定
if not np.all(np.isnan(set_azs)):
s_min, s_max = np.nanmin(set_azs), np.nanmax(set_azs)
ax4_right.set_ylim(s_min - abs(s_max-s_min)*margin, s_max + abs(s_max-s_min)*margin)
# 凡例をまとめる
h1, l1 = ax4_left.get_legend_handles_labels()
h2, l2 = ax4_right.get_legend_handles_labels()
ax4_left.legend(h1 + h2, l1 + l2, loc='upper center', ncol=2, fontsize='small')
# フォーマッタ (時刻用)
from matplotlib.ticker import FuncFormatter
time_fmt = FuncFormatter(lambda y, p: f"{int(y):02d}:{int((y-int(y))*60):02d}")
axes[0].yaxis.set_major_formatter(time_fmt)
axes[1].yaxis.set_major_formatter(time_fmt)
# X軸設定
num_ticks = 12
step = max(1, len(dates) // num_ticks)
axes[3].set_xticks(np.arange(0, len(dates), step))
axes[3].set_xticklabels([dates[i] for i in range(0, len(dates), step)], rotation=45)
axes[3].set_xlabel("月/日")
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()
# --- メイン処理 ---
if __name__ == "__main__":
print(f"--- 太陽データリスト (開始: {START_YEAR}/{START_MONTH}/{START_DAY}, {NUM_DAYS}日間) ---")
print("年月日,日の出時刻,日の出方角,日没時刻,日没方角(度),日中時間")
sun_results = calculate_sun_data(
TARGET_LATITUDE_DEG, TARGET_LONGITUDE_DEG,
START_YEAR, START_MONTH, START_DAY, NUM_DAYS
)
for r in sun_results:
r_az = f"{r['rise_az']:.1f}" if not np.isnan(r['rise_az']) else "N/A"
s_az = f"{r['set_az']:.1f}" if not np.isnan(r['set_az']) else "N/A"
print(f"{r['date']},{r['rise_t']},{r_az},{r['set_t']},{s_az},{r['duration_t']}")
print("----------------------------------------------------------------")
ans = input("グラフを表示しますか？ (y/n): ")
if ans.lower() == 'y':
plot_data(TARGET_LATITUDE_DEG, TARGET_LONGITUDE_DEG, sun_results)
