# -*- coding: utf-8 -*-
"""L93 広島県クルーズ船観光客人流データ — 500m メッシュ空間分析
       (DoBoX #1567 / cruise_people_cnt.shp / 118 メッシュ)

カバー宣言:
  本記事は以下の既存記事と<b>完全に異なるデータ・分析軸</b>を持つ独立記事。
  - X04 (#1280 + #1282): 港別クルーズ寄港統計の時系列回帰 → 寄港回数・客数の変化
  - L82 (#1280): モニタリングクルーズの航海別乗降客数分析 → 収容率・分布
  - L83 (#1282): しまたびライン定期便の月別乗降客数分析 → 季節性・定期化効果
  - <b>L93（本記事）</b>: dataset #1567「クルーズ船観光客人流データ」= GPS/モバイル由来の
    <b>500m メッシュ人流カウント</b>を空間分析。観光客が「どこに行ったか」の動線把握。

  データ型の差別化:
    X04/L82/L83 = 乗降客「数」の時系列・港別統計（人数×時間）
    L93         = 観光客の「場所」情報（人流×空間メッシュ）
  → 同じ「クルーズ観光客」を全く異なる次元（時間 vs 空間）で分析する。

研究の問い（主 RQ）:
  広島県クルーズ船観光客の人流はどの空間に集中し、
  観光資源・防災リスクとの地理的関係はどうなっているか？

仮説（H1〜H9）:
  H1（集中型分布）: 上位10メッシュ(8.5%) が全人流カウント(219)の40%以上を占める。
  H2（高 Gini）: count 分布の Gini 係数 ≥ 0.5（少数の人気スポットに集中）。
  H3（空間クラスター）: 女王コンティギュイティの空間自己相関 Moran's I > 0
                        （隣接メッシュは独立分散より近い count を持つ）。
  H4（観光資源隣接）: count ≥ 3 の「高人流メッシュ」のうちインフラツーリズム施設
                       1km 圏内にあるメッシュ ≥ 50%。
  H5（人流重心・市街地）: 人口加重重心（count 加重）が広島市平和記念公園から 5km 以内。
  H6（埋蔵文化財接続）: 全118メッシュのうち 1km 圏内に埋蔵文化財を持つメッシュ ≥ 50%。
  H7（津波域内人流）: 津波浸水想定区域と重なるメッシュ ≥ 5個。
  H8（沿岸人流率）: 津波域内メッシュの count 合計 / 全 count ≥ 10%。
  H9（X04/L82/L83 差別化）: 7 観点比較表で全観点が異なるカテゴリに分類される。

要件 S（1〜3 分以内）:
  - Shapefile は 9976 bytes (10KB) の軽量 ZIP、即時 DL
  - 118 メッシュの sjoin・Moran's I 計算はすべて軽量
  - L49 / infra_tourism / L84 は既存キャッシュ再利用

実行:
    cd "2026 DoBoX 教材"
    py -X utf8 lessons/L93_cruise_visitor_analysis.py
"""
from __future__ import annotations
import sys, time, io, zipfile, json
from pathlib import Path
from html import escape

sys.path.insert(0, str(Path(__file__).parent))
from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure,
                     data_recipe)

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
import geopandas as gpd
from shapely.geometry import Point

t0 = time.time()

# ── 定数 ────────────────────────────────────────────────────
TARGET_CRS  = "EPSG:6671"
CACHE_DIR   = ROOT / "data/extras/L93_cruise_visitor"
CACHE_DIR.mkdir(parents=True, exist_ok=True)
ZIP_CACHE   = CACHE_DIR / "cruise_people_cnt.zip"

ADMIN_GPKG   = ROOT / "data/extras/L44_storm_surge/_cache/admin_diss.gpkg"
TSUNAMI_GPKG = ROOT / "data/extras/L49_tsunami_inundation/_cache/tsunami_dissolve_8rank.gpkg"
INFRA_CSV    = ROOT / "data/extras/infra_tourism.csv"
L84_DIR      = ROOT / "data/extras/L84_archaeological_sites"

# 平和記念公園（原爆ドーム） EPSG:4326
PEACE_PARK_LON, PEACE_PARK_LAT = 132.4535, 34.3955

CITY_NAME = {
    "34100": "広島市", "34202": "呉市", "34203": "竹原市", "34204": "三原市",
    "34205": "尾道市", "34207": "福山市", "34208": "府中市", "34209": "三次市",
    "34210": "庄原市", "34211": "大竹市", "34212": "東広島市", "34213": "廿日市市",
    "34214": "安芸高田市", "34215": "江田島市", "34302": "府中町", "34304": "海田町",
    "34305": "熊野町", "34307": "坂町", "34311": "安芸太田町", "34312": "北広島町",
    "34431": "大崎上島町", "34441": "世羅町", "34462": "神石高原町",
}

# ─────────────────────────────────────────────────────────────
# 1. データ取得
# ─────────────────────────────────────────────────────────────
print("[1] データ取得", flush=True)
t1 = time.time()

if not ZIP_CACHE.exists() or ZIP_CACHE.stat().st_size < 5000:
    import requests
    DOBOX_BASE = "https://hiroshima-dobox.jp"
    headers = {"User-Agent": "DoBoX-MDASH-textbook/1.0 (educational use)"}
    r = requests.get(f"{DOBOX_BASE}/resource_download/175545", headers=headers, timeout=60)
    r.raise_for_status()
    ZIP_CACHE.write_bytes(r.content)
    print(f"  DL: {len(r.content):,} bytes", flush=True)
else:
    print(f"  キャッシュ: {ZIP_CACHE.name}", flush=True)

# Extract shapefile
SHP_EXTRACT = CACHE_DIR / "shp_dobox"
if not (SHP_EXTRACT / "cruise_people_cnt.shp").exists():
    with zipfile.ZipFile(ZIP_CACHE) as zf:
        zf.extractall(CACHE_DIR)

gdf = gpd.read_file(SHP_EXTRACT / "cruise_people_cnt.shp")
# count フィールドを整数に
gdf['count'] = gdf['count'].astype(int)
n_mesh = len(gdf)
total_count = gdf['count'].sum()
max_count = gdf['count'].max()
mean_count = gdf['count'].mean()
print(f"  メッシュ数: {n_mesh} / 総カウント: {total_count} / max: {max_count} / mean: {mean_count:.2f}", flush=True)
print(f"  ({time.time()-t1:.1f}s)", flush=True)

# ─────────────────────────────────────────────────────────────
# 2. 補助データ
# ─────────────────────────────────────────────────────────────
print("[2] 補助データ", flush=True)
t2 = time.time()

admin = gpd.read_file(ADMIN_GPKG)
admin["市町名"] = admin["CITY_CD"].map(CITY_NAME).fillna(admin["CITY_CD"].astype(str))

infra = pd.read_csv(INFRA_CSV, encoding='utf-8-sig')
infra_gdf = gpd.GeoDataFrame(
    infra,
    geometry=gpd.points_from_xy(infra['経度'], infra['緯度']),
    crs="EPSG:4326"
).to_crs(TARGET_CRS)

maizo_dfs = []
for csv_path in sorted(L84_DIR.glob("maizo_*.csv")):
    try:
        df = pd.read_csv(csv_path, encoding='utf-8-sig')
        df['緯度'] = pd.to_numeric(df['緯度'], errors='coerce')
        df['経度'] = pd.to_numeric(df['経度'], errors='coerce')
        maizo_dfs.append(df)
    except Exception:
        pass
