もう OTP に頼らない！Python と GTFS で作る公共交通の到達圏地図

はじめに

公共交通データを用いた空間解析や到達圏地図（アイソクローナ）の作成は、都市計画や交通網分析において非常に強力な手法です。

その分析に不可欠なのが、公共交通の時刻表・経路情報を提供する世界標準のオープンフォーマット、GTFS（General Transit Feed Specification）です。

通常、GTFSを用いた高度な分析にはOpenTripPlanner（OTP）のような専用ツールが使われますが、環境構築のハードルが高いケースも少なくありません。
この記事では、OTPが導入できない環境を想定し、GTFSデータから直接、Pythonコードで到達圏地図を作成し、最終的に地理情報システム（GIS）で活用する手順を解説します。

作成した到達圏地図は、オープンソースで高機能なQGISにインポートし、他の施設データや人口統計データと重ね合わせることが可能です。コードを動かして終わりではなく、その先の現実課題の解決に繋げることまでを視野に入れています。

対象読者と実行環境

対象読者

• Pythonの基本的な文法や、Pandasライブラリでのデータ操作に慣れている方
• 公共交通データや地理空間情報（GIS）の分析に興味がある方
• QGISの基本的な操作をご存知の方（必須ではありません）

実行環境：Google Colaboratory

本記事のコードは、Googleアカウントさえあれば誰でも無料で利用できる Google Colaboratory (Colab) 上で作成・実行しました。

読者の皆様も、Colab環境で実行することを強く推奨します。

なぜColabを推奨するのか？

• 環境構築が一切不要: ご自身のPCにPythonや各種ライブラリをインストールする必要がなく、ブラウザだけで実行できます。
• ライブラリの簡単インストール: 本記事で利用するgeopandasやrasterioといったライブラリも、コードの先頭で!pip installコマンドを実行するだけで簡単に導入できます。
• Google Driveとの連携: GTFSデータや出力した地図ファイルをGoogle Driveに保存し、永続化するのが容易です。

以降で提示するコードブロックに登場する !pip install ... という ! から始まるコマンドは、Colab（やJupyter Notebook）上で外部コマンドを実行するための特別な記法です。

もしローカル環境で実行する場合は、お使いのターミナルやコマンドプロンプトで!を付けずにpip install ...を実行し、ライブラリをインストールしてください。

背景：なぜOTPを使わずに「自作」するのか？

公共交通の経路探索や到達圏分析において、OpenTripPlanner (OTP) は非常に高機能で優れたオープンソースソフトウェアです。では、なぜ本記事ではあえてOTPを使わずに、Pythonで自作するアプローチをとるのでしょうか？

それは、OTPの導入や利用が難しい、あるいは自作の方が適している具体的な場面が存在するからです。

OTPの導入が難しい・不向きなケース

OTPはJavaベースの本格的なサーバーアプリケーションであり、その高機能さゆえに、時として以下のような課題に直面します。

• 環境構築の壁: Java実行環境（JDK）のセットアップや、適切なメモリ割り当てが必要となり、環境構築のハードルがやや高めです。
• 実行環境の制約: セキュリティポリシー等の理由で、任意のソフトウェアをインストールできない社内・学内のサーバーや、Pythonスクリプトだけを配置したい軽量なコンテナ環境。
• アルゴリズムの柔軟性: 「特定の路線には独自の乗り換えペナルティを加えたい」といった、研究・分析用途での細かいアルゴリズムのカスタマイズが難しい場合があります。

本記事のアプローチが提供する価値

上記のような課題に対し、PythonとGTFSだけで到達圏を算出する本記事のアプローチには、以下のようなメリットがあります。

• 圧倒的な手軽さ: Python環境さえあれば、pip installだけで必要なライブラリを揃えるだけですぐに始められます。
• アルゴリズムの透明性と学習効果: 経路探索のロジックがすべてPythonコードとして記述されているため、処理内容を完全に把握・改変できます。
• QGISとのシームレスな連携: 算出した到達圏をGeoTIFFとして出力することで、使い慣れたQGISに取り込み、高度な空間分析へとスムーズに繋げることが可能です。

---

実装概要

今回のサンプルコードでは、以下の手順で到達圏地図を作成しています。

1. 入力パラメータの設定  
      - 出発日時、歩行速度、乗換時間などをパラメータとして定義。  
      - 各施設ごとに処理を実施できるよう、施設一覧（CSV）から出発地点をインポート。

