家計調査に基づく消費重心の計算

import os
import math
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Patch

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

FIG_DIR = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

NOTE = '報告値の可視化（再計算ではない）'  # 全図共通の注記

# ══════════════════════════════════════════════════════════
# 原論文の報告値（PDFのタブル・本文を転記。単位は decimal degrees）
# ══════════════════════════════════════════════════════════

# 表1 大分類の消費重心（経度, 緯度）
DAIBUNRUI = [
    ('穀類',        137.1216, 35.41703),
    ('魚介類',      137.2183, 35.48567),
    ('肉類',        136.9404, 35.30843),
    ('乳卵類',      137.1518, 35.28413),
    ('野菜・海藻',  137.2787, 35.43841),
    ('果物',        137.2931, 35.44298),
    ('油脂・調味料', 137.0571, 35.36921),
    ('菓子類',      137.1802, 35.42998),
    ('調理食品',    137.1772, 35.32412),
    ('飲料',        137.1688, 35.39765),
    ('酒類',        137.2205, 35.55091),
    ('外食',        137.3352, 35.35948),
]
GOUKEI = ('食料（合計）', 137.1959, 35.39448)  # 表1 食料合計

# 図1・図2 小分類213品目の消費重心 経度・緯度の基本統計量（報告値）
STATS_LON = dict(平均=137.1901, 標準誤差=0.02256, 中央値=137.1637, 標準偏差=0.32928,
                 分散=0.10843, 尖度=2.17056, 歪度=-0.22126, 最小=135.7723, 最大=138.2165, n=213)
STATS_LAT = dict(平均=35.4238, 標準誤差=0.01231, 中央値=35.4069, 標準偏差=0.17970,
                 分散=0.03229, 尖度=4.04831, 歪度=0.66821, 最小=34.7792, 最大=36.1986, n=213)

# 表2 経度に特徴のある消費品目（西側=経度が小さい / 東側=経度が大きい）
WEST = [  # 西に位置する10品目（経度が小さい順）
    ('たい', 135.7723, 34.7792), ('まんじゅう', 136.2335, 35.05084),
    ('あじ', 136.2797, 34.8392), ('合いびき肉', 136.3631, 34.9684),
    ('揚げかまぼこ', 136.4847, 35.0584), ('ぶり', 136.4962, 35.0493),
    ('さば', 136.5260, 35.1315), ('そうざい材料セット', 136.5547, 35.1627),
    ('牛肉', 136.5805, 35.0234), ('いわし', 136.5829, 35.1546),
]
EAST = [  # 東に位置する10品目（経度が大きい順）
    ('まぐろ', 138.2165, 35.5758), ('グレープフルーツ', 138.0606, 35.7085),
    ('しじみ', 138.0287, 35.9265), ('魚介の漬物', 138.0014, 35.9617),
    ('しゅうまい', 137.9851, 35.4739), ('ワイン', 137.9358, 35.5691),
    ('ウイスキー', 137.8987, 35.9656), ('塩さけ', 137.8681, 35.9674),
    ('ほたて貝', 137.8290, 36.1005), ('さやまめ', 137.7646, 35.4946),
]

# 表3 緯度に特徴のある消費品目（南側=緯度が小さい / 北側=緯度が大きい）
SOUTH = [
    ('たい', 135.7723, 34.7792), ('あじ', 136.2797, 34.8392),
    ('合いびき肉', 136.3631, 34.9684), ('牛肉', 136.5805, 35.0234),
    ('ぶり', 136.4962, 35.0493), ('揚げかまぼこ', 136.4847, 35.0584),
    ('まんじゅう', 136.2335, 35.05084), ('かつお節・削り節', 136.8234, 35.0757),
    ('さば', 136.5260, 35.1315), ('いわし', 136.5829, 35.1546),
]
NORTH = [
    ('メロン', 137.681509, 36.1986493), ('ほたて貝', 137.828986, 36.1005324),
    ('塩さけ', 137.868137, 35.9674378), ('ウイスキー', 137.898738, 35.9655654),
    ('魚介の漬物', 138.001437, 35.9616854), ('かれい', 137.168282, 35.9436351),
    ('しじみ', 138.028741, 35.9265045), ('他の貝', 137.528317, 35.8431935),
    ('さんま', 137.562414, 35.7237647), ('グレープフルーツ', 138.06064, 35.7084886),
]

