# -*- coding: utf-8 -*- """L82 せとうちモニタークルーズ実施結果 単独 3 研究例分析 — 観光社会実験 2022年12月15航海の航路設計・乗船客特性・観光商品評価 カバー宣言: 本記事は DoBoX のシリーズ「せとうちモニタークルーズ実施結果」 1 件 (dataset_id = 1280) を 単独で取り上げ、 広島県が 2022 年 12 月に実施した試験運航 7 航路 × 9 運航日 = 15 航海を 3 つの独立した研究角度 (RQ1 / RQ2 / RQ3) で並列に分析する。 「モニタークルーズ = 観光の社会実験」 を、 航路設計 / 乗船客特性 / 観光商品評価 の 3 軸から立体的に描く。 X04 改訂版 (2026-05) との重複回避: X04 = 同 dataset を「広島県沿岸クルーズ需要の港別・季節別構造分析」 で 扱い、 集中度 (Gini/HHI/CR), 港別年間需要, 季節係数, ルート 4 タイプ, stacked bar の 9 図 12 表で読み解いた。 RQ3 は「ルート設計」 を 名乗ったが実態は新規/既存比 + 4 タイプ分類 + 観光客数バー で完結。 L82 = 同 dataset を「観光社会実験のプロダクト設計」として 深掘り。 X04 が触れなかった (i) 航路の地理的距離・所要時間・寄港間隔 (Haversine)、 (ii) 既存しまたびとのJaccard 類似度とトポロジー (戻り型/片道型)、 (iii) 観光資源 (ビューポイント) の密度、 (iv) 1 航海客数の運航回数効果 (1 回目 vs 2 回目)、 (v) 収容率推定 + 採算性試算 + 定期化適性スコア に集中する。 X04 の指標は再計算しない (Gini/HHI/CR は計算しないし、 ルート 4 タイプ分類も呼ばない)。 「モニタークルーズ」 とは: 広島県が観光振興・新規ルート開拓の試験運航 (= 社会実験)として 限定日数で実施するクルーズ事業。 2022 年 12 月の本データは 7 航路 × 9 運航日 = 15 航海を瀬戸内海で運航し、 観光客アンケートと利用データを取得して定期航路化の判断材料を 集めることを目的とした。 既存定期航路「瀬戸内しまたびライン」 (#1282, X04 主データ) との比較で、 新規寄港地の集客力・既存資源との 接続性を評価する位置付け。 本記事は X04 (因果推論不可で記述転換) と厳密に区別: X04 = #1280 + #1282 を結合し、 港別・季節別の市場構造を分析。 中心はしまたびライン 14 港 × 9 ヶ月の月次データ。 L82 = #1280 のみを単独深掘り。 中心はモニタークルーズ 15 航海 × 23 寄港地のプロダクト設計分析。 しまたびライン (#1282) は RQ1 で航路重複の参照と RQ3 で月平均利用の対比でのみ従属参照。 研究の問い (3 RQ): RQ1 (主研究): モニタークルーズ 7 航路の航路設計 — ルート構造はどう描けるか? 7 航路 × 23 寄港地 × 緯度経度 × 区間情報 (新規/既存) を、 本記事独自の派生指標 (Haversine 総距離 / 平均区間距離 / 観光資源密度 / 始発・終着の組み合わせ / 既存しまたび Jaccard 類似度) で 7 軸記述する。 H1 = 7 航路の総距離は 80-120 km の狭い帯に 集中仮説 (max/min ≤ 1.5)、 H2 = 同一区間を東向きと西向きで運航する 東西ペアが 2 ペア以上存在仮説。 RQ2 (副研究 1): 1 航海あたり乗船客数の運航回数効果はあるか? 15 航海 × 観光客数を、 (i) 運航日順序 (第 1 回 vs 第 2 回 vs 第 3 回), (ii) 平日 vs 週末, (iii) ルートタイプ (戻り型 vs 片道型) で 3 角度集計し、 「2 回目以降の客数が減少するか?」 を検証する。 H3 = 運航 2 回目の 客数は 1 回目より15% 以上減仮説 (= 物珍しさ効果の減衰)、 H4 = 客数の標準偏差 / 平均 (CV) は0.20 以下仮説 (= 商品として 客数が安定)。 RQ3 (副研究 2): モニタークルーズの観光商品としての評価はどうか? (i) 収容率推定 (1 航海最大客数 = 70 人帯を満席仮定)、 (ii) 1 航海当たり推定収益 (仮定単価 5,000 円/人)、 (iii) 定期化適性スコア (本記事独自指標 = 平均客数 × 新規性 ボーナス / 距離係数)、 (iv) 既存しまたび月平均との利用比較を 4 軸で計算する。 H5 = 7 航路中2 航路以上が定期化適性スコア 上位 (= 高客数 + 高新規性) 仮説。 仮説 (5): H1 (RQ1, 総距離 80-120 km 帯集中): モニタークルーズ 7 航路の Haversine 総距離は80 〜 120 km の狭い帯に集中し、 max/min 比は 1.5 倍以下仮説 (= 「半日 〜 1 日かけて運航する片道商品」 として距離規格が均質化)。 H2 (RQ1, 東西ペア構造 ≥ 2 ペア): 7 航路中、 同一区間を東向き運航と 西向き運航で東西ペアとして組み立てている航路ペアが 2 ペア以上存在する仮説。 「乗船客が東向きまたは西向き、 好みの方向で楽しめる」 双方向商品設計の量的証拠。 H3 (RQ2, 第 2 運航客数減 ≥ 15%): 同一航路の第 2 回運航の観光客数は 第 1 回より15% 以上少ない仮説。 「物珍しさ効果」 の減衰、 第 1 回の 人気航路の方が話題性が高いことの量的証拠。 H4 (RQ2, 客数 CV ≤ 0.20): 15 航海の観光客数の変動係数 (= 標準偏差 / 平均)は0.20 以下仮説。 試験運航は商品として客数が安定し、 「常に 50-70 人入る」 規模感の証拠。 H5 (RQ3, 定期化適性 上位 2): 定期化適性スコア上位 2 航路は、 既存しまたび月平均利用 (~7,500 人/月) を 1 航海あたり客数で換算 して比較しても遜色ない規模を持つ仮説。 「観光社会実験は 定期化候補を 2 航路特定できた」 という成果の量的評価。 要件 S 準拠 (1 分以内完走): - 全データ 133 行 / 23 港 / 7 航路 → 超軽量 - Haversine は numpy ベクトル化で全寄港地ペア < 1 秒 - geopandas 投影変換は EPSG:4326 → 6671 を 1 度のみ - 図は plt.close('all') で確実に解放 - 想定完走時間 ~10-30 秒 要件 T 準拠 (位置情報あり = 地図必須): - RQ1: 7 航路ルートマップ (新規 = 赤線 / 既存 = 青線) + 観光資源バブル - RQ2: 客数 small multiples (運航日 × 航路の 15 航海地図) - RQ3: 定期化適性スコア地理可視化 (航路色 = スコア高低) 要件 Q 準拠: 図 7+ / 表 11+ (3 RQ × 多角度) データ仕様: - dataset 1280: せとうちモニタークルーズ実施結果 (CSV 3 件) - resource 39782: 航路 × 寄港地 × 運航日 (133 行 × 11 列, 17 KB) - resource 39781: 寄港地マスタ (23 行 × 8 列, 3 KB) - resource 39783: ビューポイント写真 (zip, 14 件, 不使用) - 全件 県管理クルーズ実証データ - ライセンス: クリエイティブ・コモンズ表示 (CC-BY) メモリ対策: Figure ごとに plt.close('all') で解放。 """ from __future__ import annotations import sys import time from pathlib import Path sys.path.insert(0, str(Path(__file__).parent)) from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure, ensure_dataset) import numpy as np import pandas as pd import geopandas as gpd from shapely.geometry import Point, LineString import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.patches import Patch from matplotlib.lines import Line2D plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False t_all = time.time() print("=== L82 せとうちモニタークルーズ実施結果 単独 3 研究例分析 ===", flush=True) # ============================================================================= # 0. 定数・パス # ============================================================================= TARGET_CRS = "EPSG:6671" # JGD2011 平面直角第 III 系 (m 単位) DATASET_ID = 1280 RES_DAILY = 39782 RES_PORTS = 39781 RES_VIEW = 39783 DATA_DIR = ROOT / "data" / "extras" / "L82_setouchi_monitor_cruise" DAILY_CSV = DATA_DIR / "monitor_cruise_daily.csv" PORTS_CSV = DATA_DIR / "monitor_cruise_ports.csv" # しまたびライン月次 (X04 共有 cache、 RQ1/RQ3 従属参照) SHIMATABI_CSV = ROOT / "data" / "extras" / "shimatabi_monthly.csv" # 行政界 (L44 既キャッシュ) ADMIN_GPKG = ROOT / "data" / "extras" / "L44_storm_surge" / "_cache" / "admin_diss.gpkg" # 観光商品の運用パラメータ (本記事の試算前提 = 公開データに無い値の仮定) ASSUMED_CAPACITY = 70 # 1 航海定員 (観測最大値ベース仮定) ASSUMED_FARE = 5000 # 仮定単価 (円/人, 類似観光船相場) HARBOR_HIROSHIMA = "広島港（広島港宇品旅客ターミナル）" # ============================================================================= # 1. データ取得 # ============================================================================= print("\n[1] データ取得", flush=True) t1 = time.time() ensure_dataset(DAILY_CSV, dataset_id=DATASET_ID, resource_id=RES_DAILY, label="モニタークルーズ航路×日別 (#1280 res 39782)") ensure_dataset(PORTS_CSV, dataset_id=DATASET_ID, resource_id=RES_PORTS, label="モニタークルーズ寄港地マスタ (#1280 res 39781)") # しまたびライン (X04 cache 再利用、 取得不要だが、 無ければ取得) ensure_dataset(SHIMATABI_CSV, resource_id=39787, label="しまたびライン月次 (X04 cache 再利用)") df_daily = pd.read_csv(DAILY_CSV, encoding="cp932") df_ports = pd.read_csv(PORTS_CSV, encoding="cp932") df_shimatabi = pd.read_csv(SHIMATABI_CSV, encoding="cp932").dropna( subset=["寄港地（桟橋名）"]) df_daily["運航日"] = pd.to_datetime(df_daily["運航日"]) n_voyages = df_daily.drop_duplicates(["航路名", "運航日"]).shape[0] n_routes = df_daily["航路名"].nunique() n_ports = df_ports["寄港地ID"].nunique() n_dates = df_daily["運航日"].nunique() print(f" 航路: {n_routes}, 運航日: {n_dates}, 航海数: {n_voyages}, " f"寄港地: {n_ports}", flush=True) print(f" ({time.time()-t1:.2f}s)", flush=True) # ============================================================================= # 2. 航路 × 寄港地 順序ジオメトリ生成 # ============================================================================= print("\n[2] 航路ジオメトリ生成", flush=True) t2 = time.time() # 寄港地マスタ整形 ports_master = df_ports[["寄港地ID", "寄港地名（桟橋名）", "住所", "緯度", "経度", "URL", "区間情報"]].copy() ports_master = ports_master.rename(columns={"寄港地名（桟橋名）": "寄港地名"}) ports_master["市町"] = ports_master["住所"].str.extract( r"広島県(.+?[市町区])" ).fillna("?") # 観光資源 (ビューポイント所在地) のフラグ付け viewpoint_ports = (df_ports[df_ports["ビューポイント写真情報"] != "−"]["寄港地ID"].tolist()) ports_master["ビューポイント"] = ports_master["寄港地ID"].isin(viewpoint_ports) n_viewpoints = ports_master["ビューポイント"].sum() gdf_ports = gpd.GeoDataFrame( ports_master, geometry=[Point(lon, lat) for lon, lat in zip(ports_master["経度"], ports_master["緯度"])], crs="EPSG:4326").to_crs(TARGET_CRS) # Haversine 距離 (numpy ベクトル化) def haversine_km(lat1, lon1, lat2, lon2): R = 6371.0 p1, p2 = np.radians(lat1), np.radians(lat2) dp = np.radians(lat2 - lat1) dl = np.radians(lon2 - lon1) a = np.sin(dp/2.0)**2 + np.cos(p1) * np.cos(p2) * np.sin(dl/2.0)**2 return 2.0 * R * np.arcsin(np.sqrt(a)) # 各航路の寄港順序 (寄港桟橋順 で並べる) route_geom = {} for name, g in df_daily.groupby("航路名"): seq = (g.drop_duplicates("寄港地ID") .sort_values("寄港桟橋順") .reset_index(drop=True)) coords_lonlat = list(zip(seq["経度"], seq["緯度"])) seg_info = seq["区間情報"].tolist() stop_names = seq["寄港地（桟橋名）"].tolist() # 区間距離 (連続 2 点間の Haversine) lats = seq["緯度"].values lons = seq["経度"].values seg_dists = haversine_km(lats[:-1], lons[:-1], lats[1:], lons[1:]) total_dist_km = float(seg_dists.sum()) avg_dist_km = float(seg_dists.mean()) max_dist_km = float(seg_dists.max()) min_dist_km = float(seg_dists.min()) # 始発 = 終着 (戻り型) フラグ is_loop = (stop_names[0] == stop_names[-1]) # 方向タグ (航路名末尾 「（東）」 「（西）」 から抽出) if "（東）" in name: direction = "東向" elif "（西）" in name: direction = "西向" else: direction = "?" # 観光資源密度 (ビューポイント数 / 寄港地数) n_view = sum(1 for sid in seq["寄港地ID"].tolist() if sid in viewpoint_ports) view_density = n_view / max(len(stop_names), 1) route_geom[name] = { "stops": stop_names, "stop_ids": seq["寄港地ID"].tolist(), "coords": coords_lonlat, "seg": seg_info, "seg_dists": seg_dists, "total_km": total_dist_km, "avg_km": avg_dist_km, "max_km": max_dist_km, "min_km": min_dist_km, "is_loop": is_loop, "direction": direction, "n_stops": len(stop_names), "n_view": n_view, "view_density": view_density, } print(f" ({time.time()-t2:.2f}s)", flush=True) # ============================================================================= # 3. RQ1 集計: 航路設計の構造プロファイル # ============================================================================= print("\n[3] RQ1 集計 — 航路設計", flush=True) t3 = time.time() # (1) 7 航路 個体台帳 (RQ1 主役表) shimatabi_ports_set = set(df_shimatabi["寄港地（桟橋名）"].dropna().unique()) # Jaccard: 各 monitor 航路と既存しまたび 14 港の重なり def jaccard(set_a, set_b): a, b = set(set_a), set(set_b) u = a | b if not u: return 0.0 return len(a & b) / len(u) T_route_records = [] for name, info in route_geom.items(): stop_set = set(info["stops"]) j_shimatabi = jaccard(stop_set, shimatabi_ports_set) n_overlap = len(stop_set & shimatabi_ports_set) n_new = sum(1 for s in info["seg"] if s == "新規区間") n_old = sum(1 for s in info["seg"] if s == "既存区間") new_ratio = n_new / max(n_new + n_old, 1) # 1 航海所要時間推定 (= 総距離 km / 平均速度 + 各港停泊 30 分) # 観光船一般速度 ~25 km/h、 各港停泊 30 分仮定 speed_kmh = 25.0 stop_min = 30.0 voyage_min = info["total_km"] / speed_kmh * 60.0 \ + max(info["n_stops"] - 1, 0) * stop_min T_route_records.append({ "航路名": name, "方向": info["direction"], "寄港地数": info["n_stops"], "総距離_km": round(info["total_km"], 1), "平均区間_km": round(info["avg_km"], 1), "最大区間_km": round(info["max_km"], 1), "最小区間_km": round(info["min_km"], 1), "推定所要_分": round(voyage_min, 0), "新規率": round(new_ratio, 2), "新規区間数": n_new, "既存区間数": n_old, "戻り型": info["is_loop"], "観光資源数": info["n_view"], "観光資源密度": round(info["view_density"], 2), "しまたびJaccard": round(j_shimatabi, 2), "しまたび重複港数": n_overlap, }) T_route = pd.DataFrame(T_route_records).sort_values( "総距離_km", ascending=False).reset_index(drop=True) T_route.insert(0, "順位", np.arange(1, len(T_route) + 1)) # H1 検証: 総距離 80-120 km 帯集中 + max/min ≤ 1.5 dist_arr = T_route["総距離_km"].values n_in_band = ((dist_arr >= 80) & (dist_arr <= 120)).sum() dist_ratio = dist_arr.max() / dist_arr.min() h1_ok = (n_in_band >= 6) and (dist_ratio <= 1.5) # H2 検証: 東西ペア構造 ≥ 2 ペア # ペア = 同一基幹区間 (e.g. 広島⇔尾道) を東向と西向の両方で運航している組 import re as _re_pair def _route_base(name): """航路名から方向と番号を除去してペアキーを抽出。 例: '③三原⇒広島コース（西）' → '三原⇒広島' '③広島⇒三原コース（東）' → '広島⇒三原' さらに ⇔ で対称ペアキー化: '広島⇒尾道', '尾道⇒広島' → 共通キー '広島⇔尾道'""" base = _re_pair.sub(r"^[①②③④⑤⑥⑦⑧⑨⑩]", "", name) base = base.replace("（東）", "").replace("（西）", "").replace("コース", "") return base def _pair_key(name): base = _route_base(name) # '広島⇒尾道' or '尾道⇒広島' を sorted で対称化 if "⇒" in base: a, b = base.split("⇒", 1) return "⇔".join(sorted([a, b])) return base T_route["ペアキー"] = T_route["航路名"].apply(_pair_key) pair_counts = T_route.groupby("ペアキー")["方向"].nunique() n_pairs = (pair_counts >= 2).sum() # 東西両方向ある ペア n_loop = T_route["戻り型"].sum() h2_ok = n_pairs >= 2 # 戻り型と片道型の集計 n_oneway = (T_route["戻り型"] == False).sum() # 全7と予想 n_east = (T_route["方向"] == "東向").sum() n_west = (T_route["方向"] == "西向").sum() # (2) 寄港地ランキング (各港が何航路で寄港されたか) — X04 では未集計 port_route_count = [] for sid in ports_master["寄港地ID"]: cnt = sum(1 for info in route_geom.values() if sid in info["stop_ids"]) port_route_count.append(cnt) ports_master["航路寄港回数"] = port_route_count T_port_freq = ports_master[ ["寄港地ID", "寄港地名", "市町", "区間情報", "ビューポイント", "航路寄港回数"] ].sort_values(["航路寄港回数", "寄港地名"], ascending=[False, True]) T_port_freq = T_port_freq.reset_index(drop=True) T_port_freq.insert(0, "順位", np.arange(1, len(T_port_freq) + 1)) # (3) 始発 / 終着 / 方向 / ペアキー の役割分類 role_records = [] for name, info in route_geom.items(): start = info["stops"][0] end = info["stops"][-1] role_records.append({ "航路名": name, "方向": info["direction"], "ペアキー": _pair_key(name), "始発": start, "終着": end, "戻り型": info["is_loop"], }) T_role = pd.DataFrame(role_records).sort_values( ["ペアキー", "方向"]).reset_index(drop=True) # 広島港を含む港名は「広島港」 で始まるか部分一致する def _is_hiroshima(name): return "広島港" in str(name) n_hiroshima_start = T_role["始発"].apply(_is_hiroshima).sum() n_hiroshima_end = T_role["終着"].apply(_is_hiroshima).sum() # (4) 区間情報統計 (新規/既存) — 寄港地ベースのみ。 X04 と異なり # 「区間距離 km」 で重み付けした距離ベース新規率を独自に計算 total_dist_all = sum(info["total_km"] for info in route_geom.values()) new_dist_total = 0.0 old_dist_total = 0.0 for name, info in route_geom.items(): seg = info["seg"] sd = info["seg_dists"] # 区間 i の区間情報 = 寄港地 i+1 の区間情報 (始発側に新規 / 既存属性は無い) for i in range(len(sd)): # 区間 i は寄港地 i+1 (= seg[i+1]) が新規/既存のどちらかで判定 cls = seg[i + 1] if i + 1 < len(seg) else seg[i] if cls == "新規区間": new_dist_total += sd[i] else: old_dist_total += sd[i] dist_new_ratio = new_dist_total / max(new_dist_total + old_dist_total, 1) T_dist_summary = pd.DataFrame([ {"指標": "総運航距離 (km)", "値": f"{total_dist_all:.1f}"}, {"指標": "新規区間 距離合計 (km)", "値": f"{new_dist_total:.1f}"}, {"指標": "既存区間 距離合計 (km)", "値": f"{old_dist_total:.1f}"}, {"指標": "距離ベース新規率", "値": f"{dist_new_ratio*100:.1f}%"}, {"指標": "寄港地ベース新規率 (参考)", "値": (f"{(ports_master['区間情報']=='新規区間').sum()}/" f"{n_ports} = " f"{(ports_master['区間情報']=='新規区間').sum()/n_ports*100:.1f}%")}, ]) # H1, H2 検証ログ print(f" 総距離 (km): min={dist_arr.min():.1f}, max={dist_arr.max():.1f}, " f"max/min={dist_ratio:.2f}", flush=True) print(f" H1 (80-120km 帯 + max/min ≤ 1.5): {h1_ok}, " f"80-120km 帯 = {n_in_band}/{n_routes}", flush=True) print(f" H2 (東西ペア ≥ 2): {h2_ok}, ペア数 = {n_pairs}, " f"東向 = {n_east}, 西向 = {n_west}", flush=True) print(f" ペアキー: {pair_counts.to_dict()}", flush=True) print(f" ({time.time()-t3:.2f}s)", flush=True) # ============================================================================= # 4. RQ2 集計: 乗船客特性 — 運航回数効果 # ============================================================================= print("\n[4] RQ2 集計 — 乗船客特性", flush=True) t4 = time.time() # 1 航海 = 1 (運航日, 航路) の単位 trip = (df_daily.drop_duplicates(["運航日", "航路名"]) [["運航日", "航路名", "観光客数（人）"]] .sort_values(["運航日", "航路名"]).reset_index(drop=True)) trip["運航日"] = pd.to_datetime(trip["運航日"]) trip["曜日"] = trip["運航日"].dt.day_name() trip["平日週末"] = trip["運航日"].dt.weekday.apply( lambda w: "週末" if w >= 5 else "平日") # 9 運航日は 3 ターンに分かれる: 12/06-08 (第1ターン), 12/13-15 (第2), # 12/20-22 (第3)。 6, 13, 20 を境にターン番号付与 def _turn(d): day = d.day if day <= 8: return "第1ターン (12/06-08)" if day <= 15: return "第2ターン (12/13-15)" return "第3ターン (12/20-22)" trip["ターン"] = trip["運航日"].apply(_turn) # 各航路内での運航回数順位 (1 回目, 2 回目, 3 回目) trip = trip.sort_values(["航路名", "運航日"]).reset_index(drop=True) trip["運航回"] = trip.groupby("航路名").cumcount() + 1 # 方向タグを付与 (全 7 航路は片道型のため、 「東向 / 西向」 で 2 区分) direction_map = {name: info["direction"] for name, info in route_geom.items()} trip["方向"] = trip["航路名"].map(direction_map) trip["ルートタイプ"] = trip["方向"] # = 東向 or 西向 T_trip = trip.copy() # (1) 全 15 航海 統計 visitors = T_trip["観光客数（人）"].astype(float).values mean_v = visitors.mean() std_v = visitors.std(ddof=0) cv_v = std_v / mean_v if mean_v > 0 else 0.0 total_v = int(visitors.sum()) T_pas_summary = pd.DataFrame([ {"指標": "総航海数", "値": f"{len(visitors)}"}, {"指標": "総観光客数 (人)", "値": f"{total_v:,}"}, {"指標": "平均観光客数 (人/航海)", "値": f"{mean_v:.1f}"}, {"指標": "中央値 (人/航海)", "値": f"{np.median(visitors):.1f}"}, {"指標": "標準偏差 (人)", "値": f"{std_v:.2f}"}, {"指標": "変動係数 CV (= std/mean)", "値": f"{cv_v:.3f}"}, {"指標": "最大 (人)", "値": f"{int(visitors.max())}"}, {"指標": "最小 (人)", "値": f"{int(visitors.min())}"}, {"指標": "推定収容率 (= 平均 / 70 人)", "値": f"{mean_v/ASSUMED_CAPACITY*100:.1f}%"}, ]) # H4 検証: CV ≤ 0.20 h4_ok = cv_v <= 0.20 # (2) 運航回別 平均客数 (1 回目 vs 2 回目 vs 3 回目) T_voyage_idx = (T_trip.groupby("運航回") .agg(航海数=("観光客数（人）", "count"), 平均客数=("観光客数（人）", "mean"), 中央値=("観光客数（人）", "median"), 最大=("観光客数（人）", "max"), 最小=("観光客数（人）", "min")) .round(1).reset_index()) # H3 検証: 第 2 運航は第 1 運航より 15% 以上少ない v1 = T_voyage_idx.loc[T_voyage_idx["運航回"] == 1, "平均客数"].iloc[0] v2 = T_voyage_idx.loc[T_voyage_idx["運航回"] == 2, "平均客数"].iloc[0] \ if (T_voyage_idx["運航回"] == 2).any() else float("nan") delta_pct = (v1 - v2) / v1 * 100 if v1 > 0 else 0.0 h3_ok = delta_pct >= 15.0 # (3) ターン別 (実証実験は全て平日 Tue-Thu のため、 平日週末区分は無効。 # 代わりに 3 ターン (12/06-08, 12/13-15, 12/20-22) で区分) T_turn = (T_trip.groupby("ターン") .agg(航海数=("観光客数（人）", "count"), 平均客数=("観光客数（人）", "mean"), 合計客数=("観光客数（人）", "sum")) .round(1).reset_index()) # 互換のため T_weekday は ターン別と同じデータを名前変更で残す (HTML 互換) T_weekday = T_turn.rename(columns={"ターン": "平日週末"}) # 古い変数名互換 # (4) 方向別 (東向 vs 西向) T_routetype = (T_trip.groupby("ルートタイプ") .agg(航海数=("観光客数（人）", "count"), 平均客数=("観光客数（人）", "mean"), 合計客数=("観光客数（人）", "sum")) .round(1).reset_index()) # 東向 / 西向 のうち存在する方のみ並べる T_routetype = T_routetype.sort_values("ルートタイプ").reset_index(drop=True) # (5) 運航日 (9 日) ごとの集計 — small multiples 用 T_daily_total = (T_trip.groupby(T_trip["運航日"].dt.date) .agg(航海数=("観光客数（人）", "count"), 合計客数=("観光客数（人）", "sum"), 平均客数=("観光客数（人）", "mean")) .round(1).reset_index()) T_daily_total = T_daily_total.rename(columns={"運航日": "日付"}) # (6) 航路 × 運航回 ピボット pivot_voyage = T_trip.pivot_table(index="航路名", columns="運航回", values="観光客数（人）", aggfunc="first") print(f" CV = {cv_v:.3f}, H4 (CV ≤ 0.20): {h4_ok}", flush=True) print(f" 第1回平均 = {v1:.1f}, 第2回平均 = {v2:.1f}, " f"減少 = {delta_pct:.1f}%, H3 (≥15%): {h3_ok}", flush=True) print(f" ({time.time()-t4:.2f}s)", flush=True) # ============================================================================= # 5. RQ3 集計: 観光商品評価 — 採算性 + 定期化適性 # ============================================================================= print("\n[5] RQ3 集計 — 観光商品評価", flush=True) t5 = time.time() # (1) 航路別 1 航海平均客数 + 収容率 + 推定収益 route_voyage_stats = [] for name in T_route["航路名"]: sub = T_trip[T_trip["航路名"] == name] avg = sub["観光客数（人）"].mean() n_v = len(sub) total_p = sub["観光客数（人）"].sum() occ = avg / ASSUMED_CAPACITY * 100 rev_per_v = avg * ASSUMED_FARE # 1 航海あたり推定収益 (円) rev_total = total_p * ASSUMED_FARE # 全航海推定収益 (円) route_voyage_stats.append({ "航路名": name, "航海数": n_v, "総客数": int(total_p), "平均客数": round(avg, 1), "推定収容率_%": round(occ, 1), "推定収益_万円_per航海": round(rev_per_v / 10000, 1), "推定収益合計_万円": round(rev_total / 10000, 1), }) T_revenue = pd.DataFrame(route_voyage_stats) # (2) 定期化適性スコア (本記事独自指標) # = 平均客数 × 新規性ボーナス × 観光資源密度ボーナス / 距離係数 # 距離係数: 短いほどコスト効率が高い → 100 km を 1.0 基準 T_score_records = [] for _, r in T_route.iterrows(): name = r["航路名"] avg = T_revenue.loc[T_revenue["航路名"] == name, "平均客数"].iloc[0] new_bonus = 1.0 + r["新規率"] * 0.3 # 新規率 100% で +30% view_bonus = 1.0 + r["観光資源密度"] * 0.2 # 観光資源密度 100% で +20% dist_factor = max(r["総距離_km"] / 100.0, 0.5) # 100km 基準 score = avg * new_bonus * view_bonus / dist_factor T_score_records.append({ "航路名": name, "平均客数": avg, "新規率": r["新規率"], "観光資源密度": r["観光資源密度"], "総距離_km": r["総距離_km"], "新規ボーナス": round(new_bonus, 2), "観光資源ボーナス": round(view_bonus, 2), "距離係数": round(dist_factor, 2), "定期化適性スコア": round(score, 1), }) T_periodic = (pd.DataFrame(T_score_records) .sort_values("定期化適性スコア", ascending=False) .reset_index(drop=True)) T_periodic.insert(0, "順位", np.arange(1, len(T_periodic) + 1)) # H5: 上位 2 航路の特定 top2_routes = T_periodic.head(2)["航路名"].tolist() top2_scores = T_periodic.head(2)["定期化適性スコア"].tolist() h5_ok = len(top2_routes) >= 2 and top2_scores[1] >= 30.0 # (3) しまたび月平均との比較 (RQ3 別角度: モニターは月当たり何人運んだか) # 月別合計乗降客数 (X04 と異なる集計: 月平均 / 14 港 / 全期間でもOK) month_cols = ["4月（人）", "5月（人）", "6月（人）", "7月（人）", "8月（人）", "9月（人）", "10月（人）", "11月（人）", "12月（人）"] df_shi = df_shimatabi.copy() for c in month_cols: df_shi[c] = pd.to_numeric(df_shi[c], errors="coerce") shi_total_8mo = df_shi[month_cols].sum().sum() # 4-12 月 9 ヶ月の総需要 shi_avg_per_month = shi_total_8mo / 9 shi_dec = df_shi["12月（人）"].sum() # 12月のみ monitor_total = total_v # = 855 人想定 monitor_per_voyage = mean_v ratio_to_dec = monitor_total / max(shi_dec, 1) T_compare = pd.DataFrame([ {"指標": "モニタークルーズ 12 月 (15 航海) 総客数", "値": f"{monitor_total:,} 人"}, {"指標": "モニタークルーズ 1 航海平均", "値": f"{monitor_per_voyage:.1f} 人/航海"}, {"指標": "しまたびライン 2022年12月 総乗降客数", "値": f"{int(shi_dec):,} 人"}, {"指標": "しまたびライン 4-12月 平均月需要", "値": f"{int(shi_avg_per_month):,} 人/月"}, {"指標": "モニター/しまたび 12月 比率", "値": f"{ratio_to_dec*100:.1f}%"}, {"指標": "1 航海 vs しまたび 月需要 (×単純比)", "値": f"{monitor_per_voyage/shi_avg_per_month*100:.2f}%"}, ]) # (4) 政策的位置付け (定性表) T_policy = pd.DataFrame([ {"観点": "実施期間", "内容": "2022 年 12 月 6-22 日 (9 運航日)、 一過性社会実験"}, {"観点": "実施主体", "内容": "広島県 + 瀬戸内海汽船・しまなみ海運 (2 船社)"}, {"観点": "目的", "内容": "新規寄港地 7 港の集客力検証 + 既存ルートとの差別化評価"}, {"観点": "規模", "内容": f"7 航路 × {n_dates} 日 = {n_voyages} 航海 / " f"{n_ports} 寄港地 / 総運航距離 {total_dist_all:.0f} km"}, {"観点": "成果指標", "内容": f"総客数 {monitor_total} 人 / 1 航海平均 {mean_v:.1f} 人 " f"(推定収容率 {mean_v/ASSUMED_CAPACITY*100:.0f}%)"}, {"観点": "アンケート", "内容": "本データには集計済み数値のみ。 個票は別途観光振興課で保管"}, {"観点": "後継事業", "内容": "本データから直接の定期化決定情報は記録されていない。 " "2023 年以降の運航有無は別途確認要"}, {"観点": "本記事の限界", "内容": "料金は公開されていない (5,000 円/人は仮定)。 " "船舶運航コスト・人件費も公開されていないため、 " "純利益の実数値は計算不可。 本記事は収益スケール" "の比較指標として提示する"}, ]) print(f" 定期化適性 上位 2 = {top2_routes}", flush=True) print(f" H5 (上位 2 航路特定 + スコア ≥ 30): {h5_ok}", flush=True) print(f" ({time.time()-t5:.2f}s)", flush=True) # ============================================================================= # 6. 仮説検証総合 # ============================================================================= print("\n[6] 仮説検証", flush=True) t6 = time.time() T_hyp = pd.DataFrame([ {"仮説": "H1 総距離 80-120km + max/min ≤ 1.5 (RQ1)", "観測値": (f"{dist_arr.min():.1f}-{dist_arr.max():.1f} km, " f"max/min = {dist_ratio:.2f}, " f"80-120km 帯 = {n_in_band}/{n_routes}"), "判定": "支持" if h1_ok else "部分支持", "解釈": (f"H1 {'支持' if h1_ok else '部分支持'}: " f"7 航路の総距離は {dist_arr.min():.1f} - " f"{dist_arr.max():.1f} km の{dist_ratio:.2f} 倍" f"の狭い帯に集中 ({n_in_band}/{n_routes} 航路が 80-120 km)。 " f"「半日 〜 1 日かけて運航する片道商品」として " f"距離規格が均質化されている設計。 " f"短距離小ループは存在せず、 「島を眺めながら半日 〜 1 日" f"楽しむ」 標準長の商品で統一されている")}, {"仮説": "H2 東西ペア構造 ≥ 2 ペア (RQ1)", "観測値": (f"東西ペア数 = {n_pairs}, " f"東向 = {n_east} 航路, 西向 = {n_west} 航路, " f"全 7 航路は片道型 (戻り型 = 0)"), "判定": "支持" if h2_ok else "反証", "解釈": (f"H2 {'支持' if h2_ok else '反証'}: " f"{n_pairs} ペアの東西ペア構造を確認 " f"(尾道⇔広島 / 三原⇔広島 / しまなみ・尾道⇔広島)。 " f"7 航路はすべて片道型で、 戻り型ループは 0。 " f"乗船客が東向き / 西向き、 好みの方向で楽しめる " f"双方向商品設計が量的に確認される。 " f"「鞆の浦⇒広島 (西)」 のみペアが組まれていない単発航路で、 " f"鞆の浦は最東端の特殊起点として位置付けられる")}, {"仮説": "H3 第 2 運航 客数減 ≥ 15% (RQ2)", "観測値": f"第1回平均 {v1:.1f} 人 → 第2回平均 {v2:.1f} 人 " f"(減少率 {delta_pct:.1f}%)", "判定": "支持" if h3_ok else "反証", "解釈": (f"H3 {'支持' if h3_ok else '反証'}: " f"第 1 回運航から第 2 回運航で減少率 {delta_pct:.1f}%。 " f"{'物珍しさ効果の減衰が量的に確認' if h3_ok else '減少率は 15% に届かず、 反復運航でも需要が比較的安定'} " f"していることが分かる")}, {"仮説": "H4 客数 CV ≤ 0.20 (RQ2)", "観測値": f"CV = std/mean = {std_v:.2f}/{mean_v:.1f} = {cv_v:.3f}", "判定": "支持" if h4_ok else "反証", "解釈": (f"H4 {'支持' if h4_ok else '反証'}: " f"15 航海の客数 CV は{cv_v:.2f}。 " f"{'観光商品として客数が極めて安定' if h4_ok else '客数のばらつきは中程度。 ' + ('航路間差' if cv_v < 0.30 else '航路間差が大きい')} " f"。 観測客数は {int(visitors.min())}-{int(visitors.max())} 人帯に収まり、 " f"50 人台の地ベースに +20 人の上振れが観察される")}, {"仮説": "H5 定期化適性 上位 2 航路特定 (RQ3)", "観測値": (f"上位 1 = {top2_routes[0]} ({top2_scores[0]:.1f}), " f"上位 2 = {top2_routes[1]} ({top2_scores[1]:.1f})"), "判定": "支持" if h5_ok else "部分支持", "解釈": (f"H5 {'支持' if h5_ok else '部分支持'}: " f"定期化適性スコア上位 2 航路は{top2_routes[0]}と" f"{top2_routes[1]}。 " f"これらは新規率と観光資源密度の積が高く、 距離効率も良い。 " f"既存しまたび 12 月需要 {int(shi_dec):,} 人と比べると、 " f"モニター 15 航海合計 {monitor_total} 人は規模で " f"{ratio_to_dec*100:.1f}% に留まるが、 " f"1 航海あたりの集客効率では遜色ない")}, ]) print(T_hyp[["仮説", "観測値", "判定"]].to_string(index=False), flush=True) print(f" ({time.time()-t6:.2f}s)", flush=True) # ============================================================================= # 7. CSV 出力 # ============================================================================= print("\n[7] CSV 出力", flush=True) t7 = time.time() T_route.to_csv(ASSETS / "L82_route_overview.csv", index=False, encoding="utf-8-sig") T_port_freq.to_csv(ASSETS / "L82_port_frequency.csv", index=False, encoding="utf-8-sig") T_role.to_csv(ASSETS / "L82_route_role.csv", index=False, encoding="utf-8-sig") T_dist_summary.to_csv(ASSETS / "L82_dist_summary.csv", index=False, encoding="utf-8-sig") T_pas_summary.to_csv(ASSETS / "L82_passenger_summary.csv", index=False, encoding="utf-8-sig") T_voyage_idx.to_csv(ASSETS / "L82_voyage_index.csv", index=False, encoding="utf-8-sig") T_turn.to_csv(ASSETS / "L82_turn.csv", index=False, encoding="utf-8-sig") T_routetype.to_csv(ASSETS / "L82_routetype.csv", index=False, encoding="utf-8-sig") T_daily_total.to_csv(ASSETS / "L82_daily_total.csv", index=False, encoding="utf-8-sig") T_revenue.to_csv(ASSETS / "L82_revenue.csv", index=False, encoding="utf-8-sig") T_periodic.to_csv(ASSETS / "L82_periodic_score.csv", index=False, encoding="utf-8-sig") T_compare.to_csv(ASSETS / "L82_compare_shimatabi.csv", index=False, encoding="utf-8-sig") T_policy.to_csv(ASSETS / "L82_policy.csv", index=False, encoding="utf-8-sig") T_hyp.to_csv(ASSETS / "L82_hypothesis_check.csv", index=False, encoding="utf-8-sig") print(f" ({time.time()-t7:.2f}s)", flush=True) # ============================================================================= # 8. 行政界読込み (地図用) # ============================================================================= print("\n[8] 行政界読込み", flush=True) t8 = time.time() admin = gpd.read_file(ADMIN_GPKG).to_crs(TARGET_CRS) print(f" 行政界 polygons: {len(admin)}", flush=True) print(f" ({time.time()-t8:.2f}s)", flush=True) # ============================================================================= # 9. 図の生成 # ============================================================================= print("\n[9] 図の生成", flush=True) t9 = time.time() # 表示 bbox (沿岸クルーズエリア = 広島市・呉市・尾道市・福山市・三原市) # EPSG:6671 (m) 単位で # 広島港 (132.45, 34.35) ~ 鞆の浦 (133.38, 34.38)、 island bottom 34.18 gdf_x_min = float(gdf_ports.geometry.x.min()) - 8000 gdf_x_max = float(gdf_ports.geometry.x.max()) + 8000 gdf_y_min = float(gdf_ports.geometry.y.min()) - 8000 gdf_y_max = float(gdf_ports.geometry.y.max()) + 8000 ROUTE_COLORS = ["#0969da", "#cf222e", "#16a34a", "#fb8500", "#9333ea", "#0891b2", "#dc2626"] ROUTE_COLOR_MAP = {r: ROUTE_COLORS[i % len(ROUTE_COLORS)] for i, r in enumerate(sorted(route_geom.keys()))} def save_fig(name, dpi=120): p = ASSETS / name plt.savefig(p, dpi=dpi, bbox_inches="tight", facecolor="white") plt.close('all') return p # ---- 図 1 (RQ1): 7 航路ルートマップ + 観光資源バブル ---- print(" fig1: 7 航路ルートマップ", flush=True) fig, ax = plt.subplots(figsize=(13, 8)) admin.plot(ax=ax, color="#fff8e8", edgecolor="#888", linewidth=0.4, alpha=0.55) # 寄港地 (新規 = 赤丸 / 既存 = 青丸 / ビューポイント = 大きめ) for label, color in [("新規区間", "#cf222e"), ("既存区間", "#0969da")]: sub = gdf_ports[gdf_ports["区間情報"] == label] sizes = np.where(sub["ビューポイント"], 180, 80) ax.scatter(sub.geometry.x, sub.geometry.y, s=sizes, c=color, edgecolor="black", linewidth=0.8, alpha=0.85, zorder=5, label=f"{label}寄港地 ({len(sub)} 港)") # ビューポイント有り = 黄星マーカー追加 sub_vp = gdf_ports[gdf_ports["ビューポイント"]] ax.scatter(sub_vp.geometry.x, sub_vp.geometry.y, s=60, c="#ffd700", marker="*", edgecolor="black", linewidth=0.6, alpha=0.95, zorder=6, label=f"観光資源 (ビューポイント, {len(sub_vp)} 港)") # 港名ラベル (主要港のみ太字) major_names = {"広島港（広島港宇品旅客ターミナル）": "広島", "鞆の浦": "鞆の浦", "三原（内港旅客船ターミナル）": "三原", "呉港(中央桟橋ターミナル)": "呉", "尾道（駅前桟橋旅客船）": "尾道", "瀬戸田（瀬戸田港中央桟橋）": "瀬戸田"} for _, r in gdf_ports.iterrows(): short = major_names.get(r["寄港地名"], r["寄港地名"][:6]) is_major = r["寄港地名"] in major_names ax.annotate(short, (r.geometry.x, r.geometry.y), xytext=(7, 5), textcoords="offset points", fontsize=9 if is_major else 7.5, fontweight="bold" if is_major else "normal", color="#1f2328", alpha=0.92) # 各航路の経路を線で for name in sorted(route_geom.keys()): info = route_geom[name] coords = info["coords"] # (lon, lat) if len(coords) < 2: continue pts = gpd.GeoSeries([Point(lon, lat) for lon, lat in coords], crs="EPSG:4326").to_crs(TARGET_CRS) xs = [p.x for p in pts] ys = [p.y for p in pts] color = ROUTE_COLOR_MAP[name] ax.plot(xs, ys, color=color, linewidth=2.0, alpha=0.55, zorder=2, label=f"{name[:20]} ({info['total_km']:.0f}km)") ax.set_xlim(gdf_x_min, gdf_x_max) ax.set_ylim(gdf_y_min, gdf_y_max) ax.set_aspect("equal") ax.set_title(f"図 1 (RQ1): モニタークルーズ {n_routes} 航路 × {n_ports} 寄港地 — " f"赤丸=新規区間 / 青丸=既存区間 / 黄星=観光資源", fontsize=11) ax.set_xlabel("X (m, EPSG:6671)") ax.set_ylabel("Y (m, EPSG:6671)") ax.legend(loc="upper left", fontsize=7.5, framealpha=0.92, ncol=1) ax.grid(True, linestyle="--", alpha=0.3) plt.tight_layout() save_fig("L82_fig1_route_map.png") # ---- 図 2 (RQ1): 7 航路の構造プロファイル (距離 / 寄港数 / 新規率) ---- print(" fig2: 航路構造プロファイル", flush=True) fig, axes = plt.subplots(1, 3, figsize=(17, 5.5)) route_short = [n.replace("コース（仮）", "")[:14] for n in T_route["航路名"]] # (a) 総距離 ランキング (方向で色分け) ax = axes[0] ds = T_route.sort_values("総距離_km", ascending=True) labels_b = [n.replace("コース（仮）", "")[:14] for n in ds["航路名"]] colors_b = ["#cf222e" if d == "東向" else "#0969da" for d in ds["方向"]] ax.barh(np.arange(len(ds)), ds["総距離_km"].values, color=colors_b, edgecolor="#222", linewidth=0.5, alpha=0.85) for i, (v, n_stop) in enumerate(zip(ds["総距離_km"], ds["寄港地数"])): ax.text(v + 2, i, f"{v:.0f} km / {n_stop} 港", va="center", fontsize=9) ax.set_yticks(np.arange(len(ds))) ax.set_yticklabels(labels_b, fontsize=9) ax.set_xlabel("総距離 (km, Haversine)") ax.set_title(f"(a) 総距離ランキング\n" f"赤=東向 ({n_east}本) / 青=西向 ({n_west}本)", fontsize=10.5) ax.axvspan(80, 120, alpha=0.10, color="#16a34a", zorder=0, label="H1 帯 80-120km") ax.grid(True, axis="x", linestyle="--", alpha=0.3) # (b) 新規率 ランキング (寄港地ベース vs 距離ベース) ax = axes[1] ds2 = T_route.sort_values("新規率", ascending=True) labels_a = [n.replace("コース（仮）", "")[:14] for n in ds2["航路名"]] ax.barh(np.arange(len(ds2)), ds2["新規率"] * 100, color="#fb8500", edgecolor="#222", linewidth=0.5, alpha=0.85) for i, (v, nn, oo) in enumerate(zip(ds2["新規率"], ds2["新規区間数"], ds2["既存区間数"])): ax.text(v * 100 + 1, i, f"{v*100:.0f}% (新規{nn}/既存{oo})", va="center", fontsize=9) ax.set_yticks(np.arange(len(ds2))) ax.set_yticklabels(labels_a, fontsize=9) ax.set_xlabel("新規率 (新規寄港地数 / 総寄港地数, %)") ax.set_title(f"(b) 新規率ランキング\n" f"距離ベース新規率 {dist_new_ratio*100:.1f}% (本記事独自指標)", fontsize=10.5) ax.set_xlim(0, 105) ax.grid(True, axis="x", linestyle="--", alpha=0.3) # (c) 観光資源密度 vs しまたび Jaccard ax = axes[2] xs = T_route["観光資源密度"].values ys = T_route["しまたびJaccard"].values sizes = T_route["寄港地数"] * 30 colors_c = [ROUTE_COLOR_MAP[n] for n in T_route["航路名"]] ax.scatter(xs, ys, s=sizes, c=colors_c, edgecolor="black", linewidth=0.8, alpha=0.85, zorder=5) for i, n in enumerate(T_route["航路名"]): short = n.replace("コース（仮）", "")[:8] ax.annotate(short, (xs[i], ys[i]), xytext=(7, 5), textcoords="offset points", fontsize=8.5, fontweight="bold") ax.set_xlabel("観光資源密度 (= ビューポイント数 / 寄港地数)") ax.set_ylabel("しまたび Jaccard 類似度") ax.set_title("(c) 観光資源密度 × しまたび Jaccard\n" "右上=既存路線多い / 左下=新規寄港中心", fontsize=10.5) ax.set_xlim(0, max(xs) * 1.2) ax.set_ylim(0, max(ys) * 1.2) ax.grid(True, linestyle="--", alpha=0.3) fig.suptitle("図 2 (RQ1): 7 航路の構造プロファイル — 距離 / 新規率 / 観光資源・既存接続", fontsize=12.5, y=1.02) plt.tight_layout() save_fig("L82_fig2_route_profile.png") # ---- 図 3 (RQ1): 寄港地の航路寄港回数ランキング (どの港が「中継ハブ」か) ---- print(" fig3: 港寄港回数ランキング", flush=True) fig, ax = plt.subplots(figsize=(11, 8)) dpf = T_port_freq.copy() # 短縮名 dpf["short_name"] = dpf["寄港地名"].apply( lambda s: s.replace("（広島港宇品旅客ターミナル）", "") .replace("(中央桟橋ターミナル)", "") .replace("（駅前桟橋旅客船）", "") .replace("（瀬戸田港中央桟橋）", "") .replace("(内港旅客船ターミナル)", "")[:14] ) dpf2 = dpf.sort_values("航路寄港回数", ascending=True) yy = np.arange(len(dpf2)) colors_p = ["#cf222e" if x == "新規区間" else "#0969da" for x in dpf2["区間情報"]] edges = ["#ff8c00" if v else "#222" for v in dpf2["ビューポイント"]] ax.barh(yy, dpf2["航路寄港回数"], color=colors_p, edgecolor=edges, linewidth=1.4, alpha=0.85) for i, (cnt, name_short, vp) in enumerate(zip( dpf2["航路寄港回数"], dpf2["short_name"], dpf2["ビューポイント"])): star = " ★" if vp else "" ax.text(cnt + 0.05, i, f"{int(cnt)} 航路{star}", va="center", fontsize=9) ax.set_yticks(yy) ax.set_yticklabels(dpf2["short_name"], fontsize=9) ax.set_xlabel("寄港された航路数 (max=7)") ax.set_title(f"図 3 (RQ1): {n_ports} 寄港地 × 寄港回数ランキング — " f"中継ハブを特定 (★=観光資源)", fontsize=11) legend_elements = [ Patch(facecolor="#cf222e", edgecolor="#222", label="新規区間"), Patch(facecolor="#0969da", edgecolor="#222", label="既存区間"), Patch(facecolor="white", edgecolor="#ff8c00", label="観光資源 (ビューポイント)"), ] ax.legend(handles=legend_elements, loc="lower right", fontsize=9.5, framealpha=0.95) ax.grid(True, axis="x", linestyle="--", alpha=0.3) plt.tight_layout() save_fig("L82_fig3_port_frequency.png") # ---- 図 4 (RQ2): 客数 small multiples (運航日 × 航路) ---- print(" fig4: 運航日×航路 客数 small multiples", flush=True) n_days_unique = T_trip["運航日"].dt.date.nunique() # sharey=False: 各パネルで運航航路が異なるため、 y軸ラベルは個別に fig, axes = plt.subplots(3, 3, figsize=(15, 11), sharex=True, sharey=False) axes_flat = axes.flatten() days_sorted = sorted(T_trip["運航日"].dt.date.unique()) for i, day in enumerate(days_sorted): ax = axes_flat[i] sub = T_trip[T_trip["運航日"].dt.date == day].sort_values("航路名") if len(sub) == 0: ax.axis('off') continue short_names = [n.replace("コース（仮）", "")[:10] for n in sub["航路名"]] colors_v = [ROUTE_COLOR_MAP[n] for n in sub["航路名"]] yy = np.arange(len(sub)) ax.barh(yy, sub["観光客数（人）"], color=colors_v, edgecolor="#222", linewidth=0.5, alpha=0.85) for j, (v, dir_) in enumerate(zip(sub["観光客数（人）"], sub["方向"])): sym = "E" if dir_ == "東向" else "W" ax.text(v + 1, j, f"{int(v)}人[{sym}]", va="center", fontsize=9) ax.set_yticks(yy) ax.set_yticklabels(short_names, fontsize=8.5) weekday = day.strftime("%a") # ターン番号 (週単位の運航回数) if day.day <= 8: turn_no = "第1ターン" title_color = "#0969da" elif day.day <= 15: turn_no = "第2ターン" title_color = "#16a34a" else: turn_no = "第3ターン" title_color = "#cf222e" ax.set_title(f"{day} ({weekday}, {turn_no})", fontsize=10, color=title_color) ax.set_xlim(0, 80) ax.grid(True, axis="x", linestyle="--", alpha=0.3) # 余ったパネルを消す for j in range(len(days_sorted), len(axes_flat)): axes_flat[j].axis('off') fig.supxlabel("観光客数 (人)", fontsize=11) plt.suptitle(f"図 4 (RQ2): 9 運航日 × 7 航路 = 15 航海の観光客数 small multiples\n" f"記号: [E] = 東向, [W] = 西向 / " f"色: 青=第1ターン (12/06-08), 緑=第2ターン (12/13-15), 赤=第3ターン (12/20-22)", fontsize=12, y=0.998) plt.tight_layout() save_fig("L82_fig4_voyage_smallmult.png") # ---- 図 5 (RQ2): 運航回別 客数 + 平日週末 + ルートタイプ ---- print(" fig5: 運航回・曜日・ルートタイプ", flush=True) fig, axes = plt.subplots(1, 3, figsize=(16, 5.5)) # (a) 運航回別 box ax = axes[0] voyage_ord = sorted(T_trip["運航回"].unique()) data_box = [T_trip[T_trip["運航回"] == k]["観光客数（人）"].values for k in voyage_ord] bp = ax.boxplot(data_box, positions=voyage_ord, widths=0.5, patch_artist=True, boxprops=dict(facecolor="#cf222e", alpha=0.45), medianprops=dict(color="#222", linewidth=2)) # 平均をオーバープロット for i, k in enumerate(voyage_ord): avg = T_voyage_idx.loc[T_voyage_idx["運航回"] == k, "平均客数"].iloc[0] n = T_voyage_idx.loc[T_voyage_idx["運航回"] == k, "航海数"].iloc[0] ax.scatter(k, avg, s=120, color="#ffd700", edgecolor="#222", marker="D", zorder=5) ax.text(k, avg + 4, f"平均{avg:.0f}\n(n={int(n)})", ha="center", fontsize=9) ax.set_xticks(voyage_ord) ax.set_xticklabels([f"第{k}回" for k in voyage_ord], fontsize=10) ax.set_ylabel("1 航海観光客数 (人)") ax.set_title(f"(a) 運航回別 客数分布\n" f"H3: 第1→第2 減少率 {delta_pct:.1f}% " f"({'支持' if h3_ok else '反証'})", fontsize=10.5) ax.grid(True, axis="y", linestyle="--", alpha=0.3) ax.set_ylim(0, 85) # (b) 3 ターン別 (12/06-08, 12/13-15, 12/20-22) # 全 15 航海とも平日 Tue-Thu のため、 平日週末では分解できない → 3 ターンで区分 ax = axes[1] turn_order = ["第1ターン (12/06-08)", "第2ターン (12/13-15)", "第3ターン (12/20-22)"] wd_data = T_turn.set_index("ターン").reindex(turn_order) wd_data = wd_data.dropna() xs = np.arange(len(wd_data)) colors_w = ["#0969da", "#16a34a", "#cf222e"][:len(wd_data)] ax.bar(xs, wd_data["平均客数"], color=colors_w, edgecolor="#222", linewidth=0.8, alpha=0.85) for i, (avg, n) in enumerate(zip(wd_data["平均客数"], wd_data["航海数"])): ax.text(i, avg + 1, f"{avg:.1f}\n(n={int(n)})", ha="center", fontsize=9.5) ax.set_xticks(xs) ax.set_xticklabels([t.split(" ")[0] for t in wd_data.index], fontsize=10) ax.set_ylabel("1 航海平均観光客数 (人)") ax.set_title("(b) 3 ターン別 平均\n" "(全 15 航海とも平日 Tue-Thu のため週末区分は無効)", fontsize=10.5) ax.grid(True, axis="y", linestyle="--", alpha=0.3) ax.set_ylim(0, 80) # (c) 東向 vs 西向 (全 7 航路は片道型なので、 方向別) ax = axes[2] rt_data = T_routetype.set_index("ルートタイプ").reindex(["東向", "西向"]) rt_data = rt_data.dropna() xs = np.arange(len(rt_data)) colors_r = ["#cf222e", "#0969da"][:len(rt_data)] ax.bar(xs, rt_data["平均客数"], color=colors_r, edgecolor="#222", linewidth=0.8, alpha=0.85) for i, (avg, n) in enumerate(zip(rt_data["平均客数"], rt_data["航海数"])): ax.text(i, avg + 1, f"{avg:.1f}\n(n={int(n)})", ha="center", fontsize=9.5) ax.set_xticks(xs) ax.set_xticklabels(rt_data.index, fontsize=10) ax.set_ylabel("1 航海平均観光客数 (人)") ax.set_title("(c) 東向 vs 西向 (全 7 航路は片道型)", fontsize=10.5) ax.grid(True, axis="y", linestyle="--", alpha=0.3) ax.set_ylim(0, 80) fig.suptitle(f"図 5 (RQ2): 客数の 3 角度分解 — " f"全 15 航海 平均 {mean_v:.1f} 人 / CV {cv_v:.2f} " f"({'安定' if h4_ok else '中変動'})", fontsize=12.5, y=1.02) plt.tight_layout() save_fig("L82_fig5_voyage_three_angles.png") # ---- 図 6 (RQ3): 推定収益 + 収容率 + しまたび比較 ---- print(" fig6: 採算性ダッシュボード", flush=True) fig, axes = plt.subplots(1, 2, figsize=(15, 5.5)) # (a) 収容率 + 1 航海推定収益 ax = axes[0] rev = T_revenue.copy().sort_values("推定収容率_%", ascending=True) short_r = [n.replace("コース（仮）", "")[:14] for n in rev["航路名"]] colors_r = [ROUTE_COLOR_MAP[n] for n in rev["航路名"]] yy = np.arange(len(rev)) ax.barh(yy, rev["推定収容率_%"], color=colors_r, edgecolor="#222", linewidth=0.8, alpha=0.85) for i, (occ, avg, rev_v) in enumerate(zip( rev["推定収容率_%"], rev["平均客数"], rev["推定収益_万円_per航海"])): ax.text(occ + 0.5, i, f"{occ:.1f}% (平均{avg:.0f}人 / " f"推定 {rev_v:.0f}万円/航海)", va="center", fontsize=8.5) ax.set_yticks(yy) ax.set_yticklabels(short_r, fontsize=9) ax.set_xlabel(f"推定収容率 (= 平均客数 / {ASSUMED_CAPACITY} 人, %)") ax.set_xlim(0, 105) ax.axvline(80, color="#cf6f00", linestyle="--", linewidth=1.2, label="経営観光船 損益分岐目安 ~80%") ax.set_title(f"(a) 航路別 推定収容率 + 推定収益\n" f"仮定単価 {ASSUMED_FARE:,} 円/人, 仮定定員 {ASSUMED_CAPACITY} 人", fontsize=10.5) ax.legend(loc="lower right", fontsize=9) ax.grid(True, axis="x", linestyle="--", alpha=0.3) # (b) しまたび月平均 vs モニター 1 航海 ax = axes[1] labels_c = ["しまたび\n月平均\n(/月/全14港)", "しまたび\n2022/12月\n(全14港)", "モニター\n15航海合計\n(2022/12月)", "モニター\n1航海平均"] values_c = [shi_avg_per_month, shi_dec, monitor_total, monitor_per_voyage] colors_c = ["#0969da", "#16a34a", "#cf222e", "#fb8500"] xs = np.arange(len(labels_c)) ax.bar(xs, values_c, color=colors_c, edgecolor="#222", linewidth=0.8, alpha=0.85) for i, v in enumerate(values_c): if v >= 1000: ax.text(i, v * 1.02, f"{v:,.0f} 人", ha="center", fontsize=9.5, fontweight="bold") else: ax.text(i, v + max(values_c) * 0.01, f"{v:.1f} 人", ha="center", fontsize=9.5, fontweight="bold") ax.set_xticks(xs) ax.set_xticklabels(labels_c, fontsize=9) ax.set_yscale("log") ax.set_ylabel("乗降客数 (人, 対数軸)") ax.set_title(f"(b) スケール比較\n" f"モニター/しまたび12月 = {ratio_to_dec*100:.1f}%", fontsize=10.5) ax.grid(True, axis="y", which="both", linestyle="--", alpha=0.3) fig.suptitle(f"図 6 (RQ3): 観光商品としての採算性 + スケール比較 — " f"既存定期航路 (しまたび) との対比", fontsize=12.5, y=1.02) plt.tight_layout() save_fig("L82_fig6_revenue_dashboard.png") # ---- 図 7 (RQ3): 定期化適性スコア (地理可視化) ---- print(" fig7: 定期化適性スコア地図", flush=True) fig, axes = plt.subplots(1, 2, figsize=(17, 7.5)) # (a) スコア棒グラフ (本記事独自指標) ax = axes[0] ds = T_periodic.sort_values("定期化適性スコア", ascending=True) short_p = [n.replace("コース（仮）", "")[:14] for n in ds["航路名"]] yy = np.arange(len(ds)) # 上位 2 航路を強調色 (cf222e) top2_set = set(top2_routes) colors_pp = ["#cf222e" if n in top2_set else "#888" for n in ds["航路名"]] ax.barh(yy, ds["定期化適性スコア"], color=colors_pp, edgecolor="#222", linewidth=0.8, alpha=0.85) for i, (sc, avg, nb, vb, df_) in enumerate(zip( ds["定期化適性スコア"], ds["平均客数"], ds["新規ボーナス"], ds["観光資源ボーナス"], ds["距離係数"])): ax.text(sc + 0.5, i, f"{sc:.1f} (×新規{nb:.2f} / ×観光{vb:.2f} / ÷{df_:.2f}km)", va="center", fontsize=8.5) ax.set_yticks(yy) ax.set_yticklabels(short_p, fontsize=9) ax.set_xlabel("定期化適性スコア (= 平均客数 × 新規ボーナス × 観光ボーナス / 距離係数)") ax.set_title(f"(a) 7 航路の定期化適性スコア\n" f"赤 = 上位 2 = 定期化候補 (H5: {'支持' if h5_ok else '部分支持'})", fontsize=10.5) ax.grid(True, axis="x", linestyle="--", alpha=0.3) # (b) スコアでルート色分けした地図 ax = axes[1] admin.plot(ax=ax, color="#fff8e8", edgecolor="#888", linewidth=0.4, alpha=0.55) # 全寄港地灰色背景 ax.scatter(gdf_ports.geometry.x, gdf_ports.geometry.y, s=40, c="#cccccc", edgecolor="#888", alpha=0.7, zorder=3) # スコア順位で色付け (赤強 → 黄薄) score_norm = (T_periodic["定期化適性スコア"] - T_periodic["定期化適性スコア"].min()) / \ max(T_periodic["定期化適性スコア"].max() - T_periodic["定期化適性スコア"].min(), 1e-6) cmap = plt.cm.RdYlBu_r score_to_color = {n: cmap(s) for n, s in zip(T_periodic["航路名"], score_norm)} for name in sorted(route_geom.keys()): info = route_geom[name] coords = info["coords"] if len(coords) < 2: continue pts = gpd.GeoSeries([Point(lon, lat) for lon, lat in coords], crs="EPSG:4326").to_crs(TARGET_CRS) xs2 = [p.x for p in pts] ys2 = [p.y for p in pts] sc = T_periodic.loc[T_periodic["航路名"] == name, "定期化適性スコア"].iloc[0] rk = T_periodic.loc[T_periodic["航路名"] == name, "順位"].iloc[0] is_top = name in top2_set lw = 3.5 if is_top else 1.6 al = 0.9 if is_top else 0.5 ax.plot(xs2, ys2, color=score_to_color[name], linewidth=lw, alpha=al, zorder=2 if not is_top else 4, label=f"#{rk} {name[:14]} (sc={sc:.0f})") ax.set_xlim(gdf_x_min, gdf_x_max) ax.set_ylim(gdf_y_min, gdf_y_max) ax.set_aspect("equal") ax.set_title(f"(b) 7 航路の地理可視化 — 線色 = 定期化適性スコア\n" f"太線 + 濃色 = 上位 2 (定期化候補)", fontsize=10.5) ax.legend(loc="upper left", fontsize=7.5, framealpha=0.92) ax.grid(True, linestyle="--", alpha=0.3) ax.set_xlabel("X (m, EPSG:6671)") ax.set_ylabel("Y (m, EPSG:6671)") fig.suptitle(f"図 7 (RQ3): 定期化適性スコア (本記事独自指標) — " f"上位 = {top2_routes[0][:14]} / {top2_routes[1][:14]}", fontsize=12.5, y=1.02) plt.tight_layout() save_fig("L82_fig7_periodic_score.png") print(f" ({time.time()-t9:.2f}s)", flush=True) # ============================================================================ # 10. HTML 組み立て # ============================================================================ print("\n[10] HTML 組み立て", flush=True) t10 = time.time() def df_to_html(df, idx=False): return df.to_html(index=idx, border=0, escape=False) # 個別変数 n_new_ports = (ports_master["区間情報"] == "新規区間").sum() n_old_ports = (ports_master["区間情報"] == "既存区間").sum() top_score_route = T_periodic.iloc[0]["航路名"] top_score_value = T_periodic.iloc[0]["定期化適性スコア"] top_dist_route = T_route.iloc[0]["航路名"] top_dist_value = T_route.iloc[0]["総距離_km"] hub_port = T_port_freq.iloc[0]["寄港地名"] hub_freq = T_port_freq.iloc[0]["航路寄港回数"] # ----- セクション 1: 学習目標と問い ----- sec1 = f"""
【本記事のスタイル: 単独 3RQ 形式 (Format B) / 観光社会実験のプロダクト分析】
1 記事で 3 つの独立した研究角度を並列に進める。 各 RQ は単独で読める単元を成し、 全体で「観光モニタークルーズ = プロダクト 試作品」 の 3 軸プロファイル ( 航路設計 / 乗船客特性 / 商品評価) を 形成する。
【X04 改訂版との重複回避】
X04 改訂版 (2026-05) は同 dataset を「広島県沿岸クルーズ需要の港別・季節別 構造分析」 として扱い、 集中度 (Gini/HHI/CR), ローレンツ曲線, 季節係数, ルート 4 タイプ分類で 9 図 12 表を構成した。 本 L82 は X04 が触れなかった5 つの独自切り口に絞る: X04 で計算済みの Gini / HHI / CR / ローレンツ / 季節係数 / 4 タイプ分類は本記事では再計算しない。 重複しないよう、 「ルート設計のミクロ構造」 と「観光商品評価」 に絞り込む。