2. GTFS データの解凍・読み込み  
      - stops.txt, stop_times.txt, trips.txt, calendar.txt などのファイルを読み込み、必要なデータ構造に整形。  
      - カレンダー情報をもとに、対象日・曜日のサービスのみ抽出。

3. 経路探索と到達時刻の計算  
      - 各停留所間の距離は Haversine 公式を用いて計算。  
      - 出発地点から各停留所への徒歩時間、バス乗車や乗換の時間を加味し、到達可能な時間（30分、45分、60分）での到達圏を算出。  
      - ここではヒープキュー（優先度付きキュー）を用いたダイクストラ法を実装しており、各停留所までの最短到着時間を求めています。

4. アイソクローナ（到達圏）ポリゴンの生成  
      - 到着時刻に応じた「余裕時間」を使って、各停留所周辺にバッファ（円形）を作成。  
      - 複数のバッファを統合して到達圏ポリゴンとしています。  
      - 投影変換（EPSG:3857 → EPSG:4326）を行い、GeoTIFF 用の座標系に整えます。

5. Folium によるインタラクティブマップ作成  
      - 各到達圏（30分、45分、60分）を異なる色で表示したマップを作成。  
      - 停留所情報も併せてプロットすることで、経路探索の結果を視覚的に確認可能にしています。

6. GeoTIFF のエクスポート  
      - rasterio を使い、作成した到達圏ポリゴンを GeoTIFF 形式で出力。  
      - この GeoTIFF を QGIS にインポートすることで、他の空間データ（例：施設情報、人口データ）と簡単に重ね合わせることが可能です。

コード全体の流れ

以下、主要なコード部分を解説します。

1. GTFS データの準備とカレンダー処理

# GTFS データの解凍と CSV 読み込み
import zipfile, os, pandas as pd, datetime
# ...（解凍処理、CSV 読み込みのコード）

# 入力日時をもとにカレンダーから有効なサービスを抽出
input_dt = datetime.datetime.strptime(input_datetime_str, "%Y/%m/%d %H:%M")
target_date_str = input_dt.strftime("%Y%m%d")
weekday_map = {0:"monday", 1:"tuesday", 2:"wednesday", 3:"thursday", 4:"friday", 5:"saturday", 6:"sunday"}
weekday_str = weekday_map[input_dt.weekday()]
# カレンダーからサービス ID を抽出
cal_active = calendar[
    (calendar['start_date'] <= int(target_date_str)) &
    (calendar['end_date'] >= int(target_date_str)) &
    (calendar[weekday_str] == 1)
]

2. 経路探索（到達圏の計算）

def haversine(lon1, lat1, lon2, lat2):
    # Haversine 公式による直線距離の計算
    import math
    R = 6371000  # 地球の半径 [m]
    phi1 = math.radians(lat1)
    phi2 = math.radians(lat2)
    dphi = math.radians(lat2 - lat1)
    dlambda = math.radians(lon2 - lon1)
    a = math.sin(dphi/2)**2 + math.cos(phi1)*math.cos(phi2)*math.sin(dlambda/2)**2
    return 2 * R * math.asin(math.sqrt(a))

def compute_isochrones_for_dep(dep):
    # 出発地点から各停留所までの到達時刻を計算し、30, 45, 60 分圏を求める
    # ヒープキューを用いた探索アルゴリズムで最短到着時間を更新
    # 最終的に、各圏に含まれる停留所からバッファを作成しポリゴンを生成
    # ...（コード全文）
    return iso_poly_30_local, iso_poly_45_local, iso_poly_60_local

3. Folium マップ作成と GeoTIFF エクスポート

# 施設一覧 CSV から各施設についてループ
for idx, facility in facilities_df.iterrows():
    origin_lat = float(facility['lat'])
    origin_lon = float(facility['lon'])
    # 出発時刻候補（15分刻みで計算）
    departure_times = [base_departure_input, base_departure_input + 15, ...]
    
    for dep in departure_times:
        iso_poly_30_dep, iso_poly_45_dep, iso_poly_60_dep = compute_isochrones_for_dep(dep)
        
        # Folium マップ作成（各圏を異なる色で表示）
        import folium
        m_dep = folium.Map(location=[origin_lat, origin_lon], zoom_start=13)
        # GeoJSON レイヤー追加
        folium.GeoJson(data=gpd.GeoSeries(iso_poly_60_dep).__geo_interface__, ...).add_to(m_dep)
        # ...（30分、45分圏も同様に追加）
        m_dep.save(output_path_dep)
        print("→ Folium マップのHTMLを保存しました:", output_path_dep)
        
        # GeoTIFF エクスポート
        export_geotiff(iso_poly_30_dep, f"{facility_folder}/iso_poly_30_{datetime_str_dep}_{lat_lon_str}.tif")
        # ...（45分、60分圏も同様に出力）

