SpatiaLiteに新しいネットワーク解析が搭載されたので触ってみた

こんにちは。今年も、また性懲りもなくSpatiaLiteで参戦します。

---

はじめに

今年の8/3に、SpatiaLite5.0.0-betaのリリースのアナウンスがありました。
現在の安定版4.3.0の公開は2015/07/01、次の4.4.0がRC版として利用できるようになったのが2015/12/14、しかしその後一向にRCがStableになる気配もなくついにはStableにならないことが決定し、3年近くバージョンアップがされないままとなっていました。
最近では、「もはやこのプロジェクトは死んだのでは…？」的な不穏なコメントが出始めていましたが、満を持してメジャーバージョン5のリリースサイクルのアナウンスとともに、beta版が利用できるようになりました。

まだbetaではありますが、公開後継続して使っている限りでは、安定して利用できています。

メジャーバージョンアップということで、いくつか大きな変更や機能の追加がされています。

• LWGEOMライブラリを使用していた機能を全てRTTOPOに移行した
• PostgreSQLのテーブルに対してトランザクション（BEGIN/COMMIT/ROLLBACK）できるようになった
• SpatiaLite-GUIにおいて、PostgreSQL/PostGISを利用できるVirtualPGがエクステンションから標準装備になった
• ストアド処理、ストアド変数が利用できるようになった
• トポロジーに基づくジオメトリモデル（SpatiaLite Topology）の追加（Topology-GeometryとTopology-Network）
• いろいろ関数も追加された
• 詳しくはこちら

公開されたライブラリ・ツール群の構成は以下の4つで、いろいろあったのが集約された感じになっています。一番下のSpatiaLite-GUI（SpatiaLiteのGUIツール）には、上3つのライブラリ群が全てアサインされています。

• libspatialite-5.0.0
• librasterlite2-1.1.0
• virtualpg-2.0.0
• spatialite_gui-2.1.0

さて、今回のバージョンから、新しいネットワーク解析（Routing）の機能が装備されたので、ここではこれを触ってみようと思います。

---

SpatiaLiteによるネットワーク解析（VirtualRouting）

SpatiaLiteでRoutingをすることは、以前のバージョンでも可能でした。（VirtualNetwork）
Routingに使用するためのデータを生成するには、SpatiaLite-GUIに搭載されていたメニューを使用するか、spatialite_networkというCLIツールを利用して、Routing用のデータモデルを生成する必要がありました。
また、ネットワークデータとして利用できるデータ（終点と始点にユニークなIDを付与した点データと線データを用意する）に道路データを整形する手続きのハードルが高いことも難点でした。
また、以前のデータモデルでは、最短経路問題しか扱うことができなかったため、用途が限定されていました。

SpatiaLite5では、Routingエンジンが以前のVirtualNetworkから新しくVirtualRoutingというデータモデルへ移行し、データモデル生成からRoutingの実施までを全てSQL関数で処理できるようになったため、利用性が高まったと思います。

参考：SpatiaLite: VirtualRouting
（この記事は、上記Webページのチュートリアルを参考に作成しました。）

---

VirtualRoutingを触ってみる

※ここでの実証は、Windows10-64bitの環境下で行っています。

#0 ネットワークデータの準備

OSMのダウンロード

こちらのサイトから、北海道のOSMデータをshpでダウンロードしてきました。（北海道全域で約180MB）

SpatiaLite5.0.0-betaの入手

Windowsのバイナリ配布元（64bit、32bit）から、SpatiaLite-GUI（spatialite_gui-NG-win-xxxx.7z）か、昨年紹介した mod_spatialite（mod_spatialite-NG-win-xxxx.7z）をSQLiteにロードするなどの方法で利用します。（配布ファイルは7z圧縮されているので、何らかの方法で展開してください。）

以下、ここでは、SpatiaLite-GUIでの操作を載せます。
Menu > About で、バージョンが確認できます。

クエリで確認することもできます。

SELECT spatialite_version();

VirtualRoutingが利用できるかどうかは、以下で確認できます。

SELECT HasRouting(); -- 1ならTrue

空のDBの作成

Menu > Creating a New (Empty) SQLite DBから、空のSpatiaLite-DBを作成します。（data.sqliteとしました。）
mod_spatialiteやSpatiaLite-CLIを利用する場合は、SELECT InitMetaData();の実行によってメタデータの作成が必要です。（GUIの場合は不要。）