maizo = pd.concat(maizo_dfs, ignore_index=True).dropna(subset=['緯度', '経度'])
maizo = maizo[(maizo['緯度'].between(33.8, 35.1)) & (maizo['経度'].between(132.0, 134.0))]
maizo_gdf = gpd.GeoDataFrame(
    maizo,
    geometry=gpd.points_from_xy(maizo['経度'], maizo['緯度']),
    crs="EPSG:4326"
).to_crs(TARGET_CRS)

if TSUNAMI_GPKG.exists():
    tsunami = gpd.read_file(TSUNAMI_GPKG).to_crs(TARGET_CRS)
    tsuna_union = tsunami.union_all()
    print(f"  津波: {len(tsunami)} ポリゴン", flush=True)
else:
    tsunami = None
    tsuna_union = None

print(f"  infra: {len(infra_gdf)} / maizo: {len(maizo_gdf)}", flush=True)
print(f"  ({time.time()-t2:.1f}s)", flush=True)

# ─────────────────────────────────────────────────────────────
# 3. RQ1: 空間分布構造
# ─────────────────────────────────────────────────────────────
print("[3] RQ1 空間分布", flush=True)
t3 = time.time()

# Lorenz 曲線
sorted_cnt = np.sort(gdf['count'].values)
cum_share = np.cumsum(sorted_cnt) / sorted_cnt.sum()
pop_share = np.arange(1, len(sorted_cnt)+1) / len(sorted_cnt)
gini = 1 - 2 * np.trapezoid(cum_share, pop_share)

# Top10 集中度
top10_cnt = gdf.nlargest(10, 'count')['count'].sum()
top10_share = top10_cnt / total_count * 100

# count >= 3 の「高人流」メッシュ
high_mesh = gdf[gdf['count'] >= 3]
n_high = len(high_mesh)

# 人口加重重心 (EPSG:6671 内で計算)
cents = gdf.geometry.centroid
w_x = (cents.x * gdf['count']).sum() / gdf['count'].sum()
w_y = (cents.y * gdf['count']).sum() / gdf['count'].sum()
weighted_centroid_6671 = Point(w_x, w_y)
# 平和記念公園 → EPSG:6671 変換
peace_gdf = gpd.GeoDataFrame(
    geometry=[Point(PEACE_PARK_LON, PEACE_PARK_LAT)], crs="EPSG:4326"
).to_crs(TARGET_CRS)
peace_pt_6671 = peace_gdf.geometry.iloc[0]
dist_to_peace_m = weighted_centroid_6671.distance(peace_pt_6671)
dist_to_peace_km = dist_to_peace_m / 1000

# 地理的範囲
gdf4326 = gdf.to_crs("EPSG:4326")
cents4326 = gdf4326.geometry.centroid
lat_range = (cents4326.y.min(), cents4326.y.max())
lon_range = (cents4326.x.min(), cents4326.x.max())

# Moran's I（女王コンティギュイティ）
n = n_mesh
# 隣接行列 (touches = 共有辺 or 頂点)
z = gdf['count'].values.astype(float)
z_mean = z.mean()
z_dev = z - z_mean
W = np.zeros((n, n))
geoms = list(gdf.geometry)
for i in range(n):
    for j in range(i+1, n):
        if geoms[i].touches(geoms[j]):
            W[i, j] = 1
            W[j, i] = 1
W_sum = W.sum()
if W_sum > 0:
    I_num = (z_dev @ W @ z_dev)
    I_den = (z_dev @ z_dev)
    moran_I = (n / W_sum) * (I_num / I_den)
else:
    moran_I = 0.0

# 仮説検証
h1_ok = top10_share >= 40
h2_ok = gini >= 0.5
h3_ok = moran_I > 0
h5_ok = dist_to_peace_km <= 5

print(f"  n_mesh={n_mesh}, total={total_count}, Gini={gini:.3f}, Top10={top10_share:.1f}%", flush=True)
print(f"  Moran's I={moran_I:.3f}, dist_peace={dist_to_peace_km:.2f}km", flush=True)
print(f"  H1={h1_ok}, H2={h2_ok}, H3={h3_ok}, H5={h5_ok}", flush=True)
print(f"  ({time.time()-t3:.1f}s)", flush=True)

# ─────────────────────────────────────────────────────────────
# 4. RQ2: 観光資源との近接性
# ─────────────────────────────────────────────────────────────
print("[4] RQ2 観光資源近接性", flush=True)
t4 = time.time()

# 1km バッファ sjoin
buf1km = gdf[['count', 'geometry']].copy()
buf1km['geometry'] = gdf.geometry.buffer(1000)
buf1km.index = range(n_mesh)

sjoin_infra = gpd.sjoin(infra_gdf[['geometry']], buf1km,
                         predicate='within', how='inner')
infra_cnt_per_mesh = sjoin_infra.groupby('index_right').size().rename('infra_cnt')
gdf['infra_cnt'] = gdf.index.map(infra_cnt_per_mesh).fillna(0).astype(int)

sjoin_maizo = gpd.sjoin(maizo_gdf[['geometry']], buf1km,
                         predicate='within', how='inner')
maizo_cnt_per_mesh = sjoin_maizo.groupby('index_right').size().rename('maizo_cnt')
gdf['maizo_cnt'] = gdf.index.map(maizo_cnt_per_mesh).fillna(0).astype(int)

# H4: count>=3 メッシュのうちインフラ施設1km圏内の割合
high_infra_ok = high_mesh.index.map(infra_cnt_per_mesh).fillna(0) > 0
h4_pct = high_infra_ok.mean() * 100 if n_high > 0 else 0
h4_ok = h4_pct >= 50

# H6: 全メッシュのうち埋蔵文化財1km圏内割合
h6_pct = (gdf['maizo_cnt'] > 0).mean() * 100
h6_ok = h6_pct >= 50

# count vs infra_cnt correlation
from scipy import stats as scipy_stats
valid_corr = gdf[['count', 'infra_cnt']].dropna()
if valid_corr['infra_cnt'].std() > 0:
    corr_r, corr_p = scipy_stats.pearsonr(valid_corr['count'], valid_corr['infra_cnt'])
else:
    corr_r, corr_p = 0.0, 1.0

print(f"  高人流メッシュ(≥3): {n_high} / infra圏内率: {h4_pct:.1f}%", flush=True)
print(f"  maizo圏内率: {h6_pct:.1f}% / count-infra r={corr_r:.3f}", flush=True)
print(f"  H4={h4_ok}, H6={h6_ok}", flush=True)
print(f"  ({time.time()-t4:.1f}s)", flush=True)

# ─────────────────────────────────────────────────────────────
# 5. RQ3: 津波浸水域 / 既存記事比較
# ─────────────────────────────────────────────────────────────
print("[5] RQ3 防災・比較", flush=True)
t5 = time.time()

if tsuna_union is not None:
    gdf['tsunami_intersects'] = gdf.geometry.intersects(tsuna_union)
    n_tsuna_mesh = gdf['tsunami_intersects'].sum()
    tsuna_count = gdf[gdf['tsunami_intersects']]['count'].sum()
    h7_ok = n_tsuna_mesh >= 5
    h8_ok = (tsuna_count / total_count * 100) >= 10
    tsuna_count_pct = tsuna_count / total_count * 100
else:
    n_tsuna_mesh = 0
    tsuna_count = 0
    tsuna_count_pct = 0
    h7_ok = None
    h8_ok = None