実装のポイントとメリット

• OTP がなくても実現可能
  OTP のセットアップが難しい環境下でも、GTFS データから直接公共交通のアイソクローナを作成できるため、手軽に実装可能です。

• 柔軟なカスタマイズ性
  施設一覧ファイルを用いることで、複数の出発地点に対して同時に解析が可能。さらに、出発時刻や歩行速度、乗換時間などのパラメータを自由に調整できるため、実際のシナリオに合わせた検証ができます。

• QGIS との連携
  GeoTIFF 出力により、QGIS で他の空間データ（人口分布、施設配置、地価情報など）と容易に重ね合わせられるため、空間分析の幅が広がります。

まとめ

今回の例では、OTP を使わずに Python スクリプトと GTFS データを組み合わせ、バス＋徒歩の到達圏を算出しました。出力した GeoTIFF やシェープファイルを QGIS に読み込めば、街中の施設データや地形データ、人口統計データなどと容易に重ね合わせられ、より広範な空間分析に活用できます。

この手法は、公共交通の改善や都市計画、災害時の避難経路の検討など、多岐にわたる分野で応用可能です。また、今回の実装はシンプルなアルゴリズムを採用しているので、処理速度や精度、リアルタイムデータの取り込みなど、今後の拡張にも柔軟に対応できます。

ぜひ皆さんも、独自のデータセットやシナリオと組み合わせて、QGIS 上で多角的な分析に挑戦してみてください。何かご意見や改善点があれば、コメント欄でお待ちしています。皆さんの公共交通・都市計画分野での活用事例が広がることを楽しみにしています！

---

Appendix: 全コードリスト

以下のコードは、本記事中で解説した内容をすべてまとめたものです。実際の運用環境に合わせてファイルパスやパラメータを適宜変更してください。

##############################
# USER INPUT PARAMETERS
##############################

# 以下の値を必要に応じて変更してください

# ① 入力日時（"YYYY/MM/DD HH:MM" 形式）
input_datetime_str = "2024/12/20 9:00"

# ② グローバルパラメータ
walking_speed = 80       # 歩行速度 [m/分]（例：約4.8km/h）
transfer_time = 2        # 乗換時間 [分]

# ※施設の出発地点は、後で施設一覧ファイルからインポートします

##############################
# 施設一覧ファイルのパス（施設情報：id,施設名,lat,lon）
##############################
facilities_csv_path = "/content/drive/MyDrive/gtfs関係/地図作成施設一覧（愛大病院）.csv"

##############################
# ① GTFSデータの解凍・読み込み
##############################

# 必要なライブラリのインストール（Colabの場合）
!pip install folium geopandas shapely rasterio

import os, zipfile, math, heapq, datetime
import pandas as pd
import numpy as np
import geopandas as gpd
from shapely.geometry import Point
import folium
import rasterio
from rasterio.features import rasterize
from rasterio.transform import from_origin

############################################################
# GTFSデータの解凍・読み込み
############################################################

# GTFSデータのパス（適宜変更してください）
gtfs_zip_path = "/content/drive/MyDrive/gtfs関係/iyotetsu_bus_20241205134041.zip"
extract_path = "/content/drive/MyDrive/gtfs関係/gtfs_extracted"

if not os.path.exists(extract_path):
    os.makedirs(extract_path)

with zipfile.ZipFile(gtfs_zip_path, 'r') as zip_ref:
    zip_ref.extractall(extract_path)

stops = pd.read_csv(os.path.join(extract_path, 'stops.txt'))
stop_times = pd.read_csv(os.path.join(extract_path, 'stop_times.txt'))
trips = pd.read_csv(os.path.join(extract_path, 'trips.txt'))
calendar = pd.read_csv(os.path.join(extract_path, 'calendar.txt'))
# ※必要に応じて calendar_dates の補正処理を追加してください

############################################################
# 入力日時と該当曜日の決定
############################################################

