南海トラフ地震に備えて
指定避難所に3人に1人が避難できず、災害時の医療体制は本当に十分か？

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

FIG_DIR = 'html/figures'
DATA_A  = 'data/raw/SSDSE-A-2025.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ── データ読み込み（SSDSE-A：市区町村データ） ──────────────────
# 1行目=項目コード、2行目=年度、3行目=日本語の項目名。コード行をヘッダーに使う
df_a = pd.read_csv(DATA_A, encoding='cp932', header=0, skiprows=[1, 2])
df_a = df_a.rename(columns={
    'SSDSE-A-2025': '地域コード', 'Prefecture': '都道府県', 'Municipality': '市区町村',
    'A1101': '総人口', 'E2101': '小学校数', 'E3101': '中学校数',
    'E4101': '高等学校数', 'G1201': '公民館数', 'I6100': '医師数',
})
cols = ['総人口', '小学校数', '中学校数', '高等学校数', '公民館数', '医師数']
df_a[cols] = df_a[cols].apply(pd.to_numeric, errors='coerce')

# ── 都道府県庁所在市の抽出（東京は原論文にならい23区を合算） ──────
capitals = {
    '北海道': '札幌市', '青森県': '青森市', '岩手県': '盛岡市', '宮城県': '仙台市',
    '秋田県': '秋田市', '山形県': '山形市', '福島県': '福島市', '茨城県': '水戸市',
    '栃木県': '宇都宮市', '群馬県': '前橋市', '埼玉県': 'さいたま市', '千葉県': '千葉市',
    '神奈川県': '横浜市', '新潟県': '新潟市', '富山県': '富山市', '石川県': '金沢市',
    '福井県': '福井市', '山梨県': '甲府市', '長野県': '長野市', '岐阜県': '岐阜市',
    '静岡県': '静岡市', '愛知県': '名古屋市', '三重県': '津市', '滋賀県': '大津市',
    '京都府': '京都市', '大阪府': '大阪市', '兵庫県': '神戸市', '奈良県': '奈良市',
    '和歌山県': '和歌山市', '鳥取県': '鳥取市', '島根県': '松江市', '岡山県': '岡山市',
    '広島県': '広島市', '山口県': '山口市', '徳島県': '徳島市', '香川県': '高松市',
    '愛媛県': '松山市', '高知県': '高知市', '福岡県': '福岡市', '佐賀県': '佐賀市',
    '長崎県': '長崎市', '熊本県': '熊本市', '大分県': '大分市', '宮崎県': '宮崎市',
    '鹿児島県': '鹿児島市', '沖縄県': '那覇市',
}
rows = []
for pref, city in capitals.items():
    rows.append(df_a[(df_a['都道府県'] == pref) & (df_a['市区町村'] == city)][cols].iloc[0])
df_cap = pd.DataFrame(rows)
df_cap.insert(0, '市', list(capitals.values()))

# 東京23区（地域コード R13101〜R13123）を合算して1行にする
wards23 = df_a[df_a['地域コード'].astype(str).str.match(r'R131[0-2]\d')]
tokyo = wards23[cols].sum()
tokyo['市'] = '東京23区'
df_cap = pd.concat([df_cap, tokyo.to_frame().T], ignore_index=True)
df_cap[cols] = df_cap[cols].apply(pd.to_numeric)
print(f'都道府県庁所在市 {len(df_cap)} 市（東京23区を含む）')
print(f'東京23区の合算人口: {int(tokyo["総人口"]):,} 人')

# ── 指標の計算（原論文 2節の定義どおり） ─────────────────────────
# 避難所 ＝ 小学校＋中学校＋高等学校＋公民館 とみなす（原論文の定義）
df_cap['避難所数'] = (df_cap['小学校数'] + df_cap['中学校数']
                      + df_cap['高等学校数'] + df_cap['公民館数'])
df_cap['避難所1か所当たり人数'] = df_cap['総人口'] / df_cap['避難所数']
df_cap['人口1万人あたり医師数'] = df_cap['医師数'] / df_cap['総人口'] * 10000

mean_shelter = df_cap['避難所1か所当たり人数'].mean()
mean_doctor  = df_cap['人口1万人あたり医師数'].mean()
m = df_cap[df_cap['市'] == '松山市'].iloc[0]
print(f'\n【避難所1か所当たりの人数】47市平均: {mean_shelter:,.0f} 人')
print('上位5市:')
top5 = df_cap.nlargest(5, '避難所1か所当たり人数')
for _, r in top5.iterrows():
    print(f"  {r['市']}: {r['避難所1か所当たり人数']:,.0f} 人")