OSMデータのインポート

Menu > Advanced > Load Shapefile

• ジオテーブル名：gt_osm_road
• ジオメトリカラム名：geom_ln_4326
• SRID：4326（WGS84）

30秒ほど待つと、447858レコードがインサートされました。

OSMデータの修正

利用しない道路の削除

自動車が通行できる道路に限定するとともに、管理道路などを除外します。

DELETE FROM "gt_osm_road" WHERE "fclass" IN ('bridleway', 'cycleway', 'footway', 'path', 'service');
-- 45s

一方通行の方向修正とフラグカラムの追加

OSMの道路データにはonewayカラムがあり、B（双方向）、F（順方向）、T（逆方向）のフラグが入っていました。一方通行は、ジオメトリの方向にのみ通行可とするのが扱いやすいので、逆方向になっているものを反転します。

UPDATE "gt_osm_road"
SET "geom_ln_4326" = ST_Reverse("geom_ln_4326") WHERE "oneway" IS 'T';
-- 1s

後ほどRoutingデータを作成する際に、通行可能フラグがBoolean（0または1）で必要になるため、このフラグカラムflg_invを追加します。（1なら逆方向も通行可能。）
また、現状では順方向の通行可能フラグについても必ず必要であるため、全て1のフラグカラムflg_fwdも追加します。

ALTER TABLE "gt_osm_road" ADD COLUMN "flg_inv" INTEGER;
UPDATE "gt_osm_road" SET "flg_inv" = CASE WHEN "oneway" IS 'B' THEN 1 ELSE 0 END;
ALTER TABLE "gt_osm_road" ADD COLUMN "flg_fwd" INTEGER;
UPDATE "gt_osm_road" SET "flg_fwd" = 1;
-- 55s

（ジオメトリの方向を修正せず、各方向フラグで管理してもいいと思います。）

ノード間のセグメントに分解

OSMの道路データをそのままRoutingに利用しようとしたところ、交差点でラインが切れておらず、長いラインストリングになっていました。このようなデータでは、端点の交差点でしか曲がることができないため、正しいRoutingができません。さいわい、全ての交差点にノードが配置されていたので、（少し強引ですが）ラインストリングを全てのノードで切断したテーブルを生成しました。（クエリの詳細説明は割愛します。）

• ジオテーブル名：gt_road
• ジオメトリカラム名：geom_ln_4326

-- ジオメトリタイプが変更になるので、ジオメトリの登録を解除
SELECT DiscardGeometryColumn('gt_osm_road', 'geom_ln_4326');
-- セグメントに分割（マルチラインストリングになる）
UPDATE "gt_osm_road" SET "geom_ln_4326" = ST_Multi(ST_DissolveSegments("geom_ln_4326"));
-- ジオメトリの分解にはジオメトリカラムの再登録が必要
SELECT RecoverGeometryColumn('gt_osm_road', 'geom_ln_4326', 4326, 'MULTILINESTRING');
-- 要素ジオメトリへの分解
SELECT ElementaryGeometries(
  'gt_osm_road', -- 元テーブル名
  'geom_ln_4326',  -- 元テーブルのジオメトリカラム名
  'gt_road', -- 生成テーブル名
  'id_seg', -- 生成テーブルにおけるプライマリキーカラム名
  'id_src' -- 元テーブルのプライマリキー値を格納するカラム名
);
-- 1min 50s

分割されて、約341万行のデータになりました。
元のデータは不要になったため削除し、空き領域を開放します。

SELECT DropGeoTable('gt_osm_road'); -- ジオテーブルの削除
VACUUM; -- 解放
-- 1min 20s

この時点でDBファイルのサイズは、475MBになりました。
ここまでの操作は、Routing用に用意されたベクタデータであれば既に要件を満たしていて、不要かもしれません。

Routingデータの作成

やっと本題です。新たなRouting用のデータモデルであるVirtualRoutingのデータを作成するには、関数CreateRouting()を利用します。基本的な記法は以下の通りです。