# 本文で言及される代表品目（対照分析用, 報告値）
CONTRAST = {
    '牛肉':   (136.5805, 35.02335),
    '豚肉':   (137.1804, 35.43735),
    'たい':   (135.7723, 34.7792),
    'まぐろ': (138.2165, 35.5758),
}

# 人口重心（総務省統計局2017, 平成27年国勢調査, 全国）
# 原論文本文の報告値: 東経137度02分15.84秒, 北緯35度34分51.44秒
JINKO_JUSHIN_DMS = '東経137度02分15.84秒 / 北緯35度34分51.44秒'
JINKO_JUSHIN = (137 + 2/60 + 15.84/3600, 35 + 34/60 + 51.44/3600)  # ≈(137.0377, 35.5810)
REPORTED_KM_JINKO = 16.75    # 人口重心と食料全体の消費重心の距離（報告値）
REPORTED_KM_BUTA_GYU = 71.411  # 豚肉と牛肉の消費重心間距離（報告値）


def haversine_km(a, b):
    """2地点 (lon, lat) 間の大円距離[km]。報告座標の整合性確認に使用。"""
    R = 6371.0
    lon1, lat1 = math.radians(a[0]), math.radians(a[1])
    lon2, lat2 = math.radians(b[0]), math.radians(b[1])
    dlon, dlat = lon2 - lon1, lat2 - lat1
    h = math.sin(dlat/2)**2 + math.cos(lat1)*math.cos(lat2)*math.sin(dlon/2)**2
    return 2 * R * math.asin(math.sqrt(h))

# ══════════════════════════════════════════════════════════
# コンソール出力：報告値の確認
# ══════════════════════════════════════════════════════════
print('=' * 64)
print('家計調査に基づく消費重心の計算 — 原論文の報告値（可視化用）')
print('=' * 64)
print(f'食料全体の消費重心（表1 食料合計）: 経度 {GOUKEI[1]}, 緯度 {GOUKEI[2]}')
print(f'人口重心（平成27年国勢調査）      : {JINKO_JUSHIN_DMS}')
print(f'  → 消費重心は人口重心より東北東へ約 {REPORTED_KM_JINKO} km（原論文の報告値）')

print('\n── 大分類で最も西／東に位置した消費重心（表1）──')
lon_sorted = sorted(DAIBUNRUI, key=lambda t: t[1])
print(f'  最も西: {lon_sorted[0][0]}（経度 {lon_sorted[0][1]}）')
print(f'  最も東: {lon_sorted[-1][0]}（経度 {lon_sorted[-1][1]}）')

print('\n── 小分類213品目の経度・緯度の基本統計量（図1・図2, 報告値）──')
print(f'  経度: 平均{STATS_LON["平均"]} 中央値{STATS_LON["中央値"]} '
      f'標準偏差{STATS_LON["標準偏差"]} 歪度{STATS_LON["歪度"]} 尖度{STATS_LON["尖度"]}')
print(f'  緯度: 平均{STATS_LAT["平均"]} 中央値{STATS_LAT["中央値"]} '
      f'標準偏差{STATS_LAT["標準偏差"]} 歪度{STATS_LAT["歪度"]} 尖度{STATS_LAT["尖度"]}')

print('\n── 報告座標による品目間距離の整合性チェック（haversine）──')
d_bg = haversine_km(CONTRAST['牛肉'], CONTRAST['豚肉'])
d_tm = haversine_km(CONTRAST['たい'], CONTRAST['まぐろ'])
print(f'  牛肉↔豚肉: 計算 {d_bg:.3f} km（原論文の報告値 {REPORTED_KM_BUTA_GYU} km）')
print(f'  たい↔まぐろ: 計算 {d_tm:.3f} km（東西の代替的関係）')

# ══════════════════════════════════════════════════════════
# 図1: 大分類12品目の消費重心の分布（経度×緯度, 表1の報告値）
# ══════════════════════════════════════════════════════════
fig1, ax1 = plt.subplots(figsize=(10, 7))
lons = [t[1] for t in DAIBUNRUI]
lats = [t[2] for t in DAIBUNRUI]
ax1.scatter(lons, lats, s=90, color='#f0a500', edgecolors='#7c5a00',
            linewidths=0.8, zorder=3, label='大分類の消費重心（12品目）')
for name, lon, lat in DAIBUNRUI:
    ax1.annotate(name, (lon, lat), fontsize=9, ha='left', va='bottom',
                 xytext=(3, 3), textcoords='offset points')