print(f"松山市: {m['避難所1か所当たり人数']:,.0f} 人（平均の {m['避難所1か所当たり人数']/mean_shelter:.2f} 倍）")
print(f'\n【人口1万人あたりの医師数】47市平均: {mean_doctor:.1f} 人')
low5 = df_cap.nsmallest(5, '人口1万人あたり医師数')
print('下位5市:')
for _, r in low5.iterrows():
    print(f"  {r['市']}: {r['人口1万人あたり医師数']:.1f} 人")
print(f"松山市: {m['人口1万人あたり医師数']:.1f} 人（平均より "
      f"{'低い' if m['人口1万人あたり医師数'] < mean_doctor else '高い'}）")

# ── 図1：避難所1か所当たりの人数（原論文 図1 の再現） ─────────────
d1 = df_cap.sort_values('避難所1か所当たり人数', ascending=True)
colors1 = ['#C62828' if c == '松山市' else
           ('#F57C00' if c in ('東京23区', '札幌市', '名古屋市', '大阪市', '横浜市')
            else '#78909C') for c in d1['市']]
fig1, ax1 = plt.subplots(figsize=(9, 13))
ax1.barh(d1['市'], d1['避難所1か所当たり人数'], color=colors1, alpha=0.9, height=0.72)
ax1.axvline(mean_shelter, color='#333333', linewidth=2, linestyle='--',
            label=f'47市平均 {mean_shelter:,.0f} 人', zorder=5)
ax1.set_xlabel('避難所1か所当たりの人数（人）', fontsize=12)
ax1.set_title('都道府県庁所在市別 避難所1か所当たりの人数\n'
              '（避難所＝小学校＋中学校＋高等学校＋公民館、SSDSE-A-2025で再計算）',
              fontsize=13, fontweight='bold', pad=12)
ax1.tick_params(axis='y', labelsize=8.5)
ax1.grid(True, axis='x', alpha=0.3)
ax1.legend(loc='lower right', fontsize=10)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2019_H2_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print('\n図1 保存完了（避難所1か所当たりの人数）')

# ── 図2：人口1万人あたりの医師数（原論文 図2 の再現） ─────────────
d2 = df_cap.sort_values('人口1万人あたり医師数', ascending=True)
colors2 = ['#C62828' if c == '松山市' else
           ('#1565C0' if c == 'さいたま市' else '#78909C') for c in d2['市']]
fig2, ax2 = plt.subplots(figsize=(9, 13))
ax2.barh(d2['市'], d2['人口1万人あたり医師数'], color=colors2, alpha=0.9, height=0.72)
ax2.axvline(mean_doctor, color='#333333', linewidth=2, linestyle='--',
            label=f'47市平均 {mean_doctor:.1f} 人', zorder=5)
ax2.set_xlabel('人口1万人あたりの医師数（人）', fontsize=12)
ax2.set_title('都道府県庁所在市別 人口1万人あたりの医師数\n'
              '（SSDSE-A-2025で再計算。赤＝松山市、青＝さいたま市）',
              fontsize=13, fontweight='bold', pad=12)
ax2.tick_params(axis='y', labelsize=8.5)
ax2.grid(True, axis='x', alpha=0.3)
ax2.legend(loc='lower right', fontsize=10)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2019_H2_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print('図2 保存完了（人口1万人あたりの医師数）')

# ── 図3：半径3km圏の人口と収容可能人数（原論文 表1 の報告値の可視化） ──
# 以下の3つの数値は原論文 表1 の報告値（jSTAT MAP・松山市指定緊急避難場所
# 一覧表から原論文が集計したもの）。本スクリプトでの再計算ではない。
pop_3km      = 207828   # 黒線内エリアの人口総数（人）
capacity_3km = 142164   # 合計収容可能人数（人）
shelters_3km = 67       # 指定避難所の数（か所）
ratio = capacity_3km / pop_3km * 100
shortage = pop_3km - capacity_3km
print(f'\n【原論文 表1 の報告値】松山南高校 半径3km圏')
print(f'  人口総数        : {pop_3km:,} 人')
print(f'  合計収容可能人数: {capacity_3km:,} 人（避難所 {shelters_3km} か所）')
print(f'  収容率 = {capacity_3km:,} / {pop_3km:,} × 100 ≒ {ratio:.2f} %')
print(f'  → 約 {shortage:,} 人（およそ3人に1人）が指定避難所に避難できない')

