CO₂排出特性による市区町村の類型化と
地域特性の関係

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.gridspec import GridSpec
import warnings
warnings.filterwarnings('ignore')

# ── パス設定 ─────────────────────────────────────────────────────────────────
FIG_DIR = 'html/figures'
DATA_DIR = 'data/raw'
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family':         'Hiragino Sans',
    'axes.unicode_minus':  False,
    'figure.dpi':          150,
    'axes.spines.top':     False,
    'axes.spines.right':   False,
})

print("=" * 65)
print("■ SSDSE-E-2026 実データ読み込み")
print("=" * 65)

df_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'), encoding='cp932')
df_e = df_e_raw.iloc[2:].copy()
df_e.columns = df_e_raw.iloc[1].values
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)

# 数値変換
num_cols_e = [
    '総人口', 'ごみ総排出量（総量）', '1人1日当たりの排出量', 'ごみのリサイクル率',
    '農家数（販売農家）', '1人当たり県民所得（平成27年基準）',
    '総面積（北方地域及び竹島を除く）', '大型小売店数', '飲食店数', '医師数',
    '事業所数（民営）（製造業）', '事業所数（民営）（建設業）',
]
for c in num_cols_e:
    df_e[c] = pd.to_numeric(df_e[c], errors='coerce')

print(f"  読み込み完了: {len(df_e)} 都道府県")
print(f"  主要変数確認:")
print(f"    ごみ総排出量（総量）[t]: {df_e['ごみ総排出量（総量）'].describe().loc[['mean','min','max']].round(0).to_dict()}")
print(f"    ごみのリサイクル率[%]:  {df_e['ごみのリサイクル率'].describe().loc[['mean','min','max']].round(1).to_dict()}")

# ── 指標の計算 ─────────────────────────────────────────────────────────────

# 農家率（農家数 / 総人口 × 1000）
df_e['農家率'] = df_e['農家数（販売農家）'] / df_e['総人口'] * 1000

# 人口密度 (人/km²) — 面積はha単位なので÷100でkm²
df_e['人口密度'] = df_e['総人口'] / (df_e['総面積（北方地域及び竹島を除く）'] / 100)

# 商業活動指数（飲食店数 + 大型小売店数 の人口千対）
df_e['商業活動指数'] = (df_e['飲食店数'] + df_e['大型小売店数']) / df_e['総人口'] * 1000

# 医師数10万対
df_e['医師数10万対'] = df_e['医師数'] / df_e['総人口'] * 100000

# 製造業事業所率（製造業事業所 / 全事業所）
df_e['製造業率'] = df_e['事業所数（民営）（製造業）'] / (df_e['事業所数（民営）（製造業）'] + df_e['事業所数（民営）（建設業）'] + 1)

# ── ウィーバー法用の4部門シェア計算 ───────────────────────────────────────
# 廃棄物・環境の視点から4つの活動軸を定義:
#   1. ごみ排出（ごみ総排出量 → 都道府県シェア）
#   2. リサイクル（ごみのリサイクル率）
#   3. 農業活動（農家率）
#   4. 商業活動（商業活動指数）

# 各軸を全都道府県合計に対する%シェアに変換
def to_share_pct(series):
    """シリーズを全都道府県の合計に対する%シェアに変換"""
    total = series.sum()
    return (series / total * 100).values

# 4部門シェア（合計は各部門の合計シェア=100になるよう正規化）
waste_share    = to_share_pct(df_e['ごみ総排出量（総量）'])
recycle_share  = to_share_pct(df_e['ごみのリサイクル率'])
agri_share     = to_share_pct(df_e['農家率'])
commerce_share = to_share_pct(df_e['商業活動指数'])

# 各都道府県の4部門比率（行合計が100%）
shares_matrix = np.column_stack([
    waste_share, recycle_share, agri_share, commerce_share
])
# 各都道府県内で行合計を100に正規化
row_sums = shares_matrix.sum(axis=1, keepdims=True)
shares_pct = shares_matrix / row_sums * 100  # shape: (47, 4)

SECTORS = ['廃棄物排出', 'リサイクル', '農業活動', '商業活動']

