海水温からの降水量予測を目指して

# 原論文 表4：集中豪雨（1時間降水量）の年間回数と海水温の相関係数（1976–2022, 全国アメダス）
# 行 = 豪雨のしきい値, 列 = 季節。すべて «原論文の報告値» （新数値ではない）。
table4 = {
    '50mm以上':  {'年平均': 0.598, '冬': 0.525, '春': 0.365, '夏': 0.611, '秋': 0.493},
    '80mm以上':  {'年平均': 0.491, '冬': 0.451, '春': 0.319, '夏': 0.479, '秋': 0.390},
    '100mm以上': {'年平均': 0.383, '冬': 0.336, '春': 0.295, '夏': 0.440, '秋': 0.208},
}
seasons = ['年平均', '冬', '春', '夏', '秋']
print('\n【原論文 表4 の報告値】集中豪雨回数 × 海水温の相関係数（再計算ではない）')
for k, row in table4.items():
    print('  ' + k.ljust(9) + ' ' + '  '.join(f'{s}={row[s]:+.3f}' for s in seasons))
print(f'  → 年平均: 50mm(r={table4["50mm以上"]["年平均"]}) > 100mm(r={table4["100mm以上"]["年平均"]})')
print('    降水量が少ない豪雨ほど海水温との相関が強い、と原論文は報告。')

fig, ax = plt.subplots(figsize=(9.2, 5.2))
x = np.arange(len(seasons))
w = 0.26
colors = {'50mm以上': '#1565C0', '80mm以上': '#42A5F5', '100mm以上': '#90CAF9'}
for i, (k, row) in enumerate(table4.items()):
    ax.bar(x + (i - 1) * w, [row[s] for s in seasons], w, label=k, color=colors[k])
ax.axhline(0.2, color='gray', ls='--', lw=1)
ax.text(4.35, 0.205, 'r=0.2（やや弱い相関の目安）', fontsize=8, color='gray', ha='right', va='bottom')
ax.set_xticks(x)
ax.set_xticklabels(seasons)
ax.set_ylabel('相関係数 r（原論文の報告値）')
ax.set_title('図1【原論文 表4 の報告値・再計算ではない】\n集中豪雨の年間回数 × 海水温の相関（海水温は SSDSE 未収録）',
             fontsize=11)
ax.legend(fontsize=9)
ax.grid(axis='y', alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2023_H5_4_fig1_scatter2022.png'), bbox_inches='tight')
plt.close(fig)
print('  → 図1 保存: 2023_H5_4_fig1_scatter2022.png（報告値の可視化）')

df_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    header=0,            # 行0がコード名（B4101 等）
    encoding='cp932',    # 日本語Windows由来の文字コード
)

# 都道府県行だけを抽出（地域コード Code が R##### 形式）
mask_pref = df_raw['Code'].astype(str).str.match(r'^R\d{5}$', na=False)
df_b = df_raw[mask_pref].copy()
df_b = df_b.rename(columns={'SSDSE-B-2026': '年度', 'Prefecture': '都道府県'})

# 降水量（B4109）・年平均気温（B4101）・年度を数値化
for c in ['B4101', 'B4109']:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')
df_b['年度'] = pd.to_numeric(df_b['年度'], errors='coerce')

years = sorted(df_b['年度'].dropna().unique().astype(int))
print('=' * 62)
print('■ SSDSE-B-2026 読み込み完了（都道府県データ）')
print(f'  都道府県 × 年度の総行数: {len(df_b)}')
print(f'  収録年度: {years[0]}〜{years[-1]}（{len(years)}年分）')
print(f'  2023年度の都道府県数: {(df_b["年度"] == 2023).sum()}')
print('  ※ 海水温は SSDSE 未収録。降水量 B4109 のみ実再現に使用。')
print('=' * 62)

# 原論文の観測地点 → 対応する都道府県（SSDSE-B は都道府県単位）
JPSEA_PREFS = ['秋田県', '新潟県', '富山県', '石川県', '福井県', '兵庫県']  # 豊岡=兵庫
PACIFIC_PREFS = ['千葉県', '静岡県', '愛知県', '和歌山県', '高知県', '宮崎県']

def annual_mean_precip(prefs):
    sub = df_b[df_b['都道府県'].isin(prefs)]
    return sub.groupby('年度')['B4109'].mean()

jpsea = annual_mean_precip(JPSEA_PREFS)     # 日本海側6県の年降水量の平均
pacif = annual_mean_precip(PACIFIC_PREFS)   # 太平洋側6県の年降水量の平均

print('\n【SSDSE-B 実再現】日本海側/太平洋側の年降水量（B4109, mm）')
print(f'  日本海側6県  12年平均 = {jpsea.mean():.0f} mm / 変動係数CV = {jpsea.std()/jpsea.mean():.3f}')
print(f'  太平洋側6県  12年平均 = {pacif.mean():.0f} mm / 変動係数CV = {pacif.std()/pacif.mean():.3f}')

fig, ax = plt.subplots(figsize=(9.2, 5.2))
ax.plot(jpsea.index, jpsea.values, 'o-', color=C_JPSEA, lw=2.2,
        label='日本海側6県 平均（秋田・新潟・富山・石川・福井・兵庫）')
ax.plot(pacif.index, pacif.values, 's-', color=C_PACIF, lw=2.2,
        label='太平洋側6県 平均（千葉・静岡・愛知・和歌山・高知・宮崎）')