input_dt = datetime.datetime.strptime(input_datetime_str, "%Y/%m/%d %H:%M")
target_date_str = input_dt.strftime("%Y%m%d")
base_departure_input = input_dt.hour * 60 + input_dt.minute  # 分単位

weekday_map = {0:"monday", 1:"tuesday", 2:"wednesday", 3:"thursday", 4:"friday", 5:"saturday", 6:"sunday"}
weekday_str = weekday_map[input_dt.weekday()]

cal_active = calendar[
    (calendar['start_date'] <= int(target_date_str)) &
    (calendar['end_date'] >= int(target_date_str)) &
    (calendar[weekday_str] == 1)
]
active_service_ids = cal_active['service_id'].unique()

trips_active = trips[trips['service_id'].isin(active_service_ids)]
active_trip_ids = trips_active['trip_id'].unique()
stop_times_active = stop_times[stop_times['trip_id'].isin(active_trip_ids)].copy()

def time_to_minutes(t):
    try:
        h, m, s = t.split(':')
        return int(h)*60 + int(m) + int(s)/60.0
    except Exception:
        return np.nan

stop_times_active['arrival_mins'] = stop_times_active['arrival_time'].apply(time_to_minutes)
stop_times_active['departure_mins'] = stop_times_active['departure_time'].apply(time_to_minutes)
stop_times_active.sort_values(by=['trip_id', 'stop_sequence'], inplace=True)

# 各tripごとの停留所イベントリスト
trip_dict = {}
for trip_id, group in stop_times_active.groupby('trip_id'):
    events = group.sort_values('stop_sequence')[['stop_id','arrival_mins','departure_mins','stop_sequence']].to_dict('records')
    trip_dict[trip_id] = events

# 各停留所からの乗車可能なイベントリスト作成
stop_events = {}
for trip_id, events in trip_dict.items():
    for idx, event in enumerate(events):
        stop_id = event['stop_id']
        if stop_id not in stop_events:
            stop_events[stop_id] = []
        stop_events[stop_id].append({
            'trip_id': trip_id,
            'departure_mins': event['departure_mins'],
            'index': idx,
            'arrival_mins': event['arrival_mins']
        })
for stop_id in stop_events:
    stop_events[stop_id] = sorted(stop_events[stop_id], key=lambda x: x['departure_mins'])

# 停留所の座標情報
stops_dict = stops.set_index('stop_id')[['stop_lat','stop_lon']].to_dict('index')

############################################################
# Haversine の公式を用いた距離計算関数
############################################################

def haversine(lon1, lat1, lon2, lat2):
    R = 6371000  # 地球の半径 [m]
    phi1 = math.radians(lat1)
    phi2 = math.radians(lat2)
    dphi = math.radians(lat2 - lat1)
    dlambda = math.radians(lon2 - lon1)
    a = math.sin(dphi/2)**2 + math.cos(phi1)*math.cos(phi2)*math.sin(dlambda/2)**2
    return 2 * R * math.asin(math.sqrt(a))

############################################################
# 施設一覧ファイルの読み込み
############################################################

facilities_df = pd.read_csv(facilities_csv_path, encoding='utf-8-sig')
# 施設一覧のカラムは: id, 施設名, lat, lon とする

############################################################
# GeoTIFF エクスポート用関数
############################################################

def export_geotiff(geom, filename, pixel_size=0.0001, crs='EPSG:4326'):
    """GeoTIFF エクスポート用関数"""
    if geom is None:  # geomがNoneの場合の処理
        print(f"Error: Geometry is None. Skipping GeoTIFF export for {filename}")
        return  # エクスポートをスキップ

    minx, miny, maxx, maxy = geom.bounds
    width = int((maxx - minx) / pixel_size)
    height = int((maxy - miny) / pixel_size)
    transform = from_origin(minx, maxy, pixel_size, pixel_size)
    raster = rasterize(
        [(geom, 1)],
        out_shape=(height, width),
        transform=transform,
        fill=0,
        dtype='uint8'
    )
    with rasterio.open(
        filename,
        'w',
        driver='GTiff',
        height=height,
        width=width,
        count=1,
        dtype=raster.dtype,
        crs=crs,
        transform=transform,
    ) as dst:
        dst.write(raster, 1)
    print("GeoTIFF saved:", filename)

############################################################
# 指定した出発時刻での到達圏（Isochrone）を計算する関数（Haversine 版）
############################################################

