あの日みた学校の最寄駅を僕達はまだ知らない。だからSQLで探してみる。

## はじめに

MIERUNEでバックエンドを担当してる鈴木です。
どこかで見たようなタイトルですが、お気になさらず。

さて、我が家には中学校に通う子供がいますが、早いものでもう中学二年生。来シーズンの冬が終わる頃にはどこかの高校を受験してる予定です。

現在、北海道に住んでいる私ですが、出身は北海道ではありません。なので、子供が行きたい志望校がどの辺りにあるのか、またどんな交通手段で通学するのか、土地勘があまりなく、学校がありそうな場所もいまいちよくわかっていません。

そんな現実に軽い絶望感を抱いた私は、地理位置情報を取り扱う会社に所属してる身として、知りたい事を可視化してみたいと思いました。

## データ可視化のためのおおまかな工程

最終的なアウトプットがこんな感じだとして、

![スクリーンショット 2021-12-10 20.52.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/18599/e6a16f39-8cd2-add1-415d-dc7f8d7b6ff4.png)


この図を再現するために、このような工程を考えてみました。

- 学校と駅のデータをどこかから取得
- 学校と駅のデータをPostGISに格納
- PostGISで学校と駅の距離を求めるSQLを発行する
    - 二点間の距離の最も短いデータだけピックアップする
- QGISで学校・駅・距離のデータを可視化する

ということで、それぞれの詳細を書いてきます。

## 学校と駅のデータをどこかから取得

### 国土数値情報を探索する

検索元データは、国土数値情報から探ってみる事にします。

- 国土数値情報・鉄道データ :  [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v2_3.html](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v2_3.html)
- 国土数値情報・学校のデータ : [https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P29.html](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P29.html)
（北海道のデータだけ取得）

国土数値情報で得られるデータには ESRI SHAPEファイルが含まれています。
SHAPEファイルはArcGISやQGISなどのGISツールで表示できますが、今回は PostGIS を使って処理を行いたいので、SQLに変換する必要があります。

### SHAPEファイルをSQLに変換

ということで、PostGIS の付属ツール、`shp2pgsql` を使って、SHAPEファイルをSQLに変換してみます。コマンドは以下の通り。

```bash
% shp2pgsql -W cp932 -D -I -s 4326 N02-20_Station.shp staion > station.sql
% shp2pgsql -W cp932 -D -I -s 4326 P29-13_01.shp hokkaido_school > hokkaido_school.sql
```

これでSHAPEファイルから SQL を生成することができました。

## 学校と駅のデータをPostGISに格納

### SQLをPostGISに投入

引き続き `shp2pgsql` で生成したSQLファイルを PostGIS に格納してみます。

```bash
% psql -U postgres -h localhost -d sample < hokkaido_school.sql
% psql -U postgres -h localhost -d sample < station.sql
```

生成したSQL は `create table` 文も含まれてるので、データだけでなく、テーブルも同時に作成されます。

```sql
sample=# \dt
              List of relations
 Schema |      Name       | Type  |  Owner   
--------+-----------------+-------+----------
 public | hokkaido_school | table | postgres
 public | spatial_ref_sys | table | postgres
 public | station         | table | postgres
(3 rows)
```

これで、駅、および、学校のデータのテーブルができました。

それぞれのテーブルの中身を少し覗いてみます。