ax.set_xlabel('年度')
ax.set_ylabel('年降水量 B4109（mm）')
ax.set_title('図2【SSDSE-B 実再現】都道府県別・年降水量の推移（2012–2023）\n'
             '※海水温は SSDSE 未収録のため本図に含まない／原論文は市・季節別3か月降水量を使用（集計単位が異なる）',
             fontsize=11)
ax.legend(fontsize=9, loc='upper right')
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2023_H5_4_fig2_timeseries.png'), bbox_inches='tight')
plt.close(fig)
print('  → 図2 保存: 2023_H5_4_fig2_timeseries.png（実再現）')

sites = ['油津', '清水', '室戸岬', '潮岬', '伊良湖', '御前岬', '石廊岬', '館山']
# 原論文 表8（太平洋側の夏・四国沖）／表9（太平洋側の冬・四国沖）——報告値
summer_shikoku = [-0.147, -0.060, -0.289, -0.150, -0.251, -0.349, -0.379, -0.274]  # 表8 四国沖
winter_shikoku = [0.360, 0.285, 0.289, 0.303, 0.322, 0.190, 0.228, 0.077]          # 表9 四国沖
print('\n【原論文 表8・表9 の報告値】太平洋側 海水温×降水量の相関（四国沖・再計算ではない）')
print('  地点     夏(表8)   冬(表9)')
for s, a, b in zip(sites, summer_shikoku, winter_shikoku):
    print(f'  {s.ljust(4)} {a:+.3f}   {b:+.3f}')
print(f'  夏の平均 r={np.mean(summer_shikoku):+.3f}（全地点で負）／'
      f'冬の平均 r={np.mean(winter_shikoku):+.3f}（全地点で正・図4の四国沖 r=0.257）')

fig, ax = plt.subplots(figsize=(9.4, 5.2))
x = np.arange(len(sites))
w = 0.38
ax.bar(x - w / 2, summer_shikoku, w, label='夏（表8・四国沖）', color=C_NEG, alpha=0.85)
ax.bar(x + w / 2, winter_shikoku, w, label='冬（表9・四国沖）', color=C_POS, alpha=0.85)
ax.axhline(0, color='black', lw=0.8)
ax.set_xticks(x)
ax.set_xticklabels(sites, fontsize=9)
ax.set_ylabel('相関係数 r（原論文の報告値）')
ax.set_title('図3【原論文 表8・表9 の報告値・再計算ではない】\n'
             '太平洋側：海水温(四国沖)×降水量 は 夏=負・冬=正 に反転（海水温は SSDSE 未収録）',
             fontsize=11)
ax.legend(fontsize=9)
ax.grid(axis='y', alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2023_H5_4_fig3_panel_scatter.png'), bbox_inches='tight')
plt.close(fig)
print('  → 図3 保存: 2023_H5_4_fig3_panel_scatter.png（報告値の可視化）')

# 原論文 表10：8地点の平均相関（月別）——報告値
months = ['12月', '1月', '2月', '3月']
month_avg = [-0.010, 0.263, 0.193, 0.011]   # 表10「平均」列
print('\n【原論文 表10 の報告値】太平洋側の冬・月別（8地点平均・再計算ではない）')
for m, r in zip(months, month_avg):
    print(f'  {m}: r={r:+.3f}')
print('  → 1月(r=0.263)・2月(r=0.193)で正の相関。12月・3月はほぼ0。')
print('    原論文は「1〜2月に発達する南岸低気圧の影響」と解釈。')

fig, ax = plt.subplots(figsize=(8.2, 5.0))
bar_colors = [C_POS if r >= 0.2 else '#9E9E9E' for r in month_avg]
bars = ax.bar(months, month_avg, color=bar_colors, width=0.55)
for b, r in zip(bars, month_avg):
    ax.text(b.get_x() + b.get_width() / 2, r + (0.008 if r >= 0 else -0.02),
            f'{r:+.3f}', ha='center', va='bottom' if r >= 0 else 'top', fontsize=10)
ax.axhline(0, color='black', lw=0.8)
ax.axhline(0.2, color='gray', ls='--', lw=1)
ax.set_ylabel('相関係数 r（原論文の報告値・8地点平均）')
ax.set_ylim(-0.1, 0.35)
ax.set_title('図4【原論文 表10 の報告値・再計算ではない】\n'
             '太平洋側の冬・月別：1〜2月に海水温×降水量の正の相関（海水温は SSDSE 未収録）',
             fontsize=11)
ax.grid(axis='y', alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2023_H5_4_fig4_beta_ts.png'), bbox_inches='tight')
plt.close(fig)
print('  → 図4 保存: 2023_H5_4_fig4_beta_ts.png（報告値の可視化）')

print('\n' + '=' * 62)
print('全図の生成完了（4枚）')
print('  図1: 集中豪雨×海水温の相関（表4）        … 報告値の可視化')
print('  図2: 日本海側/太平洋側の年降水量 B4109    … SSDSE-B 実再現')
print('  図3: 太平洋側 夏 vs 冬（表8・表9）          … 報告値の可視化')
print('  図4: 太平洋側の冬・月別（表10）            … 報告値の可視化')
print('  ※ 海水温絡み（図1・図3・図4）は原論文の報告値であり再計算ではない。')
print('  ※ 新たな数値の捏造は一切していない。')
print('=' * 62)