def modified_weaver(shares_arr):
    """
    shares_arr: 1次元array、各部門の百分率シェア（合計100）
    戻り値   : (k_opt, dominant_idx)
      k_opt       最適なk（複合類型の次元数）
      dominant_idx 最大シェア部門のインデックス
    """
    n = len(shares_arr)
    sorted_idx = np.argsort(shares_arr)[::-1]   # 降順インデックス
    sorted_sh  = shares_arr[sorted_idx]

    best_k = 1
    best_w = np.inf
    for k in range(1, n + 1):
        top_k = sorted_sh[:k]
        equal = 100.0 / k
        w_k   = np.sum((top_k - equal) ** 2)
        if w_k < best_w:
            best_w = w_k
            best_k = k

    dominant = sorted_idx[0]
    return best_k, dominant


print("\n" + "=" * 65)
print("■ 修正ウィーバー法による都道府県類型化")
print("=" * 65)

dominant_idx = np.array([modified_weaver(shares_pct[i])[1] for i in range(47)])
dominant_k   = np.array([modified_weaver(shares_pct[i])[0] for i in range(47)])
weaver_type  = np.array([SECTORS[d] + '型' for d in dominant_idx])

df_e['ウィーバー類型'] = weaver_type
df_e['ウィーバーk']    = dominant_k

# ── 人口密度区分（地域類型の代理） ───────────────────────────────────────────
density_q33 = df_e['人口密度'].quantile(0.33)
density_q67 = df_e['人口密度'].quantile(0.67)

def density_group(d):
    if d >= density_q67:
        return '高密度（都市部）'
    elif d >= density_q33:
        return '中密度（地方中核）'
    else:
        return '低密度（農村・山間）'

df_e['密度区分'] = df_e['人口密度'].apply(density_group)

type_counts = df_e['ウィーバー類型'].value_counts()
print(f"\n総観測数: 47 都道府県")
print(f"\n【ウィーバー類型の分布】")
for t, c in type_counts.items():
    print(f"  {t:<14} {c:>3}都道府県  ({c/47*100:.1f}%)")

print(f"\n【密度区分 × ウィーバー類型のクロス集計（行%）】")
ct = pd.crosstab(df_e['密度区分'], df_e['ウィーバー類型'], normalize='index').round(3)
print(ct * 100)

print(f"\n【リサイクル率の記述統計】")
print(df_e['ごみのリサイクル率'].describe().round(2))

mfg_share_pref = df_e['事業所数（民営）（製造業）'] / (
    df_e['事業所数（民営）（製造業）'] + df_e['事業所数（民営）（建設業）']
)
mfg_nat_share = mfg_share_pref.mean()
df_e['LQ製造業'] = mfg_share_pref / mfg_nat_share

print(f"\n  製造業LQ: 平均={df_e['LQ製造業'].mean():.3f}, 最大={df_e['LQ製造業'].max():.3f} ({df_e.loc[df_e['LQ製造業'].idxmax(), '都道府県']})")

print("\n" + "=" * 65)
print("✓ 全図の生成完了（4枚）")
print("=" * 65)
print("\n【最終サマリー】")
print(f"\n  データ: SSDSE-E-2026 実公的データ（47都道府県）")
print(f"  全体リサイクル率平均: {df_e['ごみのリサイクル率'].mean():.1f}%")
print(f"  1人1日排出量平均: {df_e['1人1日当たりの排出量'].mean():.0f} g")
print(f"\n  ウィーバー類型別リサイクル率（平均）:")
for t in TYPE_ORDER:
    sub = df_e[df_e['ウィーバー類型'] == t]
    m = sub['ごみのリサイクル率'].mean()
    n = len(sub)
    print(f"    {t:<14} {m:.1f}%  (n={n})")
print(f"\n  製造業LQトップ5:")
top5 = df_e.nlargest(5, 'LQ製造業')[['都道府県', 'LQ製造業', 'ごみのリサイクル率']]
for _, row in top5.iterrows():
    print(f"    {row['都道府県']}: LQ={row['LQ製造業']:.2f}, リサイクル率={row['ごみのリサイクル率']:.1f}%")