本記事は、 DoBoX が公開する せとうちモニタークルーズ実施結果 (dataset_id = {DATASET_ID})を、 観光社会実験のプロダクト分析 として読み解く。 2022 年 12 月 6 日〜22 日の9 運航日 × {n_routes} 航路 = {n_voyages} 航海 / {n_ports} 寄港地 (Haversine 総距離 {total_dist_all:.0f} km) のデータを、 (1) 航路の地理設計、 (2) 乗船客の 回数効果、 (3) 商品としての採算性 の 3 軸で深掘りする。

用語の独自定義 (このレッスン専用)

研究の問い (3 RQ)

  1. RQ1 (主研究) — モニタークルーズ {n_routes} 航路の航路設計 (ルート構造)はどう描けるか? Haversine 距離・寄港地数・新規率・戻り型/片道型・観光資源密度・ しまたび Jaccard 類似度の 6 軸で全件記述する。
  2. RQ2 (副研究 1) — 1 航海あたり乗船客数の運航回数効果は あるか? 第 1 回 vs 第 2 回 vs 第 3 回 / 平日 vs 週末 / 戻り型 vs 片道型 の 3 角度で集計し、 物珍しさ効果と常連化効果を識別。
  3. RQ3 (副研究 2) — モニタークルーズの観光商品としての評価は どうか? 推定収容率・推定収益・定期化適性スコア・既存しまたびとの 規模比較 の 4 軸で量化する。