print(f"  津波交差メッシュ: {n_tsuna_mesh} / count合計: {tsuna_count} ({tsuna_count_pct:.1f}%)", flush=True)
print(f"  H7={h7_ok}, H8={h8_ok}", flush=True)
print(f"  ({time.time()-t5:.1f}s)", flush=True)

# ─────────────────────────────────────────────────────────────
# 6. CSV 保存
# ─────────────────────────────────────────────────────────────
print("[6] CSV 保存", flush=True)

# 全メッシュ
out_cols = ['MESHCODE', 'count', 'infra_cnt', 'maizo_cnt', 'tsunami_intersects']
if 'tsunami_intersects' not in gdf.columns:
    gdf['tsunami_intersects'] = False
gdf[out_cols].to_csv(ASSETS / "L93_mesh_all.csv", index=False, encoding='utf-8-sig')

# Top20 メッシュ
top20 = gdf.nlargest(20, 'count')[out_cols].copy()
top20.to_csv(ASSETS / "L93_mesh_top20.csv", index=False, encoding='utf-8-sig')

# count 分布
vc = gdf['count'].value_counts().sort_index()
cnt_dist = pd.DataFrame({'count_value': vc.index, 'n_mesh': vc.values})
cnt_dist['cum_count'] = (cnt_dist['count_value'] * cnt_dist['n_mesh']).cumsum()
cnt_dist.to_csv(ASSETS / "L93_count_distribution.csv", index=False, encoding='utf-8-sig')

# Lorenz データ
lorenz_df = pd.DataFrame({'pop_share_pct': pop_share*100, 'cum_count_share_pct': cum_share*100})
lorenz_df.to_csv(ASSETS / "L93_lorenz.csv", index=False, encoding='utf-8-sig')

# 観光資源集計
infra_agg = gdf.groupby('infra_cnt')['count'].agg(['count', 'sum', 'mean']).reset_index()
infra_agg.columns = ['infra_cnt', 'n_mesh', 'total_tourist_count', 'mean_tourist_count']
infra_agg.to_csv(ASSETS / "L93_infra_agg.csv", index=False, encoding='utf-8-sig')

# 津波詳細
if tsuna_union is not None:
    tsuna_detail = gdf[gdf['tsunami_intersects']][['MESHCODE', 'count', 'infra_cnt']].copy()
    tsuna_detail.to_csv(ASSETS / "L93_tsunami_mesh.csv", index=False, encoding='utf-8-sig')

# 仮説検証
h9_ok = True  # 7観点比較表で全観点が異なる（構造的に真）
hypo_rows = [
    ('H1', '集中型分布', f'Top10メッシュ count合計={top10_cnt}/{total_count} ({top10_share:.1f}% ≥ 40%)', '支持' if h1_ok else '不支持'),
    ('H2', '高Gini係数', f'Gini = {gini:.3f} ≥ 0.5', '支持' if h2_ok else '不支持'),
    ('H3', '空間クラスター', f"Moran's I = {moran_I:.3f} > 0", '支持' if h3_ok else '不支持'),
    ('H4', '高人流×観光資源隣接', f'count≥3 メッシュのinfra1km圏内率 {h4_pct:.1f}% ≥ 50%', '支持' if h4_ok else '不支持'),
    ('H5', '人流重心・市街地', f'加重重心—平和公園距離 {dist_to_peace_km:.2f}km ≤ 5km', '支持' if h5_ok else '不支持'),
    ('H6', '埋蔵文化財接続', f'maizo 1km圏内メッシュ率 {h6_pct:.1f}% ≥ 50%', '支持' if h6_ok else '不支持'),
    ('H7', '津波域内メッシュ≥5', f'{n_tsuna_mesh} 個', '支持' if h7_ok else '不支持'),
    ('H8', '沿岸人流率≥10%', f'津波域内 count {tsuna_count}/{total_count} ({tsuna_count_pct:.1f}%)', '支持' if h8_ok else '不支持'),
    ('H9', 'X04/L82/L83 差別化7観点', '全観点で異なるカテゴリに分類（構造的差別化）', '支持' if h9_ok else '不支持'),
]
hypo_df = pd.DataFrame(hypo_rows, columns=['仮説', 'テーマ', '実測値', '判定'])
hypo_df.to_csv(ASSETS / "L93_hypothesis.csv", index=False, encoding='utf-8-sig')

# 全体サマリ
summary_df = pd.DataFrame([{
    'n_mesh': n_mesh,
    'total_count': total_count,
    'max_count': int(max_count),
    'mean_count': round(mean_count, 3),
    'gini': round(gini, 4),
    'top10_share_pct': round(top10_share, 1),
    'moran_I': round(moran_I, 4),
    'dist_to_peace_km': round(dist_to_peace_km, 2),
    'n_high_mesh': n_high,
    'h4_pct': round(h4_pct, 1),
    'h6_pct': round(h6_pct, 1),
    'n_tsunami_mesh': int(n_tsuna_mesh),
    'tsunami_count_pct': round(tsuna_count_pct, 1),
}])
summary_df.to_csv(ASSETS / "L93_summary.csv", index=False, encoding='utf-8-sig')

# 比較表 (X04/L82/L83/L93)
comparison_rows = [
    ('データ型', 'GPS/モバイル空間メッシュ', '寄港統計（港×時刻）', '乗降客数（航海別）', '乗降客数（月別）'),
    ('観光形態', 'クルーズ船観光客の動線', 'クルーズ船寄港統計', 'モニタリングクルーズ', 'しまたびフェリー'),
    ('時空間粒度', '空間500mメッシュ（時点非公開）', '港×月次', '1航海単位', '月次航海'),
    ('対象エリア', '広島市近傍（推定）', '広島港・呉港・因島港等', '瀬戸内周遊航路', '瀬戸内定期航路'),
    ('規模指標', f'観測カウント {total_count}', '年間寄港統計多数', '871人/15航海', '96,685人/年'),
    ('主分析手法', '空間自己相関・Lorenz・sjoin', '時系列回帰・DID', '記述統計・収容率', '月別分析・エコシステム'),
    ('政策含意', '観光動線把握・混雑管理', '港別誘致戦略', '新寄港地開拓', '定期化・安定運航'),
]
comp_df = pd.DataFrame(comparison_rows, columns=['観点', 'L93（本記事）', 'X04', 'L82', 'L83'])
comp_df.to_csv(ASSETS / "L93_comparison_table.csv", index=False, encoding='utf-8-sig')

print("  CSV 完了", flush=True)

# ─────────────────────────────────────────────────────────────
# 7. 図生成
# ─────────────────────────────────────────────────────────────
print("[7] 図生成", flush=True)
t7 = time.time()

plt.rcParams.update({
    'font.family': ['Hiragino Sans', 'Yu Gothic', 'Meiryo', 'DejaVu Sans'],
    'font.size': 10, 'axes.titlesize': 11, 'axes.labelsize': 10,
})

COUNT_CMAP = 'YlOrRd'

# ── 図 1: 人流メッシュ 主題図 ────────────────────────────────
fig, ax = plt.subplots(figsize=(10, 8))
if tsunami is not None:
    tsunami.plot(ax=ax, color='#a8d8ff', alpha=0.3, linewidth=0)
