MySQLのGISを試してみた

この記事は[MicroAd Advent Calendar 2018](https://qiita.com/advent-calendar/2018/microad)の18日目の記事です。

#はじめに
GoogleMapで「カフェ おすすめ」っ検索したら職場からほど近い猫カフェがおすすめ度61%で表示されました。猫を飼ってる情報が筒抜けなのかそうでないのかわからないですが知らない土地でなにか探したいときってGoogleMap、というか位置情報を使ったサービスって便利ですよね。
そういえばMySQLでもGIS機能あったはずだよなぁと思ったのでちょっと試してみました。

#緯度経度から市区町村を検索

環境：MySQL 8.0.13

[s-stat](https://www.e-stat.go.jp/)から東京の境界データをとってきてMySQLにいれてみます。QGISというアプリを使うとこのシェープファイルが開けます。
![スクリーンショット 2018-12-18 11.50.27_QGIS.png](https://qiita-image-store.s3.amazonaws.com/0/21170/5a5b1a74-722c-25a3-b667-37d07fe3c1d9.png)


初めてGIS触ってみましたが、こんな境界データが提供されているんですね。
ちなみにシェープデータがShift_JISなのでこのままQGISを使ってUTF-8に変換しておきます。あとインポートにogr2ogrを使うのでGDAL2.3.2を入れておきます。

空のデータベースを作っておいて、そこを指定するとテーブル作成してインポートしてくれます。上記で変換したデータは/tmp下に置いておきます。

```shell-session
ogr2ogr -f "MySQL" MySQL:"geo_test_db,host=127.0.0.1,user=root,password=hogehoge,port=3306" /tmp/h27ka13.shp
```

問題なく終わればインポートの完了です。

```shell-session
$ mysql -uroot -p geo_test_db;
mysql> desc h27ka13;
+------------+---------------+------+-----+---------+----------------+
| Field      | Type          | Null | Key | Default | Extra          |
+------------+---------------+------+-----+---------+----------------+
| OGR_FID    | int(11)       | NO   | PRI | NULL    | auto_increment |
| SHAPE      | geometry      | NO   | MUL | NULL    |                |
| key_code   | varchar(11)   | YES  |     | NULL    |                |
| pref       | varchar(2)    | YES  |     | NULL    |                |
| city       | varchar(3)    | YES  |     | NULL    |                |
| s_area     | varchar(6)    | YES  |     | NULL    |                |
| pref_name  | varchar(12)   | YES  |     | NULL    |                |
| city_name  | varchar(15)   | YES  |     | NULL    |                |
| s_name     | varchar(96)   | YES  |     | NULL    |                |
| kigo_e     | varchar(3)    | YES  |     | NULL    |                |
| hcode      | decimal(4,0)  | YES  |     | NULL    |                |
| area       | double(19,3)  | YES  |     | NULL    |                |
| perimeter  | double(19,3)  | YES  |     | NULL    |                |
| h27kaxx_   | decimal(6,0)  | YES  |     | NULL    |                |
| h27kaxx_id | decimal(6,0)  | YES  |     | NULL    |                |
| ken        | varchar(2)    | YES  |     | NULL    |                |
| ken_name   | varchar(12)   | YES  |     | NULL    |                |
| sityo_name | varchar(22)   | YES  |     | NULL    |                |
| gst_name   | varchar(15)   | YES  |     | NULL    |                |
| css_name   | varchar(14)   | YES  |     | NULL    |                |
| kihon1     | varchar(4)    | YES  |     | NULL    |                |
| dummy1     | varchar(1)    | YES  |     | NULL    |                |
| kihon2     | varchar(2)    | YES  |     | NULL    |                |
| keycode1   | varchar(9)    | YES  |     | NULL    |                |
| keycode2   | varchar(9)    | YES  |     | NULL    |                |
| area_max_f | varchar(1)    | YES  |     | NULL    |                |
| kigo_d     | varchar(2)    | YES  |     | NULL    |                |
| n_ken      | varchar(2)    | YES  |     | NULL    |                |
| n_city     | varchar(3)    | YES  |     | NULL    |                |
| kigo_i     | varchar(1)    | YES  |     | NULL    |                |
| moji       | varchar(96)   | YES  |     | NULL    |                |
| kbsum      | decimal(4,0)  | YES  |     | NULL    |                |
| jinko      | decimal(10,0) | YES  |     | NULL    |                |
| setai      | decimal(10,0) | YES  |     | NULL    |                |
| x_code     | double(19,5)  | YES  |     | NULL    |                |
| y_code     | double(19,5)  | YES  |     | NULL    |                |
| kcode1     | varchar(7)    | YES  |     | NULL    |                |
+------------+---------------+------+-----+---------+----------------+
37 rows in set (0.00 sec)
```

こんな感じのテーブルになってます。

早速今いる場所がどこなのか検索してみましょう。
とは言っても自分の緯度経度はわからないのでgoogle先生に聞いてみました。

　**緯度: 35.65532  経度: 139.69378**

ここみたいなので境界データの中にこの緯度経度が入っているかどうかST_CONTAINSを使って検索してみます。

```shell-session
SELECT pref_name,city_name,s_name FROM h27ka13 WHERE ST_CONTAINS(SHAPE,ST_GeomFromText('POINT(35.65532 139.69378)',4612));
+-----------+-----------+-----------+
| pref_name | city_name | s_name    |
+-----------+-----------+-----------+
| 東京都    | 渋谷区    | 神泉町    |
+-----------+-----------+-----------+
1 row in set (1.80 sec)
```

ちょっとパフォーマンスが悪い気がしますが緯度経度だけでこんなに簡単に出るものなんですね。

#最寄り駅を探してみる

駅のデータが提供されているサイトがあるので下記からデータを取ってきます。
[駅データ.jp](http://www.ekidata.jp/)

取ってきたデータ用のテーブルを作成してデータをLOADします。

```shell-session
CREATE TABLE `station_location` (
  `pref_code`    int(2) NOT NULL,
  `station_code` int(9) NOT NULL,
  `station_name` varchar(255) NOT NULL,
  `postal_code`  char(8) NOT NULL,
  `town_name`    varchar(255) NOT NULL,
  `lat`          double not null,
  `lon`          double not null,
  `geo_point`       geometry ,
  `line_code`    int(6) NOT NULL,
  `line_name_h`  varchar(255) NOT NULL,
  `company_code` int(5) NOT NULL,
  `company_name` varchar(255) NOT NULL,
  PRIMARY KEY (`station_code`)
) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4;
```

ここでLOADと言いたいところですが、8.0ではオプションなしでmysqlコマンドで接続し、LOAD DATAしようとすると
`ERROR 1148 (42000): The used command is not allowed with this MySQL version`
と怒られてしまうので、--local_infile=1を付ける必要がありました。
ちなみに上記でダウンロードしたファイルは/tmp下に置いておきます。

```shell-session
$ mysql -uroot -p --local_infile=1 geo_test_db
mysql > LOAD DATA LOCAL INFILE '/tmp/station_latlon_utf8_20180611.csv'
REPLACE INTO TABLE station_location
CHARACTER SET utf8mb4
FIELDS TERMINATED BY ',' ENCLOSED BY '"'
LINES TERMINATED BY '\r\n' STARTING BY ''
IGNORE 1 LINES
(@pref_code, @station_code, @station_name, @postal_code, @town_name, @lat, @lon, @line_code, @line_name_h, @company_code, @company_name)
SET
pref_code = @pref_code,
station_code = @station_code,
station_name = @station_name,
postal_code = @postal_code,
town_name = @town_name,
lat = @lat,
lon = @lon,
geo_point = ST_GeomFromText(CONCAT('POINT(', @lon, ' ', @lat, ')')),
line_code = @line_code,
line_name_h = @line_name_h,
company_code = @company_code,
company_name = @company_name
;
```

確認してみると、

```shell-session

mysql> select * from station_location limit 10;
+-----------+--------------+-----------------+-------------+------------------------------------------------+-----------+------------+---------------------------+-----------+------------------------------------+--------------+--------------+
| pref_code | station_code | station_name    | postal_code | town_name                                      | lat       | lon        | latlon                    | line_code | line_name_h                        | company_code | company_name |
+-----------+--------------+-----------------+-------------+------------------------------------------------+-----------+------------+---------------------------+-----------+------------------------------------+--------------+--------------+
|         1 |      1110101 | 函館            | 040-0063    | 北海道函館市若松町１２-１３                    | 41.773709 | 140.726413 |        xρ�D@u�y�>�a@       |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
|         1 |      1110102 | 五稜郭          | 041-0813    | 函館市亀田本町                                 | 41.803557 | 140.733539 |        f�����D@���&y�a@          |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
|         1 |      1110103 | 桔梗            | 041-0801    | 北海道函館市桔梗３丁目４１-３６                | 41.846457 | 140.722952 |        2<�X�D@}�;l"�a@      |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
|         1 |      1110104 | 大中山          | 041-1121    | 亀田郡七飯町大字大中山                         | 41.864641 |  140.71358 |        ��h���D@�b��Ֆa@          |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
|         1 |      1110105 | 七飯            | 041-1111    | 亀田郡七飯町字本町                             | 41.886971 | 140.688556 |        �D��D@ޮ��a@        |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
|         1 |      1110106 | 新函館北斗      | 041-1242    | 北海道北斗市市渡                               |   41.9054 | 140.646525 |        2�%��D@|a2U��a@       |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
|         1 |      1110107 | 仁山            | 041-1101    | 亀田郡七飯町字仁山                             | 41.930011 | 140.635183 |        ��
�D@5LkS�a@       |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
|         1 |      1110108 | 大沼            | 041-1354    | 北海道亀田郡七飯町字大沼町４                   | 41.971954 | 140.669347 |        ���h�D@�UfJk�a@       |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
|         1 |      1110109 | 大沼公園        | 041-1354    | 北海道亀田郡七飯町字大沼町８５                 | 41.980958 | 140.669758 |        �_��D@�GT�n�a@       |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |
k E@z�sѐ�a@      |     11101 | JR函館本線(函館～長万部)           |            1 | JR北海道     |           | 42.003267 | 140.642678 |        ��
+-----------+--------------+-----------------+-------------+------------------------------------------------+-----------+------------+---------------------------+-----------+------------------------------------+--------------+--------------+
10 rows in set (0.00 sec)
```

こんな感じで入ってますね。

geo_pointカラムをgeometryデータ型にしていてLOADのときにST_GeomFromText(CONCAT('POINT(', @lon, ' ', @lat, ')'))とすることで座標空間内の単一の位置を表すデータを入れています。

INSERT文だとこんな感じで入れられます。

```shell-session
INSERT INTO latlon VALUES (ST_GeomFromText('POINT(140.726413 41.773709)'));
```

で、ここのPOINTでの指定ですが緯度・経度で入れたいのが人の常なものかどうかlatlon or lonlat、そもそもlon or lngなのかとかここはいろんな意見があると思いますが後述するST_Distance_Sphereを使う場合だと経度・緯度じゃないと使えないのでした。

先程の緯度経度を使ってみます。ちなみに弊社の最寄り駅は神泉なのでそこまでの距離を測ってみます。

```shell-session

select ST_Distance_Sphere(geo_point,ST_GeomFromText('POINT(139.69378 35.65532)'),4612) from station_location where station_name = '神泉';
+---------------------------------------------------------------------------------+
| ST_Distance_Sphere(geo_point,ST_GeomFromText('POINT(139.69378 35.65532)'),4612) |
+---------------------------------------------------------------------------------+
|                                                             0.15542853757321037 |
+---------------------------------------------------------------------------------+
1 row in set (0.02 sec)
```

155mという結果がでました。GoogleMapで見てみると196m。まあ誤差なんですかね。。
![スクリーンショット 2018-12-18 12.06.17.png](https://qiita-image-store.s3.amazonaws.com/0/21170/92d101bd-4951-5186-b0c3-426fde468446.png)

次に緯度経度から最寄り駅を出してみたいと思います。今回は駅データの件数もさほど多くないので全駅との距離を測ってみてソートしてみます。きっと他にもっといい方法があると思いますが。

```shell-session

set @here = ST_GeomFromText('POINT(139.69378 35.65532)');
select company_name,line_name_h,station_name,ST_Distance_Sphere(geo_point,@here,4612) from station_location ORDER BY ST_Distance_Sphere(geo_point,@here, 4612) limit 10;
+-----------------+-----------------------------+--------------+------------------------------------------+
| company_name    | line_name_h                 | station_name | ST_Distance_Sphere(geo_point,@here,4612) |
+-----------------+-----------------------------+--------------+------------------------------------------+
| 京王電鉄        | 京王井の頭線                | 神泉         |                      0.15542853757321037 |
| 京王電鉄        | 京王井の頭線                | 渋谷         |                       0.5377475250715124 |
| 東京メトロ      | 東京メトロ半蔵門線          | 渋谷         |                       0.5602737574619313 |
| 東京メトロ      | 東京メトロ銀座線            | 渋谷         |                       0.5602737574619313 |
| JR東日本        | JR成田エクスプレス          | 渋谷         |                       0.5653590815730241 |
| JR東日本        | JR埼京線                    | 渋谷         |                       0.5653590815730241 |
| JR東日本        | JR山手線                    | 渋谷         |                       0.5653590815730241 |
| JR東日本        | JR湘南新宿ライン            | 渋谷         |                       0.5653590815730241 |
| 東急電鉄        | 東急田園都市線              | 渋谷         |                       0.5888727700854771 |
| 東京メトロ      | 東京メトロ副都心線          | 渋谷         |                       0.6593624396346356 |
+-----------------+-----------------------------+--------------+------------------------------------------+
10 rows in set (0.15 sec)
```

まあ神泉と渋谷ですが、微妙に緯度経度が違うんですね。

さて、最寄り駅とは言ってもあくまで直線距離なので例えば緯度:35.7951367,経度:140.7409544こんな場所に置き去りにされたとして、

```shell-session

set @here = ST_GeomFromText('POINT(140.7409544 35.7951367)');
select company_name,line_name_h,station_name,ST_Distance_Sphere(geo_point,@here,4612) from station_location ORDER BY ST_Distance_Sphere(geo_point,@here, 4612) limit 10;
+--------------+----------------+--------------+------------------------------------------+
| company_name | line_name_h    | station_name | ST_Distance_Sphere(geo_point,@here,4612) |
+--------------+----------------+--------------+------------------------------------------+
| JR東日本     | JR成田線       | 下総豊里     |                       1.3145268599324935 |
| JR東日本     | JR成田線       | 椎柴         |                       2.6785846917676444 |
| JR東日本     | JR成田線       | 下総橘       |                        3.156841913995182 |
| JR東日本     | JR総武本線     | 猿田         |                       3.8509173773004597 |
| JR東日本     | JR総武本線     | 倉橋         |                        4.936536910704812 |
| JR東日本     | JR総武本線     | 松岸         |                        5.715325659747645 |
| JR東日本     | JR成田線       | 松岸         |                        5.715325659747645 |
| JR東日本     | JR総武本線     | 飯岡         |                        6.476201274210973 |
| JR東日本     | JR成田線       | 笹川         |                        6.622113877434571 |
| JR東日本     | JR総武本線     | 銚子         |                        7.743079285844181 |
+--------------+----------------+--------------+------------------------------------------+
10 rows in set (0.15 sec)
```

なるほど茨城か、下総豊里（しもうさとよさと）までは直線で1.3kmだからまあ実際は2kmくらいで2,30分歩けば駅につけそうだなぁ。などと安易に考えると大変です。
![スクリーンショット 2018-12-18 12.37.20.png](https://qiita-image-store.s3.amazonaws.com/0/21170/4742ee81-d4ae-a2bd-5441-888394268e59.png)
![スクリーンショット 2018-12-18 12.37.03.png](https://qiita-image-store.s3.amazonaws.com/0/21170/86ac441f-1b84-40a1-b092-8590dcad9ba6.png)

利根川が。。ちなみにここにはバス停があり銚子駅行のバスがあるようでした。

#まとめ

GIS系はまったく触ったことなかったんですが、この程度であればMySQLだけでも完結するんですね。（いろいろ雑な感じでしたが）おそらくここでやっていることは5.7でもできると思いますが8.0では更にGIS関連の機能が増えているとのことです。

と、いろいろ調べてたらRedisでも距離的なものは簡単にできるっぽい。機会があればパフォーマンスの比較検証してみたいです。