OSSで札幌駅に2時間以内に辿り着ける範囲を可視化してみた〜QGISプラグインも自作！〜

先日、Xにこんな投稿をしたところ大変ご好評をいただくことができました！
今日はこの地図の作り方を簡単にご説明しようかと思います。


<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">公共交通機関で札幌に2時間以内に着ける範囲を、10分ごとに可視化してみました <a href="https://t.co/Sm4tmOhbhR">pic.twitter.com/Sm4tmOhbhR</a></p>&mdash; Shota Yamamoto (@geogra_geogra) <a href="https://twitter.com/geogra_geogra/status/1849832074951393380?ref_src=twsrc%5Etfw">October 25, 2024</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>


## そもそもこれは何？
この図のような地図は、到達圏(Isochrone アイソクロン)と呼ばれる手法の表現方法です。

![1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/d6b75c4f-36e1-d541-1dcc-df7e5fa87361.png)

到達圏とは、"指定の時刻"に"指定の地点"を、"出発または到着"した場合にどこまで辿り着けるのかを図示したものです。

今回の場合は、札幌駅の南口に平日の朝9時〜夜24時までに2時間以内にたどり着くことができるエリアを対象としていて、それぞれ所要時間が何分かかるのかという設定で地図を作成しています。
従来の研究では、指定の時刻の到達圏を時刻表を手計算したり、鉄道やバスの平均速度から移動可能距離を計算することが多かったです。
しかし、これからご紹介する方法を用いることで、既存のオープンデータを用いることで簡単に到達圏を表示することができます！

:::note info
このプラグインは、2024年実施の「第33回地理情報システム学会学術研究発表大会」にて発表を行なった「[GTFSデータを使用した到達圏を解析・表示するQGISプラグインの開発](https://pub.confit.atlas.jp/ja/event/gisa2024/presentation/E5-02 )」 に合わせて作成したものです。 
発表や論文には含まれない、開発内容に特化したご説明となります。
:::

## 到達圏の作り方
### OpenTripPlannerとConveyal R5
到達圏の作成には、[OpenTripPlanner(OTP)]()というOSSを用いました。OTPは、現在Ver.2.6ですが、到達圏解析機能は、Ver.1.5までの機能となります。本来は、後継の[Conveyal R5]()を用いるべきですが、サーバーの構築の容易さや再現性の観点から、今回はOTPを用いて解析を行いました。
なお、Conveyal R5は現在も改良が進んでおり、より高精度な解析を行うことが可能です。
Conveyal社の開発するConveyalというソフトでの利用を前提としていますが、エンジン部分のR5がOSSとして公開されています。

### OpenTripPlannerの使い方
ここでは、OTPを動作させるために３つのデータを使用します
- OTP
- GTFSデータ
- OpenStreetMapの道路データ

データの加工方法や、サーバーの立て方については、私含めていくつかの記事がありますので、こちらをご参考にしてください

https://qiita.com/geogra/items/4e7feea23408a324540c

https://qiita.com/kumatira/items/5235554b12f5f2725301



## OTPを使って到達圏を作成する
実際にサーバーを動作させて、到達圏を作成してみましょう。
サーバーはREST APIの形で呼び出すので、ブラウザ上に以下のようなパラメータを入力することでジオメトリを出力できます。

例：12/1の15時に札幌駅南口へ到着するためにかかる時間を5分刻みで60分まで算出する
```
http://localhost:8080/otp/routers/default/isochrone?fromPlace=43.066887086991805,141.35168508191856&toPlace=43.06175069808622, 141.3417904693123&arriveBy=true&date=2024-12-01&time=15:00:00&mode=WALK,TRANSIT&cutoffSec=300&cutoffSec=600&cutoffSec=900&cutoffSec=1200&cutoffSec=1500&cutoffSec=1800&cutoffSec=2100&cutoffSec=2400&cutoffSec=2700&cutoffSec=3000&cutoffSec=3300&cutoffSec=3600
```
これをブラウザに打ち込むとこのようなJSONファイルが出力されました（データの大きさによってはEsri Shapefileとなることもあります。）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/c8252584-9cdd-3e06-bff4-ac9227c76c73.png)

QGIS上で表示するとこのようになりました
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/a6710d18-3828-e330-0264-2ae641662bcc.png)

レイヤプロパティからシンボロジを「カテゴリ値による定義」でtime属性を値に設定し、カラーランプSpectralにするとこのようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/70606fb7-c512-4e68-4844-e195510fd315.png)

## 到達圏を複数時間で連続して算出したい！
一つずつリクエストを作成して、QGISに表示する方法ではGTFSのスケールメリットを活かすことができません。そこで、URLを連続して準備してリクエストを行うことで、複数時間の到達圏を自動で取得することができるようになりました。
気をつけた点としては、到着・出発の両方に対応することや、計算すべき時間間隔を計算で導出すること、指示をせずに呼び出すとデータ量によってShapefileとGeoJSONの形式が統一されないことから、これを事前に設定することなどです。