```sql
sample=# select * from hokkaido_school limit 10;
 gid | p29_001 | p29_002 | p29_003 | p29_004 |        p29_005         |           p29_006            | p29_007 |                        geom                        
-----+---------+---------+---------+---------+------------------------+------------------------------+---------+----------------------------------------------------
   1 | 01101   | 16      | 16001   | 16001   | 桑園小学校ひまわり分校 | 北11条西13（市立札幌病院内） |       3 | 0101000020E61000005D6DC5FEB2AA6140537AA69718894540 
   2 | 01101   | 16      | 16001   | 16001   | 円山小学校             | 北1条西25-1-8                |       3 | 0101000020E61000006310786F39AA614055E1386765874540 
   3 | 01101   | 16      | 16001   | 16001   | 宮の森小学校           | 宮の森4条6-1-1               |       3 | 0101000020E6100000CBDCECE6AEA96140C830AE7431884540 
   4 | 01101   | 16      | 16001   | 16001   | 桑園小学校             | 北8条西17                    |       3 | 0101000020E6100000A471A8DF85AA6140A7CF0EB8AE884540 
   5 | 01101   | 16      | 16001   | 16001   | 三角山小学校           | 宮の森4条11-4-1              |       3 | 0101000020E6100000D9794BAF76A96140C353ED5F50874540 
   6 | 01101   | 16      | 16001   | 16001   | 山鼻小学校             | 南14条西10-1                 |       3 | 0101000020E6100000A4E82EB107AB614021B375ED4B854540 
   7 | 01101   | 16      | 16001   | 16001   | 山鼻南小学校           | 南29条西12-1-1               |       3 | 0101000020E6100000803F8C3903AB61407634CD09B2824540 
   8 | 01101   | 16      | 16001   | 16001   | 資生館小学校           | 南3条西7                     |       3 | 0101000020E6100000D23AAA9A20AB61407EA99F3715874540 
   9 | 01101   | 16      | 16001   | 16001   | 大倉山小学校           | 宮の森3条13-6-20             |       3 | 0101000020E61000002E1C08C982A961401650A8A78F864540 
  10 | 01101   | 16      | 16001   | 16001   | 中央小学校             | 大通東6-12                   |       3 | 0101000020E610000048C9C973B5AB61407DFBDB2FFE874540 

sample=# select * from station limit 10;
 gid  | n02_001 | n02_002 |       n02_003        |      n02_004       |          n02_005           |                                                     geom                                                     
------+---------+---------+----------------------+--------------------+----------------------------+--------------------------------------------------------------------------------------------------------------
 3687 | 11      | 2       | 境線                 | 西日本旅客鉄道     | 中浜                       | 0105000020E6100000010000000102000000020000003F912749D7A76040376C5B94D9C04140D28C45D3D9A760407094BC3AC7C04140 
 5887 | 11      | 2       | 飯田線               | 東海旅客鉄道       | 金野                       | 0105000020E61000000100000001020000000200000005FA449E243A614021CD58349DB54140F20C1AFA273A61402F51BD35B0B54140 
 6906 | 11      | 2       | 鹿島線               | 東日本旅客鉄道     | 潮来                       | 0105000020E610000001000000010200000002000000840D4FAF94916140B2D7BB3FDEF74140882EA86F999161400EA14ACD1EF84140 
    1 | 11      | 2       | 指宿枕崎線           | 九州旅客鉄道       | 二月田                     | 0105000020E6100000010000000102000000020000001361C3D32B5460407E8CB96B09413F4057EC2FBB275460407EE36BCF2C413F40 
    2 | 23      | 5       | 沖縄都市モノレール線 | 沖縄都市モノレール | 古島                       | 0105000020E6100000010000000102000000020000002C7DE882FAEC5F4034A2B437F83A3A404209336DFFEC5F4050AA7D3A1E3B3A40 
    3 | 24      | 5       | 東京臨海新交通臨海線 | ゆりかもめ         | お台場海浜公園             | 0105000020E610000001000000010200000002000000D74CBED9E6786140069E7B0F97D0414045F0BF95EC786140713D0AD7A3D04140 
    4 | 24      | 5       | 東京臨海新交通臨海線 | ゆりかもめ         | 東京国際クルーズターミナル | 0105000020E6100000010000000102000000020000004AC36214BF78614063A19DD37FCF41406A956071BB786140A24E797491CF4140 
    5 | 24      | 5       | 東京臨海新交通臨海線 | ゆりかもめ         | テレコムセンター           | 0105000020E610000001000000010200000002000000B21188D7F578614089EFC4AC17CF4140CD237F30F078614065FCFB8C0BCF4140 
    6 | 24      | 5       | 東京臨海新交通臨海線 | ゆりかもめ         | 汐留                       | 0105000020E610000001000000010200000002000000F8A57EDE54786140A20BEA5BE6D441407E3A1E335078614070B1A206D3D44140 
    7 | 24      | 5       | 東京臨海新交通臨海線 | ゆりかもめ         | 台場                       | 0105000020E61000000100000001020000000200000084471B47AC78614072E1404816D04140693524EEB1786140B9AAECBB22D04140 

(10 rows)
```