SELECT CreateRouting(
  '生成されるネットワークバイナリのテーブル名' -- 任意の文字列
  , '生成されるVirtualRoutingテーブル名' -- 任意の文字列
  , 'ネットワークのベースになるジオテーブル名' -- Spatial Viewも可
  , '始点ノードIDのカラム名'
  , '終点ノードIDのカラム名'
  , 'ラインストリングジオメトリのカラム名' -- NULLの場合は地理データを含まないネットワークとして扱われる（Logical Network）
  , 'コストのカラム名' -- NULLの場合はラインストリングの長さがコストとして扱われる*1
  , 'ネットワーク名のカラム名' -- （オプション）道路名を格納した列などを指定。NULLも可
  , 0 or 1 -- （オプション）探索に使用するアルゴリズムでDijkstra法に加えてA-star（A*）法も併用するかを指定（1がTrue、デフォルトはTrue）
  , 0 or 1 -- （オプション）0＝無向グラフ 1＝有向グラフ（デフォルトは1）
  , '正方向への通行可否を示すカラム名' -- （オプション）列は整数型で定義され、0または1のどちらかが格納されている必要がある*2
  , '逆方向への通行可否を示すカラム名' -- （オプション）同上
  , 0 or 1 -- （オプション）テーブルを上書きするか（デフォルトは0=False）
);

1. ラインストリングの長さは、ネットワークデータの空間参照で決まります。地理座標系の場合は楕円体の表面長さがメートル単位で、投影座標系の場合は投影座標上の長さが、空間参照システムで指定されている長さ単位で返されます。なお、ジオメトリを含まないネットワークとする場合は、コストを表す列は必ず必要です。
2. 無向グラフの場合は、正方向／逆方向への通行可否の指定は無効となるため、NULLのみ指定できます。

また、CreateRouting()では、ベースになるジオテーブルに「始点ノードID」と「終点ノードID」がカラムとして必要ですが、このようなカラムがない場合にこのIDを生成する関数がCreateRoutingNodes()です。自前でデータを準備するにはハードルが高かった部分だったため、この機能はとても効果的です。

SELECT CreateRoutingNodes(
  'データベースのエイリアス名' -- アタッチした外部DBのデータを利用する場合（NULLの場合はこのDB本体）
  , 'ベースになるネットワークのジオテーブル名'
  , 'ネットワークのジオメトリカラム名'
  , '始点IDを格納するカラム名'
  , '終点IDを格納するカラム名'
);

それでは、道路データにノードIDを付与してVirtualRoutingデータを生成してみましょう。

-- ノードIDの付与
SELECT CreateRoutingNodes(
  NULL -- NULLの場合はこのDB本体
  , 'gt_road' -- 道路テーブル
  , 'geom_ln_4326' -- ジオメトリカラム
  , 'node_from' -- 始点IDを格納するカラム名
  , 'node_to' -- 終点IDを格納するカラム名
);
-- 3min 13s

道路のジオテーブルには新たなカラムnode_fromとnode_toが追加されています。

-- VirtualRoutingデータの生成
SELECT CreateRouting(
  'routing_net_data'
  , 'routing_net'
  , 'gt_road'
  , 'node_from' -- from_col
  , 'node_to' -- to_col
  , 'geom_ln_4326'
  , NULL -- コストカラム（NULLの場合は長さ）
  , 'name' -- 道路名カラム
  , 1 -- A*アルゴリズムも利用可能とする
  , 1 -- 有向グラフ
  , 'flg_fwd' -- 順方向フラグ
  , 'flg_inv' -- 逆方向フラグ
  , 1 -- overwrite
);
-- 3min 43s

VirtualRoutingのデータセットである、routing_net_dataとrouting_netが生成されています。

ここまでで、データベースのサイズは987MBになりました。ずいぶん大きくなりました...
生成されたデータを覗いてみます。

VirtualRoutingテーブル

SELECT * FROM "routing_net";

routing_netには、Routingの設定が格納されています。また、探索の結果もこのテーブル（仮想テーブル）に出てきます。

Algorithm	Request	Options	Delimiter	...
Dijkstra	Shortest Path	Full	, [dec=44, hex=2c]	（他のカラムはNULL）

ネットワークバイナリデータテーブル

SELECT * FROM "routing_net_data";

こちらのテーブルには、Routingのためのバイナリデータが格納されています。通常、このテーブルに直接アクセスすることはないと思います。

id	NetworkData
0	BLOB
1	BLOB
2	BLOB
...	...

（補足）ネットワークテーブルの削除について

