"""
教育用再現コード: 2023年度 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================
論文タイトル：兵庫県の宝「いかなご」を守る
              ——兵庫県の小魚「いかなご」の漁獲量減少の要因分析

【分析概要】
  対象論文の分析テーマ：
    瀬戸内海（兵庫県沿岸）でのいかなご漁獲量の急激な減少要因を
    気候変動・環境要因から統計的に分析。
    ・気温上昇（海水温上昇の代理変数）
    ・降水量（沿岸域への淡水流入・栄養塩供給）
    ・漁業コミュニティの高齢化
    ・消費支出・食料費（需要側の変化）

  本教育用コードのアプローチ：
    SSDSE-B（都道府県パネルデータ, 2012〜2023年）を用いて、
    以下の4図を生成し、論文の分析手法（時系列分析・散布図・
    相関分析・地域比較）を教育的に再現する。

  ▶ 図1: 兵庫県の気温（年平均・最高）の時系列 + 食料費の重ね合わせ
  ▶ 図2: 全47都道府県での年平均気温 vs 食料費割合の散布図（兵庫強調）
  ▶ 図3: 食料費/消費支出比率の時系列（兵庫県 vs 全国平均, 2012〜2023）
  ▶ 図4: 沿岸主要都道府県の降水量棒グラフ（最新年度）

【変数（SSDSE-B より）】
  B4101相当: 年平均気温           — 海水温上昇の代理変数
  B4102相当: 最高気温             — 極端高温イベント
  B4109相当: 降水量（年間）        — 沿岸環境への淡水・栄養塩流入
  高齢化率  : 65歳以上人口/総人口  — 漁業コミュニティの担い手不足
  L3221相当: 消費支出             — 需要側の変化
  L322101相当: 食料費             — 水産物消費の代理変数

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計センターより公表の実データ（np.random 一切不使用）

【データサイエンス学習ポイント】
  1. 時系列データのパネル構造（都道府県×年度）の理解
  2. 環境要因の代理変数（proxy variable）の考え方
  3. 散布図による地域間比較と外れ値の視覚化
  4. 沿岸・水産業に関連する都道府県の特徴抽出
=================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2023_H5_6_shorei.py
# ============================================================


import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats
import warnings
import os

warnings.filterwarnings('ignore')

# ──────────────────────────────────────────────────────────────
# 共通設定
# ──────────────────────────────────────────────────────────────
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ================================================================
# ■ 実データ読み込み（SSDSE-B: 都道府県パネルデータ 2012〜2023）
# ================================================================
df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)

# 都道府県行のみ抽出（地域コード: R01000〜R47000）
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

# 数値変換する列（文字列列は除く）
NUM_COLS = [
    '総人口', '65歳以上人口',
    '年平均気温',
    '最高気温（日最高気温の月平均の最高値）',
    '降水量（年間）',
    '消費支出（二人以上の世帯）',
    '食料費（二人以上の世帯）',
    '月間有効求人数（一般）',
    '就職件数（一般）',
]

for col in NUM_COLS:
    if col in df_b.columns:
        df_b[col] = pd.to_numeric(df_b[col], errors='coerce')

df_b['年度'] = pd.to_numeric(df_b['年度'], errors='coerce')

# ── 派生変数の計算 ──────────────────────────────────────────────
# 高齢化率（%）
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100

# 食料費/消費支出 比率（%）— 食への支出割合
df_b['食料費割合'] = df_b['食料費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100

# 労働市場充足率（求人充足度のプロキシ）— 沿岸漁業町の雇用圧力
df_b['求人充足率'] = df_b['就職件数（一般）'] / df_b['月間有効求人数（一般）'] * 100

# ── 兵庫県データの抽出（2012〜2023）────────────────────────────
hyogo = df_b[df_b['都道府県'] == '兵庫県'].copy()
hyogo = hyogo.sort_values('年度').reset_index(drop=True)

print("=" * 60)
print("■ SSDSE-B 読み込み完了")
print(f"  全都道府県×年度行数: {len(df_b)}")
print(f"  兵庫県の年度数: {len(hyogo)}年分（{hyogo['年度'].min()}〜{hyogo['年度'].max()}）")
print("=" * 60)

print("\n兵庫県 主要変数の時系列サマリー:")
print(hyogo[['年度', '年平均気温', '最高気温（日最高気温の月平均の最高値）',
             '降水量（年間）', '高齢化率', '食料費割合']].to_string(index=False))

# ── 沿岸・水産業主要都道府県の定義 ────────────────────────────
COASTAL_PREFS = ['北海道', '青森県', '岩手県', '宮城県', '秋田県', '新潟県',
                 '富山県', '石川県', '福井県', '千葉県', '神奈川県', '静岡県',
                 '愛知県', '三重県', '兵庫県', '和歌山県', '鳥取県', '島根県',
                 '岡山県', '広島県', '山口県', '徳島県', '香川県', '愛媛県',
                 '高知県', '福岡県', '佐賀県', '長崎県', '熊本県', '大分県',
                 '宮崎県', '鹿児島県', '沖縄県']

# ── 全国平均（都道府県平均）の計算 ───────────────────────────
nat_avg = df_b.groupby('年度')[['食料費割合', '年平均気温']].mean().reset_index()
nat_avg = nat_avg.sort_values('年度')

# ── 最新年度（2023年）の全都道府県スナップショット ──────────────
df_2023 = df_b[df_b['年度'] == 2023].copy().dropna(subset=['年平均気温', '食料費割合'])

# 沿岸都道府県の2023年スナップ
df_coastal_2023 = df_2023[df_2023['都道府県'].isin(COASTAL_PREFS)].copy()
df_coastal_2023 = df_coastal_2023.dropna(subset=['降水量（年間）'])
df_coastal_2023 = df_coastal_2023.sort_values('降水量（年間）', ascending=False).head(10)

print(f"\n■ 2023年 降水量上位10沿岸都道府県:")
print(df_coastal_2023[['都道府県', '降水量（年間）', '年平均気温', '食料費割合']].to_string(index=False))

# ================================================================
# ■ 図1: 兵庫県の気温時系列 + 食料費の重ね合わせ（2012〜2023）
# ================================================================
fig1, ax1a = plt.subplots(figsize=(10, 5))

years = hyogo['年度'].values
temp_avg = hyogo['年平均気温'].values
temp_max = hyogo['最高気温（日最高気温の月平均の最高値）'].values
food_ratio = hyogo['食料費割合'].values

color_temp = '#E65100'
color_max  = '#C62828'
color_food = '#1565C0'

# 左軸: 気温
ax1a.plot(years, temp_avg, 'o-', color=color_temp, linewidth=2.0,
          markersize=6, label='年平均気温 (°C)')
ax1a.plot(years, temp_max, 's--', color=color_max, linewidth=1.5,
          markersize=5, alpha=0.7, label='最高気温 (°C)')
ax1a.fill_between(years, temp_avg, alpha=0.10, color=color_temp)

# トレンド線（年平均気温）
z1 = stats.linregress(years, temp_avg)
ax1a.plot(years, z1.intercept + z1.slope * years, '-',
          color=color_temp, linewidth=1.0, alpha=0.5, linestyle=':')

ax1a.set_xlabel('年度', fontsize=11)
ax1a.set_ylabel('気温（°C）', fontsize=11, color=color_temp)
ax1a.tick_params(axis='y', labelcolor=color_temp)
ax1a.set_xticks(years)
ax1a.set_xticklabels([str(y) for y in years], rotation=45, fontsize=9)
ax1a.grid(True, alpha=0.25)

# 右軸: 食料費割合
ax1b = ax1a.twinx()
ax1b.plot(years, food_ratio, 'D-', color=color_food, linewidth=2.0,
          markersize=5, label='食料費/消費支出 (%)', alpha=0.8)
ax1b.set_ylabel('食料費/消費支出 (%)', fontsize=11, color=color_food)
ax1b.tick_params(axis='y', labelcolor=color_food)

# 凡例を統合
lines1, labels1 = ax1a.get_legend_handles_labels()
lines2, labels2 = ax1b.get_legend_handles_labels()
ax1a.legend(lines1 + lines2, labels1 + labels2,
            loc='upper left', fontsize=9, framealpha=0.85)

r1, p1 = stats.pearsonr(temp_avg, food_ratio)
ax1a.set_title(
    f'兵庫県の気温・食料費の推移（2012〜2023年）\n'
    f'気温×食料費割合 相関係数 r={r1:.3f}（p={p1:.3f}）',
    fontsize=12, fontweight='bold'
)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2023_H5_6_fig1_timeseries.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2023_H5_6_fig1_timeseries.png")

# ================================================================
# ■ 図2: 散布図 — 年平均気温 vs 食料費割合（全47都道府県, 2023年）
#         兵庫県をハイライト
# ================================================================
fig2, ax2 = plt.subplots(figsize=(10, 6))

# 全都道府県プロット
mask_not_hyogo = df_2023['都道府県'] != '兵庫県'
ax2.scatter(
    df_2023.loc[mask_not_hyogo, '年平均気温'],
    df_2023.loc[mask_not_hyogo, '食料費割合'],
    color='#78909C', s=45, alpha=0.65, label='他都道府県', zorder=2
)

# 兵庫県ハイライト
hyogo_2023 = df_2023[df_2023['都道府県'] == '兵庫県']
ax2.scatter(
    hyogo_2023['年平均気温'],
    hyogo_2023['食料費割合'],
    color='#C62828', s=180, zorder=5, marker='*', label='兵庫県'
)
if not hyogo_2023.empty:
    hx = float(hyogo_2023['年平均気温'].values[0])
    hy = float(hyogo_2023['食料費割合'].values[0])
    ax2.annotate(
        '兵庫県',
        xy=(hx, hy), xytext=(hx + 0.5, hy + 0.3),
        fontsize=10, color='#C62828', fontweight='bold',
        arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.2)
    )

# 回帰線
x2 = df_2023['年平均気温'].dropna()
y2 = df_2023['食料費割合'].dropna()
common2 = df_2023[['年平均気温', '食料費割合']].dropna()
x2v = common2['年平均気温'].values
y2v = common2['食料費割合'].values
z2 = stats.linregress(x2v, y2v)
xs2 = [x2v.min(), x2v.max()]
ax2.plot(xs2, [z2.intercept + z2.slope * x for x in xs2],
         'r-', linewidth=1.5, alpha=0.55, label=f'回帰直線 (r={z2.rvalue:.3f})')

r2, p2 = stats.pearsonr(x2v, y2v)

ax2.set_xlabel('年平均気温（°C）', fontsize=12)
ax2.set_ylabel('食料費/消費支出 (%)', fontsize=12)
ax2.set_title(
    f'年平均気温 vs 食料費割合（2023年, 47都道府県）\n'
    f'r={r2:.3f}（p={p2:.3f}）　気温が高い地域ほど食料費割合の傾向を確認',
    fontsize=12, fontweight='bold'
)
ax2.legend(fontsize=9, framealpha=0.85)
ax2.grid(True, alpha=0.25)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2023_H5_6_fig2_scatter.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2023_H5_6_fig2_scatter.png")

# ================================================================
# ■ 図3: 食料費/消費支出 比率の時系列
#         兵庫県 vs 全国平均（2012〜2023）
# ================================================================
fig3, ax3 = plt.subplots(figsize=(10, 5))

# 兵庫県
ax3.plot(years, food_ratio, 'o-', color='#C62828', linewidth=2.5,
         markersize=7, label='兵庫県', zorder=3)

# 全国平均
nat_yr = nat_avg['年度'].values
nat_fr = nat_avg['食料費割合'].values
ax3.plot(nat_yr, nat_fr, 's--', color='#1565C0', linewidth=2.0,
         markersize=5, alpha=0.8, label='全国平均（47都道府県平均）', zorder=2)

# 差分の塗りつぶし
common_yrs = sorted(set(years) & set(nat_yr))
hy_dict = dict(zip(years, food_ratio))
na_dict = dict(zip(nat_yr, nat_fr))
c_y = [hy_dict[y] for y in common_yrs]
c_n = [na_dict[y] for y in common_yrs]
ax3.fill_between(common_yrs, c_y, c_n, alpha=0.08, color='#C62828',
                 label='兵庫と全国の差')

# データラベル（兵庫県）
for yr, fr in zip(years, food_ratio):
    if yr % 2 == 0:
        ax3.annotate(f'{fr:.1f}%', xy=(yr, fr), xytext=(0, 8),
                     textcoords='offset points', ha='center', fontsize=8,
                     color='#C62828')

ax3.set_xlabel('年度', fontsize=11)
ax3.set_ylabel('食料費/消費支出 (%)', fontsize=11)
ax3.set_xticks(common_yrs)
ax3.set_xticklabels([str(y) for y in common_yrs], rotation=45, fontsize=9)
ax3.set_title(
    '食料費/消費支出 比率の推移（2012〜2023年）\n兵庫県 vs 全国平均',
    fontsize=12, fontweight='bold'
)
ax3.legend(fontsize=9, framealpha=0.85)
ax3.grid(True, alpha=0.25)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2023_H5_6_fig3_foodratio.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2023_H5_6_fig3_foodratio.png")

# ================================================================
# ■ 図4: 沿岸主要都道府県の降水量比較（2023年, 上位10道府県）
# ================================================================
fig4, ax4 = plt.subplots(figsize=(11, 6))

pref_names = df_coastal_2023['都道府県'].values
rain_vals  = df_coastal_2023['降水量（年間）'].values
colors4 = ['#C62828' if p == '兵庫県' else '#1565C0' for p in pref_names]

bars = ax4.bar(range(len(pref_names)), rain_vals,
               color=colors4, alpha=0.78, edgecolor='white', width=0.65)

# 数値ラベル
for i, (bar, val) in enumerate(zip(bars, rain_vals)):
    ax4.text(bar.get_x() + bar.get_width() / 2, val + 20,
             f'{val:.0f} mm', ha='center', va='bottom', fontsize=9)

ax4.set_xticks(range(len(pref_names)))
ax4.set_xticklabels(pref_names, rotation=30, ha='right', fontsize=10)
ax4.set_ylabel('年間降水量（mm）', fontsize=11)
ax4.set_title(
    '沿岸主要都道府県の年間降水量（2023年, 上位10道府県）\n'
    '降水量は沿岸生態系への栄養塩供給・海水塩分濃度に影響',
    fontsize=12, fontweight='bold'
)
ax4.grid(axis='y', alpha=0.25)

from matplotlib.patches import Patch
ax4.legend(handles=[
    Patch(color='#C62828', alpha=0.78, label='兵庫県'),
    Patch(color='#1565C0', alpha=0.78, label='その他沿岸都道府県'),
], fontsize=9, loc='upper right')

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2023_H5_6_fig4_rainfall.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2023_H5_6_fig4_rainfall.png")

# ================================================================
# ■ 相関サマリー（コンソール出力）
# ================================================================
print("\n" + "=" * 60)
print("■ 兵庫県 環境要因の相関分析サマリー（2012〜2023年, n=12）")
print("=" * 60)
hyogo_clean = hyogo.dropna(subset=['年平均気温', '降水量（年間）',
                                    '高齢化率', '食料費割合'])

pairs = [
    ('年平均気温', '食料費割合', '気温 × 食料費割合'),
    ('年度',       '年平均気温', '年度 × 年平均気温（温暖化トレンド）'),
    ('年度',       '高齢化率',   '年度 × 高齢化率（高齢化進行）'),
    ('年度',       '食料費割合', '年度 × 食料費割合'),
    ('降水量（年間）', '食料費割合', '降水量 × 食料費割合'),
]
print(f"\n  {'分析':<35} {'r':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 63)
for x_col, y_col, label in pairs:
    tmp = hyogo_clean[[x_col, y_col]].dropna()
    if len(tmp) >= 4:
        r, p = stats.pearsonr(tmp[x_col].values, tmp[y_col].values)
        sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
        print(f"  {label:<35} {r:>8.4f} {p:>10.4f} {sig:>6}")

print("\n全図の生成完了（4枚）")
print("  2023_H5_6_fig1_timeseries.png : 気温・食料費の時系列（兵庫県）")
print("  2023_H5_6_fig2_scatter.png    : 気温vs食料費割合 散布図（全47都道府県）")
print("  2023_H5_6_fig3_foodratio.png  : 食料費割合 時系列（兵庫 vs 全国）")
print("  2023_H5_6_fig4_rainfall.png   : 降水量比較（沿岸上位10道府県）")