def compute_isochrones_for_dep(dep):
    """
    dep: 出発時刻（分単位、絶対時刻）
    返り値: iso_poly_30, iso_poly_45, iso_poly_60（EPSG:4326）
    ※徒歩移動はHaversine公式による直線距離を用います。
    """

    earliest_local = {stop_id: float('inf') for stop_id in stops['stop_id'].unique()}
    pq = []

    # 出発地点から各停留所までの徒歩時間（直線距離）
    for stop_id, coords in stops_dict.items():
        dist = haversine(origin_lon, origin_lat, coords['stop_lon'], coords['stop_lat'])
        walk_time = dist / walking_speed
        if walk_time <= 60 and walk_time < earliest_local[stop_id]:
            earliest_local[stop_id] = walk_time
            heapq.heappush(pq, (walk_time, stop_id))

    while pq:
        current_time, current_stop = heapq.heappop(pq)
        if current_time > 60:
            continue

        if current_stop in stop_events:
            for event in stop_events[current_stop]:
                if event['departure_mins'] >= dep + current_time + transfer_time:
                    trip_id = event['trip_id']
                    idx = event['index']
                    events = trip_dict[trip_id]
                    for next_event in events[idx+1:]:
                        arrival_offset = next_event['arrival_mins'] - dep
                        if arrival_offset < earliest_local[next_event['stop_id']]:
                            earliest_local[next_event['stop_id']] = arrival_offset
                            heapq.heappush(pq, (arrival_offset, next_event['stop_id']))

        for other_stop, coords in stops_dict.items():
            if other_stop == current_stop:
                continue
            current_coords = stops_dict[current_stop]
            d = haversine(
                current_coords['stop_lon'], current_coords['stop_lat'],
                coords['stop_lon'], coords['stop_lat']
            )
            walk_time = d / walking_speed
            new_time = current_time + walk_time
            if new_time < earliest_local[other_stop] and new_time <= 60:
                earliest_local[other_stop] = new_time
                heapq.heappush(pq, (new_time, other_stop))

    iso_stops_30 = [stop_id for stop_id, t in earliest_local.items() if t <= 30]
    iso_stops_45 = [stop_id for stop_id, t in earliest_local.items() if t <= 45]
    iso_stops_60 = [stop_id for stop_id, t in earliest_local.items() if t <= 60]

    gdf_stops_local = gpd.GeoDataFrame(
        stops,
        geometry=stops.apply(lambda row: Point(row['stop_lon'], row['stop_lat']), axis=1),
        crs="EPSG:4326"
    )
    gdf_30_local = gdf_stops_local[gdf_stops_local['stop_id'].isin(iso_stops_30)]
    gdf_45_local = gdf_stops_local[gdf_stops_local['stop_id'].isin(iso_stops_45)]
    gdf_60_local = gdf_stops_local[gdf_stops_local['stop_id'].isin(iso_stops_60)]

    gdf_30_3857_local = gdf_30_local.to_crs(epsg=3857)
    gdf_45_3857_local = gdf_45_local.to_crs(epsg=3857)
    gdf_60_3857_local = gdf_60_local.to_crs(epsg=3857)

    def create_isochrone_polygon_local(gdf, time_limit, earliest_local):
        buffers = []
        for idx, row in gdf.iterrows():
            stop_id = row['stop_id']
            t_arrival = earliest_local[stop_id]
            remaining = max(time_limit - t_arrival, 0)
            buffer_dist = remaining * walking_speed  # [m]
            buffers.append(row.geometry.buffer(buffer_dist))
        if buffers:
            return gpd.GeoSeries(buffers).unary_union
        else:
            return None

    iso_poly_30_3857_local = create_isochrone_polygon_local(gdf_30_3857_local, 30, earliest_local)
    iso_poly_45_3857_local = create_isochrone_polygon_local(gdf_45_3857_local, 45, earliest_local)
    iso_poly_60_3857_local = create_isochrone_polygon_local(gdf_60_3857_local, 60, earliest_local)

    def to_epsg4326_geom_local(geom_3857):
        gdf_temp = gpd.GeoDataFrame(geometry=[geom_3857], crs="EPSG:3857")
        return gdf_temp.to_crs(epsg=4326).iloc[0].geometry

    iso_poly_30_local = to_epsg4326_geom_local(iso_poly_30_3857_local)
    iso_poly_45_local = to_epsg4326_geom_local(iso_poly_45_3857_local)
    iso_poly_60_local = to_epsg4326_geom_local(iso_poly_60_3857_local)

    return iso_poly_30_local, iso_poly_45_local, iso_poly_60_local