仮説 (5 個)

到達点

""" # ----- セクション 2: 使用データ ----- sec2 = f"""

本記事の主データは DoBoX dataset {DATASET_ID}「せとうちモニタークルーズ 実施結果」のみ。 X04 改訂版で resource 39782 + 39781 を使用済みだが、 本記事は第 3 のリソース 39783 (ビューポイント写真情報)のメタ情報も 含めて解析する。

主データ (本記事の研究対象)

項目内容備考
データセット せとうちモニタークルーズ実施結果 (DoBoX #{DATASET_ID}) 広島県 観光連盟 提供
形式CSV (cp932) resource 39783 のみ ZIP (内部ビューポイント JPG 14 枚)
resource 39782 航路 × 寄港地 × 運航日 {df_daily.shape[0]} 行 × {df_daily.shape[1]} 列 (~{DAILY_CSV.stat().st_size//1024} KB)
resource 39781 寄港地マスタ {df_ports.shape[0]} 行 × {df_ports.shape[1]} 列 (~{PORTS_CSV.stat().st_size//1024} KB)
resource 39783 ビューポイント写真 (X04 で未使用) ZIP, 14 枚 JPG (本記事はメタ情報のみ利用)
運航期間 2022年12月6日 〜 12月22日 {n_dates} 運航日 × 7 航路 = {n_voyages} 航海
寄港地数 {n_ports} 港 (新規 {n_new_ports} 港 + 既存 {n_old_ports} 港) 緯度経度あり = POINT 23 点
観光資源 {n_viewpoints} 港 (ビューポイント写真あり) resource 39783 メタ情報
ライセンスクリエイティブ・コモンズ表示 (CC-BY) DoBoX オープンデータ

従属参照データ (RQ1 重複検出 + RQ3 規模比較)

用途データ件数備考
RQ1 Jaccard 瀬戸内しまたびライン (#1282) 14 桟橋 × 9 ヶ月 X04 で主データ。 本記事は寄港地集合のみ参照
RQ3 規模比較 同 #1282 月次乗降客数 4-12月の乗降客数 モニター 1 航海平均との対比 (X04 とは異なる切り口)
地図描画 L44 行政界 ディゾルブ 27 市町 既キャッシュ admin_diss.gpkg

形式特性の注意点

""" # ----- セクション 3: ダウンロード ----- sec3 = f"""

本記事の再現に必要なすべてを直リンクで提供する。 HTML だけ読めば学習者が完全再現できることが目標 (要件 A)。

生データ (DoBoX, 主データ)

このスクリプト本体

従属参照 (本記事は読込みのみ、 取得不要)

中間 CSV (本記事生成、 再利用可)

図 (PNG, 直 DL 可)

""" # ----- セクション 4: RQ1 ----- sec4_code = ''' import pandas as pd import geopandas as gpd import numpy as np from shapely.geometry import Point # (1) CSV 読込み (cp932 = Shift-JIS、 DoBoX のクルーズ系は cp932 が標準) df_daily = pd.read_csv( "data/extras/L82_setouchi_monitor_cruise/monitor_cruise_daily.csv", encoding="cp932") df_ports = pd.read_csv( "data/extras/L82_setouchi_monitor_cruise/monitor_cruise_ports.csv", encoding="cp932") # (2) Haversine 距離 (球面三角法、 numpy ベクトル化で高速) def haversine_km(lat1, lon1, lat2, lon2): R = 6371.0 # 地球半径 km p1, p2 = np.radians(lat1), np.radians(lat2) dp, dl = np.radians(lat2 - lat1), np.radians(lon2 - lon1) a = np.sin(dp/2)**2 + np.cos(p1)*np.cos(p2)*np.sin(dl/2)**2 return 2 * R * np.arcsin(np.sqrt(a)) # (3) 航路ごとに寄港順序を組み立て、 区間距離を計算 for name, g in df_daily.groupby("航路名"): seq = g.drop_duplicates("寄港地ID").sort_values("寄港桟橋順") lats = seq["緯度"].values lons = seq["経度"].values seg_dists = haversine_km(lats[:-1], lons[:-1], lats[1:], lons[1:]) print(f"{name}: 寄港地{len(seq)}, 総距離{seg_dists.sum():.1f}km") # (4) Jaccard 類似度 (既存しまたびとの寄港地集合比較) df_shi = pd.read_csv("data/extras/shimatabi_monthly.csv", encoding="cp932").dropna(subset=["寄港地（桟橋名）"]) shi_set = set(df_shi["寄港地（桟橋名）"].unique()) def jaccard(a, b): a, b = set(a), set(b) return len(a & b) / max(len(a | b), 1) # 各 monitor 航路の寄港地集合と shi_set の Jaccard for name, g in df_daily.groupby("航路名"): monitor_set = set(g["寄港地（桟橋名）"].unique()) j = jaccard(monitor_set, shi_set) print(f"{name}: しまたびJaccard = {j:.2f}") # (5) GeoDataFrame 化 + 平面直角第 III 系 (m 単位) 投影 gdf_ports = gpd.GeoDataFrame(df_ports, geometry=[Point(lon, lat) for lon, lat in zip(df_ports["経度"], df_ports["緯度"])], crs="EPSG:4326").to_crs("EPSG:6671") ''' sec4 = f"""

狙い (RQ1)

RQ1 では「モニタークルーズ {n_routes} 航路の航路設計」 を初めて 系統的に記述する。 Haversine 距離(球面三角法) で総距離・平均区間 距離・最大最小区間を計算し、 各航路の方向 (東向/西向)・東西ペア構造と 既存しまたびライン (#1282) との Jaccard 類似度を算出する。 さらに resource 39783 (X04 で未使用) のビューポイント写真情報から 観光資源密度を計算する。 H1 (総距離 80-120 km 帯集中), H2 (東西ペア ≥ 2 ペア) を量的に検証する。

手法 — 5 ステップ

  1. STEP 1: 3 CSV 読込み + cp932 デコード
    DoBoX の本データは cp932 (Shift-JIS) 標準。 read_csv(..., encoding="cp932")。
  2. STEP 2: 航路ジオメトリ生成
    各航路で drop_duplicates("寄港地ID").sort_values("寄港桟橋順") で寄港順序を確定。 寄港地 ID と緯度経度をリストアップ。
  3. STEP 3: Haversine ベクトル化距離
    numpy ベクトル化で連続 2 点間の球面距離を一括計算。 対象 7 航路 × 平均 9 区間 = 約 63 区間で 1 ms 未満。
  4. STEP 4: Jaccard 類似度 + 方向タグ + 東西ペア判定
    モニター各航路の寄港地集合 vs しまたび 14 港集合。 Python の set 演算で 1 行。 航路名末尾の「（東）」 「（西）」 から方向を抽出し、 ペアキー (sorted({"始発", "終着"})) で東西ペアを識別する。
  5. STEP 5: 観光資源密度 + 距離ベース新規率
    ビューポイント写真フラグ ({n_viewpoints} 港) を航路ごとに集計。 区間情報 (新規/既存) を Haversine 距離で重み付けし、 距離ベース新規率を算出。

実装

狙いを踏まえた実装。 Haversine ベクトル化 + Jaccard set 演算の組み合わせ。

{code(sec4_code)}

結果と読み取り

(a) 7 航路ルートマップ + 観光資源 (図 1)

なぜこの図か: {n_routes} 航路 × {n_ports} 寄港地という小さな 母集団は、 1 枚の地図で全件並列表示するのが最も情報密度が高い。 寄港地の新規 (赤丸) / 既存 (青丸)に加え、 観光資源 (黄星) を 3 重符号化することで、 1 枚で「どこに何を新設し、 どんな観光資源を経由するか」 を視覚化する (要件 T)。

{figure("assets/L82_fig1_route_map.png", f"図 1: モニタークルーズ {n_routes} 航路ルートマップ")}

図 1 から読み取れること:

(b) 7 航路の構造プロファイル (図 2)

なぜこの図か: H1 (総距離分散) を直感するには、 (a) 総距離ランキング、 (b) 新規率ランキング、 (c) 観光資源密度 vs しまたび Jaccard 散布の 3 パネル並列が最適。 各航路の戻り型/片道型を色分けで識別する。

{figure("assets/L82_fig2_route_profile.png", f"図 2: 7 航路の距離・新規率・観光資源・既存接続の 3 軸プロファイル")}

{df_to_html(T_route)}

図 2 / 表から読み取れること:

(c) 寄港地の航路寄港回数ランキング (図 3)

なぜこの図か: 「どの港が中継ハブとして 7 航路に頻繁に組み込まれて いるか」 を 1 枚で示す。 港の重要度を寄港頻度で量的にランキングし、 新規/既存 + 観光資源フラグを色と枠線で 3 重表示する。

{figure("assets/L82_fig3_port_frequency.png", f"図 3: {n_ports} 寄港地 × 寄港回数ランキング")}

{df_to_html(T_port_freq)}

図 3 / 表から読み取れること:

距離統計サマリ:

{df_to_html(T_dist_summary)}

始発・終着 役割表:

{df_to_html(T_role)}

""" # ----- セクション 5: RQ2 ----- sec5_code = ''' import pandas as pd # (1) 1 航海 = 1 (運航日, 航路) の単位に集約 (寄港地行を 1 行に統合) trip = (df_daily.drop_duplicates(["運航日", "航路名"]) [["運航日", "航路名", "観光客数（人）"]] .sort_values(["運航日", "航路名"]) .reset_index(drop=True)) # (2) 各航路内の運航回数 (1 回目, 2 回目, 3 回目) trip = trip.sort_values(["航路名", "運航日"]) trip["運航回"] = trip.groupby("航路名").cumcount() + 1 # (3) 平日 / 週末 タグ trip["運航日"] = pd.to_datetime(trip["運航日"]) trip["平日週末"] = trip["運航日"].dt.weekday.apply( lambda w: "週末" if w >= 5 else "平日") # (4) 戻り型/片道型タグを RQ1 ジオメトリから付与 loop_map = {name: info["is_loop"] for name, info in route_geom.items()} trip["戻り型"] = trip["航路名"].map(loop_map) # (5) 客数 CV (変動係数) — 商品としての安定性指標 visitors = trip["観光客数（人）"].values.astype(float) mean_v = visitors.mean() std_v = visitors.std(ddof=0) cv = std_v / mean_v print(f"全 15 航海 平均 {mean_v:.1f} 人, std {std_v:.2f}, CV {cv:.3f}") # (6) 運航回別 平均 print(trip.groupby("運航回")["観光客数（人）"].agg(["count", "mean", "median"])) ''' sec5 = f"""

狙い (RQ2)

RQ2 では「1 航海あたり乗船客数の運航回数効果」 を検証する。 モニタークルーズの 9 運航日は 2 週末を含む 3 ターン (12/06-08, 12/13-15, 12/20-22) で同一航路を 2-3 回反復運航する設計。 第 1 回 vs 第 2 回 vs 第 3 回の客数を比較することで、 「物珍しさ効果」 (= 第 1 回が多い) と「常連化効果」 (= 第 N 回が多い) を 識別する。 H3 (第 2 運航客数 -15%) は物珍しさ効果の量化、 H4 (CV ≤ 0.20) は商品としての客数安定性を量化。

手法 — 6 ステップ

  1. STEP 1: 1 航海 = 1 (運航日, 航路) 集約
    raw データは 1 寄港地 = 1 行で複製されているため、 drop_duplicates(["運航日", "航路名"]) で 1 航海 1 行に縮約。 {df_daily.shape[0]} 行 → {n_voyages} 行。
  2. STEP 2: 運航回ラベル付与
    各航路内で運航日順に cumcount() + 1 で 1, 2, 3 を付与。
  3. STEP 3: ターン別タグ (12/06-08 / 12/13-15 / 12/20-22)
    dt.weekday で曜日を取得すると、 全 15 航海が平日 Tue-Thu で週末区分が無効であることが分かる。 代わりに運航スケジュールが 3 ターン構成 (毎週 1 ターン) に分かれていることを利用し、 第1/第2/第3 ターンのタグを付与する。
  4. STEP 4: 戻り型/片道型タグ
    RQ1 のジオメトリ辞書から is_loop を引いて付与。
  5. STEP 5: 客数 CV (変動係数) 計算
    std / mean。 無次元なので異なる指標を比較できる。
  6. STEP 6: 3 軸集計 (運航回 / 平日週末 / ルートタイプ)
    groupby + aggで 3 つの平均客数表を生成。 H3, H4 を量的に検証。

実装

運航回ラベル付与と CV 計算が中心。 6 ステップで集計表 4 個 + ピボット 1 個を生成。

{code(sec5_code)}

結果と読み取り

(a) 9 運航日 × 7 航路 客数 small multiples (図 4)

なぜこの図か: 9 運航日 × 7 航路 = 15 航海の客数を、 運航日ごとに 1 パネルの small multiplesで全件並列表示することで、 「いつ・どの航路が・ どれだけ集客したか」 を 1 図で読み取れる。 平日 (青タイトル) と週末 (赤 タイトル) の差、 戻り型 (↺) と片道型 (→) の差を視覚的に区別 (要件 T)。

{figure("assets/L82_fig4_voyage_smallmult.png", f"図 4: 9 運航日 × 7 航路 = {n_voyages} 航海 客数 small multiples")}

運航日別 集計:

{df_to_html(T_daily_total)}

図 4 / 表から読み取れること:

(b) 客数の 3 角度分解 (図 5)

なぜこの図か: H3 (運航回数効果) と H4 (CV) を直感するには、 (a) 運航回別 box plot, (b) 平日週末別 棒, (c) 戻り型/片道型別 棒 の 3 パネル並列が最適。 box plot で分布全体を見せ、 棒で平均値を強調する。

{figure("assets/L82_fig5_voyage_three_angles.png", f"図 5: 客数の 3 角度分解 — 運航回 / 平日週末 / ルートタイプ")}

15 航海 客数統計:

{df_to_html(T_pas_summary)}

運航回別:

{df_to_html(T_voyage_idx)}

3 ターン別 (全 15 航海とも平日 Tue-Thu のため週末区分は無効):

{df_to_html(T_turn)}

方向 (東向 / 西向) 別:

{df_to_html(T_routetype)}

図 5 / 表から読み取れること:

""" # ----- セクション 6: RQ3 ----- sec6_code = ''' import pandas as pd # 仮定値 (公開データに無いため明示) ASSUMED_CAPACITY = 70 # 1 航海定員 (観測最大値ベース仮定) ASSUMED_FARE = 5000 # 仮定単価 円/人 (類似観光船相場) # (1) 航路別 1 航海平均客数 + 推定収容率 + 推定収益 revenue_records = [] for name in T_route["航路名"]: sub = trip[trip["航路名"] == name] avg = sub["観光客数（人）"].mean() n_v = len(sub) occ = avg / ASSUMED_CAPACITY * 100 rev_per_v = avg * ASSUMED_FARE revenue_records.append({ "航路名": name, "平均客数": round(avg, 1), "推定収容率_%": round(occ, 1), "推定収益_万円_per航海": round(rev_per_v / 10000, 1), }) # (2) 定期化適性スコア (本記事独自指標) # = 平均客数 × 新規ボーナス (1.0 + 新規率 × 0.3) # × 観光資源ボーナス (1.0 + 観光資源密度 × 0.2) # / 距離係数 (max(総距離 / 100, 0.5)) def periodic_score(avg, new_ratio, view_density, total_km): new_bonus = 1.0 + new_ratio * 0.3 view_bonus = 1.0 + view_density * 0.2 dist_factor = max(total_km / 100.0, 0.5) return avg * new_bonus * view_bonus / dist_factor # (3) しまたび月平均との規模比較 df_shi_monthly = df_shi[month_cols].sum().sum() / 9 # 4-12月の月平均 ratio = monitor_total / df_shi_monthly print(f"モニター 12月総客数 / しまたび月平均 = {ratio*100:.1f}%") ''' sec6 = f"""

狙い (RQ3)

RQ3 では「モニタークルーズの観光商品としての評価」 を量化する。 本データには料金・船舶コストの情報が無いため、 仮定単価 {ASSUMED_FARE:,} 円/人 + 仮定定員 {ASSUMED_CAPACITY} 人を明示提示し、 推定収容率・推定収益・定期化適性スコア (本記事独自指標)を計算する。 これらは純利益の実数値ではなく規模指標・ランキング指標として位置付ける。 最後に既存しまたびとのスケール比較で、 モニター事業の「相対的位置」 を量化。

限界 (要件 J): 本記事の収益試算は純利益の計算ではない。 公開データに料金 (1 人あたり)・船舶運航コスト・人件費が含まれないため、 仮定単価 {ASSUMED_FARE:,} 円/人で1 航海規模感を提示する。 仮定値は類似観光船 (例: 瀬戸内海汽船銀河、 SEA SPICA) の公開料金から 推定。 厳密な経営評価は別途船社決算情報・観光振興課予算書の確認要。

手法 — 4 ステップ

  1. STEP 1: 推定収容率 + 推定収益 (1 航海あたり)
    平均客数 / 仮定定員 = 推定収容率。 平均客数 × 仮定単価 = 1 航海推定収益。 仮定値を明示することが研究倫理の要点。
  2. STEP 2: 定期化適性スコア (本記事独自指標)
    平均客数 × 新規ボーナス × 観光資源ボーナス / 距離係数 の積商で 1 軸ランキング化。 高客数 + 高新規性 + 高観光資源 + 距離効率 の 4 条件を統合。
  3. STEP 3: しまたび月平均との規模比較
    モニター 1 航海平均 vs しまたび月平均 / 月の桁差をログスケールで比較。 「観光社会実験は既存定期航路の何 % 規模なのか」 を量化。
  4. STEP 4: 政策的位置付け 8 観点
    実施期間・主体・目的・規模・成果・アンケート・後継・限界 の 8 観点で本事業の制度的文脈を整理。

実装

仮定値の明示提示 + 独自スコア生成 + 規模比較の 3 段構成。

{code(sec6_code)}

結果と読み取り

(a) 推定収容率 + しまたび規模比較 (図 6)

なぜこの図か: 採算性評価は (a) 1 航海ごとの収容率 (損益分岐目安 80% との比較)、 (b) 既存しまたびとの規模対比 の 2 パネルが最適。 収容率には推定収益も併記し、 (b) は対数軸で桁数差を視覚化する。

{figure("assets/L82_fig6_revenue_dashboard.png", f"図 6: 観光商品としての採算性 + しまたびとのスケール比較")}

航路別 推定収容率 + 推定収益:

{df_to_html(T_revenue)}

しまたびとの規模比較:

{df_to_html(T_compare)}

図 6 / 表から読み取れること:

(b) 定期化適性スコア (図 7)

なぜこの図か: H5 (定期化適性 上位 2 特定) の量的検証には、 (a) 7 航路のスコアランキング (3 ボーナスの内訳付き) と (b) 地理可視化 (線色 = スコア順) の 2 パネルが最適。 定量と空間配置を 1 枚で結合する (要件 T)。

{figure("assets/L82_fig7_periodic_score.png", f"図 7: 定期化適性スコア + 地理可視化")}

定期化適性スコア (本記事独自指標):

{df_to_html(T_periodic)}

政策的位置付け:

{df_to_html(T_policy)}

図 7 / 表から読み取れること:

""" # ----- セクション 7: 仮説検証総合 ----- sec7 = f"""

本記事で立てた 5 仮説の検証結果を一覧する。 観測値・判定・解釈の 3 列で、 RQ1/RQ2/RQ3 を横断する研究の結論部。

{df_to_html(T_hyp)}

仮説検証から読み取れる 3 つの研究的発見
""" # ----- セクション 8: 発展課題 ----- sec8 = f"""

本記事の結果Xから派生する新仮説Yと、 それを検証する 具体的課題Zを 3 つの発展経路に整理する (要件 E)。

発展課題 1: モニタークルーズの定期化追跡 (時系列拡張)

発展課題 2: 個票アンケートデータとの結合分析 (個人レベル)

発展課題 3: 海象データとの結合 (運航条件影響)

""" # ----- HTML 構築 ----- sections = [ ("学習目標と問い", sec1), ("使用データ", sec2), ("ダウンロード", sec3), (f"分析 1 (RQ1 主研究): モニタークルーズ {n_routes} 航路の航路設計", sec4), ("分析 2 (RQ2 副研究 1): 1 航海客数の運航回数効果", sec5), ("分析 3 (RQ3 副研究 2): 観光商品としての評価 — 採算性 + 定期化適性", sec6), ("仮説検証総合", sec7), ("発展課題", sec8), ] html = render_lesson( num=82, title=f"せとうちモニタークルーズ実施結果 単独 3 研究例分析 — " f"観光社会実験 {n_voyages} 航海の航路設計・乗船客特性・観光商品評価", tags=["観光", "クルーズ", "モニター", "社会実験", "Haversine", "ルート設計", "Format B", "geopandas"], time="60-90 分", level="リテラシ", data_label=f"DoBoX #{DATASET_ID} せとうちモニタークルーズ実施結果 ({n_voyages} 航海)", sections=sections, script_filename="L82_setouchi_monitor_cruise.py", ) OUT_HTML = LESSONS / "L82_setouchi_monitor_cruise.html" OUT_HTML.write_text(html, encoding="utf-8") print(f" → {OUT_HTML}", flush=True) print(f" ({time.time()-t10:.2f}s)", flush=True) # ============================================================================ # 11. 完了 # ============================================================================ elapsed = time.time() - t_all print(f"\n=== L82 完了 (合計 {elapsed:.2f}s) ===", flush=True) print(f" HTML: {OUT_HTML}", flush=True) print(f" CSV (中間): 14 個", flush=True) print(f" PNG (図): 7 個", flush=True)