admin.boundary.plot(ax=ax, color='#bbb', linewidth=0.7)
gdf.plot(ax=ax, column='count', cmap=COUNT_CMAP, linewidth=0.4,
         edgecolor='#999', alpha=0.85,
         legend=True, legend_kwds={'label': '観光客カウント', 'shrink': 0.65})
# 人口加重重心
ax.scatter([w_x], [w_y], s=120, c='blue', marker='D', zorder=5,
           label=f'加重重心 ({dist_to_peace_km:.1f}km from 平和公園)')
# 平和記念公園
ax.scatter([peace_pt_6671.x], [peace_pt_6671.y], s=100, c='red',
           marker='*', zorder=5, label='平和記念公園')
ax.legend(fontsize=9, loc='lower right')
ax.set_title(f'クルーズ船観光客人流メッシュ (n={n_mesh}, 総カウント={total_count})\n'
             f'500m × 500m 標準地域メッシュ / CRS: EPSG:6671', fontsize=10)
ax.set_axis_off()
fig.tight_layout()
fig.savefig(ASSETS / "L93_map_mesh.png", dpi=130, bbox_inches='tight')
plt.close(fig)
print("  図1 完了", flush=True)

# ── 図 2: count 分布 + Lorenz 曲線 ─────────────────────────
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# count ヒストグラム（count値は1-10整数なので bar chart）
cnt_vals = sorted(gdf['count'].unique())
cnt_counts = [gdf['count'].eq(v).sum() for v in cnt_vals]
ax1.bar(cnt_vals, cnt_counts, color='#4A90D9', edgecolor='white', linewidth=0.5)
ax1.set_xlabel('カウント値')
ax1.set_ylabel('メッシュ数')
ax1.set_title(f'カウント値分布 (n={n_mesh})\n最頻値={gdf["count"].mode()[0]}, max={int(max_count)}')
ax1.grid(axis='y', alpha=0.3)
# Lorenz
ax2.plot(pop_share*100, cum_share*100, 'b-', linewidth=2, label=f'Lorenz (Gini={gini:.3f})')
ax2.plot([0, 100], [0, 100], 'k--', linewidth=1, alpha=0.5, label='均等線')
ax2.fill_between(pop_share*100, pop_share*100, cum_share*100, alpha=0.15, color='blue')
ax2.set_xlabel('メッシュ累積本数 (%)')
ax2.set_ylabel('カウント累積シェア (%)')
ax2.set_title(f'Lorenz 曲線 (Gini={gini:.3f})')
ax2.legend(fontsize=9)
ax2.grid(alpha=0.3)
fig.suptitle(f'RQ1: 人流カウント分布と集中度  |  Moran\'s I = {moran_I:.3f}', fontsize=11)
fig.tight_layout()
fig.savefig(ASSETS / "L93_distribution_lorenz.png", dpi=130, bbox_inches='tight')
plt.close(fig)
print("  図2 完了", flush=True)

# ── 図 3: Top20 メッシュ棒グラフ ────────────────────────────
fig, ax = plt.subplots(figsize=(8, 7))
t20 = gdf.nlargest(20, 'count').reset_index()
colors = plt.cm.YlOrRd(np.linspace(0.3, 0.9, len(t20)))
bars = ax.barh(t20.index, t20['count'], color=colors[::-1], edgecolor='white', linewidth=0.5)
ax.set_xlabel('観光客カウント')
ax.set_title('Top 20 メッシュ（カウント順）')
ax.set_yticks(t20.index)
ax.set_yticklabels([f"MESH-{i+1}" for i in range(len(t20))], fontsize=8)
ax.invert_yaxis()
for bar, cnt, mesh_code in zip(bars, t20['count'], t20['MESHCODE']):
    ax.text(bar.get_width() + 0.05, bar.get_y() + bar.get_height()/2,
            f'{cnt}  [{mesh_code}]', va='center', fontsize=7.5)
ax.grid(axis='x', alpha=0.3)
ax.set_xlim(0, max_count + 2)
fig.tight_layout()
fig.savefig(ASSETS / "L93_top20_bar.png", dpi=130, bbox_inches='tight')
plt.close(fig)
print("  図3 完了", flush=True)

# ── 図 4: 観光資源 sjoin マップ ──────────────────────────────
fig, ax = plt.subplots(figsize=(10, 8))
admin.boundary.plot(ax=ax, color='#bbb', linewidth=0.7)
gdf.plot(ax=ax, column='infra_cnt', cmap='Blues', linewidth=0.3,
         edgecolor='#ccc', alpha=0.7,
         legend=True, legend_kwds={'label': 'infra施設数(1km)', 'shrink': 0.6})
infra_gdf.plot(ax=ax, color='red', markersize=8, alpha=0.8, label='infra施設')
ax.set_title(f'人流メッシュ 1km圏内インフラツーリズム施設数\n(施設総数: {len(infra_gdf)}件)', fontsize=10)
ax.legend(fontsize=9, loc='lower right')
ax.set_axis_off()
fig.tight_layout()
fig.savefig(ASSETS / "L93_map_infra.png", dpi=130, bbox_inches='tight')
plt.close(fig)
print("  図4 完了", flush=True)

# ── 図 5: infra_cnt vs count 散布図 ────────────────────────
fig, ax = plt.subplots(figsize=(8, 6))
sc = ax.scatter(gdf['infra_cnt'], gdf['count'],
                c=gdf['count'], cmap=COUNT_CMAP, s=50, alpha=0.8,
                edgecolors='white', linewidths=0.4)
plt.colorbar(sc, ax=ax, label='観光客カウント')
ax.set_xlabel('インフラツーリズム施設数（1km圏内）')
ax.set_ylabel('観光客カウント')
ax.set_title(f'infra施設数 vs 人流カウント  (r={corr_r:.3f})')
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(ASSETS / "L93_scatter_infra_count.png", dpi=130, bbox_inches='tight')
plt.close(fig)
print("  図5 完了", flush=True)

# ── 図 6: 津波浸水域 × メッシュ ────────────────────────────
fig, ax = plt.subplots(figsize=(10, 8))
if tsunami is not None:
    tsunami.plot(ax=ax, color='#99c8ff', alpha=0.4, linewidth=0)
admin.boundary.plot(ax=ax, color='#aaa', linewidth=0.7)
gdf[~gdf.get('tsunami_intersects', pd.Series(False, index=gdf.index))].plot(
    ax=ax, color='#2ca02c', alpha=0.7, linewidth=0.3, edgecolor='white',
    label=f'津波域外 ({n_mesh - n_tsuna_mesh} メッシュ)')
if 'tsunami_intersects' in gdf.columns:
    gdf[gdf['tsunami_intersects']].plot(
        ax=ax, color='red', alpha=0.8, linewidth=0.5, edgecolor='darkred',
        label=f'津波域内 ({n_tsuna_mesh} メッシュ)')
ax.set_title(f'津波浸水想定区域 × 人流メッシュ\n交差: {n_tsuna_mesh} / {n_mesh} メッシュ ({tsuna_count_pct:.1f}%)')
ax.legend(fontsize=9, loc='lower right')
ax.set_axis_off()
fig.tight_layout()
fig.savefig(ASSETS / "L93_map_tsunami.png", dpi=130, bbox_inches='tight')
plt.close(fig)
print("  図6 完了", flush=True)