ひとまず大丈夫そうですね。

### 距離を求めるためテーブルに Geography 列を追加

テーブルに定義されてる geom 列の型は Geometry 型のデータですが、二点間の距離を求めるためには、Geography 型の列を用いるほうが計算が速くなり、求めている結果が得られます。

そこで、駅と学校それぞれのテーブルに、距離計算用の Geography 列を追加します。

なお、駅テーブルの geom 列は MULTILINESTRING としてデータが格納されてるので、LINESTRING の中心座標である POINT のデータを求めて Geography 列に追加することにします。

```sql
-- 距離の計算を行いたいので、geography 型の列をテーブルに追加。
sample=# alter table hokkaido_school add column geog geography(POINT, 4326);
sample=# alter table station add column geog geography(POINT, 4326);

-- 登録済みの geometry 型の列より、geography 型の列にデータを追加。
-- station については、MULTILINESTRING の中心座標だけ欲しいので、
-- ST_Centroid 関数により求めた中心座標を格納。
sample=# update hokkaido_school set geog = geom::geography;
sample=# update station set geog = ST_Centroid(geom)::geography;

-- 一応インデックスを作成する
sample=# CREATE INDEX idx_hokkaido_school_geo ON hokkaido_school USING gist (geog);
sample=# CREATE INDEX idx_station_geo ON station USING gist (geog);
```

## PostGISで学校と駅の距離を求めるSQLを発行する

前準備が終わったので、学校と駅の最短距離を求める SQL を書いてみます。
（二点間の直線距離だけ。道路の経路については考慮しません）

学校はひとまず私が見つけたい「高校」だけの検索とします。

```sql
-- p29_003 列を '16004' のコードで検索すると高校の検索結果だけ出てくる。
SELECT 
a.gid, 
a.p29_005 AS school_name, 
closest_pt.line_name, closest_pt.station_name, closest_pt.distance
FROM hokkaido_school AS a 
CROSS JOIN LATERAL (
    SELECT 
    n02_003 as line_name, 
    n02_005 as station_name, 
    a.geog <-> b.geog as distance,
    b.geom,
    b.geog
    FROM station AS b 
    ORDER BY a.geog <-> b.geog LIMIT 1
) AS closest_pt
WHERE a.p29_003 = '16004'
ORDER BY a.gid
```

### 二点間の距離の最も短いデータだけピックアップする

SQL は学校を求めるメインクエリ、駅の情報を求めるサブクエリを、LATERAL でクロスジョインを行います。

サブクエリ中の`a.geog <-> b.geog` の部分、ここでメインクエリで求めた各学校と駅テーブルの行、それぞれで持っている Geography 列のデータ同士の距離を求めます。

その後、`ORDER BY a.geog <-> b.geog` の句により、学校ごと駅ごとに二点間の距離を求めた集合のうち最も距離が近い順番でソートを行い、`LIMIT 1` の句で求めた集合の先頭の1行、すなわち、ある学校から最も近い駅を示す行を抽出します。

### 検索結果

このSQLの結果はこのようになります。