# 食料合計
ax1.scatter([GOUKEI[1]], [GOUKEI[2]], s=180, marker='*', color='#c62828',
            edgecolors='black', linewidths=0.8, zorder=4, label='食料（合計）の消費重心')
ax1.annotate('食料（合計）', (GOUKEI[1], GOUKEI[2]), fontsize=10, fontweight='bold',
             color='#c62828', ha='right', va='top', xytext=(-4, -4), textcoords='offset points')
# 最西・最東の強調（重ね塗り）
ax1.scatter([lon_sorted[0][1]], [lon_sorted[0][2]], s=90, color='#1565C0',
            edgecolors='black', linewidths=0.8, zorder=5)
ax1.scatter([lon_sorted[-1][1]], [lon_sorted[-1][2]], s=90, color='#2E7D32',
            edgecolors='black', linewidths=0.8, zorder=5)
ax1.text(0.02, 0.97, '青＝最も西（肉類）　緑＝最も東（外食）',
         transform=ax1.transAxes, fontsize=10, va='top', ha='left',
         bbox=dict(boxstyle='round,pad=0.3', fc='#FFF8E1', ec='#F9A825', alpha=0.9))

ax1.set_xlabel('経度（東経・度）→ 東', fontsize=12)
ax1.set_ylabel('緯度（北緯・度）→ 北', fontsize=12)
ax1.set_title(f'図1  大分類の消費重心の分布（経度×緯度）\n原論文 表1・{NOTE}',
              fontsize=13, fontweight='bold')
ax1.legend(loc='lower right', fontsize=10, framealpha=0.9)
ax1.grid(True, alpha=0.3, linestyle='--')
ax1.spines[['top', 'right']].set_visible(False)
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_U5_5_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print('\n図1 保存完了: 2020_U5_5_fig1.png')

# ══════════════════════════════════════════════════════════
# 図2: 経度に特徴のある消費品目（西10・東10, 表2の報告値）
# ══════════════════════════════════════════════════════════
items = list(reversed(WEST)) + list(reversed(EAST))  # 下=西, 上=東 になるよう配置
names = [t[0] for t in items]
vals = [t[1] for t in items]
colors = ['#1565C0'] * len(WEST) + ['#2E7D32'] * len(EAST)

fig2, ax2 = plt.subplots(figsize=(10, 9))
ax2.barh(range(len(items)), vals, color=colors, alpha=0.85, edgecolor='white', linewidth=0.5)
ax2.axvline(GOUKEI[1], color='#c62828', linewidth=1.5, linestyle='--',
            label=f'食料（合計）の経度 {GOUKEI[1]}')
ax2.set_yticks(range(len(items)))
ax2.set_yticklabels(names, fontsize=10)
for i, v in enumerate(vals):
    ax2.text(v + 0.006, i, f'{v:.4f}', va='center', fontsize=8, color='#333')
ax2.set_xlabel('経度（東経・度）　←西　　東→', fontsize=12)
ax2.set_title(f'図2  経度に特徴のある消費品目（西10・東10）\n原論文 表2・{NOTE}',
              fontsize=13, fontweight='bold')
ax2.set_xlim(135.5, 138.5)
handles = [Patch(facecolor='#1565C0', alpha=0.85, label='西に位置（暖流・西日本の漁業等）'),
           Patch(facecolor='#2E7D32', alpha=0.85, label='東に位置（寒流・東日本の産地等）'),
           plt.Line2D([0], [0], color='#c62828', linestyle='--', label=f'食料（合計）{GOUKEI[1]}')]
ax2.legend(handles=handles, loc='lower right', fontsize=9, framealpha=0.9)
ax2.grid(True, axis='x', alpha=0.3, linestyle='--')
ax2.spines[['top', 'right']].set_visible(False)
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_U5_5_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print('図2 保存完了: 2020_U5_5_fig2.png')

# ══════════════════════════════════════════════════════════
# 図3: 緯度に特徴のある消費品目（南10・北10, 表3の報告値）
# ══════════════════════════════════════════════════════════
items3 = list(reversed(SOUTH)) + list(reversed(NORTH))  # 下=南, 上=北
names3 = [t[0] for t in items3]
vals3 = [t[2] for t in items3]
colors3 = ['#E65100'] * len(SOUTH) + ['#00695C'] * len(NORTH)