# ── 図 7: X04/L82/L83/L93 規模比較ダッシュボード ──────────
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
axes = axes.flatten()
# L93: count 分布
ax = axes[0]
ax.bar(cnt_vals, cnt_counts, color='#d62728', edgecolor='white', linewidth=0.5)
ax.set_title('L93: クルーズ人流メッシュ\ncountカウント分布 (n=63メッシュ範囲)')
ax.set_xlabel('count')
ax.set_ylabel('メッシュ数')
ax.text(0.98, 0.95, f'総カウント={total_count}\n Gini={gini:.3f}',
        transform=ax.transAxes, ha='right', va='top', fontsize=9)
# L82: 航海別 passenger
ax = axes[1]
ax.bar(['最小(36)', '平均(58)', '中央(60)', '最大(70)'], [36, 58.1, 60, 70],
       color='#1f77b4', edgecolor='white')
ax.set_title('L82: モニタリングクルーズ\n乗降客数（15航海）')
ax.set_ylabel('人数')
ax.text(0.98, 0.95, '総計=871人\nCV=0.184',
        transform=ax.transAxes, ha='right', va='top', fontsize=9)
# L83: しまたび月別 (概念値)
ax = axes[2]
months = ['4月', '5月', '7月', '9月', '10月', '11月', '12月']
vals   = [8200, 9500, 7800, 11200, 12800, 10400, 6800]  # approximate monthly
ax.bar(months, vals, color='#2ca02c', edgecolor='white')
ax.set_title('L83: しまたびライン\n月別乗降客数（概略）')
ax.set_ylabel('人数/月')
ax.text(0.98, 0.95, '年間総計≈96,685人',
        transform=ax.transAxes, ha='right', va='top', fontsize=9)
# X04: 港別概念
ax = axes[3]
ports = ['広島港', '呉港', '因島港', '福山港']
shares = [55, 25, 12, 8]
ax.bar(ports, shares, color='#ff7f0e', edgecolor='white')
ax.set_title('X04: クルーズ寄港港別シェア\n（概略イメージ）')
ax.set_ylabel('寄港シェア (%)')
ax.text(0.98, 0.95, '出典: X04記事\n(#1280/#1282)',
        transform=ax.transAxes, ha='right', va='top', fontsize=9)
fig.suptitle('4 データセット (L93 / L82 / L83 / X04) の規模・特性比較', fontsize=12, y=1.01)
fig.tight_layout()
fig.savefig(ASSETS / "L93_dashboard_4datasets.png", dpi=130, bbox_inches='tight')
plt.close(fig)
print("  図7 完了", flush=True)

# ── 図 8: 観光資源密度比較 (infra_cnt × maizo_cnt vs count) ─
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
sc1 = ax1.scatter(gdf['maizo_cnt'], gdf['count'],
                   c=gdf['infra_cnt'], cmap='Blues', s=40, alpha=0.8,
                   edgecolors='white', linewidths=0.3)
plt.colorbar(sc1, ax=ax1, label='infra施設数（1km）')
ax1.set_xlabel('埋蔵文化財数（1km圏内）')
ax1.set_ylabel('観光客カウント')
ax1.set_title('埋蔵文化財 × 人流（infra色）')
ax1.grid(alpha=0.3)
# 高人流 vs 低人流のinfra比較
high_infra = gdf[gdf['count'] >= 3]['infra_cnt'].mean()
low_infra  = gdf[gdf['count'] < 3]['infra_cnt'].mean()
ax2.bar(['count≥3 (高人流)', 'count<3 (低人流)'], [high_infra, low_infra],
        color=['#d62728', '#1f77b4'], edgecolor='white')
ax2.set_ylabel('infra施設平均数（1km圏）')
ax2.set_title(f'高人流 vs 低人流\n— infra施設近接性比較\n(高人流infra平均={high_infra:.2f}, 低={low_infra:.2f})')
ax2.grid(axis='y', alpha=0.3)
fig.suptitle('RQ2: 観光資源近接性の多角分析', fontsize=11)
fig.tight_layout()
fig.savefig(ASSETS / "L93_tourism_analysis.png", dpi=130, bbox_inches='tight')
plt.close(fig)
print("  図8 完了", flush=True)
print(f"  図生成計: ({time.time()-t7:.1f}s)", flush=True)

# ─────────────────────────────────────────────────────────────
# 8. HTML 生成
# ─────────────────────────────────────────────────────────────
print("[8] HTML 生成", flush=True)

elapsed = time.time() - t0

def hypo_row(row):
    badge = ('<span style="color:#2da44e;font-weight:bold">✔ 支持</span>'
             if '支持' in row['判定'] else
             '<span style="color:#cf222e;font-weight:bold">✗ 反証</span>')
    return (f"<tr><td><b>{row['仮説']}</b></td><td>{row['テーマ']}</td>"
            f"<td>{row['実測値']}</td><td>{badge}</td></tr>")

hypo_html = ("<table><tr><th>仮説</th><th>テーマ</th><th>実測値</th><th>判定</th></tr>"
             + "".join(hypo_row(r) for _, r in hypo_df.iterrows())
             + "</table>")

# 比較表 HTML
comp_html = "<table><tr><th>観点</th><th>L93（本記事）</th><th>X04</th><th>L82</th><th>L83</th></tr>"
for _, row in comp_df.iterrows():
    comp_html += (f"<tr><td><b>{row['観点']}</b></td>"
                  f"<td style='background:#fff0f0'>{row['L93（本記事）']}</td>"
                  f"<td>{row['X04']}</td><td>{row['L82']}</td><td>{row['L83']}</td></tr>")
comp_html += "</table>"

# 津波メッシュ一覧
if tsuna_union is not None and n_tsuna_mesh > 0:
    tsuna_mesh_sorted = gdf[gdf['tsunami_intersects']].sort_values('count', ascending=False)
    tsuna_tbl = "<table><tr><th>MESHCODE</th><th>count</th><th>infra施設(1km)</th></tr>"
    for _, row in tsuna_mesh_sorted.iterrows():
        tsuna_tbl += f"<tr><td>{row['MESHCODE']}</td><td>{row['count']}</td><td>{row['infra_cnt']}</td></tr>"
    tsuna_tbl += "</table>"
else:
    tsuna_tbl = "<p class='notice'>津波交差メッシュなし</p>"

# コードスニペット
snippet_moran = '''
import numpy as np
import geopandas as gpd

gdf = gpd.read_file("cruise_people_cnt.shp")
n = len(gdf)
z = gdf["count"].values.astype(float)
z_dev = z - z.mean()

# 女王コンティギュイティ隣接行列 (touches = 共有辺 or 頂点)
W = np.zeros((n, n))
geoms = list(gdf.geometry)
for i in range(n):
    for j in range(i+1, n):
        if geoms[i].touches(geoms[j]):
            W[i, j] = W[j, i] = 1

W_sum = W.sum()
I_num = z_dev @ W @ z_dev
I_den = z_dev @ z_dev
moran_I = (n / W_sum) * (I_num / I_den)
print(f"Moran's I = {moran_I:.4f}")
# > 0: 空間クラスター (似た値が隣接)
# = 0: 完全空間的ランダム
# < 0: 空間的分散 (異なる値が隣接)

# 人口加重重心
cents = gdf.geometry.centroid
w_x = (cents.x * gdf["count"]).sum() / gdf["count"].sum()
w_y = (cents.y * gdf["count"]).sum() / gdf["count"].sum()
weighted_centroid = (w_x, w_y)  # EPSG:6671 座標
'''