```sql
gid  |         school_name          |       line_name        |  station_name  |     distance     
------+------------------------------+------------------------+----------------+------------------
   29 | 札幌旭丘高等学校             | 山鼻西線               | 西線16条       | 1434.53066067268
   30 | 札幌聖心女子学院高等学校     | 東西線                 | 西28丁目       | 2476.45511720429
   31 | 札幌西高等学校               | 東西線                 | 西28丁目       | 1247.31568312826
   32 | 札幌静修高等学校             | 山鼻線                 | 静修学園前     | 69.8133121630928
   33 | 札幌南高等学校               | 山鼻線                 | 山鼻19条       | 194.644220875434
   34 | 札幌龍谷学園高等学校         | 東西線                 | 西18丁目       | 649.425190390717
   35 | 北星学園女子高等学校         | 山鼻西線               | 西線6条        | 391.189208366334
   88 | 札幌工業高等学校             | 札沼線                 | 八軒           | 1053.42633862242
   89 | 札幌国際情報高等学校         | 函館線                 | 発寒           | 2791.86360601168
   90 | 札幌篠路高等学校             | 札沼線                 | 拓北           | 1010.99863059946
   91 | 札幌新川高等学校             | 札沼線                 | 新川           | 2438.85888628389

...(略)...
```

Geography 型のデータ同士を `<->` オペレーターで比較することで求めた距離は、単位がメートルとして表現されます。

検索結果のパッと見は、なんか良さそうな感じですね。

ということで、ここまで得られたデータを QGIS を用いて可視化してみます。

## QGISで学校・駅・距離のデータを可視化する

QGISを起動し、

- 駅のデータ
- 路線のデータ
- 「高校」で絞り込みをおこなった学校のデータ
- 背景地図にOpenStreetMap

これらをそれぞれレイヤに貼り付けます。

そして、先に発行した SQL を少し加工し、各学校の最短距離にある駅までを結ぶ LINE 情報を追加。

この SQL を QGIS の DBマネージャー から発行し、得られた結果を QGIS のレイヤとして表示してみます。

```sql
SELECT 
a.gid, a.p29_001, a.p29_003, 
a.p29_005 AS school_name, 
closest_pt.line_name, closest_pt.station_name, closest_pt.distance,
ST_MakeLine(a.geom, closest_pt.geog::geometry)
FROM hokkaido_school AS a 
CROSS JOIN LATERAL (
    SELECT 
    n02_003 as line_name, 
    n02_005 as station_name, 
    a.geog <-> b.geog as distance,
    b.geom,
    b.geog
    FROM station AS b 
    ORDER BY a.geog <-> b.geog LIMIT 1
) AS closest_pt
WHERE a.p29_003 = '16004'
```

すると、OpenStreetMap のベース地図の上に、このようなポリゴンが描写されます。

![スクリーンショット 2021-12-10 18.45.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/18599/6fce2db3-ad9b-1e13-275f-3ab70a0c8f3b.png)


イメージ通りの絵が出てきました。

今回行った可視化では、学校と駅を直線で結んだだけなので、各学校が最寄駅としている場所とは一致しない所もあります。妙に線が長いところについてはスクールバスで通ったり、近所に寮があるのかもしれません。バスのデータを入手してデータベースに混ぜれば最寄りのバス停もわかると思います。

ただ、学校の土地勘がない私にとっては、どの辺りにどんな学校があり、近そうな駅はどこか、おおよそわかったので、ひとまずここまでにしたいと思います。

### おまけ

今回求めたデータは札幌近辺のみでなく、北海道全体の高校を対象としてるので、地図のズームを引いてみると、こんな感じの絵が出てきます。北海道はでっかいどう。

![スクリーンショット 2021-12-10 18.45.46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/18599/14dc27f3-98e7-e947-54ce-6a7880a9ae34.png)


ということで、地理位置情報に関するいろいろな悩みを解決する会社・MIERUNE を今後とも御贔屓いただけますようよろしくお願いいたします。

## 参考


- SQLはこちらの記事を参考にさせていただきました。[https://qiita.com/n_slender/items/efe3ccad0021334d828a](https://qiita.com/n_slender/items/efe3ccad0021334d828a)
- [PostGIS Distance Operators](https://postgis.net/docs/geometry_distance_knn.html)