fig3, ax3 = plt.subplots(figsize=(9, 4.6))
ax3.barh(['合計収容可能人数', '人口総数'], [capacity_3km, pop_3km],
         color=['#1565C0', '#78909C'], alpha=0.9, height=0.55)
ax3.barh(['合計収容可能人数'], [pop_3km - capacity_3km], left=[capacity_3km],
         color='#C62828', alpha=0.35, height=0.55,
         label=f'収容できない人数 約{shortage:,}人')
for y, v in zip([0, 1], [capacity_3km, pop_3km]):
    ax3.annotate(f'{v:,} 人', (v, y), xytext=(6, 0), textcoords='offset points',
                 va='center', fontsize=11, fontweight='bold')
ax3.annotate(f'収容率 {ratio:.1f}%\n＝ 3人に1人が避難できない',
             (capacity_3km * 0.5, 0), va='center', ha='center',
             fontsize=11, color='white', fontweight='bold')
ax3.set_xlim(0, 235000)
ax3.set_title('松山南高校 半径3km圏の人口総数と指定避難所の収容可能人数\n'
              '（原論文 表1 の報告値の可視化。再計算ではない）',
              fontsize=12.5, fontweight='bold', pad=10)
ax3.set_xlabel('人数（人）', fontsize=11)
ax3.grid(True, axis='x', alpha=0.3)
ax3.legend(loc='lower right', fontsize=10)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2019_H2_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print('図3 保存完了（表1の報告値の可視化）')

# ── 図4：エリア①〜④の報告値の可視化（原論文 4.1節・6節） ─────────
# いずれも原論文の報告値。徒歩10分＝分速80m×10分＝800mの到達圏に入らない
# 4つの空白エリアの人数と、提案場所の導入で新たに避難可能となる人数。
areas = ['エリア①\n（私立幼稚園2園＋\n大学グラウンド）',
         'エリア②\n（ショッピングセンター＋\n認定こども園）',
         'エリア③\n（自動車教習所＋\n私立保育所・幼稚園）',
         'エリア④\n（松山市総合公園＋\n南江戸公園）']
unreachable = [1400, 7000, 6000, 4000]      # 徒歩10分で避難所に到達できない人数
new_capacity = [30000, 16000, 6600, 20000]  # 提案導入で新たに避難可能となる人数

print('\n【原論文の報告値】徒歩10分で指定避難所に到達できないエリア')
for a, u, n in zip(['①', '②', '③', '④'], unreachable, new_capacity):
    print(f'  エリア{a}: 到達困難 約{u:,}人 → 提案場所の導入で新たに約{n:,}人が避難可能')
print(f'  合計    : 到達困難 約{sum(unreachable):,}人 / 新規避難可能 約{sum(new_capacity):,}人')

y = np.arange(len(areas))
fig4, ax4 = plt.subplots(figsize=(10, 6.2))
ax4.barh(y + 0.21, unreachable, height=0.38, color='#C62828', alpha=0.9,
         label='徒歩10分以内に指定避難所へ到達できない人数')
ax4.barh(y - 0.21, new_capacity, height=0.38, color='#2E7D32', alpha=0.9,
         label='提案場所の導入で新たに避難可能となる人数')
for yi, v in zip(y + 0.21, unreachable):
    ax4.annotate(f'約{v:,}人', (v, yi), xytext=(5, 0), textcoords='offset points',
                 va='center', fontsize=10)
for yi, v in zip(y - 0.21, new_capacity):
    ax4.annotate(f'約{v:,}人', (v, yi), xytext=(5, 0), textcoords='offset points',
                 va='center', fontsize=10)
ax4.set_yticks(y)
ax4.set_yticklabels(areas, fontsize=9.5)
ax4.invert_yaxis()
ax4.set_xlim(0, 34000)
ax4.set_xlabel('人数（人）', fontsize=11)
ax4.set_title('避難所空白エリア①〜④の到達困難人数と提案の効果\n'
              '（原論文 4.1節・6節の報告値の可視化。再計算ではない）',
              fontsize=12.5, fontweight='bold', pad=10)
ax4.grid(True, axis='x', alpha=0.3)
ax4.legend(loc='lower right', fontsize=10)
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2019_H2_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print('図4 保存完了（エリア①〜④の報告値の可視化）')