```py
def request_isochrone(self, lat, lon):
    start_time = self.ui.startTime.dateTime()
    finish_time = self.ui.finishTime.dateTime()

    # maxtime: 処理時間の上限（分）
    cutoff_query = self.generate_cutoff_secs(
        self.ui.maxTime.value(),
        self.ui.outputPolygonInterval.value(),
    )

    # デフォルトの始点と終点を設定
    # 到着地点を指定しないとエラーが出るため、一応デフォルトを用意してある
    lat_dep, lon_dep = ("43.0815", "141.3074")
    lat_ari, lon_ari = ("43.0815", "141.3074")
    self.arrive_boolean = False

    # buttonGroupの選択に基づいて到着なのか出発なのかを判定
    if self.ui.arriveby.checkedButton() == self.ui.button_dep:
        lat_dep, lon_dep = (lat, lon)
        self.arrive_boolean = False
    elif self.ui.arriveby.checkedButton() == self.ui.button_ari:
        lat_ari, lon_ari = (lat, lon)
        self.arrive_boolean = True

    arrive_boolean_str = "true" if self.arrive_boolean else "false"

    # サーバーURLを設定
    # 例：server_url = "http://localhost:8080/"
    if server_url.endswith("/"):
        server_url = server_url[:-1]  # 最後の「/」を取り除く

    # start_timeとfinish_timeの間を指定された間隔（分）ごとに処理
    current_time = start_time

    # リクエストループ
    while current_time <= finish_time and not self.error_occurred:
        current_date_str = current_time.toString("yyyy-MM-dd")
        current_time_str = current_time.toString("HH:mm:ss")  # 秒を含めた時間形式

        # フルURLの生成
        full_url = (
            f"{server_url}/otp/routers/default/isochrone"
            f"?fromPlace={lat_dep},{lon_dep}&toPlace={lat_ari},{lon_ari}"
            f"&date={current_date_str}&time={current_time_str}&mode=WALK,TRANSIT"
            f"&arriveBy={arrive_boolean_str}&{cutoff_query}"
        )

        # リクエスト送信とレスポンス処理
        request = requests.get(full_url, headers=headers)
        # データ量が多くなると勝手にShapefileが選ばれるので、ヘッダーでgeojsonを指定
        request.setRawHeader(b"Accept", b"application/json")

        # リクエストを送信して次の時刻に進める
        self.send_request_and_wait(request, current_date_str, current_time_str)
        current_progress += progress_increment
        self.ui.progressBar.setValue(int(current_progress))

        # 指定された分の間隔で current_time を進める
        # 例：time_interval_minutes=15 の場合、15分ごとに処理を行う
        current_time = current_time.addSecs(time_interval_minutes * 60)
```

## QGISプラグインの作成
上記のコードを用いることで、連続してGeoJSON形式のジオメトリを出力することができましたが、このままQGISにインプットすると、全てのデータにシンボロジを適用する必要があります。
そこで、シンボロジの適用に加えて、各種解析までを済ませるプラグインを作成しました。

プラグインは[こちら](https://plugins.qgis.org/plugins/qgis_otp_multi_isochrone_plugin/)

https://plugins.qgis.org/plugins/qgis_otp_multi_isochrone_plugin/

https://github.com/geogra-geogra/qgis-otp-multi-isochrone-plugin

### 到達圏にシンボロジを適用する
このように、出力が表示された時点でシンボロジが適用されます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/8c5856b2-5fa5-1f88-99c8-afa07888d101.png)
*地理院地図（白地図）を背景に利用*

### グリッドを作成する
グリッドを作成する場合、このように作成したすべてのジオメトリが包含されるグリッドを作成します。
グリッドの幅は指定することができる他、標準地域メッシュなどの既存のメッシュデータを利用することも可能です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/3ac63686-5acc-b341-9a5f-481ba8272fad.png)
*地理院地図（白地図）を背景に利用*

### ラスタを作る
作成したグリッドを基準にラスタデータを作成します。
データ自体は保存されますが、QGIS上には表示されません。

### 統計値を計算する
統計値として、すべての時刻データのラスタを重ね合わせて、平均値・中央値・最大値・最小値、そして最大値と最小値の差分をラスタで表示します。

**最小値**
1日で札幌駅に最短で到着する所要時間。「地域のポテンシャル」をより正確に表示することができます。
![min.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/e0f272c0-c352-fb29-6921-305c4278e371.png)

**最大値**
最大の所要時間(2時間以内)
この場合、朝イチには電車がないため地下鉄のエリアを中心に境界があります。
![max.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/59af4432-304d-8341-dfc4-4e0caad949af.png)

**中央値**
1日での所要時間の中央値なので、これが概ねその地域の通常の所要時間と言えるでしょう。
ただし、今回の計測方法では2時間以上かかる場合の計測を一切していません。最短の所要時間がちょうど2時間の場合、統計値がすべて2時間になることに注意が必要です。
![median.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2747613/6c8958f1-ff0c-54ea-6994-179ccf668b5e.png)

このようにして、今までにない方法で到達圏を計測することが可能になります。

## まとめ
このように、自動化を行うことで、今まで手計算で求められたり、平均速度でまとめられることが多かった到達圏を細かく解析することができるようになりました。
また別の記事で、詳しい活用方法などをご紹介していこうと思っています！


弊社新オフィス移転について分析したきj

https://qiita.com/geogra/items/dbec185059e2882008df