生成したVirtualRoutingデータを削除する際は、はじめにVirtualRoutingテーブル（上の例ではrouting_net）を削除し、次にBinaryDataテーブル（上の例ではrouting_net_data）を削除します。先にBinaryDataテーブルを削除してしまうと、VirtualRoutingテーブルを削除できなくなるので注意が必要です。

端点ノードのジオテーブルの作成

ここまでに作成したVirtualRoutingデータセットがあればRoutingはできるのですが、「じゃあこの端点のIDは何なのよ」というのがちょっとわかりにくいので、端点にIDをもたせたポイントデータを別に生成することにします。

-- 端点ノードのリストを作成
CREATE TABLE "gt_road_nodes" ("id" INTEGER PRIMARY KEY);
SELECT AddGeometryColumn('gt_road_nodes', 'geom_pt_4326', 4326, 'POINT');
INSERT INTO "gt_road_nodes"
  SELECT * FROM (
    SELECT
      "node_from" AS "id",
      ST_StartPoint("geom_ln_4326")
    FROM "gt_road"
    UNION
    SELECT
      "node_to",
      ST_EndPoint("geom_ln_4326")
    FROM "gt_road"
  )
  ORDER BY id
;
-- 空間インデックスも付与
SELECT CreateSpatialIndex('gt_road_nodes', 'geom_pt_4326');
-- 3min 21s

ネットワークのデータと、ノードのデータが利用できるようになりました。

以上で下準備は終わりです（長
ファイルサイズは1.3GBまで大きくなりました。

---

#1 最短経路（Node to Node）

最短経路の探索は、次のクエリで得られます。

SELECT * FROM "routing_net" -- VirtualRoutingテーブル名
WHERE "NodeFrom" == 100000 -- 始点ノードID（整数）
AND "NodeTo" == 200000 -- 終点ノードID（整数）
;

クエリ結果の一行目に、経路全体のコストと経路のジオメトリが示されています。（Route）
二行目以下に、個別の道路区間の情報が示されています。（Link）
この結果をジオテーブルに挿入し、QGISで確認してみます。

-- ジオテーブルの生成
CREATE TABLE "gt_path" (
  "PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT
  , "Algorithm" TEXT -- 使用したアルゴリズム
  , "Options" TEXT -- オプション
  , "NodeFrom" INTEGER -- 始点ノードID
  , "NodeTo" INTEGER -- 終点ノードID
  , "Cost" REAL -- コスト
);
SELECT AddGeometryColumn('gt_path', 'geom_lnm_4326', 4326, 'LINESTRINGM');
-- Routing
UPDATE "routing_net" SET "Options" = 'NO LINKS'; -- ルートだけが返る設定に変更
INSERT INTO "gt_path"
SELECT 1, "Algorithm", "Options", "NodeFrom", "NodeTo", "Cost", "Geometry"
FROM "routing_net"
WHERE "NodeFrom" == 100000
AND "NodeTo" == 200000
;
-- 2s

なお、Routingクエリを投げる前に、Optionsカラムを"NO LINKS"としておくことで、ジオメトリを持つRouteだけを取得しています。
なお、VirtualRoutingで得られるRouteのジオメトリはLINESTRING Mというジオメトリ形式で得られます。（ネットワークデータが3次元の場合はLINESTRING ZM。）M次元を持つラインストリングは線形参照（Linear Referencing）をサポートするジオメトリタイプで、これをサポートする関数によって、Routeの一部（出発後xx[m]からyy[m]までとか）を抽出したり、特定の位置を抽出したりすることができます。（詳細は本家チュートリアルページを参照。）

函館市の海岸から出発して有珠山の南側まで、167kmの道程でした。

出発地と目的地を設定するにあたっては、近傍のノードIDをGISで確認して、そのノード間の経路を探索するということになります。

始点と終点をグラフィカルに設定する

せっかくなので（？）始点と終点をポイントで設定し、その近傍にある端点間のルートを得られるようにしてみましょう。
新しく、ジオテーブルgt_from_toを作成し、ポイントのジオメトリカラムを登録します。