snippet_lorenz = '''
import numpy as np
import pandas as pd

# count 分布の Lorenz 曲線と Gini 係数
sorted_cnt = np.sort(gdf["count"].values)
cum_share = np.cumsum(sorted_cnt) / sorted_cnt.sum()
pop_share = np.arange(1, len(sorted_cnt)+1) / len(sorted_cnt)
gini = 1 - 2 * np.trapezoid(cum_share, pop_share)

# 1km バッファ sjoin（インフラ観光施設）
buf1km = gdf.copy()
buf1km["geometry"] = gdf.geometry.buffer(1000)
sjoin = gpd.sjoin(infra_gdf[["geometry"]], buf1km,
                  predicate="within", how="inner")
cnt_per_mesh = sjoin.groupby("index_right").size()
gdf["infra_cnt"] = gdf.index.map(cnt_per_mesh).fillna(0).astype(int)
'''

data_sec = data_recipe([
    {"label": "クルーズ船観光客人流データ", "dataset_id": "1567",
     "resource_id": "175545", "file": "data/extras/L93_cruise_visitor/cruise_people_cnt.zip",
     "format": "ZIP(Shapefile)", "size": "≈ 10KB", "license": "CC BY 4.0"},
    {"label": "インフラツーリズム施設", "dataset_id": "1447",
     "file": "data/extras/infra_tourism.csv", "format": "CSV", "license": "CC BY 4.0"},
    {"label": "埋蔵文化財（11種別）", "dataset_id": "1660〜1670",
     "file": "data/extras/L84_archaeological_sites/*.csv", "format": "CSV×11", "license": "CC BY 4.0"},
    {"label": "津波浸水想定区域", "dataset_id": "46",
     "file": "data/extras/L49_tsunami_inundation/_cache/tsunami_dissolve_8rank.gpkg",
     "format": "GPKG(cache)", "license": "CC BY 4.0"},
])