fig3, ax3 = plt.subplots(figsize=(10, 9))
ax3.barh(range(len(items3)), vals3, color=colors3, alpha=0.85, edgecolor='white', linewidth=0.5)
ax3.axvline(GOUKEI[2], color='#c62828', linewidth=1.5, linestyle='--',
            label=f'食料（合計）の緯度 {GOUKEI[2]}')
ax3.set_yticks(range(len(items3)))
ax3.set_yticklabels(names3, fontsize=10)
for i, v in enumerate(vals3):
    ax3.text(v + 0.004, i, f'{v:.4f}', va='center', fontsize=8, color='#333')
ax3.set_xlabel('緯度（北緯・度）　←南　　北→', fontsize=12)
ax3.set_title(f'図3  緯度に特徴のある消費品目（南10・北10）\n原論文 表3・{NOTE}',
              fontsize=13, fontweight='bold')
ax3.set_xlim(34.5, 36.4)
handles3 = [Patch(facecolor='#E65100', alpha=0.85, label='南に位置（九州・四国の影響）'),
            Patch(facecolor='#00695C', alpha=0.85, label='北に位置（北海道の名産品）'),
            plt.Line2D([0], [0], color='#c62828', linestyle='--', label=f'食料（合計）{GOUKEI[2]}')]
ax3.legend(handles=handles3, loc='lower right', fontsize=9, framealpha=0.9)
ax3.grid(True, axis='x', alpha=0.3, linestyle='--')
ax3.spines[['top', 'right']].set_visible(False)
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_U5_5_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print('図3 保存完了: 2020_U5_5_fig3.png')

# ══════════════════════════════════════════════════════════
# 図4: 品目間の対照（牛肉↔豚肉, たい↔まぐろ）と距離（報告値）
# ══════════════════════════════════════════════════════════
fig4, ax4 = plt.subplots(figsize=(10, 7))
pair_style = {
    '牛肉':   ('#C62828', 'o'), '豚肉': ('#C62828', 's'),
    'たい':   ('#1565C0', 'o'), 'まぐろ': ('#1565C0', 's'),
}
for name, (lon, lat) in CONTRAST.items():
    col, mk = pair_style[name]
    ax4.scatter([lon], [lat], s=140, color=col, marker=mk, edgecolors='black',
                linewidths=0.7, zorder=3)
    ax4.annotate(name, (lon, lat), fontsize=11, fontweight='bold', color=col,
                 ha='left', va='bottom', xytext=(4, 4), textcoords='offset points')
# 対照ペアを線で結ぶ
for pair, kmtext in [(('牛肉', '豚肉'), f'牛肉と豚肉 {REPORTED_KM_BUTA_GYU} km'),
                     (('たい', 'まぐろ'), 'たいとまぐろ（東西の代替）')]:
    a, b = CONTRAST[pair[0]], CONTRAST[pair[1]]
    ax4.plot([a[0], b[0]], [a[1], b[1]], color='#888', linestyle='--', linewidth=1.2, zorder=2)
    mx, my = (a[0]+b[0])/2, (a[1]+b[1])/2
    ax4.annotate(kmtext, (mx, my), fontsize=9, color='#555', ha='center', va='center',
                 bbox=dict(boxstyle='round,pad=0.2', fc='white', ec='#bbb', alpha=0.9))
# 食料合計
ax4.scatter([GOUKEI[1]], [GOUKEI[2]], s=200, marker='*', color='#f0a500',
            edgecolors='black', linewidths=0.8, zorder=4, label='食料（合計）')
ax4.annotate('食料（合計）', (GOUKEI[1], GOUKEI[2]), fontsize=10,
             ha='right', va='top', xytext=(-4, -4), textcoords='offset points')

ax4.set_xlabel('経度（東経・度）→ 東', fontsize=12)
ax4.set_ylabel('緯度（北緯・度）→ 北', fontsize=12)
ax4.set_title(f'図4  品目間の対照：牛肉と豚肉・たいとまぐろの消費重心\n原論文本文・{NOTE}',
              fontsize=13, fontweight='bold')
ax4.legend(loc='lower right', fontsize=10, framealpha=0.9)
ax4.grid(True, alpha=0.3, linestyle='--')
ax4.spines[['top', 'right']].set_visible(False)
fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_U5_5_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print('図4 保存完了: 2020_U5_5_fig4.png')

print('\n全図生成完了。')