print("\n図1: ウィーバー類型別環境指標グラフを作成中...")

fig1, axes1 = plt.subplots(1, 2, figsize=(13, 5))
fig1.suptitle('修正ウィーバー法による類型化と廃棄物・環境特性\n（実データ：SSDSE-E-2026）',
              fontsize=13, fontweight='bold')

# (a) 類型別リサイクル率（平均±SD）
ax = axes1[0]
means_rc = [df_e[df_e['ウィーバー類型'] == t]['ごみのリサイクル率'].mean() for t in TYPE_ORDER]
stds_rc  = [df_e[df_e['ウィーバー類型'] == t]['ごみのリサイクル率'].std()  for t in TYPE_ORDER]
colors   = [TYPE_COLORS.get(t, '#888') for t in TYPE_ORDER]
x = np.arange(len(TYPE_ORDER))
bars = ax.bar(x, means_rc, yerr=stds_rc, capsize=6, color=colors, alpha=0.85,
              edgecolor='white', error_kw={'elinewidth': 1.5})
ax.axhline(df_e['ごみのリサイクル率'].mean(), color='gray', linestyle='--',
           linewidth=1.2, label=f'全体平均 ({df_e["ごみのリサイクル率"].mean():.1f}%)')
ax.set_xticks(x)
ax.set_xticklabels(TYPE_ORDER, fontsize=10)
ax.set_ylabel('ごみのリサイクル率（%）', fontsize=11)
ax.set_title('ウィーバー類型別 リサイクル率\n（平均±1SD）', fontsize=11, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(axis='y', alpha=0.3)
for i, (m, s) in enumerate(zip(means_rc, stds_rc)):
    n_type = (df_e['ウィーバー類型'] == TYPE_ORDER[i]).sum()
    ax.text(i, m + (s if not np.isnan(s) else 0) + 0.5,
            f'{m:.1f}%\n(n={n_type})', ha='center', fontsize=8.5, fontweight='bold')

# (b) 密度区分別の1人1日当たり排出量（ボックスプロット）
ax2 = axes1[1]
data_by_density = [df_e[df_e['密度区分'] == d]['1人1日当たりの排出量'].values
                   for d in DENSITY_ORDER]
bp = ax2.boxplot(data_by_density, patch_artist=True, notch=False, widths=0.55)
dp_colors = [DENSITY_COLORS[d] for d in DENSITY_ORDER]
for patch, color in zip(bp['boxes'], dp_colors):
    patch.set_facecolor(color)
    patch.set_alpha(0.75)
for med in bp['medians']:
    med.set(color='white', linewidth=2)
ax2.set_xticklabels(['高密度\n（都市部）', '中密度\n（地方中核）', '低密度\n（農村・山間）'],
                    fontsize=9)
ax2.set_ylabel('1人1日当たりの排出量（g）', fontsize=11)
ax2.set_title('人口密度区分別\n1人1日当たりごみ排出量分布', fontsize=11, fontweight='bold')
ax2.grid(axis='y', alpha=0.3)
ax2.axhline(df_e['1人1日当たりの排出量'].mean(), color='gray', linestyle='--', linewidth=1.0)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_U5_1_fig1_co2_bar.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  → 2025_U5_1_fig1_co2_bar.png 保存完了")

print("図2: クロス集計ヒートマップを作成中...")

ct_abs = pd.crosstab(df_e['密度区分'], df_e['ウィーバー類型'])
# 存在する列のみ使う
type_order_exist   = [t for t in TYPE_ORDER if t in ct_abs.columns]
density_order_exist = [d for d in DENSITY_ORDER if d in ct_abs.index]
ct_abs  = ct_abs.reindex(index=density_order_exist, columns=type_order_exist, fill_value=0)
ct_row  = ct_abs.div(ct_abs.sum(axis=1), axis=0) * 100

fig2, axes2 = plt.subplots(1, 2, figsize=(13, 5))
fig2.suptitle('人口密度区分 × ウィーバー類型 クロス集計\n（実データ：SSDSE-E-2026, 47都道府県）',
              fontsize=13, fontweight='bold')

for ax, data, title, fmt in zip(
    axes2,
    [ct_abs, ct_row],
    ['実数（都道府県数）', '行%（密度区分別の構成比）'],
    ['.0f', '.1f'],
):
    im = ax.imshow(data.values, cmap='Blues', aspect='auto')
    ax.set_xticks(range(len(type_order_exist)))
    ax.set_yticks(range(len(density_order_exist)))
    ax.set_xticklabels(type_order_exist, fontsize=10, rotation=20, ha='right')
    ax.set_yticklabels(density_order_exist, fontsize=10)
    ax.set_title(title, fontsize=11, fontweight='bold')
    for i in range(len(density_order_exist)):
        for j in range(len(type_order_exist)):
            val = data.values[i, j]
            max_val = data.values.max()
            text_color = 'white' if (max_val > 0 and val > max_val * 0.6) else 'black'
            ax.text(j, i, f'{val:{fmt}}', ha='center', va='center',
                    fontsize=9, fontweight='bold', color=text_color)
    plt.colorbar(im, ax=ax, shrink=0.85)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_U5_1_fig2_crosstab.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  → 2025_U5_1_fig2_crosstab.png 保存完了")

print("図3: 産業特化係数×環境指標散布図を作成中...")

fig3, axes3 = plt.subplots(1, 2, figsize=(13, 5))
fig3.suptitle('産業特化係数（製造業）と廃棄物・環境指標\n（実データ：SSDSE-E-2026）',
              fontsize=13, fontweight='bold')

# (a) LQ vs リサイクル率（密度区分で色分け）
ax3a = axes3[0]
for d, col in DENSITY_COLORS.items():
    mask = df_e['密度区分'] == d
    ax3a.scatter(df_e.loc[mask, 'LQ製造業'],
                 df_e.loc[mask, 'ごみのリサイクル率'],
                 color=col, alpha=0.75, s=60, label=d, zorder=3)
    # 都道府県名をラベル（LQ極値のみ）
    for _, row in df_e[mask].iterrows():
        if row['LQ製造業'] > 1.3 or row['LQ製造業'] < 0.6:
            ax3a.annotate(row['都道府県'][:2], (row['LQ製造業'], row['ごみのリサイクル率']),
                          fontsize=7, ha='left', va='bottom', color='#333')

ax3a.axvline(1.0, color='gray', linestyle='--', linewidth=1.0, label='LQ=1（全国平均）')
ax3a.set_xlabel('製造業の産業特化係数（LQ）', fontsize=11)
ax3a.set_ylabel('ごみのリサイクル率（%）', fontsize=11)
ax3a.set_title('LQ（製造業）vs リサイクル率\n（人口密度区分別）', fontsize=11, fontweight='bold')
ax3a.legend(fontsize=8, markerscale=1.5)
ax3a.grid(True, alpha=0.2)

# 回帰直線
lq_vals = df_e['LQ製造業'].values
rc_vals = df_e['ごみのリサイクル率'].values
coef = np.polyfit(lq_vals, rc_vals, 1)
x_range = np.linspace(lq_vals.min(), lq_vals.max(), 100)
ax3a.plot(x_range, np.polyval(coef, x_range), color='black', linewidth=1.5,
          linestyle='-', label=f'回帰直線 (β={coef[0]:+.2f})', zorder=2)
ax3a.legend(fontsize=8, markerscale=1.5)

# (b) ウィーバー類型別の4部門平均シェア（積み上げ棒グラフ）
ax3b = axes3[1]
sector_colors_list = ['#1565C0', '#2E7D32', '#E65100', '#6A1B9A']
bottom = np.zeros(len(TYPE_ORDER))
for i, s in enumerate(SECTORS):
    means_s = [shares_pct[df_e['ウィーバー類型'] == t, i].mean()
               for t in TYPE_ORDER]
    ax3b.bar(np.arange(len(TYPE_ORDER)), means_s, bottom=bottom,
             color=sector_colors_list[i], label=s, alpha=0.85, edgecolor='white')
    bottom += np.array(means_s)

ax3b.set_xticks(np.arange(len(TYPE_ORDER)))
ax3b.set_xticklabels(TYPE_ORDER, fontsize=9, rotation=15, ha='right')
ax3b.set_ylabel('4部門平均シェア（%）', fontsize=11)
ax3b.set_title('ウィーバー類型別 部門構成比\n（廃棄物・農業・商業活動軸）',
               fontsize=11, fontweight='bold')
ax3b.legend(fontsize=9, loc='upper right')
ax3b.grid(axis='y', alpha=0.3)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_U5_1_fig3_lq.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  → 2025_U5_1_fig3_lq.png 保存完了")

print("図4: ウィーバー法説明図を作成中...")

# 代表的な都道府県を選ぶ（LQ高い＝製造業特化）
top_mfg_pref_idx = df_e['LQ製造業'].idxmax()
top_mfg_pref = df_e.loc[top_mfg_pref_idx, '都道府県']
example_shares = shares_pct[top_mfg_pref_idx]  # 実際の4部門シェア

fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('修正ウィーバー法（Modified Weaver Method）の解説\n（実データ例示）',
              fontsize=13, fontweight='bold')

# (a) W_k の変化
ks   = range(1, len(SECTORS) + 1)
wks  = []
sorted_sh_ex = np.sort(example_shares)[::-1]
for k in ks:
    top_k = sorted_sh_ex[:k]
    w = np.sum((top_k - 100.0 / k) ** 2)
    wks.append(w)

ax4a = axes4[0]
ax4a.plot(list(ks), wks, 'o-', color='#1565C0', linewidth=2.5, markersize=9,
          markerfacecolor='white', markeredgewidth=2.5)
opt_k = int(np.argmin(wks)) + 1
ax4a.scatter([opt_k], [wks[opt_k - 1]], color='#E53935', s=130, zorder=5,
             label=f'最適 k={opt_k}（W_k 最小）')
ax4a.set_xlabel('上位k部門数', fontsize=11)
ax4a.set_ylabel('W_k（SSD）', fontsize=11)
ax4a.set_title(f'W_k の変化（{top_mfg_pref} の例）\n部門シェア: {", ".join([f"{v:.0f}" for v in example_shares])}%',
               fontsize=11, fontweight='bold')
ax4a.legend(fontsize=10)
ax4a.grid(True, alpha=0.3)
for k, w in zip(ks, wks):
    ax4a.annotate(f'W_{k}={w:.0f}', (k, w), textcoords='offset points',
                  xytext=(8, 5), fontsize=9,
                  color='#E53935' if k == opt_k else '#555')

# (b) 実際のシェア vs 均等分配（k=opt_k）の比較棒グラフ
ax4b = axes4[1]
x_pos = np.arange(len(SECTORS))
width = 0.35
# 降順に並べ直して表示
sorted_sector_idx = np.argsort(example_shares)[::-1]
sorted_labels = [SECTORS[i] for i in sorted_sector_idx]
sorted_shares = example_shares[sorted_sector_idx]

ax4b.bar(x_pos - width / 2, sorted_shares, width, label='実際のシェア（実データ）',
         color='#1565C0', alpha=0.8, edgecolor='white')
eq = np.zeros(len(SECTORS))
eq[:opt_k] = 100.0 / opt_k
ax4b.bar(x_pos + width / 2, eq, width, label=f'均等分配（k={opt_k}）',
         color='#E65100', alpha=0.7, edgecolor='white')
ax4b.set_xticks(x_pos)
ax4b.set_xticklabels([f'{s}' for s in sorted_labels], fontsize=9, rotation=15, ha='right')
ax4b.set_ylabel('部門シェア（%）', fontsize=11)
dominant_sector = SECTORS[np.argmax(example_shares)]
ax4b.set_title(f'実際のシェア vs 均等分配（k={opt_k}）\n→ {top_mfg_pref}は「{dominant_sector}型」に分類',
               fontsize=11, fontweight='bold')
ax4b.legend(fontsize=9)
ax4b.grid(axis='y', alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_U5_1_fig4_weaver.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  → 2025_U5_1_fig4_weaver.png 保存完了")