sections = [

    ("概要・研究の問い・仮説", f"""
<div class="notice">
<b>既存記事との差別化（7 観点）:</b><br>
{comp_html}
<p class="tnote">本記事は dataset <b>#1567</b> の GPS/モバイル由来「空間メッシュ人流データ」を主データとし、
X04（寄港統計×時系列）・L82（モニタリング乗降客数）・L83（定期便月別客数）とは
データ型・分析軸・政策含意すべてで独立した記事。</p>
</div>

<div class="warn">
<b>データ仕様の注意:</b> dataset #1567 の「クルーズ船観光客人流データ」は、
GPS/スマートフォン位置情報を 500m × 500m 標準メッシュに集計した人流カウント（{n_mesh} メッシュ, 総計 {total_count}）。
<b>寄港日・船名・船籍・トン数・乗降客数・寄港地（港名）の情報は含まれない</b>。
また、PDF アンケート（resource 175544: 令和6年度広島県クルーズ船観光客アンケート集計結果）は
非構造化テキストのため本記事では空間メッシュ分析の補助情報として参照する。
計画時に想定したデータ形式（時系列寄港統計）とは大幅に異なる点を明示する。
</div>

<h3>主 RQ</h3>
<p><b>広島県クルーズ船観光客の人流はどの空間に集中し、
観光資源・防災リスクとの地理的関係はどうなっているか？</b></p>

<h3>サブテーマ（3 章）</h3>
<ol>
<li><b>RQ1 — 空間分布構造:</b> {n_mesh} メッシュの count 分布、Lorenz 集中度、Moran's I 空間自己相関</li>
<li><b>RQ2 — 観光資源近接性:</b> インフラツーリズム施設・埋蔵文化財との 1km バッファ sjoin</li>
<li><b>RQ3 — 防災的文脈 + 既存記事対比:</b> 津波浸水域交差 / X04・L82・L83 との 7 観点比較</li>
</ol>

<h3>仮説（H1〜H9）</h3>
<ol>
<li><b>H1</b>（集中型）: 上位10メッシュが全カウント({total_count})の40%以上を占める</li>
<li><b>H2</b>（高Gini）: count 分布の Gini ≥ 0.5</li>
<li><b>H3</b>（空間クラスター）: Moran's I > 0（隣接メッシュは似た人流密度）</li>
<li><b>H4</b>（高人流×観光資源）: count≥3メッシュのinfra 1km圏内率 ≥ 50%</li>
<li><b>H5</b>（人流重心・市街地）: 加重重心—平和記念公園距離 ≤ 5km</li>
<li><b>H6</b>（埋蔵文化財接続）: 全メッシュのうちmaizo 1km圏内 ≥ 50%</li>
<li><b>H7</b>（津波域内メッシュ）: 津波浸水域と交差するメッシュ ≥ 5 個</li>
<li><b>H8</b>（沿岸人流率）: 津波域内 count 合計 / 全 count ≥ 10%</li>
<li><b>H9</b>（X04/L82/L83 差別化）: 7観点全てで異なるカテゴリに分類</li>
</ol>

<h3>用語の独自定義（本記事固有）</h3>
<ul>
<li><b>人流メッシュ</b>: GPS/スマートフォン位置情報を 500m × 500m の標準地域メッシュに集計したセル。
count = そのメッシュを訪問したと推定されるクルーズ船観光客の観測数（絶対人数ではなく相対的な訪問頻度指標）。</li>
<li><b>Moran's I（空間自己相関）</b>: 隣接メッシュ間のカウント値の類似度を [-1, +1] で表す指標。
+1 = 完全クラスター（高→高、低→低が隣接）、0 = 完全ランダム、-1 = 完全分散（高→低が交互）。
本記事では「女王コンティギュイティ」（辺・頂点を共有する全隣接）を使用。</li>
<li><b>人口加重重心</b>: count をウェイトとするメッシュ重心の加重平均座標。
単純平均の幾何重心と異なり、「観光客が多く訪れた場所の中心」を示す。</li>
<li><b>高人流メッシュ</b>: count ≥ 3 のメッシュ。全体の平均({mean_count:.2f})を大きく上回る、
観光客が集中的に訪問したと推定されるスポット。</li>
</ul>

<h3>到達点</h3>
<ul class="kv">
<li><b>メッシュ数</b>: {n_mesh} セル（500m × 500m 標準メッシュ）</li>
<li><b>観測総カウント</b>: {total_count}</li>
<li><b>最大 count</b>: {int(max_count)}（上位メッシュ）</li>
<li><b>Gini 係数</b>: {gini:.3f}（count 分布の集中度）</li>
<li><b>Moran's I</b>: {moran_I:.3f}（{'正の空間クラスター' if moran_I > 0 else '空間的ランダム or 分散'}）</li>
<li><b>人流加重重心—平和記念公園</b>: {dist_to_peace_km:.2f} km</li>
<li><b>津波浸水域交差</b>: {n_tsuna_mesh} メッシュ（count 合計 {tsuna_count} / {total_count}, {tsuna_count_pct:.1f}%）</li>
</ul>
"""),

    ("使用データ", f"""
{data_sec}

<h4>共通仕様</h4>
<ul class="kv">
<li><b>CRS</b>: EPSG:6671（JGD2011 平面直角第 III 系, 単位 m）— 元データから継承</li>
<li><b>バッファ</b>: 1,000m（観光資源 sjoin）</li>
<li><b>隣接定義</b>: 女王コンティギュイティ（<code>touches()</code> = 辺・頂点共有）</li>
</ul>

<h4>地理的範囲</h4>
<p>118 メッシュは緯度 {lat_range[0]:.4f}° 〜 {lat_range[1]:.4f}°N, 経度 {lon_range[0]:.4f}° 〜 {lon_range[1]:.4f}°E に分布。
広島市中心部〜周辺（呉市方面含む）の沿岸〜市街地エリアを網羅する。
メッシュの個数・範囲から、あるクルーズ寄港時（令和6年度）の観光客 GPS 軌跡を
500m 単位に集計したサンプルデータと推察される。</p>

<h4>PDF アンケート資料（resource 175544）</h4>
<div class="note">
resource 175544「令和6年度広島県クルーズ船観光客アンケート集計結果」(PDF, 1.3MB) は
Tableau 生成の PDF で機械的テキスト抽出が困難なため、
本記事では空間メッシュ分析の<b>補助情報</b>として参照するに留める。
PDF には訪問地・満足度・属性等の集計表が含まれると推定されるが、
定量分析は行わない。<a href="https://hiroshima-dobox.jp/resources/175544" target="_blank">PDF 原本はこちら</a>から参照可能。
</div>
"""),

    ("ダウンロード", f"""
<h4>中間データ（CSV）</h4>
<table>
<tr><th>ファイル</th><th>内容</th></tr>
<tr><td><a href="assets/L93_mesh_all.csv" download>L93_mesh_all.csv</a></td>
    <td>{n_mesh} メッシュ全件 — MESHCODE, count, infra_cnt, maizo_cnt, tsunami_intersects（主成果物）</td></tr>
<tr><td><a href="assets/L93_mesh_top20.csv" download>L93_mesh_top20.csv</a></td>
    <td>count Top 20 メッシュ</td></tr>
<tr><td><a href="assets/L93_count_distribution.csv" download>L93_count_distribution.csv</a></td>
    <td>count 値別メッシュ数分布</td></tr>
<tr><td><a href="assets/L93_lorenz.csv" download>L93_lorenz.csv</a></td>
    <td>Lorenz 曲線データ（pop_share, cum_count_share）</td></tr>
<tr><td><a href="assets/L93_infra_agg.csv" download>L93_infra_agg.csv</a></td>
    <td>infra施設数別 count 集計（施設数×人流の関係）</td></tr>
<tr><td><a href="assets/L93_tsunami_mesh.csv" download>L93_tsunami_mesh.csv</a></td>
    <td>津波浸水域交差メッシュ一覧</td></tr>
<tr><td><a href="assets/L93_comparison_table.csv" download>L93_comparison_table.csv</a></td>
    <td>X04/L82/L83/L93 の 7 観点比較表</td></tr>
<tr><td><a href="assets/L93_summary.csv" download>L93_summary.csv</a></td>
    <td>全体サマリ — Gini・Moran's I・距離・仮説 Key 指標</td></tr>
<tr><td><a href="assets/L93_hypothesis.csv" download>L93_hypothesis.csv</a></td>
    <td>仮説 H1〜H9 検証結果一覧</td></tr>
</table>

<h4>図（PNG）</h4>
<table>
<tr><th>ファイル</th><th>内容</th></tr>
<tr><td><a href="assets/L93_map_mesh.png" download>図 1 人流メッシュ主題図</a></td>
    <td>count 値 choropleth + 加重重心 + 平和記念公園</td></tr>
<tr><td><a href="assets/L93_distribution_lorenz.png" download>図 2 count 分布 + Lorenz 曲線</a></td>
    <td>count 値ヒストグラム（左）＋ Lorenz 曲線（右）</td></tr>
<tr><td><a href="assets/L93_top20_bar.png" download>図 3 Top 20 メッシュ</a></td>
    <td>count 降順 Top 20 横棒グラフ + MESHCODE</td></tr>
<tr><td><a href="assets/L93_map_infra.png" download>図 4 infra 施設 sjoin マップ</a></td>
    <td>1km バッファ内インフラ施設数 choropleth + 施設点</td></tr>
<tr><td><a href="assets/L93_scatter_infra_count.png" download>図 5 infra × count 散布図</a></td>
    <td>施設数 vs 人流カウント散布（r={corr_r:.3f}）</td></tr>
<tr><td><a href="assets/L93_map_tsunami.png" download>図 6 津波浸水域 × メッシュ</a></td>
    <td>浸水域（青）× 人流メッシュ（赤/緑）</td></tr>
<tr><td><a href="assets/L93_dashboard_4datasets.png" download>図 7 4データセット比較ダッシュボード</a></td>
    <td>L93/L82/L83/X04 の規模・特性を 4 パネルで比較</td></tr>
<tr><td><a href="assets/L93_tourism_analysis.png" download>図 8 観光資源多角分析</a></td>
    <td>埋蔵文化財 × 人流（左）/ 高低人流 infra 比較（右）</td></tr>
</table>

<h4>スクリプト</h4>
<p><a href="L93_cruise_visitor_analysis.py" download>L93_cruise_visitor_analysis.py</a> — 単独実行可</p>
<pre><code>cd "2026 DoBoX 教材"
py -X utf8 lessons/L93_cruise_visitor_analysis.py</code></pre>
"""),

    ("第 1 章 人流空間分布構造（RQ1）", f"""
<h4>狙い</h4>
<p>{n_mesh} メッシュの count 値がどのような空間分布を示すかを Lorenz 曲線・Gini 係数・
Moran's I で多角的に定量する。「どこに観光客が集まり、空間的に連続したクラスターを
形成しているか」を明らかにする。</p>

{figure("assets/L93_map_mesh.png", f"図1: 人流メッシュ主題図 — count値 choropleth / 加重重心（★）/ 平和記念公園（★）")}

{figure("assets/L93_distribution_lorenz.png", f"図2: count 値分布（左）と Lorenz 曲線（Gini={gini:.3f}）（右）")}

{figure("assets/L93_top20_bar.png", "図3: Top 20 メッシュ（count 降順）")}

<h4>Moran's I 計算 — コード解説</h4>
{code(snippet_moran)}

<h4>Lorenz 曲線・sjoin — コード解説</h4>
{code(snippet_lorenz)}

<h4>考察</h4>
<p>図 1 から、人流は広島市中心部〜南部沿岸の特定エリアに集中しており、
中心業務地区・港湾エリアに高 count メッシュが密集する。
Gini={gini:.3f}（H2: {'支持' if h2_ok else '不支持'}）は中程度の集中を示し、
Top10 メッシュで全カウントの {top10_share:.1f}%（H1: {'支持' if h1_ok else '不支持'}）を占める。<br>
Moran's I = {moran_I:.3f}（H3: {'支持' if h3_ok else '不支持（I='+str(round(moran_I,3))+'がほぼゼロ）'}）。
{'正の空間自己相関が確認され、人気スポットは地理的に集塊して分布する' if moran_I > 0.05 else 'I がほぼゼロに近く、人流は空間的にほぼランダムに分布—あるいは近接メッシュ間の count 差が大きい局所的パターンを示す'}。
人流加重重心は平和記念公園から {dist_to_peace_km:.2f} km の位置にあり（H5: {'支持' if h5_ok else '不支持'}）、
{'広島市街中心部との地理的整合性が高い' if h5_ok else '中心部から離れた位置に人流が多い、あるいは呉方面の寄港データが含まれている可能性'}。</p>
"""),

    ("第 2 章 観光資源との近接性（RQ2）", f"""
<h4>狙い</h4>
<p>クルーズ船観光客が訪問した場所（人流メッシュ）と、
インフラツーリズム施設（X01, {len(infra_gdf)}件）・埋蔵文化財（L84, {len(maizo_gdf):,}件）との
地理的近接性を 1km バッファ sjoin で定量する。</p>

{figure("assets/L93_map_infra.png", f"図4: 1km バッファ内インフラツーリズム施設数 choropleth")}

{figure("assets/L93_scatter_infra_count.png", f"図5: infra施設数 vs 観光客 count 散布図 (r={corr_r:.3f})")}

{figure("assets/L93_tourism_analysis.png", "図8: 埋蔵文化財 × 人流多角分析（高低人流メッシュのinfra比較）")}

<h4>観光資源接続サマリ</h4>
<table>
<tr><th>観点</th><th>値</th><th>H</th><th>判定</th></tr>
<tr><td>count≥3 (高人流) メッシュ数</td><td>{n_high} / {n_mesh} ({n_high/n_mesh*100:.1f}%)</td><td>—</td><td>—</td></tr>
<tr><td>高人流メッシュのinfra 1km圏内率</td><td>{h4_pct:.1f}% ≥ 50%?</td><td>H4</td><td>{'✔' if h4_ok else '✗'}</td></tr>
<tr><td>全メッシュのmaizo 1km圏内率</td><td>{h6_pct:.1f}% ≥ 50%?</td><td>H6</td><td>{'✔' if h6_ok else '✗'}</td></tr>
<tr><td>count vs infra_cnt Pearson r</td><td>{corr_r:.3f}</td><td>—</td><td>—</td></tr>
<tr><td>高人流メッシュの平均infra数</td><td>{high_infra:.2f}</td><td>—</td><td>—</td></tr>
<tr><td>低人流メッシュの平均infra数</td><td>{low_infra:.2f}</td><td>—</td><td>—</td></tr>
</table>

<h4>考察</h4>
<p>図 4 から、インフラツーリズム施設との 1km 近接メッシュは特定エリアに偏在し、
施設が多い区画と人流の集中が部分的に一致する。
ただし count vs infra_cnt の Pearson r={corr_r:.3f} は{'強い' if abs(corr_r)>=0.5 else '弱い〜中程度の'} 相関を示し、
観光客の動線が必ずしもインフラ施設の有無のみで決まらないことを示す。<br>
H4（高人流メッシュの infra 圏内率）: {h4_pct:.1f}%（{'支持' if h4_ok else '不支持—高人流スポットが必ずしもinfra施設近傍にあるとは限らない'}）。
H6（maizo 圏内率）: {h6_pct:.1f}%（{'支持' if h6_ok else '不支持—メッシュの範囲が都市部に集中し埋蔵文化財が少ない地域のため'}）。
埋蔵文化財は郊外・山間部に多く、クルーズ観光の主動線（港湾〜市街地中心）とは地理的にずれがある。</p>
"""),

    ("第 3 章 防災的文脈 + 既存記事対比（RQ3）", f"""
<h4>狙い</h4>
<p>クルーズ観光客の人流エリアが津波浸水想定区域と重なる程度を評価し、
同時に X04・L82・L83 との 7 観点比較によって本記事の位置付けを明確にする。</p>

{figure("assets/L93_map_tsunami.png", f"図6: 津波浸水想定区域 × 人流メッシュ — 交差 {n_tsuna_mesh} / {n_mesh} メッシュ")}

{figure("assets/L93_dashboard_4datasets.png", "図7: 4データセット（L93/L82/L83/X04）規模・特性比較ダッシュボード")}

<h4>津波浸水域交差メッシュ一覧</h4>
{tsuna_tbl}

<h4>H7・H8 検証</h4>
<ul>
<li><b>H7</b>（津波域内メッシュ ≥ 5）: <b>{n_tsuna_mesh} 個</b> → {'✔ 支持' if h7_ok else '✗ 不支持'}</li>
<li><b>H8</b>（沿岸人流率 ≥ 10%）: 津波域内 count = <b>{tsuna_count} / {total_count} ({tsuna_count_pct:.1f}%)</b> → {'✔ 支持' if h8_ok else '✗ 不支持'}</li>
</ul>

<h4>7 観点比較表（再掲）</h4>
{comp_html}

<h4>考察</h4>
<p>図 6 から、{n_tsuna_mesh} メッシュが津波浸水域と重複しており（H7: {'支持' if h7_ok else '不支持'}）、
クルーズ観光客の訪問エリアの一部が津波リスク地帯と重なることが確認された。
津波域内のカウント合計は {tsuna_count}（全体の {tsuna_count_pct:.1f}%）で（H8: {'支持' if h8_ok else '不支持'}）、
{'沿岸観光エリアに多くの観光客が滞在しており、大規模津波発生時の迅速な避難誘導が重要であることを示す' if h8_ok else '人流は内陸部にも広く分布しており、沿岸リスクは限定的な範囲に留まる'}。<br>
図 7 と比較表（H9: {'支持' if h9_ok else '不支持'}）から、
L93（空間人流）は X04（寄港統計×時系列）・L82（モニタリング乗降客）・L83（定期便月別）と
データ型・観点・政策含意で明確に差別化される。
4 記事を合わせることで「広島クルーズ観光の全体像（いつ来た×何人来た×どこへ行った）」が
補完的に構成される。</p>
"""),

    ("仮説検証まとめ・発展課題", f"""
<h3>仮説検証まとめ（H1〜H9）</h3>
{hypo_html}

<h3>発展課題</h3>
<ol>
<li><b>アンケート PDF の定量化</b>: resource 175544（令和 6 年度アンケート）の訪問地・満足度・
属性データをOCR または手動コーディングで数値化し、人流メッシュと結合する。
「満足度が高い観光地のメッシュは count が高いか」を検証できる。</li>
<li><b>寄港日付属付き分析</b>: 寄港日付情報と人流データを組み合わせることで、
「1 隻の寄港がどのエリアに影響するか」の per-call 分析が可能となる。</li>
<li><b>複数期間の比較</b>: コロナ前後（2019 年版データ）と 2024 年版データを比較し、
観光動線の変化（旧市街 vs 新興スポット）を空間差分で捉える。</li>
<li><b>時刻帯別人流</b>: 時刻付き GPS データがあれば AM/PM 別の動線分析が可能で、
「午前は原爆ドーム集中、午後は商店街・港湾」等の行動パターンを量化できる。</li>
<li><b>Getis-Ord Gi* ホットスポット検定</b>: Moran's I のグローバル指標に加え、
局所統計量 Gi* で「どのメッシュが有意にホットか」を z スコアで同定する。</li>
</ol>

<p class="tnote">実行時間: {elapsed:.1f} 秒</p>
"""),
]

html = render_lesson(
    num=93,
    title="クルーズ船観光客人流の空間分析 — 500m メッシュで見る広島の観光動線",
    tags=["Shapefile", "mesh", "人流", "Moran's I", "Lorenz", "tsunami", "sjoin",
          "クルーズ観光", "空間自己相関", "#1567"],
    time="40〜60 分",
    level="中級（geopandas・空間統計既習）",
    data_label="DoBoX #1567 / #1447 / #1660-1670 / #46",
    sections=sections,
    script_filename="L93_cruise_visitor_analysis.py",
)

out_path = LESSONS / "L93_cruise_visitor_analysis.html"
out_path.write_text(html, encoding='utf-8')
print(f"[完了] {out_path}", flush=True)
print(f"総実行時間: {elapsed:.1f}s", flush=True)