############################################################
# ⑦ 各施設についてマップ作成＆エクスポート（複数施設対応）
############################################################

# 施設一覧ファイルから施設情報を読み込み（カラム: id, 施設名, lat, lon）
facilities_df = pd.read_csv(facilities_csv_path)

# 施設一覧の各行についてループ
for idx, facility in facilities_df.iterrows():
    facility_name = str(facility['施設名']).strip()  # 施設名
    # ※必要に応じてファイル名用にサニタイズ（記号除去など）してください
    facility_folder = f"/content/drive/MyDrive/gtfs関係/{facility_name}"
    os.makedirs(facility_folder, exist_ok=True)

    # 施設の出発地点として設定（グローバル変数として上書き）
    origin_lat = float(facility['lat'])
    origin_lon = float(facility['lon'])

    # 施設ごとの出力用緯度経度文字列
    lat_lon_str = f"{origin_lat:.5f}_{origin_lon:.5f}"

    print(f"【施設: {facility_name} (lat,lon: {lat_lon_str})】")

    # 出発時刻候補（分単位）：入力時刻そのもの（＋0）に加えて、15分刻みで1時間後まで
    departure_times = [
        base_departure_input,
        base_departure_input + 15,
        base_departure_input + 30,
        base_departure_input + 45,
        base_departure_input + 60
    ]

    # 施設毎に、各出発時刻のマップ・GeoTIFFを作成
    for dep in departure_times:
        # ① 各出発時刻で到達圏ポリゴンを計算
        iso_poly_30_dep, iso_poly_45_dep, iso_poly_60_dep = compute_isochrones_for_dep(dep)

        # ② Foliumマップの作成
        m_dep = folium.Map(location=[origin_lat, origin_lon], zoom_start=13)

        folium.GeoJson(
            data=gpd.GeoSeries(iso_poly_60_dep).__geo_interface__,
            name='60分圏',
            style_function=lambda feature: {
                'fillColor': 'lightblue',
                'color': 'blue',
                'weight': 2,
                'fillOpacity': 0.5
            }
        ).add_to(m_dep)

        folium.GeoJson(
            data=gpd.GeoSeries(iso_poly_45_dep).__geo_interface__,
            name='45分圏',
            style_function=lambda feature: {
                'fillColor': 'lightgreen',
                'color': 'green',
                'weight': 2,
                'fillOpacity': 0.5
            }
        ).add_to(m_dep)

        folium.GeoJson(
            data=gpd.GeoSeries(iso_poly_30_dep).__geo_interface__,
            name='30分圏',
            style_function=lambda feature: {
                'fillColor': 'salmon',
                'color': 'red',
                'weight': 2,
                'fillOpacity': 0.5
            }
        ).add_to(m_dep)

        # 補助：各停留所をCircleMarkerとして追加
        gdf_stops = gpd.GeoDataFrame(
            stops,
            geometry=stops.apply(lambda row: Point(row['stop_lon'], row['stop_lat']), axis=1),
            crs="EPSG:4326"
        )
        for i, row in gdf_stops.iterrows():
            folium.CircleMarker(
                location=[row.geometry.y, row.geometry.x],
                radius=2,
                color='black',
                fill=True,
                fill_opacity=1
            ).add_to(m_dep)

        folium.LayerControl().add_to(m_dep)

        # ③ 出力ファイル名作成：日時（target_date_str + 出発時刻）と緯度経度を含む
        hour_dep = dep // 60
        minute_dep = dep % 60
        datetime_str_dep = f"{target_date_str}_{hour_dep:02d}{minute_dep:02d}"

        # HTML出力
        output_path_dep = f"{facility_folder}/interactive_map_{datetime_str_dep}_{lat_lon_str}.html"
        m_dep.save(output_path_dep)
        print("  → FoliumマップのHTMLを保存しました:", output_path_dep)

        # ④ GeoTIFF出力（各到達圏ごと）
        export_geotiff(iso_poly_30_dep, f"{facility_folder}/iso_poly_30_{datetime_str_dep}_{lat_lon_str}.tif")
        export_geotiff(iso_poly_45_dep, f"{facility_folder}/iso_poly_45_{datetime_str_dep}_{lat_lon_str}.tif")
        export_geotiff(iso_poly_60_dep, f"{facility_folder}/iso_poly_60_{datetime_str_dep}_{lat_lon_str}.tif")

``