CREATE TABLE "gt_from_to" ("PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT , "PointName" TEXT);
SELECT AddGeometryColumn('gt_from_to', 'geom_pt_4326', 4326, 'POINT');

レイヤをQGISに読み込み、StartとEndのノードを追加します。

出発点を札幌駅南口の広場、目的地を札幌ドームにしてみました。

このレイヤの編集を確定した上で（必ずしも編集モードをオフにする必要はありません）、以下のSQLでデータを更新します。

INSERT OR REPLACE INTO "gt_path"
SELECT 1, "Algorithm", "Options", "NodeFrom", "NodeTo", "Cost", "Geometry"
FROM "routing_net"
WHERE "NodeFrom" == (
  SELECT "id"
  FROM (
    -- 始点に最も近いネットワークノードを検索
    SELECT
      "gt_road_nodes"."id"
      , Min(ST_Distance("gt_from_to"."geom_pt_4326", "gt_road_nodes"."geom_pt_4326"))
    FROM "gt_from_to", "gt_road_nodes"
    WHERE "gt_from_to"."PointName" IS 'Start'
  )
)
AND "NodeTo" == (
  SELECT "id"
  FROM (
    -- 終点に最も近いネットワークノードを検索
    SELECT
      "gt_road_nodes"."id"
      , Min(ST_Distance("gt_from_to"."geom_pt_4326", "gt_road_nodes"."geom_pt_4326"))
    FROM "gt_from_to", "gt_road_nodes"
    WHERE "gt_from_to"."PointName" IS 'End'
  )
)
;
-- 20s

少し待つと、クエリが完了します。
QGISの表示をリフレッシュすると、新たな経路が正しく表示されました。

出発地を拡大して見てみると、近傍点として15081が選択されていました。

トリガーや空間ビューによる即時反映について

ここで紹介している例では、始点と終点の位置を確定した後、上記のINSERTを手動で実行する必要があります。SQLiteの機能であるトリガーや、SpatiaLiteの機能である空間ビューを利用することで、始点と終点を設定したら即、経路が反映されるようにできないかと思い試してみましたが、現在のQGISにマウントされているSpatiaLiteが4.3.0であり、VirtualRoutingに対応していないことから、現時点では実現できませんでした。
どちらの方法も、SpatiaLite5がマウントされている環境では可能になると考えられます。

KNN（k近傍法）による近接点の抽出

少し脱線しますが、前記の経路探索は20秒と、やや時間がかかりました。これに対して、IDを直接入力して経路を得る場合は2秒以下で終了していました。これは、出発点に指定した地点から最近傍となるノードを選択するのに時間がかかっていることが原因です。（いや、数百万もの点から20秒で抽出できるのもすごいですけど^^;）
上記のクエリでは、Min()関数を使用して二点間の距離が最小となる点を抽出していますが、例えば点からバッファを生成して、これに掛かる点を先に抽出するなどの方法を取っても、近傍点の抽出コストはほとんど改善しません。

この近傍点の抽出はよくある問題であり、SpatiaLiteではここにかかるコストを改善するためにVirtualKNNという機能が実装されています。この機能を利用することで、コストは顕著に改善します。（SpatiaLite4.4.0以降の機能のため、現行のQGISに搭載されている4.3.0aでは利用できません。）
SpatiaLite: KNN

※なお、VirtualKNNはジオメトリにSpatialIndexが設定されている必要があります -- 2021-04-21 追記

それでは試してみましょう。

-- Min(ST_Distance())による抽出
SELECT "gt_road_nodes"."id", Min(ST_Distance("gt_from_to"."geom_pt_4326", "gt_road_nodes"."geom_pt_4326"))
FROM "gt_from_to", "gt_road_nodes"
WHERE "gt_from_to"."PointName" IS 'Start'
;
-- 10s

-- VirtualKNNの利用
SELECT * FROM "KNN"
WHERE "f_table_name" = 'gt_road_nodes'
AND "ref_geometry" = (SELECT "geom_pt_4326" FROM "gt_from_to" WHERE "PointName" IS 'Start')
AND "max_items" = 1
;
-- 0.3s

一点あたり10秒だったのが、0.3秒に改善しました！
先ほどの経路探索をVirtualKNN対応にしたものが以下のクエリです。

INSERT OR REPLACE INTO "gt_path"
SELECT 1, "Algorithm", "Options", "NodeFrom", "NodeTo", "Cost", "Geometry"
FROM "routing_net"
WHERE "NodeFrom" == (
  -- 始点に最も近いネットワークノードを検索（KNNを使用）
  SELECT  "fid" FROM "KNN"
  WHERE "f_table_name" = 'gt_road_nodes'
  AND "ref_geometry" = (SELECT "geom_pt_4326" FROM "gt_from_to" WHERE "PointName" IS 'Start')
  AND "max_items" = 1 -- 最近傍のみを取得
)
AND "NodeTo" == (
  -- 終点に最も近いネットワークノードを検索（KNNを使用）
  SELECT "fid" FROM "KNN"
  WHERE "f_table_name" = 'gt_road_nodes'
  AND "ref_geometry" = (SELECT "geom_pt_4326" FROM "gt_from_to" WHERE "PointName" IS 'End')
  AND "max_items" = 1 -- 最近傍のみを取得
)
;
-- 1.4s

2秒未満で結果が得られました。

経路探索のメソッドの変更

経路探索メソッドはデフォルトではDijkstra（ダイクストラ法）ですが、CreateRouting()を実行する際に有効にしていればA*（A-star）も利用できます。経路探索のメソッドを変更するには、以下のUPDATEを行います。

UPDATE "routing_net" SET "Algorithm" = 'A*';

比較的直線的な経路の場合パフォーマンスに差が出るようですが、違いはよくわかりません。

探索結果の取得オプション

VirtualRoutingテーブルのOptionsカラムは、以下のオプションが指定可能です。

UPDATE "routing_net" SET "Options" = 'FULL'; -- 全部返る
UPDATE "routing_net" SET "Options" = 'NO LINKS'; -- リンクが返らない（ルートだけが返る）
UPDATE "routing_net" SET "Options" = 'NO GEOMETRIES'; -- ジオメトリが返らない
UPDATE "routing_net" SET "Options" = 'SIMPLE'; -- リンクもジオメトリも返らない

#2 複数の目的地

以下のクエリで、特定の出発地点から複数の目的地への最短経路の集合を得ることができます。

UPDATE "routing_net" SET "Options" = 'FULL';
SELECT *
FROM "routing_net"
WHERE "NodeFrom" == 500000
AND "NodeTo" == '300000, 400000, 600000, 700000'  -- TEXTでリストを指定
;
-- 35s

最短経路の探索との違いは、NodeToを一つのID（整数）とするか、複数のIDをリスト文字列で指定するかの違いだけです。
先ほどと同様に、ジオテーブルとしてQGISで可視化してみました。

UPDATE "routing_net" SET "Options" = 'NO LINKS';
DELETE FROM "gt_path";
INSERT INTO "gt_path"
SELECT NULL, "Algorithm", "Options", "NodeFrom", "NodeTo", "Cost", "Geometry"
FROM "routing_net"
WHERE "NodeFrom" == 500000
AND "NodeTo" == '300000, 400000, 600000, 700000'  -- TEXTでリストを指定
;

稚内市から、岩内町、留寿都村、遠軽町、弟子屈町へ向かう経路の集合でした。

せっかくなので（？）、こちらもグラフィカルに解を得られるようにしたいと思います。先ほどは、目的地点（ラベル：End）を一点だけ作成しましたが、さらにいくつもの目的地を追加します。札幌ドームに加えて、円山球場、札幌コミュニティドーム（つどーむ）、サッポロテイネスキー場に点を追加しました。

複数のIDリストを生成するには、KNNによる近傍探索を繰り返す必要があるため、再帰クエリを使って実装しました。

UPDATE "routing_net" SET "Options" = 'NO LINKS';
DELETE FROM "gt_path"; -- いったん全て削除
WITH x(n) AS (
  SELECT 1
  UNION
  SELECT n + 1 FROM x
  LIMIT (SELECT COUNT(*) FROM "gt_from_to" WHERE "PointName" IS 'End')
)
INSERT OR REPLACE INTO "gt_path"
SELECT NULL, "Algorithm", "Options", "NodeFrom", "NodeTo", "Cost", "Geometry"
FROM "routing_net"
WHERE "NodeFrom" == (
  -- 始点に最も近いネットワークノードを検索（KNNを使用）
  SELECT  "fid" FROM "KNN"
  WHERE "f_table_name" = 'gt_road_nodes'
  AND "ref_geometry" = (SELECT "geom_pt_4326" FROM "gt_from_to" WHERE "PointName" IS 'Start')
  AND "max_items" = 1 -- 最近傍のみを取得
)
AND "NodeTo" == (
  SELECT GROUP_CONCAT("fid") FROM "KNN", x
  WHERE "f_table_name" = 'gt_road_nodes'
  AND "ref_geometry" = ST_GeometryN(
    (SELECT ST_Collect("geom_pt_4326") FROM "gt_from_to" WHERE "PointName" IS 'End')
    , n)
  AND "max_items" = 1 -- 最近傍のみを取得
)
;
-- 3s

ちゃんとできてますね。

#3 到達可能範囲の把握（Isochrone map）

次のようなクエリで、始点からの距離（コスト）内にあるノードを抽出することができます。

SELECT * FROM "routing_net"
WHERE "NodeFrom" == 10000
AND "Cost" <= 10000
;

結果を抽出すると以下のようになります。

Role	NodeFrom	NodeTo	Cost	Geometry
Solution	10000	1816	8108.150948	BLOB sz=60 GEOMETRY
Solution	10000	1817	8196.418481	BLOB sz=60 GEOMETRY
Solution	10000	1818	8279.085923	BLOB sz=60 GEOMETRY
...	...	...	...	...

サイズ60のジオメトリはポイントジオメトリで、つまりコストの値を持った点群が返ります。
この点群を使って、ConcaveHull（凹包）を使って到達可能範囲を可視化してみます。

-- ジオテーブルの作成
CREATE TABLE "gt_isochrone_map" (
  "PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT, "Cost_within" REAL
);
-- マルチポリゴンになる場合があるので
SELECT AddGeometryColumn('gt_isochrone_map', 'geom_pgm_4326', 4326, 'MULTIPOLYGON');
-- 結果の挿入
INSERT INTO "gt_isochrone_map"
SELECT NULL, 2000, ST_Multi(ST_ConcaveHull(ST_Collect("Geometry"))) -- 点群を内包する凹包ポリゴンを生成
FROM "routing_net"
WHERE "NodeFrom" == (
  -- 始点に最も近いネットワークノードを検索（KNNを使用）
  SELECT "fid" FROM "KNN"
  WHERE "f_table_name" = 'gt_road_nodes'
  AND "ref_geometry" = (SELECT "geom_pt_4326" FROM "gt_from_to" WHERE "PointName" IS 'Start')
  AND "max_items" = 1 -- 最近傍のみを取得
)
AND "Cost" <= 2000
;

札幌駅から2km圏内の凹包の結果です。北大構内の道路は削除されているので、到達可能範囲がポコンとへこんでいます。

せっかくなので（？）、1km毎の凹包ポリゴンを再帰クエリで生成してみます。

DELETE FROM "gt_isochrone_map"; -- いったんレコードを削除
-- 1kmずつの等値線を書く
WITH x(n) AS (SELECT 1 UNION SELECT n + 1 FROM x LIMIT 10)
INSERT INTO "gt_isochrone_map"
SELECT NULL, n * 1000, ST_Multi(ST_ConcaveHull(ST_Collect("Geometry")))
FROM "routing_net", x
WHERE "NodeFrom" == (
  -- 始点に最も近いネットワークノードを検索（KNNを使用）
  SELECT "fid" FROM "KNN"
  WHERE "f_table_name" = 'gt_road_nodes'
  AND "ref_geometry" = (SELECT "geom_pt_4326" FROM "gt_from_to" WHERE "PointName" IS 'Start')
  AND "max_items" = 1 -- 最近傍のみを取得
)
AND "Cost" <= n * 1000
GROUP BY n
;
-- 1min 38s

なんだかカオスな図になりました...

点群を取得する範囲が広がると、クエリにも時間がかかります。
シンプルに点群をコストで色分けするか、点群からラスタを生成する方が傾向がわかりやすいかもしれません。

#4 巡回セールスマン問題（TSP）

VirtualRoutingでは、巡回セールスマン問題も扱うことができます。
巡回セールスマン問題の探索にあたっては、以下のUPDATEでTSPモードへの設定をしておきます。

UPDATE "routing_net" SET "Request" = 'TSP';

そのうえで複数の行き先探索のときと同じクエリを渡すと、複数の目的地を巡回先とした経路探索の結果が得られます。

稚内市を出発し、日本海沿いの経路で留寿都村へ向かい、その後岩内町、内陸を経由して弟子屈町、遠軽町と巡り、オホーツク海を北上して稚内へ向かいます。（ちなみに、距離をコストにしているので、あまり高速道路は利用しないようです笑）

これまで同様に、出発地と目的地をポイントから取得するようにしてみます。

UPDATE "routing_net" SET "Options" = 'NO LINKS'; -- ルートだけを返す
DELETE FROM  "gt_path";
WITH x(n) AS (
  SELECT 1
  UNION
  SELECT n + 1 FROM x
  LIMIT (SELECT COUNT(*) FROM "gt_from_to" WHERE "PointName" IS 'End')
)
INSERT OR REPLACE INTO "gt_path"
SELECT NULL, "Algorithm", "Options", "NodeFrom", "NodeTo", "Cost", "Geometry"
FROM "routing_net"
WHERE "NodeFrom" == (
  -- 始点に最も近いネットワークノードを検索（KNNを使用）
  SELECT  "fid" FROM "KNN"
  WHERE "f_table_name" = 'gt_road_nodes'
  AND "ref_geometry" = (SELECT "geom_pt_4326" FROM "gt_from_to" WHERE "PointName" IS 'Start')
  AND "max_items" = 1 -- 最近傍のみを取得
)
AND "NodeTo" == (
  SELECT GROUP_CONCAT("fid") FROM "KNN", x
  WHERE "f_table_name" = 'gt_road_nodes'
  AND "ref_geometry" = ST_GeometryN(
    (SELECT ST_Collect("geom_pt_4326") FROM "gt_from_to" WHERE "PointName" IS 'End')
    , n)
  AND "max_items" = 1 -- 最近傍のみを取得
)
;
-- 3s

札幌駅から、円山球場、つどーむ、サッポロテイネスキー場、札幌ドームと巡って戻ってくるコースをおすすめされました。個人的には、誰かに観光プランとして聞かれたらチョットマテとつっこみたくなりそうなコース　です^^;

なお、巡回セールスマン問題の探索では、最近傍法の他に遺伝的アルゴリズムを利用することもできます。結果が違うかどうかは... 試してみてください^^

UPDATE "routing_net" SET "Request" = 'TSP'; -- 最近傍法（TSP NNオプションと同じ）
UPDATE "routing_net" SET "Request" = 'TSP NN'; -- 最近傍法（TSPオプションと同じ）
UPDATE "routing_net" SET "Request" = 'TSP GA'; -- 遺伝的アルゴリズム

#5 最短経路（Point to Point）

さて、ここにきて再び最短経路です。

これまでは、始点と終点を端点ノードIDで指定する方法でした。しかし、始点と終点にポイントジオメトリそのものを指定することができるようになったそうです。しかもこの場合、始点および終点に最近接となるリンクに（端点ではなくリンクの途中に）接続し、リンクまでの接続部分＆最初と最後のリンクは（端点からスタートしない場合は）部分コストを計算してくれるという機能です。

クエリの文法は以下のように、最短経路探索の際に"NodeFrom"、"NodeTo"としていたのが"PointFrom"、"PointTo"となり、ポイントジオメトリを直接指定します。

SELECT *
FROM "routing_net"
WHERE "PointFrom" = 始点のポイントジオメトリ
AND "PointTo" = 終点のポイントジオメトリ
;

しかし、、、試してみましたがうまくいきません。（クエリが終わらない。）
「適切に空間インデックスを設定してないと結果が得られないよ」との記述もあり、リンクやノードに空間インデックスを適用してみたのですが、、うーん。。どうもうまくいきません。

ということで、このモードの検証はおあずけにします。うまくできたら追記するか、別に記事を上げたいと思います。（解決できた方がいれば、教えていただけると喜びます。）

おわりに

以前のRoutingモデル（VirtualNetwork）から比べると、ずいぶん幅が広がったなあと思う今回のアップグレードでした。
RDBMSでRoutingといえば、PostGISのpgRoutingが圧倒的な存在感を発揮していますが、ローカル・サーバーレス環境で動作するVirtualRoutingも活躍できる場所があるのではないでしょうか。

SpatiaLite、某所でもロゴが微妙とか情報探せないとか（ホントそうだよね…）いろいろタイムリーに話題に上がっていますが^^;　せいぜい40MBくらいで機能一式持ち歩ける素敵なやつでもあるので、試しに触ってみようと思ってくれる人がいるといいなと思います。

ストアド・プロシージャも利用できるようになったらしいので、今年こそはSpatiaLiteで函館から稚内まで行ける真っ直ぐ道を探したい！…と思っていたのですが、学習コストが間に合いませんでした（泣

ということで、今年もSQL勢として参加してみました。
ではではみなさま、よいお年を！